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_REUSE_MATRIX && A == *B && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place"); 148 149 if (reuse == MAT_INITIAL_MATRIX || A == *B) { 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_REUSE_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 || A != *B) { 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) { 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 MatGetSubMatrix_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 #undef __FUNCT__ 1402 #define __FUNCT__ "MatConvert_Nest_SeqAIJ_fast" 1403 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1404 { 1405 Mat_Nest *nest = (Mat_Nest*)A->data; 1406 Mat *trans; 1407 PetscScalar **avv; 1408 PetscScalar *vv; 1409 PetscInt **aii,**ajj; 1410 PetscInt *ii,*jj,*ci; 1411 PetscInt nr,nc,nnz,i,j; 1412 PetscBool done; 1413 PetscErrorCode ierr; 1414 1415 PetscFunctionBegin; 1416 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 1417 if (reuse == MAT_REUSE_MATRIX) { 1418 PetscInt rnr; 1419 1420 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr); 1421 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ"); 1422 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1423 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr); 1424 } 1425 /* extract CSR for nested SeqAIJ matrices */ 1426 nnz = 0; 1427 ierr = PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);CHKERRQ(ierr); 1428 for (i=0; i<nest->nr; ++i) { 1429 for (j=0; j<nest->nc; ++j) { 1430 Mat B = nest->m[i][j]; 1431 if (B) { 1432 PetscScalar *naa; 1433 PetscInt *nii,*njj,nnr; 1434 PetscBool istrans; 1435 1436 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1437 if (istrans) { 1438 Mat Bt; 1439 1440 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1441 ierr = MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);CHKERRQ(ierr); 1442 B = trans[i*nest->nc+j]; 1443 } 1444 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr); 1445 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ"); 1446 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr); 1447 nnz += nii[nnr]; 1448 1449 aii[i*nest->nc+j] = nii; 1450 ajj[i*nest->nc+j] = njj; 1451 avv[i*nest->nc+j] = naa; 1452 } 1453 } 1454 } 1455 if (reuse != MAT_REUSE_MATRIX) { 1456 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr); 1457 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr); 1458 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr); 1459 } else { 1460 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1461 } 1462 1463 /* new row pointer */ 1464 ierr = PetscMemzero(ii,(nr+1)*sizeof(PetscInt));CHKERRQ(ierr); 1465 for (i=0; i<nest->nr; ++i) { 1466 PetscInt ncr,rst; 1467 1468 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1469 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1470 for (j=0; j<nest->nc; ++j) { 1471 if (aii[i*nest->nc+j]) { 1472 PetscInt *nii = aii[i*nest->nc+j]; 1473 PetscInt ir; 1474 1475 for (ir=rst; ir<ncr+rst; ++ir) { 1476 ii[ir+1] += nii[1]-nii[0]; 1477 nii++; 1478 } 1479 } 1480 } 1481 } 1482 for (i=0; i<nr; i++) ii[i+1] += ii[i]; 1483 1484 /* construct CSR for the new matrix */ 1485 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr); 1486 for (i=0; i<nest->nr; ++i) { 1487 PetscInt ncr,rst; 1488 1489 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1490 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1491 for (j=0; j<nest->nc; ++j) { 1492 if (aii[i*nest->nc+j]) { 1493 PetscScalar *nvv = avv[i*nest->nc+j]; 1494 PetscInt *nii = aii[i*nest->nc+j]; 1495 PetscInt *njj = ajj[i*nest->nc+j]; 1496 PetscInt ir,cst; 1497 1498 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr); 1499 for (ir=rst; ir<ncr+rst; ++ir) { 1500 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir]; 1501 1502 for (ij=0;ij<rsize;ij++) { 1503 jj[ist+ij] = *njj+cst; 1504 vv[ist+ij] = *nvv; 1505 njj++; 1506 nvv++; 1507 } 1508 ci[ir] += rsize; 1509 nii++; 1510 } 1511 } 1512 } 1513 } 1514 ierr = PetscFree(ci);CHKERRQ(ierr); 1515 1516 /* restore info */ 1517 for (i=0; i<nest->nr; ++i) { 1518 for (j=0; j<nest->nc; ++j) { 1519 Mat B = nest->m[i][j]; 1520 if (B) { 1521 PetscInt nnr = 0, k = i*nest->nc+j; 1522 1523 B = (trans[k] ? trans[k] : B); 1524 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr); 1525 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ"); 1526 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr); 1527 ierr = MatDestroy(&trans[k]);CHKERRQ(ierr); 1528 } 1529 } 1530 } 1531 ierr = PetscFree4(aii,ajj,avv,trans);CHKERRQ(ierr); 1532 1533 /* finalize newmat */ 1534 if (reuse == MAT_INITIAL_MATRIX) { 1535 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr); 1536 } else if (reuse == MAT_INPLACE_MATRIX) { 1537 Mat B; 1538 1539 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr); 1540 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1541 } 1542 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1543 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1544 { 1545 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data); 1546 a->free_a = PETSC_TRUE; 1547 a->free_ij = PETSC_TRUE; 1548 } 1549 PetscFunctionReturn(0); 1550 } 1551 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "MatConvert_Nest_AIJ" 1554 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1555 { 1556 PetscErrorCode ierr; 1557 Mat_Nest *nest = (Mat_Nest*)A->data; 1558 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1559 PetscInt cstart,cend; 1560 PetscMPIInt size; 1561 Mat C; 1562 1563 PetscFunctionBegin; 1564 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1565 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1566 PetscInt nf; 1567 PetscBool fast; 1568 1569 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1570 if (!fast) { 1571 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1572 } 1573 for (i=0; i<nest->nr && fast; ++i) { 1574 for (j=0; j<nest->nc && fast; ++j) { 1575 Mat B = nest->m[i][j]; 1576 if (B) { 1577 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1578 if (!fast) { 1579 PetscBool istrans; 1580 1581 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1582 if (istrans) { 1583 Mat Bt; 1584 1585 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1586 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr); 1587 } 1588 } 1589 } 1590 } 1591 } 1592 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1593 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1594 if (fast) { 1595 PetscInt f,s; 1596 1597 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 1598 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1599 else { 1600 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 1601 nf += f; 1602 } 1603 } 1604 } 1605 for (i=0, nf=0; i<nest->nc && fast; ++i) { 1606 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1607 if (fast) { 1608 PetscInt f,s; 1609 1610 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 1611 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1612 else { 1613 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 1614 nf += f; 1615 } 1616 } 1617 } 1618 if (fast) { 1619 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 1620 PetscFunctionReturn(0); 1621 } 1622 } 1623 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1624 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1625 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1626 switch (reuse) { 1627 case MAT_INITIAL_MATRIX: 1628 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1629 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1630 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1631 *newmat = C; 1632 break; 1633 case MAT_REUSE_MATRIX: 1634 C = *newmat; 1635 break; 1636 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1637 } 1638 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1639 onnz = dnnz + m; 1640 for (k=0; k<m; k++) { 1641 dnnz[k] = 0; 1642 onnz[k] = 0; 1643 } 1644 for (j=0; j<nest->nc; ++j) { 1645 IS bNis; 1646 PetscInt bN; 1647 const PetscInt *bNindices; 1648 /* Using global column indices and ISAllGather() is not scalable. */ 1649 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1650 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1651 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1652 for (i=0; i<nest->nr; ++i) { 1653 PetscSF bmsf; 1654 PetscSFNode *iremote; 1655 Mat B; 1656 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1657 const PetscInt *bmindices; 1658 B = nest->m[i][j]; 1659 if (!B) continue; 1660 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1661 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1662 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1663 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1664 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1665 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1666 for (k = 0; k < bm; ++k){ 1667 sub_dnnz[k] = 0; 1668 sub_onnz[k] = 0; 1669 } 1670 /* 1671 Locate the owners for all of the locally-owned global row indices for this row block. 1672 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1673 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1674 */ 1675 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1676 for (br = 0; br < bm; ++br) { 1677 PetscInt row = bmindices[br], rowowner = 0, brncols, col; 1678 const PetscInt *brcols; 1679 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1680 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1681 /* how many roots */ 1682 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1683 /* get nonzero pattern */ 1684 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1685 for (k=0; k<brncols; k++) { 1686 col = bNindices[brcols[k]]; 1687 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 1688 sub_dnnz[br]++; 1689 } else { 1690 sub_onnz[br]++; 1691 } 1692 } 1693 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1694 } 1695 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1696 /* bsf will have to take care of disposing of bedges. */ 1697 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1698 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1699 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1700 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1701 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1702 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1703 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1704 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1705 } 1706 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1707 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1708 } 1709 /* Resize preallocation if overestimated */ 1710 for (i=0;i<m;i++) { 1711 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 1712 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 1713 } 1714 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1715 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1716 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1717 1718 /* Fill by row */ 1719 for (j=0; j<nest->nc; ++j) { 1720 /* Using global column indices and ISAllGather() is not scalable. */ 1721 IS bNis; 1722 PetscInt bN; 1723 const PetscInt *bNindices; 1724 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1725 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1726 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1727 for (i=0; i<nest->nr; ++i) { 1728 Mat B; 1729 PetscInt bm, br; 1730 const PetscInt *bmindices; 1731 B = nest->m[i][j]; 1732 if (!B) continue; 1733 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1734 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1735 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1736 for (br = 0; br < bm; ++br) { 1737 PetscInt row = bmindices[br], brncols, *cols; 1738 const PetscInt *brcols; 1739 const PetscScalar *brcoldata; 1740 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1741 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1742 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1743 /* 1744 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1745 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1746 */ 1747 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1748 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1749 ierr = PetscFree(cols);CHKERRQ(ierr); 1750 } 1751 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1752 } 1753 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1754 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1755 } 1756 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1757 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1758 PetscFunctionReturn(0); 1759 } 1760 1761 /*MC 1762 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1763 1764 Level: intermediate 1765 1766 Notes: 1767 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1768 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1769 It is usually used with DMComposite and DMCreateMatrix() 1770 1771 .seealso: MatCreate(), MatType, MatCreateNest() 1772 M*/ 1773 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1774 { 1775 Mat_Nest *s; 1776 PetscErrorCode ierr; 1777 1778 PetscFunctionBegin; 1779 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1780 A->data = (void*)s; 1781 1782 s->nr = -1; 1783 s->nc = -1; 1784 s->m = NULL; 1785 s->splitassembly = PETSC_FALSE; 1786 1787 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1788 1789 A->ops->mult = MatMult_Nest; 1790 A->ops->multadd = MatMultAdd_Nest; 1791 A->ops->multtranspose = MatMultTranspose_Nest; 1792 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1793 A->ops->transpose = MatTranspose_Nest; 1794 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1795 A->ops->assemblyend = MatAssemblyEnd_Nest; 1796 A->ops->zeroentries = MatZeroEntries_Nest; 1797 A->ops->copy = MatCopy_Nest; 1798 A->ops->duplicate = MatDuplicate_Nest; 1799 A->ops->getsubmatrix = MatGetSubMatrix_Nest; 1800 A->ops->destroy = MatDestroy_Nest; 1801 A->ops->view = MatView_Nest; 1802 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1803 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1804 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1805 A->ops->getdiagonal = MatGetDiagonal_Nest; 1806 A->ops->diagonalscale = MatDiagonalScale_Nest; 1807 A->ops->scale = MatScale_Nest; 1808 A->ops->shift = MatShift_Nest; 1809 A->ops->diagonalset = MatDiagonalSet_Nest; 1810 A->ops->setrandom = MatSetRandom_Nest; 1811 1812 A->spptr = 0; 1813 A->assembled = PETSC_FALSE; 1814 1815 /* expose Nest api's */ 1816 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1817 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1818 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1819 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1820 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1821 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1822 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1823 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1824 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1825 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 1826 1827 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1828 PetscFunctionReturn(0); 1829 } 1830