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 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 724 PetscFunctionReturn(0); 725 } 726 727 static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 728 { 729 Mat_Nest *bA = (Mat_Nest*)A->data; 730 Mat *b; 731 PetscInt i,j,nr = bA->nr,nc = bA->nc; 732 PetscErrorCode ierr; 733 734 PetscFunctionBegin; 735 ierr = PetscMalloc1(nr*nc,&b);CHKERRQ(ierr); 736 for (i=0; i<nr; i++) { 737 for (j=0; j<nc; j++) { 738 if (bA->m[i][j]) { 739 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 740 } else { 741 b[i*nc+j] = NULL; 742 } 743 } 744 } 745 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 746 /* Give the new MatNest exclusive ownership */ 747 for (i=0; i<nr*nc; i++) { 748 ierr = MatDestroy(&b[i]);CHKERRQ(ierr); 749 } 750 ierr = PetscFree(b);CHKERRQ(ierr); 751 752 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 753 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 754 PetscFunctionReturn(0); 755 } 756 757 /* nest api */ 758 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 759 { 760 Mat_Nest *bA = (Mat_Nest*)A->data; 761 762 PetscFunctionBegin; 763 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 764 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 765 *mat = bA->m[idxm][jdxm]; 766 PetscFunctionReturn(0); 767 } 768 769 /*@ 770 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 771 772 Not collective 773 774 Input Parameters: 775 + A - nest matrix 776 . idxm - index of the matrix within the nest matrix 777 - jdxm - index of the matrix within the nest matrix 778 779 Output Parameter: 780 . sub - matrix at index idxm,jdxm within the nest matrix 781 782 Level: developer 783 784 .seealso: MatNestGetSize(), MatNestGetSubMats() 785 @*/ 786 PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 787 { 788 PetscErrorCode ierr; 789 790 PetscFunctionBegin; 791 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 792 PetscFunctionReturn(0); 793 } 794 795 PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat) 796 { 797 Mat_Nest *bA = (Mat_Nest*)A->data; 798 PetscInt m,n,M,N,mi,ni,Mi,Ni; 799 PetscErrorCode ierr; 800 801 PetscFunctionBegin; 802 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 803 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 804 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 805 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 806 ierr = ISGetLocalSize(bA->isglobal.row[idxm],&mi);CHKERRQ(ierr); 807 ierr = ISGetSize(bA->isglobal.row[idxm],&Mi);CHKERRQ(ierr); 808 ierr = ISGetLocalSize(bA->isglobal.col[jdxm],&ni);CHKERRQ(ierr); 809 ierr = ISGetSize(bA->isglobal.col[jdxm],&Ni);CHKERRQ(ierr); 810 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); 811 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); 812 813 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 814 ierr = MatDestroy(&bA->m[idxm][jdxm]);CHKERRQ(ierr); 815 bA->m[idxm][jdxm] = mat; 816 PetscFunctionReturn(0); 817 } 818 819 /*@ 820 MatNestSetSubMat - Set a single submatrix in the nest matrix. 821 822 Logically collective on the submatrix communicator 823 824 Input Parameters: 825 + A - nest matrix 826 . idxm - index of the matrix within the nest matrix 827 . jdxm - index of the matrix within the nest matrix 828 - sub - matrix at index idxm,jdxm within the nest matrix 829 830 Notes: 831 The new submatrix must have the same size and communicator as that block of the nest. 832 833 This increments the reference count of the submatrix. 834 835 Level: developer 836 837 .seealso: MatNestSetSubMats(), MatNestGetSubMats() 838 @*/ 839 PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub) 840 { 841 PetscErrorCode ierr; 842 843 PetscFunctionBegin; 844 ierr = PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 848 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 849 { 850 Mat_Nest *bA = (Mat_Nest*)A->data; 851 852 PetscFunctionBegin; 853 if (M) *M = bA->nr; 854 if (N) *N = bA->nc; 855 if (mat) *mat = bA->m; 856 PetscFunctionReturn(0); 857 } 858 859 /*@C 860 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 861 862 Not collective 863 864 Input Parameters: 865 . A - nest matrix 866 867 Output Parameter: 868 + M - number of rows in the nest matrix 869 . N - number of cols in the nest matrix 870 - mat - 2d array of matrices 871 872 Notes: 873 874 The user should not free the array mat. 875 876 In Fortran, this routine has a calling sequence 877 $ call MatNestGetSubMats(A, M, N, mat, ierr) 878 where the space allocated for the optional argument mat is assumed large enough (if provided). 879 880 Level: developer 881 882 .seealso: MatNestGetSize(), MatNestGetSubMat() 883 @*/ 884 PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 885 { 886 PetscErrorCode ierr; 887 888 PetscFunctionBegin; 889 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 890 PetscFunctionReturn(0); 891 } 892 893 PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 894 { 895 Mat_Nest *bA = (Mat_Nest*)A->data; 896 897 PetscFunctionBegin; 898 if (M) *M = bA->nr; 899 if (N) *N = bA->nc; 900 PetscFunctionReturn(0); 901 } 902 903 /*@ 904 MatNestGetSize - Returns the size of the nest matrix. 905 906 Not collective 907 908 Input Parameters: 909 . A - nest matrix 910 911 Output Parameter: 912 + M - number of rows in the nested mat 913 - N - number of cols in the nested mat 914 915 Notes: 916 917 Level: developer 918 919 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 920 @*/ 921 PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 922 { 923 PetscErrorCode ierr; 924 925 PetscFunctionBegin; 926 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 927 PetscFunctionReturn(0); 928 } 929 930 static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[]) 931 { 932 Mat_Nest *vs = (Mat_Nest*)A->data; 933 PetscInt i; 934 935 PetscFunctionBegin; 936 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i]; 937 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i]; 938 PetscFunctionReturn(0); 939 } 940 941 /*@C 942 MatNestGetISs - Returns the index sets partitioning the row and column spaces 943 944 Not collective 945 946 Input Parameters: 947 . A - nest matrix 948 949 Output Parameter: 950 + rows - array of row index sets 951 - cols - array of column index sets 952 953 Level: advanced 954 955 Notes: 956 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 957 958 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs() 959 @*/ 960 PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[]) 961 { 962 PetscErrorCode ierr; 963 964 PetscFunctionBegin; 965 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 966 ierr = PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 967 PetscFunctionReturn(0); 968 } 969 970 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[]) 971 { 972 Mat_Nest *vs = (Mat_Nest*)A->data; 973 PetscInt i; 974 975 PetscFunctionBegin; 976 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i]; 977 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i]; 978 PetscFunctionReturn(0); 979 } 980 981 /*@C 982 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces 983 984 Not collective 985 986 Input Parameters: 987 . A - nest matrix 988 989 Output Parameter: 990 + rows - array of row index sets (or NULL to ignore) 991 - cols - array of column index sets (or NULL to ignore) 992 993 Level: advanced 994 995 Notes: 996 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 997 998 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs() 999 @*/ 1000 PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[]) 1001 { 1002 PetscErrorCode ierr; 1003 1004 PetscFunctionBegin; 1005 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1006 ierr = PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 1007 PetscFunctionReturn(0); 1008 } 1009 1010 PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype) 1011 { 1012 PetscErrorCode ierr; 1013 PetscBool flg; 1014 1015 PetscFunctionBegin; 1016 ierr = PetscStrcmp(vtype,VECNEST,&flg);CHKERRQ(ierr); 1017 /* In reality, this only distinguishes VECNEST and "other" */ 1018 if (flg) A->ops->getvecs = MatCreateVecs_Nest; 1019 else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0; 1020 PetscFunctionReturn(0); 1021 } 1022 1023 /*@C 1024 MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs() 1025 1026 Not collective 1027 1028 Input Parameters: 1029 + A - nest matrix 1030 - vtype - type to use for creating vectors 1031 1032 Notes: 1033 1034 Level: developer 1035 1036 .seealso: MatCreateVecs() 1037 @*/ 1038 PetscErrorCode MatNestSetVecType(Mat A,VecType vtype) 1039 { 1040 PetscErrorCode ierr; 1041 1042 PetscFunctionBegin; 1043 ierr = PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));CHKERRQ(ierr); 1044 PetscFunctionReturn(0); 1045 } 1046 1047 PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1048 { 1049 Mat_Nest *s = (Mat_Nest*)A->data; 1050 PetscInt i,j,m,n,M,N; 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 s->nr = nr; 1055 s->nc = nc; 1056 1057 /* Create space for submatrices */ 1058 ierr = PetscMalloc1(nr,&s->m);CHKERRQ(ierr); 1059 for (i=0; i<nr; i++) { 1060 ierr = PetscMalloc1(nc,&s->m[i]);CHKERRQ(ierr); 1061 } 1062 for (i=0; i<nr; i++) { 1063 for (j=0; j<nc; j++) { 1064 s->m[i][j] = a[i*nc+j]; 1065 if (a[i*nc+j]) { 1066 ierr = PetscObjectReference((PetscObject)a[i*nc+j]);CHKERRQ(ierr); 1067 } 1068 } 1069 } 1070 1071 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1072 1073 ierr = PetscMalloc1(nr,&s->row_len);CHKERRQ(ierr); 1074 ierr = PetscMalloc1(nc,&s->col_len);CHKERRQ(ierr); 1075 for (i=0; i<nr; i++) s->row_len[i]=-1; 1076 for (j=0; j<nc; j++) s->col_len[j]=-1; 1077 1078 ierr = MatNestGetSizes_Private(A,&m,&n,&M,&N);CHKERRQ(ierr); 1079 1080 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1081 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1082 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1083 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1084 1085 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1086 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1087 1088 ierr = PetscCalloc2(nr,&s->left,nc,&s->right);CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 /*@ 1093 MatNestSetSubMats - Sets the nested submatrices 1094 1095 Collective on Mat 1096 1097 Input Parameter: 1098 + N - nested matrix 1099 . nr - number of nested row blocks 1100 . is_row - index sets for each nested row block, or NULL to make contiguous 1101 . nc - number of nested column blocks 1102 . is_col - index sets for each nested column block, or NULL to make contiguous 1103 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1104 1105 Level: advanced 1106 1107 .seealso: MatCreateNest(), MATNEST 1108 @*/ 1109 PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1110 { 1111 PetscErrorCode ierr; 1112 PetscInt i,nr_nc; 1113 1114 PetscFunctionBegin; 1115 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1116 if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1117 if (nr && is_row) { 1118 PetscValidPointer(is_row,3); 1119 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1120 } 1121 if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1122 if (nc && is_col) { 1123 PetscValidPointer(is_col,5); 1124 for (i=0; i<nc; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1125 } 1126 nr_nc=nr*nc; 1127 if (nr_nc) PetscValidPointer(a,6); 1128 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog) 1133 { 1134 PetscErrorCode ierr; 1135 PetscBool flg; 1136 PetscInt i,j,m,mi,*ix; 1137 1138 PetscFunctionBegin; 1139 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1140 if (islocal[i]) { 1141 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1142 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1143 } else { 1144 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1145 } 1146 m += mi; 1147 } 1148 if (flg) { 1149 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1150 for (i=0,n=0; i<n; i++) { 1151 ISLocalToGlobalMapping smap = NULL; 1152 VecScatter scat; 1153 IS isreq; 1154 Vec lvec,gvec; 1155 union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x; 1156 Mat sub; 1157 1158 if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers"); 1159 if (colflg) { 1160 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1161 } else { 1162 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1163 } 1164 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr);} 1165 if (islocal[i]) { 1166 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1167 } else { 1168 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1169 } 1170 for (j=0; j<mi; j++) ix[m+j] = j; 1171 if (smap) {ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);} 1172 /* 1173 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1174 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1175 The approach here is ugly because it uses VecScatter to move indices. 1176 */ 1177 ierr = VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);CHKERRQ(ierr); 1178 ierr = VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);CHKERRQ(ierr); 1179 ierr = ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);CHKERRQ(ierr); 1180 ierr = VecScatterCreate(gvec,isreq,lvec,NULL,&scat);CHKERRQ(ierr); 1181 ierr = VecGetArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1182 for (j=0; j<mi; j++) x[j].integer = ix[m+j]; 1183 ierr = VecRestoreArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1184 ierr = VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1185 ierr = VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1186 ierr = VecGetArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1187 for (j=0; j<mi; j++) ix[m+j] = x[j].integer; 1188 ierr = VecRestoreArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1189 ierr = VecDestroy(&lvec);CHKERRQ(ierr); 1190 ierr = VecDestroy(&gvec);CHKERRQ(ierr); 1191 ierr = ISDestroy(&isreq);CHKERRQ(ierr); 1192 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 1193 m += mi; 1194 } 1195 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1196 } else { 1197 *ltog = NULL; 1198 } 1199 PetscFunctionReturn(0); 1200 } 1201 1202 1203 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1204 /* 1205 nprocessors = NP 1206 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1207 proc 0: => (g_0,h_0,) 1208 proc 1: => (g_1,h_1,) 1209 ... 1210 proc nprocs-1: => (g_NP-1,h_NP-1,) 1211 1212 proc 0: proc 1: proc nprocs-1: 1213 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1214 1215 proc 0: 1216 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1217 proc 1: 1218 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1219 1220 proc NP-1: 1221 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1222 */ 1223 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1224 { 1225 Mat_Nest *vs = (Mat_Nest*)A->data; 1226 PetscInt i,j,offset,n,nsum,bs; 1227 PetscErrorCode ierr; 1228 Mat sub = NULL; 1229 1230 PetscFunctionBegin; 1231 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1232 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1233 if (is_row) { /* valid IS is passed in */ 1234 /* refs on is[] are incremeneted */ 1235 for (i=0; i<vs->nr; i++) { 1236 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1237 1238 vs->isglobal.row[i] = is_row[i]; 1239 } 1240 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1241 nsum = 0; 1242 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1243 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1244 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1245 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1246 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1247 nsum += n; 1248 } 1249 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1250 offset -= nsum; 1251 for (i=0; i<vs->nr; i++) { 1252 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1253 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1254 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1255 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1256 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1257 offset += n; 1258 } 1259 } 1260 1261 if (is_col) { /* valid IS is passed in */ 1262 /* refs on is[] are incremeneted */ 1263 for (j=0; j<vs->nc; j++) { 1264 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1265 1266 vs->isglobal.col[j] = is_col[j]; 1267 } 1268 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1269 offset = A->cmap->rstart; 1270 nsum = 0; 1271 for (j=0; j<vs->nc; j++) { 1272 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1273 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1274 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1275 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1276 nsum += n; 1277 } 1278 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1279 offset -= nsum; 1280 for (j=0; j<vs->nc; j++) { 1281 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1282 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1283 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1284 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1285 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1286 offset += n; 1287 } 1288 } 1289 1290 /* Set up the local ISs */ 1291 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1292 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1293 for (i=0,offset=0; i<vs->nr; i++) { 1294 IS isloc; 1295 ISLocalToGlobalMapping rmap = NULL; 1296 PetscInt nlocal,bs; 1297 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1298 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1299 if (rmap) { 1300 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1301 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1302 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1303 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1304 } else { 1305 nlocal = 0; 1306 isloc = NULL; 1307 } 1308 vs->islocal.row[i] = isloc; 1309 offset += nlocal; 1310 } 1311 for (i=0,offset=0; i<vs->nc; i++) { 1312 IS isloc; 1313 ISLocalToGlobalMapping cmap = NULL; 1314 PetscInt nlocal,bs; 1315 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1316 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1317 if (cmap) { 1318 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1319 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1320 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1321 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1322 } else { 1323 nlocal = 0; 1324 isloc = NULL; 1325 } 1326 vs->islocal.col[i] = isloc; 1327 offset += nlocal; 1328 } 1329 1330 /* Set up the aggregate ISLocalToGlobalMapping */ 1331 { 1332 ISLocalToGlobalMapping rmap,cmap; 1333 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1334 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1335 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1336 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1337 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1338 } 1339 1340 #if defined(PETSC_USE_DEBUG) 1341 for (i=0; i<vs->nr; i++) { 1342 for (j=0; j<vs->nc; j++) { 1343 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1344 Mat B = vs->m[i][j]; 1345 if (!B) continue; 1346 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1347 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1348 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1349 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1350 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1351 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1352 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); 1353 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); 1354 } 1355 } 1356 #endif 1357 1358 /* Set A->assembled if all non-null blocks are currently assembled */ 1359 for (i=0; i<vs->nr; i++) { 1360 for (j=0; j<vs->nc; j++) { 1361 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1362 } 1363 } 1364 A->assembled = PETSC_TRUE; 1365 PetscFunctionReturn(0); 1366 } 1367 1368 /*@C 1369 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1370 1371 Collective on Mat 1372 1373 Input Parameter: 1374 + comm - Communicator for the new Mat 1375 . nr - number of nested row blocks 1376 . is_row - index sets for each nested row block, or NULL to make contiguous 1377 . nc - number of nested column blocks 1378 . is_col - index sets for each nested column block, or NULL to make contiguous 1379 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1380 1381 Output Parameter: 1382 . B - new matrix 1383 1384 Level: advanced 1385 1386 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST 1387 @*/ 1388 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1389 { 1390 Mat A; 1391 PetscErrorCode ierr; 1392 1393 PetscFunctionBegin; 1394 *B = 0; 1395 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1396 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1397 A->preallocated = PETSC_TRUE; 1398 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1399 *B = A; 1400 PetscFunctionReturn(0); 1401 } 1402 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 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1553 { 1554 PetscErrorCode ierr; 1555 Mat_Nest *nest = (Mat_Nest*)A->data; 1556 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1557 PetscInt cstart,cend; 1558 PetscMPIInt size; 1559 Mat C; 1560 1561 PetscFunctionBegin; 1562 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1563 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1564 PetscInt nf; 1565 PetscBool fast; 1566 1567 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1568 if (!fast) { 1569 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1570 } 1571 for (i=0; i<nest->nr && fast; ++i) { 1572 for (j=0; j<nest->nc && fast; ++j) { 1573 Mat B = nest->m[i][j]; 1574 if (B) { 1575 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1576 if (!fast) { 1577 PetscBool istrans; 1578 1579 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1580 if (istrans) { 1581 Mat Bt; 1582 1583 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1584 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr); 1585 } 1586 } 1587 } 1588 } 1589 } 1590 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1591 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1592 if (fast) { 1593 PetscInt f,s; 1594 1595 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 1596 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1597 else { 1598 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 1599 nf += f; 1600 } 1601 } 1602 } 1603 for (i=0, nf=0; i<nest->nc && fast; ++i) { 1604 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1605 if (fast) { 1606 PetscInt f,s; 1607 1608 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 1609 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1610 else { 1611 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 1612 nf += f; 1613 } 1614 } 1615 } 1616 if (fast) { 1617 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 1618 PetscFunctionReturn(0); 1619 } 1620 } 1621 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1622 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1623 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1624 switch (reuse) { 1625 case MAT_INITIAL_MATRIX: 1626 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1627 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1628 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1629 *newmat = C; 1630 break; 1631 case MAT_REUSE_MATRIX: 1632 C = *newmat; 1633 break; 1634 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1635 } 1636 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1637 onnz = dnnz + m; 1638 for (k=0; k<m; k++) { 1639 dnnz[k] = 0; 1640 onnz[k] = 0; 1641 } 1642 for (j=0; j<nest->nc; ++j) { 1643 IS bNis; 1644 PetscInt bN; 1645 const PetscInt *bNindices; 1646 /* Using global column indices and ISAllGather() is not scalable. */ 1647 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1648 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1649 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1650 for (i=0; i<nest->nr; ++i) { 1651 PetscSF bmsf; 1652 PetscSFNode *iremote; 1653 Mat B; 1654 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1655 const PetscInt *bmindices; 1656 B = nest->m[i][j]; 1657 if (!B) continue; 1658 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1659 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1660 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1661 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1662 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1663 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1664 for (k = 0; k < bm; ++k){ 1665 sub_dnnz[k] = 0; 1666 sub_onnz[k] = 0; 1667 } 1668 /* 1669 Locate the owners for all of the locally-owned global row indices for this row block. 1670 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1671 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1672 */ 1673 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1674 for (br = 0; br < bm; ++br) { 1675 PetscInt row = bmindices[br], rowowner = 0, brncols, col; 1676 const PetscInt *brcols; 1677 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1678 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1679 /* how many roots */ 1680 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1681 /* get nonzero pattern */ 1682 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1683 for (k=0; k<brncols; k++) { 1684 col = bNindices[brcols[k]]; 1685 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 1686 sub_dnnz[br]++; 1687 } else { 1688 sub_onnz[br]++; 1689 } 1690 } 1691 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1692 } 1693 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1694 /* bsf will have to take care of disposing of bedges. */ 1695 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1696 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1697 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1698 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1699 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1700 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1701 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1702 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1703 } 1704 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1705 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1706 } 1707 /* Resize preallocation if overestimated */ 1708 for (i=0;i<m;i++) { 1709 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 1710 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 1711 } 1712 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1713 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1714 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1715 1716 /* Fill by row */ 1717 for (j=0; j<nest->nc; ++j) { 1718 /* Using global column indices and ISAllGather() is not scalable. */ 1719 IS bNis; 1720 PetscInt bN; 1721 const PetscInt *bNindices; 1722 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1723 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1724 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1725 for (i=0; i<nest->nr; ++i) { 1726 Mat B; 1727 PetscInt bm, br; 1728 const PetscInt *bmindices; 1729 B = nest->m[i][j]; 1730 if (!B) continue; 1731 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1732 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1733 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1734 for (br = 0; br < bm; ++br) { 1735 PetscInt row = bmindices[br], brncols, *cols; 1736 const PetscInt *brcols; 1737 const PetscScalar *brcoldata; 1738 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1739 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1740 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1741 /* 1742 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1743 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1744 */ 1745 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1746 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1747 ierr = PetscFree(cols);CHKERRQ(ierr); 1748 } 1749 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1750 } 1751 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1752 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1753 } 1754 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1755 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1756 PetscFunctionReturn(0); 1757 } 1758 1759 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has) 1760 { 1761 PetscFunctionBegin; 1762 *has = PETSC_FALSE; 1763 if (op == MATOP_MULT_TRANSPOSE) { 1764 Mat_Nest *bA = (Mat_Nest*)mat->data; 1765 PetscInt i,j,nr = bA->nr,nc = bA->nc; 1766 PetscErrorCode ierr; 1767 PetscBool flg; 1768 1769 for (j=0; j<nc; j++) { 1770 for (i=0; i<nr; i++) { 1771 if (!bA->m[i][j]) continue; 1772 ierr = MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);CHKERRQ(ierr); 1773 if (!flg) PetscFunctionReturn(0); 1774 } 1775 } 1776 } 1777 if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 1778 PetscFunctionReturn(0); 1779 } 1780 1781 /*MC 1782 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1783 1784 Level: intermediate 1785 1786 Notes: 1787 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1788 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1789 It is usually used with DMComposite and DMCreateMatrix() 1790 1791 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 1792 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 1793 than the nest matrix. 1794 1795 .seealso: MatCreate(), MatType, MatCreateNest() 1796 M*/ 1797 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1798 { 1799 Mat_Nest *s; 1800 PetscErrorCode ierr; 1801 1802 PetscFunctionBegin; 1803 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1804 A->data = (void*)s; 1805 1806 s->nr = -1; 1807 s->nc = -1; 1808 s->m = NULL; 1809 s->splitassembly = PETSC_FALSE; 1810 1811 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1812 1813 A->ops->mult = MatMult_Nest; 1814 A->ops->multadd = MatMultAdd_Nest; 1815 A->ops->multtranspose = MatMultTranspose_Nest; 1816 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1817 A->ops->transpose = MatTranspose_Nest; 1818 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1819 A->ops->assemblyend = MatAssemblyEnd_Nest; 1820 A->ops->zeroentries = MatZeroEntries_Nest; 1821 A->ops->copy = MatCopy_Nest; 1822 A->ops->duplicate = MatDuplicate_Nest; 1823 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 1824 A->ops->destroy = MatDestroy_Nest; 1825 A->ops->view = MatView_Nest; 1826 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1827 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1828 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1829 A->ops->getdiagonal = MatGetDiagonal_Nest; 1830 A->ops->diagonalscale = MatDiagonalScale_Nest; 1831 A->ops->scale = MatScale_Nest; 1832 A->ops->shift = MatShift_Nest; 1833 A->ops->diagonalset = MatDiagonalSet_Nest; 1834 A->ops->setrandom = MatSetRandom_Nest; 1835 A->ops->hasoperation = MatHasOperation_Nest; 1836 1837 A->spptr = 0; 1838 A->assembled = PETSC_FALSE; 1839 1840 /* expose Nest api's */ 1841 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1842 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1843 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1844 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1845 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1846 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1847 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1848 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1849 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1850 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 1851 1852 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1853 PetscFunctionReturn(0); 1854 } 1855