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