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