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