1 2 #include <../src/mat/impls/nest/matnestimpl.h> /*I "petscmat.h" I*/ 3 #include <../src/mat/impls/aij/seq/aij.h> 4 #include <petscsf.h> 5 6 static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]); 7 static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left); 8 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*); 9 10 /* private functions */ 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatNestGetSizes_Private" 13 static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N) 14 { 15 Mat_Nest *bA = (Mat_Nest*)A->data; 16 PetscInt i,j; 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 *m = *n = *M = *N = 0; 21 for (i=0; i<bA->nr; i++) { /* rows */ 22 PetscInt sm,sM; 23 ierr = ISGetLocalSize(bA->isglobal.row[i],&sm);CHKERRQ(ierr); 24 ierr = ISGetSize(bA->isglobal.row[i],&sM);CHKERRQ(ierr); 25 *m += sm; 26 *M += sM; 27 } 28 for (j=0; j<bA->nc; j++) { /* cols */ 29 PetscInt sn,sN; 30 ierr = ISGetLocalSize(bA->isglobal.col[j],&sn);CHKERRQ(ierr); 31 ierr = ISGetSize(bA->isglobal.col[j],&sN);CHKERRQ(ierr); 32 *n += sn; 33 *N += sN; 34 } 35 PetscFunctionReturn(0); 36 } 37 38 /* operations */ 39 #undef __FUNCT__ 40 #define __FUNCT__ "MatMult_Nest" 41 static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y) 42 { 43 Mat_Nest *bA = (Mat_Nest*)A->data; 44 Vec *bx = bA->right,*by = bA->left; 45 PetscInt i,j,nr = bA->nr,nc = bA->nc; 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 for (i=0; i<nr; i++) {ierr = VecGetSubVector(y,bA->isglobal.row[i],&by[i]);CHKERRQ(ierr);} 50 for (i=0; i<nc; i++) {ierr = VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 51 for (i=0; i<nr; i++) { 52 ierr = VecZeroEntries(by[i]);CHKERRQ(ierr); 53 for (j=0; j<nc; j++) { 54 if (!bA->m[i][j]) continue; 55 /* y[i] <- y[i] + A[i][j] * x[j] */ 56 ierr = MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);CHKERRQ(ierr); 57 } 58 } 59 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);CHKERRQ(ierr);} 60 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "MatMultAdd_Nest" 66 static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z) 67 { 68 Mat_Nest *bA = (Mat_Nest*)A->data; 69 Vec *bx = bA->right,*bz = bA->left; 70 PetscInt i,j,nr = bA->nr,nc = bA->nc; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 for (i=0; i<nr; i++) {ierr = VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);CHKERRQ(ierr);} 75 for (i=0; i<nc; i++) {ierr = VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 76 for (i=0; i<nr; i++) { 77 if (y != z) { 78 Vec by; 79 ierr = VecGetSubVector(y,bA->isglobal.row[i],&by);CHKERRQ(ierr); 80 ierr = VecCopy(by,bz[i]);CHKERRQ(ierr); 81 ierr = VecRestoreSubVector(y,bA->isglobal.row[i],&by);CHKERRQ(ierr); 82 } 83 for (j=0; j<nc; j++) { 84 if (!bA->m[i][j]) continue; 85 /* y[i] <- y[i] + A[i][j] * x[j] */ 86 ierr = MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);CHKERRQ(ierr); 87 } 88 } 89 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);CHKERRQ(ierr);} 90 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatMultTranspose_Nest" 96 static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y) 97 { 98 Mat_Nest *bA = (Mat_Nest*)A->data; 99 Vec *bx = bA->left,*by = bA->right; 100 PetscInt i,j,nr = bA->nr,nc = bA->nc; 101 PetscErrorCode ierr; 102 103 PetscFunctionBegin; 104 for (i=0; i<nr; i++) {ierr = VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 105 for (i=0; i<nc; i++) {ierr = VecGetSubVector(y,bA->isglobal.col[i],&by[i]);CHKERRQ(ierr);} 106 for (j=0; j<nc; j++) { 107 ierr = VecZeroEntries(by[j]);CHKERRQ(ierr); 108 for (i=0; i<nr; i++) { 109 if (!bA->m[i][j]) continue; 110 /* y[j] <- y[j] + (A[i][j])^T * x[i] */ 111 ierr = MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);CHKERRQ(ierr); 112 } 113 } 114 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 115 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);CHKERRQ(ierr);} 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "MatMultTransposeAdd_Nest" 121 static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z) 122 { 123 Mat_Nest *bA = (Mat_Nest*)A->data; 124 Vec *bx = bA->left,*bz = bA->right; 125 PetscInt i,j,nr = bA->nr,nc = bA->nc; 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 for (i=0; i<nr; i++) {ierr = VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 130 for (i=0; i<nc; i++) {ierr = VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);CHKERRQ(ierr);} 131 for (j=0; j<nc; j++) { 132 if (y != z) { 133 Vec by; 134 ierr = VecGetSubVector(y,bA->isglobal.col[j],&by);CHKERRQ(ierr); 135 ierr = VecCopy(by,bz[j]);CHKERRQ(ierr); 136 ierr = VecRestoreSubVector(y,bA->isglobal.col[j],&by);CHKERRQ(ierr); 137 } 138 for (i=0; i<nr; i++) { 139 if (!bA->m[i][j]) continue; 140 /* z[j] <- y[j] + (A[i][j])^T * x[i] */ 141 ierr = MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);CHKERRQ(ierr); 142 } 143 } 144 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 145 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);CHKERRQ(ierr);} 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "MatTranspose_Nest" 151 static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B) 152 { 153 Mat_Nest *bA = (Mat_Nest*)A->data, *bC; 154 Mat C; 155 PetscInt i,j,nr = bA->nr,nc = bA->nc; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 if (reuse == MAT_REUSE_MATRIX && A == *B && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place"); 160 161 if (reuse == MAT_INITIAL_MATRIX || A == *B) { 162 Mat *subs; 163 IS *is_row,*is_col; 164 165 ierr = PetscCalloc1(nr * nc,&subs);CHKERRQ(ierr); 166 ierr = PetscMalloc2(nr,&is_row,nc,&is_col);CHKERRQ(ierr); 167 ierr = MatNestGetISs(A,is_row,is_col);CHKERRQ(ierr); 168 if (reuse == MAT_REUSE_MATRIX) { 169 for (i=0; i<nr; i++) { 170 for (j=0; j<nc; j++) { 171 subs[i + nr * j] = bA->m[i][j]; 172 } 173 } 174 } 175 176 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);CHKERRQ(ierr); 177 ierr = PetscFree(subs);CHKERRQ(ierr); 178 ierr = PetscFree2(is_row,is_col);CHKERRQ(ierr); 179 } else { 180 C = *B; 181 } 182 183 bC = (Mat_Nest*)C->data; 184 for (i=0; i<nr; i++) { 185 for (j=0; j<nc; j++) { 186 if (bA->m[i][j]) { 187 ierr = MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));CHKERRQ(ierr); 188 } else { 189 bC->m[j][i] = NULL; 190 } 191 } 192 } 193 194 if (reuse == MAT_INITIAL_MATRIX || A != *B) { 195 *B = C; 196 } else { 197 ierr = MatHeaderMerge(A, &C);CHKERRQ(ierr); 198 } 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatNestDestroyISList" 204 static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list) 205 { 206 PetscErrorCode ierr; 207 IS *lst = *list; 208 PetscInt i; 209 210 PetscFunctionBegin; 211 if (!lst) PetscFunctionReturn(0); 212 for (i=0; i<n; i++) if (lst[i]) {ierr = ISDestroy(&lst[i]);CHKERRQ(ierr);} 213 ierr = PetscFree(lst);CHKERRQ(ierr); 214 *list = NULL; 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "MatDestroy_Nest" 220 static PetscErrorCode MatDestroy_Nest(Mat A) 221 { 222 Mat_Nest *vs = (Mat_Nest*)A->data; 223 PetscInt i,j; 224 PetscErrorCode ierr; 225 226 PetscFunctionBegin; 227 /* release the matrices and the place holders */ 228 ierr = MatNestDestroyISList(vs->nr,&vs->isglobal.row);CHKERRQ(ierr); 229 ierr = MatNestDestroyISList(vs->nc,&vs->isglobal.col);CHKERRQ(ierr); 230 ierr = MatNestDestroyISList(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 231 ierr = MatNestDestroyISList(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 232 233 ierr = PetscFree(vs->row_len);CHKERRQ(ierr); 234 ierr = PetscFree(vs->col_len);CHKERRQ(ierr); 235 236 ierr = PetscFree2(vs->left,vs->right);CHKERRQ(ierr); 237 238 /* release the matrices and the place holders */ 239 if (vs->m) { 240 for (i=0; i<vs->nr; i++) { 241 for (j=0; j<vs->nc; j++) { 242 ierr = MatDestroy(&vs->m[i][j]);CHKERRQ(ierr); 243 } 244 ierr = PetscFree(vs->m[i]);CHKERRQ(ierr); 245 } 246 ierr = PetscFree(vs->m);CHKERRQ(ierr); 247 } 248 ierr = PetscFree(A->data);CHKERRQ(ierr); 249 250 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);CHKERRQ(ierr); 251 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);CHKERRQ(ierr); 252 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);CHKERRQ(ierr); 253 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);CHKERRQ(ierr); 254 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);CHKERRQ(ierr); 255 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);CHKERRQ(ierr); 256 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);CHKERRQ(ierr); 257 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);CHKERRQ(ierr); 258 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);CHKERRQ(ierr); 259 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 263 #undef __FUNCT__ 264 #define __FUNCT__ "MatAssemblyBegin_Nest" 265 static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type) 266 { 267 Mat_Nest *vs = (Mat_Nest*)A->data; 268 PetscInt i,j; 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 for (i=0; i<vs->nr; i++) { 273 for (j=0; j<vs->nc; j++) { 274 if (vs->m[i][j]) { 275 ierr = MatAssemblyBegin(vs->m[i][j],type);CHKERRQ(ierr); 276 if (!vs->splitassembly) { 277 /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested 278 * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was 279 * already performing an assembly, but the result would by more complicated and appears to offer less 280 * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an 281 * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives". 282 */ 283 ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr); 284 } 285 } 286 } 287 } 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatAssemblyEnd_Nest" 293 static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type) 294 { 295 Mat_Nest *vs = (Mat_Nest*)A->data; 296 PetscInt i,j; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 for (i=0; i<vs->nr; i++) { 301 for (j=0; j<vs->nc; j++) { 302 if (vs->m[i][j]) { 303 if (vs->splitassembly) { 304 ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr); 305 } 306 } 307 } 308 } 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "MatNestFindNonzeroSubMatRow" 314 static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B) 315 { 316 PetscErrorCode ierr; 317 Mat_Nest *vs = (Mat_Nest*)A->data; 318 PetscInt j; 319 Mat sub; 320 321 PetscFunctionBegin; 322 sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */ 323 for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j]; 324 if (sub) {ierr = MatSetUp(sub);CHKERRQ(ierr);} /* Ensure that the sizes are available */ 325 *B = sub; 326 PetscFunctionReturn(0); 327 } 328 329 #undef __FUNCT__ 330 #define __FUNCT__ "MatNestFindNonzeroSubMatCol" 331 static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B) 332 { 333 PetscErrorCode ierr; 334 Mat_Nest *vs = (Mat_Nest*)A->data; 335 PetscInt i; 336 Mat sub; 337 338 PetscFunctionBegin; 339 sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */ 340 for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col]; 341 if (sub) {ierr = MatSetUp(sub);CHKERRQ(ierr);} /* Ensure that the sizes are available */ 342 *B = sub; 343 PetscFunctionReturn(0); 344 } 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "MatNestFindIS" 348 static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found) 349 { 350 PetscErrorCode ierr; 351 PetscInt i; 352 PetscBool flg; 353 354 PetscFunctionBegin; 355 PetscValidPointer(list,3); 356 PetscValidHeaderSpecific(is,IS_CLASSID,4); 357 PetscValidIntPointer(found,5); 358 *found = -1; 359 for (i=0; i<n; i++) { 360 if (!list[i]) continue; 361 ierr = ISEqual(list[i],is,&flg);CHKERRQ(ierr); 362 if (flg) { 363 *found = i; 364 PetscFunctionReturn(0); 365 } 366 } 367 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set"); 368 PetscFunctionReturn(0); 369 } 370 371 #undef __FUNCT__ 372 #define __FUNCT__ "MatNestGetRow" 373 /* Get a block row as a new MatNest */ 374 static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B) 375 { 376 Mat_Nest *vs = (Mat_Nest*)A->data; 377 char keyname[256]; 378 PetscErrorCode ierr; 379 380 PetscFunctionBegin; 381 *B = NULL; 382 ierr = PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);CHKERRQ(ierr); 383 ierr = PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);CHKERRQ(ierr); 384 if (*B) PetscFunctionReturn(0); 385 386 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);CHKERRQ(ierr); 387 388 (*B)->assembled = A->assembled; 389 390 ierr = PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);CHKERRQ(ierr); 391 ierr = PetscObjectDereference((PetscObject)*B);CHKERRQ(ierr); /* Leave the only remaining reference in the composition */ 392 PetscFunctionReturn(0); 393 } 394 395 #undef __FUNCT__ 396 #define __FUNCT__ "MatNestFindSubMat" 397 static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B) 398 { 399 Mat_Nest *vs = (Mat_Nest*)A->data; 400 PetscErrorCode ierr; 401 PetscInt row,col; 402 PetscBool same,isFullCol,isFullColGlobal; 403 404 PetscFunctionBegin; 405 /* Check if full column space. This is a hack */ 406 isFullCol = PETSC_FALSE; 407 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);CHKERRQ(ierr); 408 if (same) { 409 PetscInt n,first,step,i,an,am,afirst,astep; 410 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 411 ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr); 412 isFullCol = PETSC_TRUE; 413 for (i=0,an=A->cmap->rstart; i<vs->nc; i++) { 414 ierr = ISStrideGetInfo(is->col[i],&afirst,&astep);CHKERRQ(ierr); 415 ierr = ISGetLocalSize(is->col[i],&am);CHKERRQ(ierr); 416 if (afirst != an || astep != step) isFullCol = PETSC_FALSE; 417 an += am; 418 } 419 if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE; 420 } 421 ierr = MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));CHKERRQ(ierr); 422 423 if (isFullColGlobal) { 424 PetscInt row; 425 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr); 426 ierr = MatNestGetRow(A,row,B);CHKERRQ(ierr); 427 } else { 428 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr); 429 ierr = MatNestFindIS(A,vs->nc,is->col,iscol,&col);CHKERRQ(ierr); 430 if (!vs->m[row][col]) { 431 PetscInt lr,lc; 432 433 ierr = MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);CHKERRQ(ierr); 434 ierr = ISGetLocalSize(vs->isglobal.row[row],&lr);CHKERRQ(ierr); 435 ierr = ISGetLocalSize(vs->isglobal.col[col],&lc);CHKERRQ(ierr); 436 ierr = MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 437 ierr = MatSetUp(vs->m[row][col]);CHKERRQ(ierr); 438 ierr = MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 439 ierr = MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 440 } 441 *B = vs->m[row][col]; 442 } 443 PetscFunctionReturn(0); 444 } 445 446 #undef __FUNCT__ 447 #define __FUNCT__ "MatGetSubMatrix_Nest" 448 static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B) 449 { 450 PetscErrorCode ierr; 451 Mat_Nest *vs = (Mat_Nest*)A->data; 452 Mat sub; 453 454 PetscFunctionBegin; 455 ierr = MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);CHKERRQ(ierr); 456 switch (reuse) { 457 case MAT_INITIAL_MATRIX: 458 if (sub) { ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr); } 459 *B = sub; 460 break; 461 case MAT_REUSE_MATRIX: 462 if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call"); 463 break; 464 case MAT_IGNORE_MATRIX: /* Nothing to do */ 465 break; 466 case MAT_INPLACE_MATRIX: /* Nothing to do */ 467 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet"); 468 break; 469 } 470 PetscFunctionReturn(0); 471 } 472 473 #undef __FUNCT__ 474 #define __FUNCT__ "MatGetLocalSubMatrix_Nest" 475 PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 476 { 477 PetscErrorCode ierr; 478 Mat_Nest *vs = (Mat_Nest*)A->data; 479 Mat sub; 480 481 PetscFunctionBegin; 482 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 483 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */ 484 if (sub) {ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr);} 485 *B = sub; 486 PetscFunctionReturn(0); 487 } 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "MatRestoreLocalSubMatrix_Nest" 491 static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 492 { 493 PetscErrorCode ierr; 494 Mat_Nest *vs = (Mat_Nest*)A->data; 495 Mat sub; 496 497 PetscFunctionBegin; 498 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 499 if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten"); 500 if (sub) { 501 if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times"); 502 ierr = MatDestroy(B);CHKERRQ(ierr); 503 } 504 PetscFunctionReturn(0); 505 } 506 507 #undef __FUNCT__ 508 #define __FUNCT__ "MatGetDiagonal_Nest" 509 static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v) 510 { 511 Mat_Nest *bA = (Mat_Nest*)A->data; 512 PetscInt i; 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 for (i=0; i<bA->nr; i++) { 517 Vec bv; 518 ierr = VecGetSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 519 if (bA->m[i][i]) { 520 ierr = MatGetDiagonal(bA->m[i][i],bv);CHKERRQ(ierr); 521 } else { 522 ierr = VecSet(bv,0.0);CHKERRQ(ierr); 523 } 524 ierr = VecRestoreSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 525 } 526 PetscFunctionReturn(0); 527 } 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "MatDiagonalScale_Nest" 531 static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r) 532 { 533 Mat_Nest *bA = (Mat_Nest*)A->data; 534 Vec bl,*br; 535 PetscInt i,j; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 ierr = PetscCalloc1(bA->nc,&br);CHKERRQ(ierr); 540 if (r) { 541 for (j=0; j<bA->nc; j++) {ierr = VecGetSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);} 542 } 543 bl = NULL; 544 for (i=0; i<bA->nr; i++) { 545 if (l) { 546 ierr = VecGetSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr); 547 } 548 for (j=0; j<bA->nc; j++) { 549 if (bA->m[i][j]) { 550 ierr = MatDiagonalScale(bA->m[i][j],bl,br[j]);CHKERRQ(ierr); 551 } 552 } 553 if (l) { 554 ierr = VecRestoreSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr); 555 } 556 } 557 if (r) { 558 for (j=0; j<bA->nc; j++) {ierr = VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);} 559 } 560 ierr = PetscFree(br);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563 564 #undef __FUNCT__ 565 #define __FUNCT__ "MatScale_Nest" 566 static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a) 567 { 568 Mat_Nest *bA = (Mat_Nest*)A->data; 569 PetscInt i,j; 570 PetscErrorCode ierr; 571 572 PetscFunctionBegin; 573 for (i=0; i<bA->nr; i++) { 574 for (j=0; j<bA->nc; j++) { 575 if (bA->m[i][j]) { 576 ierr = MatScale(bA->m[i][j],a);CHKERRQ(ierr); 577 } 578 } 579 } 580 PetscFunctionReturn(0); 581 } 582 583 #undef __FUNCT__ 584 #define __FUNCT__ "MatShift_Nest" 585 static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a) 586 { 587 Mat_Nest *bA = (Mat_Nest*)A->data; 588 PetscInt i; 589 PetscErrorCode ierr; 590 591 PetscFunctionBegin; 592 for (i=0; i<bA->nr; i++) { 593 if (!bA->m[i][i]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for shifting an empty diagonal block, insert a matrix in block (%D,%D)",i,i); 594 ierr = MatShift(bA->m[i][i],a);CHKERRQ(ierr); 595 } 596 PetscFunctionReturn(0); 597 } 598 599 #undef __FUNCT__ 600 #define __FUNCT__ "MatDiagonalSet_Nest" 601 static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is) 602 { 603 Mat_Nest *bA = (Mat_Nest*)A->data; 604 PetscInt i; 605 PetscErrorCode ierr; 606 607 PetscFunctionBegin; 608 for (i=0; i<bA->nr; i++) { 609 Vec bv; 610 ierr = VecGetSubVector(D,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 611 if (bA->m[i][i]) { 612 ierr = MatDiagonalSet(bA->m[i][i],bv,is);CHKERRQ(ierr); 613 } 614 ierr = VecRestoreSubVector(D,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 615 } 616 PetscFunctionReturn(0); 617 } 618 619 #undef __FUNCT__ 620 #define __FUNCT__ "MatSetRandom_Nest" 621 static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx) 622 { 623 Mat_Nest *bA = (Mat_Nest*)A->data; 624 PetscInt i,j; 625 PetscErrorCode ierr; 626 627 PetscFunctionBegin; 628 for (i=0; i<bA->nr; i++) { 629 for (j=0; j<bA->nc; j++) { 630 if (bA->m[i][j]) { 631 ierr = MatSetRandom(bA->m[i][j],rctx);CHKERRQ(ierr); 632 } 633 } 634 } 635 PetscFunctionReturn(0); 636 } 637 638 #undef __FUNCT__ 639 #define __FUNCT__ "MatCreateVecs_Nest" 640 static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left) 641 { 642 Mat_Nest *bA = (Mat_Nest*)A->data; 643 Vec *L,*R; 644 MPI_Comm comm; 645 PetscInt i,j; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 650 if (right) { 651 /* allocate R */ 652 ierr = PetscMalloc1(bA->nc, &R);CHKERRQ(ierr); 653 /* Create the right vectors */ 654 for (j=0; j<bA->nc; j++) { 655 for (i=0; i<bA->nr; i++) { 656 if (bA->m[i][j]) { 657 ierr = MatCreateVecs(bA->m[i][j],&R[j],NULL);CHKERRQ(ierr); 658 break; 659 } 660 } 661 if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 662 } 663 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr); 664 /* hand back control to the nest vector */ 665 for (j=0; j<bA->nc; j++) { 666 ierr = VecDestroy(&R[j]);CHKERRQ(ierr); 667 } 668 ierr = PetscFree(R);CHKERRQ(ierr); 669 } 670 671 if (left) { 672 /* allocate L */ 673 ierr = PetscMalloc1(bA->nr, &L);CHKERRQ(ierr); 674 /* Create the left vectors */ 675 for (i=0; i<bA->nr; i++) { 676 for (j=0; j<bA->nc; j++) { 677 if (bA->m[i][j]) { 678 ierr = MatCreateVecs(bA->m[i][j],NULL,&L[i]);CHKERRQ(ierr); 679 break; 680 } 681 } 682 if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 683 } 684 685 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr); 686 for (i=0; i<bA->nr; i++) { 687 ierr = VecDestroy(&L[i]);CHKERRQ(ierr); 688 } 689 690 ierr = PetscFree(L);CHKERRQ(ierr); 691 } 692 PetscFunctionReturn(0); 693 } 694 695 #undef __FUNCT__ 696 #define __FUNCT__ "MatView_Nest" 697 static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer) 698 { 699 Mat_Nest *bA = (Mat_Nest*)A->data; 700 PetscBool isascii; 701 PetscInt i,j; 702 PetscErrorCode ierr; 703 704 PetscFunctionBegin; 705 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 706 if (isascii) { 707 708 ierr = PetscViewerASCIIPrintf(viewer,"Matrix object: \n");CHKERRQ(ierr); 709 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 710 ierr = PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);CHKERRQ(ierr); 711 712 ierr = PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");CHKERRQ(ierr); 713 for (i=0; i<bA->nr; i++) { 714 for (j=0; j<bA->nc; j++) { 715 MatType type; 716 char name[256] = "",prefix[256] = ""; 717 PetscInt NR,NC; 718 PetscBool isNest = PETSC_FALSE; 719 720 if (!bA->m[i][j]) { 721 CHKERRQ(ierr);PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);CHKERRQ(ierr); 722 continue; 723 } 724 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr); 725 ierr = MatGetType(bA->m[i][j], &type);CHKERRQ(ierr); 726 if (((PetscObject)bA->m[i][j])->name) {ierr = PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);CHKERRQ(ierr);} 727 if (((PetscObject)bA->m[i][j])->prefix) {ierr = PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);CHKERRQ(ierr);} 728 ierr = PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr); 729 730 ierr = PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);CHKERRQ(ierr); 731 732 if (isNest) { 733 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); /* push1 */ 734 ierr = MatView(bA->m[i][j],viewer);CHKERRQ(ierr); 735 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop1 */ 736 } 737 } 738 } 739 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop0 */ 740 } 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "MatZeroEntries_Nest" 746 static PetscErrorCode MatZeroEntries_Nest(Mat A) 747 { 748 Mat_Nest *bA = (Mat_Nest*)A->data; 749 PetscInt i,j; 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 for (i=0; i<bA->nr; i++) { 754 for (j=0; j<bA->nc; j++) { 755 if (!bA->m[i][j]) continue; 756 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr); 757 } 758 } 759 PetscFunctionReturn(0); 760 } 761 762 #undef __FUNCT__ 763 #define __FUNCT__ "MatCopy_Nest" 764 static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str) 765 { 766 Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data; 767 PetscInt i,j,nr = bA->nr,nc = bA->nc; 768 PetscErrorCode ierr; 769 770 PetscFunctionBegin; 771 if (nr != bB->nr || nc != bB->nc) SETERRQ4(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot copy a Mat_Nest of block size (%D,%D) to a Mat_Nest of block size (%D,%D)",bB->nr,bB->nc,nr,nc); 772 for (i=0; i<nr; i++) { 773 for (j=0; j<nc; j++) { 774 if (bA->m[i][j] && bB->m[i][j]) { 775 ierr = MatCopy(bA->m[i][j],bB->m[i][j],str);CHKERRQ(ierr); 776 } else if (bA->m[i][j] || bB->m[i][j]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j); 777 } 778 } 779 PetscFunctionReturn(0); 780 } 781 782 #undef __FUNCT__ 783 #define __FUNCT__ "MatDuplicate_Nest" 784 static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 785 { 786 Mat_Nest *bA = (Mat_Nest*)A->data; 787 Mat *b; 788 PetscInt i,j,nr = bA->nr,nc = bA->nc; 789 PetscErrorCode ierr; 790 791 PetscFunctionBegin; 792 ierr = PetscMalloc1(nr*nc,&b);CHKERRQ(ierr); 793 for (i=0; i<nr; i++) { 794 for (j=0; j<nc; j++) { 795 if (bA->m[i][j]) { 796 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 797 } else { 798 b[i*nc+j] = NULL; 799 } 800 } 801 } 802 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 803 /* Give the new MatNest exclusive ownership */ 804 for (i=0; i<nr*nc; i++) { 805 ierr = MatDestroy(&b[i]);CHKERRQ(ierr); 806 } 807 ierr = PetscFree(b);CHKERRQ(ierr); 808 809 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 810 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 811 PetscFunctionReturn(0); 812 } 813 814 /* nest api */ 815 #undef __FUNCT__ 816 #define __FUNCT__ "MatNestGetSubMat_Nest" 817 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 818 { 819 Mat_Nest *bA = (Mat_Nest*)A->data; 820 821 PetscFunctionBegin; 822 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 823 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 824 *mat = bA->m[idxm][jdxm]; 825 PetscFunctionReturn(0); 826 } 827 828 #undef __FUNCT__ 829 #define __FUNCT__ "MatNestGetSubMat" 830 /*@ 831 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 832 833 Not collective 834 835 Input Parameters: 836 + A - nest matrix 837 . idxm - index of the matrix within the nest matrix 838 - jdxm - index of the matrix within the nest matrix 839 840 Output Parameter: 841 . sub - matrix at index idxm,jdxm within the nest matrix 842 843 Level: developer 844 845 .seealso: MatNestGetSize(), MatNestGetSubMats() 846 @*/ 847 PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 848 { 849 PetscErrorCode ierr; 850 851 PetscFunctionBegin; 852 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 853 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "MatNestSetSubMat_Nest" 858 PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat) 859 { 860 Mat_Nest *bA = (Mat_Nest*)A->data; 861 PetscInt m,n,M,N,mi,ni,Mi,Ni; 862 PetscErrorCode ierr; 863 864 PetscFunctionBegin; 865 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 866 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 867 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 868 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 869 ierr = ISGetLocalSize(bA->isglobal.row[idxm],&mi);CHKERRQ(ierr); 870 ierr = ISGetSize(bA->isglobal.row[idxm],&Mi);CHKERRQ(ierr); 871 ierr = ISGetLocalSize(bA->isglobal.col[jdxm],&ni);CHKERRQ(ierr); 872 ierr = ISGetSize(bA->isglobal.col[jdxm],&Ni);CHKERRQ(ierr); 873 if (M != Mi || N != Ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix dimension (%D,%D) incompatible with nest block (%D,%D)",M,N,Mi,Ni); 874 if (m != mi || n != ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix local dimension (%D,%D) incompatible with nest block (%D,%D)",m,n,mi,ni); 875 876 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 877 ierr = MatDestroy(&bA->m[idxm][jdxm]);CHKERRQ(ierr); 878 bA->m[idxm][jdxm] = mat; 879 PetscFunctionReturn(0); 880 } 881 882 #undef __FUNCT__ 883 #define __FUNCT__ "MatNestSetSubMat" 884 /*@ 885 MatNestSetSubMat - Set a single submatrix in the nest matrix. 886 887 Logically collective on the submatrix communicator 888 889 Input Parameters: 890 + A - nest matrix 891 . idxm - index of the matrix within the nest matrix 892 . jdxm - index of the matrix within the nest matrix 893 - sub - matrix at index idxm,jdxm within the nest matrix 894 895 Notes: 896 The new submatrix must have the same size and communicator as that block of the nest. 897 898 This increments the reference count of the submatrix. 899 900 Level: developer 901 902 .seealso: MatNestSetSubMats(), MatNestGetSubMat() 903 @*/ 904 PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub) 905 { 906 PetscErrorCode ierr; 907 908 PetscFunctionBegin; 909 ierr = PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));CHKERRQ(ierr); 910 PetscFunctionReturn(0); 911 } 912 913 #undef __FUNCT__ 914 #define __FUNCT__ "MatNestGetSubMats_Nest" 915 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 916 { 917 Mat_Nest *bA = (Mat_Nest*)A->data; 918 919 PetscFunctionBegin; 920 if (M) *M = bA->nr; 921 if (N) *N = bA->nc; 922 if (mat) *mat = bA->m; 923 PetscFunctionReturn(0); 924 } 925 926 #undef __FUNCT__ 927 #define __FUNCT__ "MatNestGetSubMats" 928 /*@C 929 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 930 931 Not collective 932 933 Input Parameters: 934 . A - nest matrix 935 936 Output Parameter: 937 + M - number of rows in the nest matrix 938 . N - number of cols in the nest matrix 939 - mat - 2d array of matrices 940 941 Notes: 942 943 The user should not free the array mat. 944 945 In Fortran, this routine has a calling sequence 946 $ call MatNestGetSubMats(A, M, N, mat, ierr) 947 where the space allocated for the optional argument mat is assumed large enough (if provided). 948 949 Level: developer 950 951 .seealso: MatNestGetSize(), MatNestGetSubMat() 952 @*/ 953 PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 954 { 955 PetscErrorCode ierr; 956 957 PetscFunctionBegin; 958 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 959 PetscFunctionReturn(0); 960 } 961 962 #undef __FUNCT__ 963 #define __FUNCT__ "MatNestGetSize_Nest" 964 PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 965 { 966 Mat_Nest *bA = (Mat_Nest*)A->data; 967 968 PetscFunctionBegin; 969 if (M) *M = bA->nr; 970 if (N) *N = bA->nc; 971 PetscFunctionReturn(0); 972 } 973 974 #undef __FUNCT__ 975 #define __FUNCT__ "MatNestGetSize" 976 /*@ 977 MatNestGetSize - Returns the size of the nest matrix. 978 979 Not collective 980 981 Input Parameters: 982 . A - nest matrix 983 984 Output Parameter: 985 + M - number of rows in the nested mat 986 - N - number of cols in the nested mat 987 988 Notes: 989 990 Level: developer 991 992 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 993 @*/ 994 PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 995 { 996 PetscErrorCode ierr; 997 998 PetscFunctionBegin; 999 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 1000 PetscFunctionReturn(0); 1001 } 1002 1003 #undef __FUNCT__ 1004 #define __FUNCT__ "MatNestGetISs_Nest" 1005 static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[]) 1006 { 1007 Mat_Nest *vs = (Mat_Nest*)A->data; 1008 PetscInt i; 1009 1010 PetscFunctionBegin; 1011 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i]; 1012 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i]; 1013 PetscFunctionReturn(0); 1014 } 1015 1016 #undef __FUNCT__ 1017 #define __FUNCT__ "MatNestGetISs" 1018 /*@C 1019 MatNestGetISs - Returns the index sets partitioning the row and column spaces 1020 1021 Not collective 1022 1023 Input Parameters: 1024 . A - nest matrix 1025 1026 Output Parameter: 1027 + rows - array of row index sets 1028 - cols - array of column index sets 1029 1030 Level: advanced 1031 1032 Notes: 1033 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 1034 1035 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs() 1036 @*/ 1037 PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[]) 1038 { 1039 PetscErrorCode ierr; 1040 1041 PetscFunctionBegin; 1042 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1043 ierr = PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 1044 PetscFunctionReturn(0); 1045 } 1046 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "MatNestGetLocalISs_Nest" 1049 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[]) 1050 { 1051 Mat_Nest *vs = (Mat_Nest*)A->data; 1052 PetscInt i; 1053 1054 PetscFunctionBegin; 1055 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i]; 1056 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i]; 1057 PetscFunctionReturn(0); 1058 } 1059 1060 #undef __FUNCT__ 1061 #define __FUNCT__ "MatNestGetLocalISs" 1062 /*@C 1063 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces 1064 1065 Not collective 1066 1067 Input Parameters: 1068 . A - nest matrix 1069 1070 Output Parameter: 1071 + rows - array of row index sets (or NULL to ignore) 1072 - cols - array of column index sets (or NULL to ignore) 1073 1074 Level: advanced 1075 1076 Notes: 1077 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 1078 1079 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs() 1080 @*/ 1081 PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[]) 1082 { 1083 PetscErrorCode ierr; 1084 1085 PetscFunctionBegin; 1086 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1087 ierr = PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 1088 PetscFunctionReturn(0); 1089 } 1090 1091 #undef __FUNCT__ 1092 #define __FUNCT__ "MatNestSetVecType_Nest" 1093 PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype) 1094 { 1095 PetscErrorCode ierr; 1096 PetscBool flg; 1097 1098 PetscFunctionBegin; 1099 ierr = PetscStrcmp(vtype,VECNEST,&flg);CHKERRQ(ierr); 1100 /* In reality, this only distinguishes VECNEST and "other" */ 1101 if (flg) A->ops->getvecs = MatCreateVecs_Nest; 1102 else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0; 1103 PetscFunctionReturn(0); 1104 } 1105 1106 #undef __FUNCT__ 1107 #define __FUNCT__ "MatNestSetVecType" 1108 /*@C 1109 MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs() 1110 1111 Not collective 1112 1113 Input Parameters: 1114 + A - nest matrix 1115 - vtype - type to use for creating vectors 1116 1117 Notes: 1118 1119 Level: developer 1120 1121 .seealso: MatCreateVecs() 1122 @*/ 1123 PetscErrorCode MatNestSetVecType(Mat A,VecType vtype) 1124 { 1125 PetscErrorCode ierr; 1126 1127 PetscFunctionBegin; 1128 ierr = PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 #undef __FUNCT__ 1133 #define __FUNCT__ "MatNestSetSubMats_Nest" 1134 PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1135 { 1136 Mat_Nest *s = (Mat_Nest*)A->data; 1137 PetscInt i,j,m,n,M,N; 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 s->nr = nr; 1142 s->nc = nc; 1143 1144 /* Create space for submatrices */ 1145 ierr = PetscMalloc1(nr,&s->m);CHKERRQ(ierr); 1146 for (i=0; i<nr; i++) { 1147 ierr = PetscMalloc1(nc,&s->m[i]);CHKERRQ(ierr); 1148 } 1149 for (i=0; i<nr; i++) { 1150 for (j=0; j<nc; j++) { 1151 s->m[i][j] = a[i*nc+j]; 1152 if (a[i*nc+j]) { 1153 ierr = PetscObjectReference((PetscObject)a[i*nc+j]);CHKERRQ(ierr); 1154 } 1155 } 1156 } 1157 1158 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1159 1160 ierr = PetscMalloc1(nr,&s->row_len);CHKERRQ(ierr); 1161 ierr = PetscMalloc1(nc,&s->col_len);CHKERRQ(ierr); 1162 for (i=0; i<nr; i++) s->row_len[i]=-1; 1163 for (j=0; j<nc; j++) s->col_len[j]=-1; 1164 1165 ierr = MatNestGetSizes_Private(A,&m,&n,&M,&N);CHKERRQ(ierr); 1166 1167 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1168 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1169 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1170 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1171 1172 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1173 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1174 1175 ierr = PetscCalloc2(nr,&s->left,nc,&s->right);CHKERRQ(ierr); 1176 PetscFunctionReturn(0); 1177 } 1178 1179 #undef __FUNCT__ 1180 #define __FUNCT__ "MatNestSetSubMats" 1181 /*@ 1182 MatNestSetSubMats - Sets the nested submatrices 1183 1184 Collective on Mat 1185 1186 Input Parameter: 1187 + N - nested matrix 1188 . nr - number of nested row blocks 1189 . is_row - index sets for each nested row block, or NULL to make contiguous 1190 . nc - number of nested column blocks 1191 . is_col - index sets for each nested column block, or NULL to make contiguous 1192 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1193 1194 Level: advanced 1195 1196 .seealso: MatCreateNest(), MATNEST 1197 @*/ 1198 PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1199 { 1200 PetscErrorCode ierr; 1201 PetscInt i; 1202 1203 PetscFunctionBegin; 1204 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1205 if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1206 if (nr && is_row) { 1207 PetscValidPointer(is_row,3); 1208 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1209 } 1210 if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1211 if (nc && is_col) { 1212 PetscValidPointer(is_col,5); 1213 for (i=0; i<nc; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1214 } 1215 if (nr*nc) PetscValidPointer(a,6); 1216 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr); 1217 PetscFunctionReturn(0); 1218 } 1219 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "MatNestCreateAggregateL2G_Private" 1222 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog) 1223 { 1224 PetscErrorCode ierr; 1225 PetscBool flg; 1226 PetscInt i,j,m,mi,*ix; 1227 1228 PetscFunctionBegin; 1229 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1230 if (islocal[i]) { 1231 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1232 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1233 } else { 1234 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1235 } 1236 m += mi; 1237 } 1238 if (flg) { 1239 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1240 for (i=0,n=0; i<n; i++) { 1241 ISLocalToGlobalMapping smap = NULL; 1242 VecScatter scat; 1243 IS isreq; 1244 Vec lvec,gvec; 1245 union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x; 1246 Mat sub; 1247 1248 if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers"); 1249 if (colflg) { 1250 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1251 } else { 1252 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1253 } 1254 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr);} 1255 if (islocal[i]) { 1256 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1257 } else { 1258 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1259 } 1260 for (j=0; j<mi; j++) ix[m+j] = j; 1261 if (smap) {ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);} 1262 /* 1263 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1264 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1265 The approach here is ugly because it uses VecScatter to move indices. 1266 */ 1267 ierr = VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);CHKERRQ(ierr); 1268 ierr = VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);CHKERRQ(ierr); 1269 ierr = ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);CHKERRQ(ierr); 1270 ierr = VecScatterCreate(gvec,isreq,lvec,NULL,&scat);CHKERRQ(ierr); 1271 ierr = VecGetArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1272 for (j=0; j<mi; j++) x[j].integer = ix[m+j]; 1273 ierr = VecRestoreArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1274 ierr = VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1275 ierr = VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1276 ierr = VecGetArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1277 for (j=0; j<mi; j++) ix[m+j] = x[j].integer; 1278 ierr = VecRestoreArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1279 ierr = VecDestroy(&lvec);CHKERRQ(ierr); 1280 ierr = VecDestroy(&gvec);CHKERRQ(ierr); 1281 ierr = ISDestroy(&isreq);CHKERRQ(ierr); 1282 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 1283 m += mi; 1284 } 1285 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1286 } else { 1287 *ltog = NULL; 1288 } 1289 PetscFunctionReturn(0); 1290 } 1291 1292 1293 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1294 /* 1295 nprocessors = NP 1296 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1297 proc 0: => (g_0,h_0,) 1298 proc 1: => (g_1,h_1,) 1299 ... 1300 proc nprocs-1: => (g_NP-1,h_NP-1,) 1301 1302 proc 0: proc 1: proc nprocs-1: 1303 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1304 1305 proc 0: 1306 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1307 proc 1: 1308 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1309 1310 proc NP-1: 1311 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1312 */ 1313 #undef __FUNCT__ 1314 #define __FUNCT__ "MatSetUp_NestIS_Private" 1315 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1316 { 1317 Mat_Nest *vs = (Mat_Nest*)A->data; 1318 PetscInt i,j,offset,n,nsum,bs; 1319 PetscErrorCode ierr; 1320 Mat sub = NULL; 1321 1322 PetscFunctionBegin; 1323 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1324 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1325 if (is_row) { /* valid IS is passed in */ 1326 /* refs on is[] are incremeneted */ 1327 for (i=0; i<vs->nr; i++) { 1328 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1329 1330 vs->isglobal.row[i] = is_row[i]; 1331 } 1332 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1333 nsum = 0; 1334 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1335 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1336 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1337 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1338 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1339 nsum += n; 1340 } 1341 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1342 offset -= nsum; 1343 for (i=0; i<vs->nr; i++) { 1344 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1345 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1346 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1347 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1348 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1349 offset += n; 1350 } 1351 } 1352 1353 if (is_col) { /* valid IS is passed in */ 1354 /* refs on is[] are incremeneted */ 1355 for (j=0; j<vs->nc; j++) { 1356 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1357 1358 vs->isglobal.col[j] = is_col[j]; 1359 } 1360 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1361 offset = A->cmap->rstart; 1362 nsum = 0; 1363 for (j=0; j<vs->nc; j++) { 1364 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1365 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1366 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1367 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1368 nsum += n; 1369 } 1370 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1371 offset -= nsum; 1372 for (j=0; j<vs->nc; j++) { 1373 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1374 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1375 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1376 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1377 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1378 offset += n; 1379 } 1380 } 1381 1382 /* Set up the local ISs */ 1383 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1384 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1385 for (i=0,offset=0; i<vs->nr; i++) { 1386 IS isloc; 1387 ISLocalToGlobalMapping rmap = NULL; 1388 PetscInt nlocal,bs; 1389 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1390 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1391 if (rmap) { 1392 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1393 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1394 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1395 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1396 } else { 1397 nlocal = 0; 1398 isloc = NULL; 1399 } 1400 vs->islocal.row[i] = isloc; 1401 offset += nlocal; 1402 } 1403 for (i=0,offset=0; i<vs->nc; i++) { 1404 IS isloc; 1405 ISLocalToGlobalMapping cmap = NULL; 1406 PetscInt nlocal,bs; 1407 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1408 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1409 if (cmap) { 1410 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1411 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1412 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1413 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1414 } else { 1415 nlocal = 0; 1416 isloc = NULL; 1417 } 1418 vs->islocal.col[i] = isloc; 1419 offset += nlocal; 1420 } 1421 1422 /* Set up the aggregate ISLocalToGlobalMapping */ 1423 { 1424 ISLocalToGlobalMapping rmap,cmap; 1425 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1426 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1427 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1428 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1429 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1430 } 1431 1432 #if defined(PETSC_USE_DEBUG) 1433 for (i=0; i<vs->nr; i++) { 1434 for (j=0; j<vs->nc; j++) { 1435 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1436 Mat B = vs->m[i][j]; 1437 if (!B) continue; 1438 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1439 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1440 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1441 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1442 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1443 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1444 if (M != Mi || N != Ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Global sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",M,N,i,j,Mi,Ni); 1445 if (m != mi || n != ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Local sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",m,n,i,j,mi,ni); 1446 } 1447 } 1448 #endif 1449 1450 /* Set A->assembled if all non-null blocks are currently assembled */ 1451 for (i=0; i<vs->nr; i++) { 1452 for (j=0; j<vs->nc; j++) { 1453 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1454 } 1455 } 1456 A->assembled = PETSC_TRUE; 1457 PetscFunctionReturn(0); 1458 } 1459 1460 #undef __FUNCT__ 1461 #define __FUNCT__ "MatCreateNest" 1462 /*@C 1463 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1464 1465 Collective on Mat 1466 1467 Input Parameter: 1468 + comm - Communicator for the new Mat 1469 . nr - number of nested row blocks 1470 . is_row - index sets for each nested row block, or NULL to make contiguous 1471 . nc - number of nested column blocks 1472 . is_col - index sets for each nested column block, or NULL to make contiguous 1473 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1474 1475 Output Parameter: 1476 . B - new matrix 1477 1478 Level: advanced 1479 1480 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST 1481 @*/ 1482 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1483 { 1484 Mat A; 1485 PetscErrorCode ierr; 1486 1487 PetscFunctionBegin; 1488 *B = 0; 1489 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1490 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1491 A->preallocated = PETSC_TRUE; 1492 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1493 *B = A; 1494 PetscFunctionReturn(0); 1495 } 1496 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "MatConvert_Nest_SeqAIJ_fast" 1499 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1500 { 1501 Mat_Nest *nest = (Mat_Nest*)A->data; 1502 PetscScalar **avv; 1503 PetscScalar *vv; 1504 PetscInt **aii,**ajj; 1505 PetscInt *ii,*jj,*ci; 1506 PetscInt nr,nc,nnz,i,j; 1507 PetscBool done; 1508 PetscErrorCode ierr; 1509 1510 PetscFunctionBegin; 1511 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 1512 if (reuse == MAT_REUSE_MATRIX) { 1513 PetscInt rnr; 1514 1515 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr); 1516 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ"); 1517 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1518 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr); 1519 } 1520 /* extract CSR for nested SeqAIJ matrices */ 1521 nnz = 0; 1522 ierr = PetscCalloc3(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv);CHKERRQ(ierr); 1523 for (i=0; i<nest->nr; ++i) { 1524 for (j=0; j<nest->nc; ++j) { 1525 Mat B = nest->m[i][j]; 1526 if (B) { 1527 PetscScalar *naa; 1528 PetscInt *nii,*njj,nnr; 1529 1530 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr); 1531 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ"); 1532 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr); 1533 nnz += nii[nnr]; 1534 1535 aii[i*nest->nc+j] = nii; 1536 ajj[i*nest->nc+j] = njj; 1537 avv[i*nest->nc+j] = naa; 1538 } 1539 } 1540 } 1541 if (reuse != MAT_REUSE_MATRIX) { 1542 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr); 1543 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr); 1544 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr); 1545 } else { 1546 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1547 } 1548 1549 /* new row pointer */ 1550 ierr = PetscMemzero(ii,(nr+1)*sizeof(PetscInt));CHKERRQ(ierr); 1551 for (i=0; i<nest->nr; ++i) { 1552 PetscInt ncr,rst; 1553 1554 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1555 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1556 for (j=0; j<nest->nc; ++j) { 1557 if (aii[i*nest->nc+j]) { 1558 PetscInt *nii = aii[i*nest->nc+j]; 1559 PetscInt ir; 1560 1561 for (ir=rst; ir<ncr+rst; ++ir) { 1562 ii[ir+1] += nii[1]-nii[0]; 1563 nii++; 1564 } 1565 } 1566 } 1567 } 1568 for (i=0; i<nr; i++) ii[i+1] += ii[i]; 1569 1570 /* construct CSR for the new matrix */ 1571 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr); 1572 for (i=0; i<nest->nr; ++i) { 1573 PetscInt ncr,rst; 1574 1575 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1576 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1577 for (j=0; j<nest->nc; ++j) { 1578 if (aii[i*nest->nc+j]) { 1579 PetscScalar *nvv = avv[i*nest->nc+j]; 1580 PetscInt *nii = aii[i*nest->nc+j]; 1581 PetscInt *njj = ajj[i*nest->nc+j]; 1582 PetscInt ir,cst; 1583 1584 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr); 1585 for (ir=rst; ir<ncr+rst; ++ir) { 1586 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir]; 1587 1588 for (ij=0;ij<rsize;ij++) { 1589 jj[ist+ij] = *njj+cst; 1590 vv[ist+ij] = *nvv; 1591 njj++; 1592 nvv++; 1593 } 1594 ci[ir] += rsize; 1595 nii++; 1596 } 1597 } 1598 } 1599 } 1600 ierr = PetscFree(ci);CHKERRQ(ierr); 1601 1602 /* restore info */ 1603 for (i=0; i<nest->nr; ++i) { 1604 for (j=0; j<nest->nc; ++j) { 1605 Mat B = nest->m[i][j]; 1606 if (B) { 1607 PetscInt nnr = 0, k = i*nest->nc+j; 1608 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr); 1609 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ"); 1610 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr); 1611 } 1612 } 1613 } 1614 ierr = PetscFree3(aii,ajj,avv);CHKERRQ(ierr); 1615 1616 /* finalize newmat */ 1617 if (reuse == MAT_INITIAL_MATRIX) { 1618 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr); 1619 } else if (reuse == MAT_INPLACE_MATRIX) { 1620 Mat B; 1621 1622 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr); 1623 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1624 } 1625 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1626 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1627 { 1628 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data); 1629 a->free_a = PETSC_TRUE; 1630 a->free_ij = PETSC_TRUE; 1631 } 1632 PetscFunctionReturn(0); 1633 } 1634 1635 #undef __FUNCT__ 1636 #define __FUNCT__ "MatConvert_Nest_AIJ" 1637 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1638 { 1639 PetscErrorCode ierr; 1640 Mat_Nest *nest = (Mat_Nest*)A->data; 1641 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1642 PetscInt cstart,cend; 1643 PetscMPIInt size; 1644 Mat C; 1645 1646 PetscFunctionBegin; 1647 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1648 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1649 PetscInt nf; 1650 PetscBool fast; 1651 1652 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1653 if (!fast) { 1654 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1655 } 1656 for (i=0; i<nest->nr && fast; ++i) { 1657 for (j=0; j<nest->nc && fast; ++j) { 1658 Mat B = nest->m[i][j]; 1659 if (B) { 1660 ierr = PetscObjectTypeCompare((PetscObject)B,MATAIJ,&fast);CHKERRQ(ierr); 1661 if (!fast) { 1662 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1663 } 1664 } 1665 } 1666 } 1667 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1668 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1669 if (fast) { 1670 PetscInt f,s; 1671 1672 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 1673 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1674 else { 1675 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 1676 nf += f; 1677 } 1678 } 1679 } 1680 for (i=0, nf=0; i<nest->nc && fast; ++i) { 1681 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1682 if (fast) { 1683 PetscInt f,s; 1684 1685 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 1686 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1687 else { 1688 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 1689 nf += f; 1690 } 1691 } 1692 } 1693 if (fast) { 1694 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 1695 PetscFunctionReturn(0); 1696 } 1697 } 1698 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1699 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1700 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1701 switch (reuse) { 1702 case MAT_INITIAL_MATRIX: 1703 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1704 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1705 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1706 *newmat = C; 1707 break; 1708 case MAT_REUSE_MATRIX: 1709 C = *newmat; 1710 break; 1711 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1712 } 1713 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1714 onnz = dnnz + m; 1715 for (k=0; k<m; k++) { 1716 dnnz[k] = 0; 1717 onnz[k] = 0; 1718 } 1719 for (j=0; j<nest->nc; ++j) { 1720 IS bNis; 1721 PetscInt bN; 1722 const PetscInt *bNindices; 1723 /* Using global column indices and ISAllGather() is not scalable. */ 1724 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1725 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1726 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1727 for (i=0; i<nest->nr; ++i) { 1728 PetscSF bmsf; 1729 PetscSFNode *iremote; 1730 Mat B; 1731 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1732 const PetscInt *bmindices; 1733 B = nest->m[i][j]; 1734 if (!B) continue; 1735 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1736 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1737 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1738 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1739 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1740 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1741 for (k = 0; k < bm; ++k){ 1742 sub_dnnz[k] = 0; 1743 sub_onnz[k] = 0; 1744 } 1745 /* 1746 Locate the owners for all of the locally-owned global row indices for this row block. 1747 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1748 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1749 */ 1750 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1751 for (br = 0; br < bm; ++br) { 1752 PetscInt row = bmindices[br], rowowner = 0, brncols, col; 1753 const PetscInt *brcols; 1754 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1755 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1756 /* how many roots */ 1757 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1758 /* get nonzero pattern */ 1759 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1760 for (k=0; k<brncols; k++) { 1761 col = bNindices[brcols[k]]; 1762 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 1763 sub_dnnz[br]++; 1764 } else { 1765 sub_onnz[br]++; 1766 } 1767 } 1768 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1769 } 1770 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1771 /* bsf will have to take care of disposing of bedges. */ 1772 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1773 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1774 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1775 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1776 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1777 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1778 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1779 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1780 } 1781 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1782 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1783 } 1784 /* Resize preallocation if overestimated */ 1785 for (i=0;i<m;i++) { 1786 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 1787 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 1788 } 1789 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1790 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1791 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1792 1793 /* Fill by row */ 1794 for (j=0; j<nest->nc; ++j) { 1795 /* Using global column indices and ISAllGather() is not scalable. */ 1796 IS bNis; 1797 PetscInt bN; 1798 const PetscInt *bNindices; 1799 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1800 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1801 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1802 for (i=0; i<nest->nr; ++i) { 1803 Mat B; 1804 PetscInt bm, br; 1805 const PetscInt *bmindices; 1806 B = nest->m[i][j]; 1807 if (!B) continue; 1808 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1809 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1810 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1811 for (br = 0; br < bm; ++br) { 1812 PetscInt row = bmindices[br], brncols, *cols; 1813 const PetscInt *brcols; 1814 const PetscScalar *brcoldata; 1815 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1816 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1817 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1818 /* 1819 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1820 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1821 */ 1822 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1823 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1824 ierr = PetscFree(cols);CHKERRQ(ierr); 1825 } 1826 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1827 } 1828 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1829 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1830 } 1831 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1832 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1833 PetscFunctionReturn(0); 1834 } 1835 1836 /*MC 1837 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1838 1839 Level: intermediate 1840 1841 Notes: 1842 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1843 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1844 It is usually used with DMComposite and DMCreateMatrix() 1845 1846 .seealso: MatCreate(), MatType, MatCreateNest() 1847 M*/ 1848 #undef __FUNCT__ 1849 #define __FUNCT__ "MatCreate_Nest" 1850 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1851 { 1852 Mat_Nest *s; 1853 PetscErrorCode ierr; 1854 1855 PetscFunctionBegin; 1856 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1857 A->data = (void*)s; 1858 1859 s->nr = -1; 1860 s->nc = -1; 1861 s->m = NULL; 1862 s->splitassembly = PETSC_FALSE; 1863 1864 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1865 1866 A->ops->mult = MatMult_Nest; 1867 A->ops->multadd = MatMultAdd_Nest; 1868 A->ops->multtranspose = MatMultTranspose_Nest; 1869 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1870 A->ops->transpose = MatTranspose_Nest; 1871 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1872 A->ops->assemblyend = MatAssemblyEnd_Nest; 1873 A->ops->zeroentries = MatZeroEntries_Nest; 1874 A->ops->copy = MatCopy_Nest; 1875 A->ops->duplicate = MatDuplicate_Nest; 1876 A->ops->getsubmatrix = MatGetSubMatrix_Nest; 1877 A->ops->destroy = MatDestroy_Nest; 1878 A->ops->view = MatView_Nest; 1879 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1880 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1881 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1882 A->ops->getdiagonal = MatGetDiagonal_Nest; 1883 A->ops->diagonalscale = MatDiagonalScale_Nest; 1884 A->ops->scale = MatScale_Nest; 1885 A->ops->shift = MatShift_Nest; 1886 A->ops->diagonalset = MatDiagonalSet_Nest; 1887 A->ops->setrandom = MatSetRandom_Nest; 1888 1889 A->spptr = 0; 1890 A->assembled = PETSC_FALSE; 1891 1892 /* expose Nest api's */ 1893 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1894 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1895 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1896 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1897 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1898 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1899 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1900 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1901 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1902 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 1903 1904 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1905 PetscFunctionReturn(0); 1906 } 1907