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; 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 if (nr*nc) PetscValidPointer(a,6); 1126 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 1130 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog) 1131 { 1132 PetscErrorCode ierr; 1133 PetscBool flg; 1134 PetscInt i,j,m,mi,*ix; 1135 1136 PetscFunctionBegin; 1137 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1138 if (islocal[i]) { 1139 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1140 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1141 } else { 1142 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1143 } 1144 m += mi; 1145 } 1146 if (flg) { 1147 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1148 for (i=0,n=0; i<n; i++) { 1149 ISLocalToGlobalMapping smap = NULL; 1150 VecScatter scat; 1151 IS isreq; 1152 Vec lvec,gvec; 1153 union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x; 1154 Mat sub; 1155 1156 if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers"); 1157 if (colflg) { 1158 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1159 } else { 1160 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1161 } 1162 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr);} 1163 if (islocal[i]) { 1164 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1165 } else { 1166 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1167 } 1168 for (j=0; j<mi; j++) ix[m+j] = j; 1169 if (smap) {ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);} 1170 /* 1171 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1172 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1173 The approach here is ugly because it uses VecScatter to move indices. 1174 */ 1175 ierr = VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);CHKERRQ(ierr); 1176 ierr = VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);CHKERRQ(ierr); 1177 ierr = ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);CHKERRQ(ierr); 1178 ierr = VecScatterCreate(gvec,isreq,lvec,NULL,&scat);CHKERRQ(ierr); 1179 ierr = VecGetArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1180 for (j=0; j<mi; j++) x[j].integer = ix[m+j]; 1181 ierr = VecRestoreArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1182 ierr = VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1183 ierr = VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1184 ierr = VecGetArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1185 for (j=0; j<mi; j++) ix[m+j] = x[j].integer; 1186 ierr = VecRestoreArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1187 ierr = VecDestroy(&lvec);CHKERRQ(ierr); 1188 ierr = VecDestroy(&gvec);CHKERRQ(ierr); 1189 ierr = ISDestroy(&isreq);CHKERRQ(ierr); 1190 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 1191 m += mi; 1192 } 1193 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1194 } else { 1195 *ltog = NULL; 1196 } 1197 PetscFunctionReturn(0); 1198 } 1199 1200 1201 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1202 /* 1203 nprocessors = NP 1204 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1205 proc 0: => (g_0,h_0,) 1206 proc 1: => (g_1,h_1,) 1207 ... 1208 proc nprocs-1: => (g_NP-1,h_NP-1,) 1209 1210 proc 0: proc 1: proc nprocs-1: 1211 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1212 1213 proc 0: 1214 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1215 proc 1: 1216 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1217 1218 proc NP-1: 1219 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1220 */ 1221 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1222 { 1223 Mat_Nest *vs = (Mat_Nest*)A->data; 1224 PetscInt i,j,offset,n,nsum,bs; 1225 PetscErrorCode ierr; 1226 Mat sub = NULL; 1227 1228 PetscFunctionBegin; 1229 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1230 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1231 if (is_row) { /* valid IS is passed in */ 1232 /* refs on is[] are incremeneted */ 1233 for (i=0; i<vs->nr; i++) { 1234 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1235 1236 vs->isglobal.row[i] = is_row[i]; 1237 } 1238 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1239 nsum = 0; 1240 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1241 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1242 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1243 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1244 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1245 nsum += n; 1246 } 1247 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1248 offset -= nsum; 1249 for (i=0; i<vs->nr; i++) { 1250 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1251 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1252 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1253 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1254 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1255 offset += n; 1256 } 1257 } 1258 1259 if (is_col) { /* valid IS is passed in */ 1260 /* refs on is[] are incremeneted */ 1261 for (j=0; j<vs->nc; j++) { 1262 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1263 1264 vs->isglobal.col[j] = is_col[j]; 1265 } 1266 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1267 offset = A->cmap->rstart; 1268 nsum = 0; 1269 for (j=0; j<vs->nc; j++) { 1270 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1271 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1272 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1273 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1274 nsum += n; 1275 } 1276 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1277 offset -= nsum; 1278 for (j=0; j<vs->nc; j++) { 1279 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1280 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1281 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1282 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1283 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1284 offset += n; 1285 } 1286 } 1287 1288 /* Set up the local ISs */ 1289 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1290 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1291 for (i=0,offset=0; i<vs->nr; i++) { 1292 IS isloc; 1293 ISLocalToGlobalMapping rmap = NULL; 1294 PetscInt nlocal,bs; 1295 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1296 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1297 if (rmap) { 1298 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1299 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1300 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1301 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1302 } else { 1303 nlocal = 0; 1304 isloc = NULL; 1305 } 1306 vs->islocal.row[i] = isloc; 1307 offset += nlocal; 1308 } 1309 for (i=0,offset=0; i<vs->nc; i++) { 1310 IS isloc; 1311 ISLocalToGlobalMapping cmap = NULL; 1312 PetscInt nlocal,bs; 1313 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1314 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1315 if (cmap) { 1316 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1317 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1318 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1319 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1320 } else { 1321 nlocal = 0; 1322 isloc = NULL; 1323 } 1324 vs->islocal.col[i] = isloc; 1325 offset += nlocal; 1326 } 1327 1328 /* Set up the aggregate ISLocalToGlobalMapping */ 1329 { 1330 ISLocalToGlobalMapping rmap,cmap; 1331 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1332 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1333 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1334 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1335 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1336 } 1337 1338 #if defined(PETSC_USE_DEBUG) 1339 for (i=0; i<vs->nr; i++) { 1340 for (j=0; j<vs->nc; j++) { 1341 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1342 Mat B = vs->m[i][j]; 1343 if (!B) continue; 1344 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1345 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1346 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1347 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1348 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1349 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1350 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); 1351 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); 1352 } 1353 } 1354 #endif 1355 1356 /* Set A->assembled if all non-null blocks are currently assembled */ 1357 for (i=0; i<vs->nr; i++) { 1358 for (j=0; j<vs->nc; j++) { 1359 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1360 } 1361 } 1362 A->assembled = PETSC_TRUE; 1363 PetscFunctionReturn(0); 1364 } 1365 1366 /*@C 1367 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1368 1369 Collective on Mat 1370 1371 Input Parameter: 1372 + comm - Communicator for the new Mat 1373 . nr - number of nested row blocks 1374 . is_row - index sets for each nested row block, or NULL to make contiguous 1375 . nc - number of nested column blocks 1376 . is_col - index sets for each nested column block, or NULL to make contiguous 1377 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1378 1379 Output Parameter: 1380 . B - new matrix 1381 1382 Level: advanced 1383 1384 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST 1385 @*/ 1386 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1387 { 1388 Mat A; 1389 PetscErrorCode ierr; 1390 1391 PetscFunctionBegin; 1392 *B = 0; 1393 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1394 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1395 A->preallocated = PETSC_TRUE; 1396 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1397 *B = A; 1398 PetscFunctionReturn(0); 1399 } 1400 1401 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1402 { 1403 Mat_Nest *nest = (Mat_Nest*)A->data; 1404 Mat *trans; 1405 PetscScalar **avv; 1406 PetscScalar *vv; 1407 PetscInt **aii,**ajj; 1408 PetscInt *ii,*jj,*ci; 1409 PetscInt nr,nc,nnz,i,j; 1410 PetscBool done; 1411 PetscErrorCode ierr; 1412 1413 PetscFunctionBegin; 1414 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 1415 if (reuse == MAT_REUSE_MATRIX) { 1416 PetscInt rnr; 1417 1418 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr); 1419 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ"); 1420 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1421 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr); 1422 } 1423 /* extract CSR for nested SeqAIJ matrices */ 1424 nnz = 0; 1425 ierr = PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);CHKERRQ(ierr); 1426 for (i=0; i<nest->nr; ++i) { 1427 for (j=0; j<nest->nc; ++j) { 1428 Mat B = nest->m[i][j]; 1429 if (B) { 1430 PetscScalar *naa; 1431 PetscInt *nii,*njj,nnr; 1432 PetscBool istrans; 1433 1434 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1435 if (istrans) { 1436 Mat Bt; 1437 1438 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1439 ierr = MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);CHKERRQ(ierr); 1440 B = trans[i*nest->nc+j]; 1441 } 1442 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr); 1443 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ"); 1444 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr); 1445 nnz += nii[nnr]; 1446 1447 aii[i*nest->nc+j] = nii; 1448 ajj[i*nest->nc+j] = njj; 1449 avv[i*nest->nc+j] = naa; 1450 } 1451 } 1452 } 1453 if (reuse != MAT_REUSE_MATRIX) { 1454 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr); 1455 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr); 1456 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr); 1457 } else { 1458 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1459 } 1460 1461 /* new row pointer */ 1462 ierr = PetscMemzero(ii,(nr+1)*sizeof(PetscInt));CHKERRQ(ierr); 1463 for (i=0; i<nest->nr; ++i) { 1464 PetscInt ncr,rst; 1465 1466 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1467 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1468 for (j=0; j<nest->nc; ++j) { 1469 if (aii[i*nest->nc+j]) { 1470 PetscInt *nii = aii[i*nest->nc+j]; 1471 PetscInt ir; 1472 1473 for (ir=rst; ir<ncr+rst; ++ir) { 1474 ii[ir+1] += nii[1]-nii[0]; 1475 nii++; 1476 } 1477 } 1478 } 1479 } 1480 for (i=0; i<nr; i++) ii[i+1] += ii[i]; 1481 1482 /* construct CSR for the new matrix */ 1483 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr); 1484 for (i=0; i<nest->nr; ++i) { 1485 PetscInt ncr,rst; 1486 1487 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1488 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1489 for (j=0; j<nest->nc; ++j) { 1490 if (aii[i*nest->nc+j]) { 1491 PetscScalar *nvv = avv[i*nest->nc+j]; 1492 PetscInt *nii = aii[i*nest->nc+j]; 1493 PetscInt *njj = ajj[i*nest->nc+j]; 1494 PetscInt ir,cst; 1495 1496 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr); 1497 for (ir=rst; ir<ncr+rst; ++ir) { 1498 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir]; 1499 1500 for (ij=0;ij<rsize;ij++) { 1501 jj[ist+ij] = *njj+cst; 1502 vv[ist+ij] = *nvv; 1503 njj++; 1504 nvv++; 1505 } 1506 ci[ir] += rsize; 1507 nii++; 1508 } 1509 } 1510 } 1511 } 1512 ierr = PetscFree(ci);CHKERRQ(ierr); 1513 1514 /* restore info */ 1515 for (i=0; i<nest->nr; ++i) { 1516 for (j=0; j<nest->nc; ++j) { 1517 Mat B = nest->m[i][j]; 1518 if (B) { 1519 PetscInt nnr = 0, k = i*nest->nc+j; 1520 1521 B = (trans[k] ? trans[k] : B); 1522 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr); 1523 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ"); 1524 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr); 1525 ierr = MatDestroy(&trans[k]);CHKERRQ(ierr); 1526 } 1527 } 1528 } 1529 ierr = PetscFree4(aii,ajj,avv,trans);CHKERRQ(ierr); 1530 1531 /* finalize newmat */ 1532 if (reuse == MAT_INITIAL_MATRIX) { 1533 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr); 1534 } else if (reuse == MAT_INPLACE_MATRIX) { 1535 Mat B; 1536 1537 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr); 1538 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1539 } 1540 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1541 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1542 { 1543 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data); 1544 a->free_a = PETSC_TRUE; 1545 a->free_ij = PETSC_TRUE; 1546 } 1547 PetscFunctionReturn(0); 1548 } 1549 1550 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1551 { 1552 PetscErrorCode ierr; 1553 Mat_Nest *nest = (Mat_Nest*)A->data; 1554 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1555 PetscInt cstart,cend; 1556 PetscMPIInt size; 1557 Mat C; 1558 1559 PetscFunctionBegin; 1560 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1561 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1562 PetscInt nf; 1563 PetscBool fast; 1564 1565 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1566 if (!fast) { 1567 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1568 } 1569 for (i=0; i<nest->nr && fast; ++i) { 1570 for (j=0; j<nest->nc && fast; ++j) { 1571 Mat B = nest->m[i][j]; 1572 if (B) { 1573 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1574 if (!fast) { 1575 PetscBool istrans; 1576 1577 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1578 if (istrans) { 1579 Mat Bt; 1580 1581 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1582 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr); 1583 } 1584 } 1585 } 1586 } 1587 } 1588 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1589 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1590 if (fast) { 1591 PetscInt f,s; 1592 1593 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 1594 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1595 else { 1596 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 1597 nf += f; 1598 } 1599 } 1600 } 1601 for (i=0, nf=0; i<nest->nc && fast; ++i) { 1602 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1603 if (fast) { 1604 PetscInt f,s; 1605 1606 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 1607 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1608 else { 1609 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 1610 nf += f; 1611 } 1612 } 1613 } 1614 if (fast) { 1615 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 1616 PetscFunctionReturn(0); 1617 } 1618 } 1619 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1620 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1621 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1622 switch (reuse) { 1623 case MAT_INITIAL_MATRIX: 1624 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1625 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1626 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1627 *newmat = C; 1628 break; 1629 case MAT_REUSE_MATRIX: 1630 C = *newmat; 1631 break; 1632 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1633 } 1634 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1635 onnz = dnnz + m; 1636 for (k=0; k<m; k++) { 1637 dnnz[k] = 0; 1638 onnz[k] = 0; 1639 } 1640 for (j=0; j<nest->nc; ++j) { 1641 IS bNis; 1642 PetscInt bN; 1643 const PetscInt *bNindices; 1644 /* Using global column indices and ISAllGather() is not scalable. */ 1645 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1646 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1647 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1648 for (i=0; i<nest->nr; ++i) { 1649 PetscSF bmsf; 1650 PetscSFNode *iremote; 1651 Mat B; 1652 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1653 const PetscInt *bmindices; 1654 B = nest->m[i][j]; 1655 if (!B) continue; 1656 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1657 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1658 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1659 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1660 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1661 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1662 for (k = 0; k < bm; ++k){ 1663 sub_dnnz[k] = 0; 1664 sub_onnz[k] = 0; 1665 } 1666 /* 1667 Locate the owners for all of the locally-owned global row indices for this row block. 1668 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1669 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1670 */ 1671 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1672 for (br = 0; br < bm; ++br) { 1673 PetscInt row = bmindices[br], rowowner = 0, brncols, col; 1674 const PetscInt *brcols; 1675 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1676 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1677 /* how many roots */ 1678 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1679 /* get nonzero pattern */ 1680 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1681 for (k=0; k<brncols; k++) { 1682 col = bNindices[brcols[k]]; 1683 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 1684 sub_dnnz[br]++; 1685 } else { 1686 sub_onnz[br]++; 1687 } 1688 } 1689 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1690 } 1691 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1692 /* bsf will have to take care of disposing of bedges. */ 1693 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1694 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1695 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1696 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1697 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1698 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1699 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1700 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1701 } 1702 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1703 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1704 } 1705 /* Resize preallocation if overestimated */ 1706 for (i=0;i<m;i++) { 1707 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 1708 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 1709 } 1710 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1711 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1712 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1713 1714 /* Fill by row */ 1715 for (j=0; j<nest->nc; ++j) { 1716 /* Using global column indices and ISAllGather() is not scalable. */ 1717 IS bNis; 1718 PetscInt bN; 1719 const PetscInt *bNindices; 1720 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1721 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1722 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1723 for (i=0; i<nest->nr; ++i) { 1724 Mat B; 1725 PetscInt bm, br; 1726 const PetscInt *bmindices; 1727 B = nest->m[i][j]; 1728 if (!B) continue; 1729 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1730 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1731 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1732 for (br = 0; br < bm; ++br) { 1733 PetscInt row = bmindices[br], brncols, *cols; 1734 const PetscInt *brcols; 1735 const PetscScalar *brcoldata; 1736 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1737 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1738 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1739 /* 1740 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1741 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1742 */ 1743 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1744 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1745 ierr = PetscFree(cols);CHKERRQ(ierr); 1746 } 1747 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1748 } 1749 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1750 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1751 } 1752 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1753 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1754 PetscFunctionReturn(0); 1755 } 1756 1757 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has) 1758 { 1759 PetscFunctionBegin; 1760 *has = PETSC_FALSE; 1761 if (op == MATOP_MULT_TRANSPOSE) { 1762 Mat_Nest *bA = (Mat_Nest*)mat->data; 1763 PetscInt i,j,nr = bA->nr,nc = bA->nc; 1764 PetscErrorCode ierr; 1765 PetscBool flg; 1766 1767 for (j=0; j<nc; j++) { 1768 for (i=0; i<nr; i++) { 1769 if (!bA->m[i][j]) continue; 1770 ierr = MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);CHKERRQ(ierr); 1771 if (!flg) PetscFunctionReturn(0); 1772 } 1773 } 1774 } 1775 if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 1776 PetscFunctionReturn(0); 1777 } 1778 1779 /*MC 1780 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1781 1782 Level: intermediate 1783 1784 Notes: 1785 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1786 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1787 It is usually used with DMComposite and DMCreateMatrix() 1788 1789 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 1790 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 1791 than the nest matrix. 1792 1793 .seealso: MatCreate(), MatType, MatCreateNest() 1794 M*/ 1795 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1796 { 1797 Mat_Nest *s; 1798 PetscErrorCode ierr; 1799 1800 PetscFunctionBegin; 1801 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1802 A->data = (void*)s; 1803 1804 s->nr = -1; 1805 s->nc = -1; 1806 s->m = NULL; 1807 s->splitassembly = PETSC_FALSE; 1808 1809 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1810 1811 A->ops->mult = MatMult_Nest; 1812 A->ops->multadd = MatMultAdd_Nest; 1813 A->ops->multtranspose = MatMultTranspose_Nest; 1814 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1815 A->ops->transpose = MatTranspose_Nest; 1816 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1817 A->ops->assemblyend = MatAssemblyEnd_Nest; 1818 A->ops->zeroentries = MatZeroEntries_Nest; 1819 A->ops->copy = MatCopy_Nest; 1820 A->ops->duplicate = MatDuplicate_Nest; 1821 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 1822 A->ops->destroy = MatDestroy_Nest; 1823 A->ops->view = MatView_Nest; 1824 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1825 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1826 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1827 A->ops->getdiagonal = MatGetDiagonal_Nest; 1828 A->ops->diagonalscale = MatDiagonalScale_Nest; 1829 A->ops->scale = MatScale_Nest; 1830 A->ops->shift = MatShift_Nest; 1831 A->ops->diagonalset = MatDiagonalSet_Nest; 1832 A->ops->setrandom = MatSetRandom_Nest; 1833 A->ops->hasoperation = MatHasOperation_Nest; 1834 1835 A->spptr = 0; 1836 A->assembled = PETSC_FALSE; 1837 1838 /* expose Nest api's */ 1839 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1840 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1841 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1842 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1843 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1844 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1845 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1846 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1847 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1848 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 1849 1850 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1851 PetscFunctionReturn(0); 1852 } 1853