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