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