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