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