1 #define PETSCMAT_DLL 2 3 #include "matnestimpl.h" 4 5 /* private functions */ 6 #undef __FUNCT__ 7 #define __FUNCT__ "MatNestGetSize_Recursive" 8 static PetscErrorCode MatNestGetSize_Recursive(Mat A,PetscBool globalsize,PetscInt *M,PetscInt *N) 9 { 10 Mat_Nest *bA = (Mat_Nest*)A->data; 11 PetscInt sizeM,sizeN,i,j,index; 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 /* rows */ 16 /* now descend recursively */ 17 for (i=0; i<bA->nr; i++) { 18 for (j=0; j<bA->nc; j++) { 19 if (bA->m[i][j]) { index = j; break; } 20 } 21 if (bA->m[i][index]) { 22 if (globalsize) { ierr = MatGetSize(bA->m[i][index],&sizeM,&sizeN);CHKERRQ(ierr); } 23 else { ierr = MatGetLocalSize(bA->m[i][index],&sizeM,&sizeN);CHKERRQ(ierr); } 24 *M = *M + sizeM; 25 } 26 } 27 28 /* cols */ 29 /* now descend recursively */ 30 for (j=0; j<bA->nc; j++) { 31 for (i=0; i<bA->nr; i++) { 32 if (bA->m[i][j]) { index = i; break; } 33 } 34 if (bA->m[index][j]) { 35 if (globalsize) { ierr = MatGetSize(bA->m[index][j],&sizeM,&sizeN);CHKERRQ(ierr); } 36 else { ierr = MatGetLocalSize(bA->m[index][j],&sizeM,&sizeN);CHKERRQ(ierr); } 37 *N = *N + sizeN; 38 } 39 } 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "PETSc_MatNest_CheckMatVecCompatibility2" 45 static PetscErrorCode PETSc_MatNest_CheckMatVecCompatibility2(Mat A,Vec x,Vec y) 46 { 47 PetscBool isANest, isxNest, isyNest; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 isANest = isxNest = isyNest = PETSC_FALSE; 52 ierr = PetscTypeCompare( (PetscObject)A, MATNEST, &isANest );CHKERRQ(ierr); 53 ierr = PetscTypeCompare( (PetscObject)x, VECNEST, &isxNest );CHKERRQ(ierr); 54 ierr = PetscTypeCompare( (PetscObject)y, VECNEST, &isyNest );CHKERRQ(ierr); 55 56 if (isANest && isxNest && isyNest) { 57 PetscFunctionReturn(0); 58 } else { 59 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_SUP, "Operation only valid for Mat(Nest)-Vec(Nest) combinations"); 60 PetscFunctionReturn(0); 61 } 62 PetscFunctionReturn(0); 63 } 64 65 /* 66 This function is called every time we insert a sub matrix the Nest. 67 It ensures that every Nest along row r, and col c has the same dimensions 68 as the submat being inserted. 69 */ 70 #undef __FUNCT__ 71 #define __FUNCT__ "PETSc_MatNest_CheckConsistency" 72 PETSC_UNUSED 73 static PetscErrorCode PETSc_MatNest_CheckConsistency(Mat A,Mat submat,PetscInt r,PetscInt c) 74 { 75 Mat_Nest *b = (Mat_Nest*)A->data; 76 PetscInt i,j; 77 PetscInt nr,nc; 78 PetscInt sM,sN,M,N; 79 Mat mat; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 nr = b->nr; 84 nc = b->nc; 85 ierr = MatGetSize(submat,&sM,&sN);CHKERRQ(ierr); 86 for (i=0; i<nr; i++) { 87 mat = b->m[i][c]; 88 if (mat) { 89 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 90 /* Check columns match */ 91 if (sN != N) { 92 SETERRQ3(((PetscObject)A)->comm,PETSC_ERR_SUP,"Inserted incompatible submatrix into Nest at (%D,%D). Submatrix must have %D rows",r,c,N ); 93 } 94 } 95 } 96 97 for (j=0; j<nc; j++) { 98 mat = b->m[r][j]; 99 if (mat) { 100 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 101 /* Check rows match */ 102 if (sM != M) { 103 SETERRQ3(((PetscObject)A)->comm,PETSC_ERR_SUP,"Inserted incompatible submatrix into Nest at (%D,%D). Submatrix must have %D cols",r,c,M ); 104 } 105 } 106 } 107 PetscFunctionReturn(0); 108 } 109 110 /* 111 Checks if the row/col's contain a complete set of IS's. 112 Once they are filled, the offsets are computed. This saves having to update 113 the index set entries each time we insert something new. 114 */ 115 #undef __FUNCT__ 116 #define __FUNCT__ "PETSc_MatNest_UpdateStructure" 117 PETSC_UNUSED 118 static PetscErrorCode PETSc_MatNest_UpdateStructure(Mat A,PetscInt ridx,PetscInt cidx) 119 { 120 Mat_Nest *b = (Mat_Nest*)A->data; 121 PetscInt m,n,i,j; 122 PetscBool fullrow,fullcol; 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 ierr = MatGetLocalSize(b->m[ridx][cidx],&m,&n);CHKERRQ(ierr); 127 b->row_len[ridx] = m; 128 b->col_len[cidx] = n; 129 130 fullrow = PETSC_TRUE; 131 for (i=0; i<b->nr; i++) { 132 if (b->row_len[i]==-1) { fullrow = PETSC_FALSE; break; } 133 } 134 if ( (fullrow) && (!b->isglobal.row[0]) ){ 135 PetscInt cnt = 0; 136 for (i=0; i<b->nr; i++) { 137 ierr = ISCreateStride(PETSC_COMM_SELF,b->row_len[i],cnt,1,&b->isglobal.row[i]);CHKERRQ(ierr); 138 cnt = cnt + b->row_len[i]; 139 } 140 } 141 142 fullcol = PETSC_TRUE; 143 for (j=0; j<b->nc; j++) { 144 if (b->col_len[j]==-1) { fullcol = PETSC_FALSE; break; } 145 } 146 if( (fullcol) && (!b->isglobal.col[0]) ){ 147 PetscInt cnt = 0; 148 for (j=0; j<b->nc; j++) { 149 ierr = ISCreateStride(PETSC_COMM_SELF,b->col_len[j],cnt,1,&b->isglobal.col[j]);CHKERRQ(ierr); 150 cnt = cnt + b->col_len[j]; 151 } 152 } 153 PetscFunctionReturn(0); 154 } 155 156 /* operations */ 157 #undef __FUNCT__ 158 #define __FUNCT__ "MatMult_Nest" 159 PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y) 160 { 161 Mat_Nest *bA = (Mat_Nest*)A->data; 162 Vec *bx; 163 Vec *by; 164 PetscInt i,j; 165 PetscErrorCode ierr; 166 167 PetscFunctionBegin; 168 ierr = PETSc_MatNest_CheckMatVecCompatibility2(A,x,y);CHKERRQ(ierr); 169 ierr = VecNestGetSubVecs(x,PETSC_NULL,&bx);CHKERRQ(ierr); 170 ierr = VecNestGetSubVecs(y,PETSC_NULL,&by);CHKERRQ(ierr); 171 172 for (i=0; i<bA->nr; i++) { 173 ierr = VecZeroEntries(by[i]);CHKERRQ(ierr); 174 for (j=0; j<bA->nc; j++) { 175 if (!bA->m[i][j]) { 176 continue; 177 } 178 /* y[i] <- y[i] + A[i][j] * x[j] */ 179 ierr = MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);CHKERRQ(ierr); 180 } 181 } 182 PetscFunctionReturn(0); 183 } 184 185 #undef __FUNCT__ 186 #define __FUNCT__ "MatMultTranspose_Nest" 187 PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y) 188 { 189 Mat_Nest *bA = (Mat_Nest*)A->data; 190 Vec *bx; 191 Vec *by; 192 PetscInt i,j; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 ierr = PETSc_MatNest_CheckMatVecCompatibility2(A,x,y);CHKERRQ(ierr); 197 ierr = VecNestGetSubVecs(x,PETSC_NULL,&bx);CHKERRQ(ierr); 198 ierr = VecNestGetSubVecs(y,PETSC_NULL,&by);CHKERRQ(ierr); 199 if (A->symmetric) { 200 ierr = MatMult_Nest(A,x,y);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 for (i=0; i<bA->nr; i++) { 204 ierr = VecZeroEntries(by[i]);CHKERRQ(ierr); 205 for (j=0; j<bA->nc; j++) { 206 if (!bA->m[j][i]) { 207 continue; 208 } 209 /* y[i] <- y[i] + A^T[i][j] * x[j], so we swap i,j in mat[][] */ 210 ierr = MatMultTransposeAdd(bA->m[j][i],bx[j],by[i],by[i]);CHKERRQ(ierr); 211 } 212 } 213 PetscFunctionReturn(0); 214 } 215 216 /* Returns the sum of the global size of all the consituent vectors in the nest */ 217 #undef __FUNCT__ 218 #define __FUNCT__ "MatGetSize_Nest" 219 PetscErrorCode MatGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 220 { 221 PetscFunctionBegin; 222 if (M) { *M = A->rmap->N; } 223 if (N) { *N = A->cmap->N; } 224 PetscFunctionReturn(0); 225 } 226 227 /* Returns the sum of the local size of all the consituent vectors in the nest */ 228 #undef __FUNCT__ 229 #define __FUNCT__ "MatGetLocalSize_Nest" 230 PetscErrorCode MatGetLocalSize_Nest(Mat A,PetscInt *m,PetscInt *n) 231 { 232 PetscFunctionBegin; 233 if (m) { *m = A->rmap->n; } 234 if (n) { *n = A->cmap->n; } 235 PetscFunctionReturn(0); 236 } 237 238 #undef __FUNCT__ 239 #define __FUNCT__ "MatNestDestroyISList" 240 static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list) 241 { 242 PetscErrorCode ierr; 243 IS *lst = *list; 244 PetscInt i; 245 246 PetscFunctionBegin; 247 if (!lst) PetscFunctionReturn(0); 248 for (i=0; i<n; i++) if (lst[i]) {ierr = ISDestroy(lst[i]);CHKERRQ(ierr);} 249 ierr = PetscFree(lst);CHKERRQ(ierr); 250 *list = PETSC_NULL; 251 PetscFunctionReturn(0); 252 } 253 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatDestroy_Nest" 257 PetscErrorCode MatDestroy_Nest(Mat A) 258 { 259 Mat_Nest *vs = (Mat_Nest*)A->data; 260 PetscInt i,j; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 /* release the matrices and the place holders */ 265 ierr = MatNestDestroyISList(vs->nr,&vs->isglobal.row);CHKERRQ(ierr); 266 ierr = MatNestDestroyISList(vs->nc,&vs->isglobal.col);CHKERRQ(ierr); 267 ierr = MatNestDestroyISList(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 268 ierr = MatNestDestroyISList(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 269 270 ierr = PetscFree(vs->row_len);CHKERRQ(ierr); 271 ierr = PetscFree(vs->col_len);CHKERRQ(ierr); 272 273 /* release the matrices and the place holders */ 274 if (vs->m) { 275 for (i=0; i<vs->nr; i++) { 276 for (j=0; j<vs->nc; j++) { 277 278 if (vs->m[i][j]) { 279 ierr = MatDestroy(vs->m[i][j]);CHKERRQ(ierr); 280 vs->m[i][j] = PETSC_NULL; 281 } 282 } 283 ierr = PetscFree( vs->m[i] );CHKERRQ(ierr); 284 vs->m[i] = PETSC_NULL; 285 } 286 ierr = PetscFree(vs->m);CHKERRQ(ierr); 287 vs->m = PETSC_NULL; 288 } 289 ierr = PetscFree(vs);CHKERRQ(ierr); 290 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "MatAssemblyBegin_Nest" 296 PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type) 297 { 298 Mat_Nest *vs = (Mat_Nest*)A->data; 299 PetscInt i,j; 300 PetscErrorCode ierr; 301 302 PetscFunctionBegin; 303 for (i=0; i<vs->nr; i++) { 304 for (j=0; j<vs->nc; j++) { 305 if (vs->m[i][j]) { ierr = MatAssemblyBegin(vs->m[i][j],type);CHKERRQ(ierr); } 306 } 307 } 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "MatAssemblyEnd_Nest" 313 PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type) 314 { 315 Mat_Nest *vs = (Mat_Nest*)A->data; 316 PetscInt i,j; 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 for (i=0; i<vs->nr; i++) { 321 for (j=0; j<vs->nc; j++) { 322 if (vs->m[i][j]) { ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr); } 323 } 324 } 325 PetscFunctionReturn(0); 326 } 327 328 #undef __FUNCT__ 329 #define __FUNCT__ "MatNestFindNonzeroSubMatRow" 330 static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B) 331 { 332 Mat_Nest *vs = (Mat_Nest*)A->data; 333 PetscInt j; 334 Mat sub; 335 336 PetscFunctionBegin; 337 sub = vs->m[row][row]; /* Prefer to find on the diagonal */ 338 for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j]; 339 if (!sub) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",row); 340 *B = sub; 341 PetscFunctionReturn(0); 342 } 343 344 #undef __FUNCT__ 345 #define __FUNCT__ "MatNestFindNonzeroSubMatCol" 346 static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B) 347 { 348 Mat_Nest *vs = (Mat_Nest*)A->data; 349 PetscInt i; 350 Mat sub; 351 352 PetscFunctionBegin; 353 sub = vs->m[col][col]; /* Prefer to find on the diagonal */ 354 for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col]; 355 if (!sub) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",col); 356 *B = sub; 357 PetscFunctionReturn(0); 358 } 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "MatNestFindIS" 362 static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found) 363 { 364 PetscErrorCode ierr; 365 PetscInt i; 366 PetscBool flg; 367 368 PetscFunctionBegin; 369 PetscValidPointer(list,3); 370 PetscValidHeaderSpecific(is,IS_CLASSID,4); 371 PetscValidIntPointer(found,5); 372 *found = -1; 373 for (i=0; i<n; i++) { 374 ierr = ISEqual(list[i],is,&flg);CHKERRQ(ierr); 375 if (flg) { 376 *found = i; 377 PetscFunctionReturn(0); 378 } 379 } 380 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_INCOMP,"Could not find index set"); 381 PetscFunctionReturn(0); 382 } 383 384 #undef __FUNCT__ 385 #define __FUNCT__ "MatNestFindSubMat" 386 static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B) 387 { 388 PetscErrorCode ierr; 389 Mat_Nest *vs = (Mat_Nest*)A->data; 390 PetscInt row,col; 391 392 PetscFunctionBegin; 393 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr); 394 ierr = MatNestFindIS(A,vs->nc,is->col,iscol,&col);CHKERRQ(ierr); 395 *B = vs->m[row][col]; 396 PetscFunctionReturn(0); 397 } 398 399 #undef __FUNCT__ 400 #define __FUNCT__ "MatGetSubMatrix_Nest" 401 static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B) 402 { 403 PetscErrorCode ierr; 404 Mat_Nest *vs = (Mat_Nest*)A->data; 405 Mat sub; 406 407 PetscFunctionBegin; 408 ierr = MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);CHKERRQ(ierr); 409 switch (reuse) { 410 case MAT_INITIAL_MATRIX: 411 ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr); 412 *B = sub; 413 break; 414 case MAT_REUSE_MATRIX: 415 if (sub != *B) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call"); 416 break; 417 case MAT_IGNORE_MATRIX: /* Nothing to do */ 418 break; 419 } 420 PetscFunctionReturn(0); 421 } 422 423 #undef __FUNCT__ 424 #define __FUNCT__ "MatGetLocalSubMatrix_Nest" 425 PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 426 { 427 PetscErrorCode ierr; 428 Mat_Nest *vs = (Mat_Nest*)A->data; 429 Mat sub; 430 431 PetscFunctionBegin; 432 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 433 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */ 434 if (sub) {ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr);} 435 *B = sub; 436 PetscFunctionReturn(0); 437 } 438 439 #undef __FUNCT__ 440 #define __FUNCT__ "MatRestoreLocalSubMatrix_Nest" 441 PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 442 { 443 PetscErrorCode ierr; 444 Mat_Nest *vs = (Mat_Nest*)A->data; 445 Mat sub; 446 447 PetscFunctionBegin; 448 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 449 if (*B != sub) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten"); 450 if (sub) { 451 if (((PetscObject)sub)->refct <= 1) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times"); 452 ierr = MatDestroy(*B);CHKERRQ(ierr); 453 } 454 *B = PETSC_NULL; 455 PetscFunctionReturn(0); 456 } 457 458 #undef __FUNCT__ 459 #define __FUNCT__ "MatGetVecs_Nest" 460 PetscErrorCode MatGetVecs_Nest(Mat A,Vec *right,Vec *left) 461 { 462 Mat_Nest *bA = (Mat_Nest*)A->data; 463 Vec *L,*R; 464 MPI_Comm comm; 465 PetscInt i,j; 466 PetscErrorCode ierr; 467 468 PetscFunctionBegin; 469 comm = ((PetscObject)A)->comm; 470 if (right) { 471 /* allocate R */ 472 ierr = PetscMalloc( sizeof(Vec) * bA->nc, &R );CHKERRQ(ierr); 473 /* Create the right vectors */ 474 for (j=0; j<bA->nc; j++) { 475 for (i=0; i<bA->nr; i++) { 476 if (bA->m[i][j]) { 477 ierr = MatGetVecs(bA->m[i][j],&R[j],PETSC_NULL);CHKERRQ(ierr); 478 break; 479 } 480 } 481 if (i==bA->nr) { 482 /* have an empty column */ 483 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 484 } 485 } 486 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr); 487 /* hand back control to the nest vector */ 488 for (j=0; j<bA->nc; j++) { 489 ierr = VecDestroy(R[j]);CHKERRQ(ierr); 490 } 491 ierr = PetscFree(R);CHKERRQ(ierr); 492 } 493 494 if (left) { 495 /* allocate L */ 496 ierr = PetscMalloc( sizeof(Vec) * bA->nr, &L );CHKERRQ(ierr); 497 /* Create the left vectors */ 498 for (i=0; i<bA->nr; i++) { 499 for (j=0; j<bA->nc; j++) { 500 if (bA->m[i][j]) { 501 ierr = MatGetVecs(bA->m[i][j],PETSC_NULL,&L[i]);CHKERRQ(ierr); 502 break; 503 } 504 } 505 if (j==bA->nc) { 506 /* have an empty row */ 507 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 508 } 509 } 510 511 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr); 512 for (i=0; i<bA->nr; i++) { 513 ierr = VecDestroy(L[i]);CHKERRQ(ierr); 514 } 515 516 ierr = PetscFree(L);CHKERRQ(ierr); 517 } 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNCT__ 522 #define __FUNCT__ "MatView_Nest" 523 PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer) 524 { 525 Mat_Nest *bA = (Mat_Nest*)A->data; 526 PetscBool isascii; 527 PetscInt i,j; 528 PetscErrorCode ierr; 529 530 PetscFunctionBegin; 531 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 532 if (isascii) { 533 534 PetscViewerASCIIPrintf(viewer,"Matrix object: \n" ); 535 PetscViewerASCIIPushTab( viewer ); /* push0 */ 536 PetscViewerASCIIPrintf( viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc); 537 538 PetscViewerASCIIPrintf(viewer,"MatNest structure: \n" ); 539 for (i=0; i<bA->nr; i++) { 540 for (j=0; j<bA->nc; j++) { 541 const MatType type; 542 const char *name; 543 PetscInt NR,NC; 544 PetscBool isNest = PETSC_FALSE; 545 546 if (!bA->m[i][j]) { 547 PetscViewerASCIIPrintf( viewer, "(%D,%D) : PETSC_NULL \n",i,j); 548 continue; 549 } 550 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr); 551 ierr = MatGetType( bA->m[i][j], &type );CHKERRQ(ierr); 552 name = ((PetscObject)bA->m[i][j])->prefix; 553 ierr = PetscTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr); 554 555 PetscViewerASCIIPrintf( viewer, "(%D,%D) : name=\"%s\", type=%s, rows=%D, cols=%D \n",i,j,name,type,NR,NC); 556 557 if (isNest) { 558 PetscViewerASCIIPushTab( viewer ); /* push1 */ 559 ierr = MatView( bA->m[i][j], viewer );CHKERRQ(ierr); 560 PetscViewerASCIIPopTab(viewer); /* pop1 */ 561 } 562 } 563 } 564 PetscViewerASCIIPopTab(viewer); /* pop0 */ 565 } 566 PetscFunctionReturn(0); 567 } 568 569 #undef __FUNCT__ 570 #define __FUNCT__ "MatZeroEntries_Nest" 571 PetscErrorCode MatZeroEntries_Nest(Mat A) 572 { 573 Mat_Nest *bA = (Mat_Nest*)A->data; 574 PetscInt i,j; 575 PetscErrorCode ierr; 576 577 PetscFunctionBegin; 578 for (i=0; i<bA->nr; i++) { 579 for (j=0; j<bA->nc; j++) { 580 if (!bA->m[i][j]) continue; 581 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr); 582 } 583 } 584 PetscFunctionReturn(0); 585 } 586 587 #undef __FUNCT__ 588 #define __FUNCT__ "MatDuplicate_Nest" 589 PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 590 { 591 Mat_Nest *bA = (Mat_Nest*)A->data; 592 Mat *b; 593 PetscInt i,j,nr = bA->nr,nc = bA->nc; 594 PetscErrorCode ierr; 595 596 PetscFunctionBegin; 597 ierr = PetscMalloc(nr*nc*sizeof(Mat),&b);CHKERRQ(ierr); 598 for (i=0; i<nr; i++) { 599 for (j=0; j<nc; j++) { 600 if (bA->m[i][j]) { 601 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 602 } else { 603 b[i*nc+j] = PETSC_NULL; 604 } 605 } 606 } 607 ierr = MatCreateNest(((PetscObject)A)->comm,nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 608 /* Give the new MatNest exclusive ownership */ 609 for (i=0; i<nr*nc; i++) { 610 if (b[i]) {ierr = MatDestroy(b[i]);CHKERRQ(ierr);} 611 } 612 ierr = PetscFree(b);CHKERRQ(ierr); 613 614 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 615 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 616 PetscFunctionReturn(0); 617 } 618 619 /* nest api */ 620 EXTERN_C_BEGIN 621 #undef __FUNCT__ 622 #define __FUNCT__ "MatNestGetSubMat_Nest" 623 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 624 { 625 Mat_Nest *bA = (Mat_Nest*)A->data; 626 PetscFunctionBegin; 627 if (idxm >= bA->nr) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 628 if (jdxm >= bA->nc) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 629 *mat = bA->m[idxm][jdxm]; 630 PetscFunctionReturn(0); 631 } 632 EXTERN_C_END 633 634 #undef __FUNCT__ 635 #define __FUNCT__ "MatNestGetSubMat" 636 /*@C 637 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 638 639 Not collective 640 641 Input Parameters: 642 . A - nest matrix 643 . idxm - index of the matrix within the nest matrix 644 . jdxm - index of the matrix within the nest matrix 645 646 Output Parameter: 647 . sub - matrix at index idxm,jdxm within the nest matrix 648 649 Notes: 650 651 Level: developer 652 653 .seealso: MatNestGetSize(), MatNestGetSubMats() 654 @*/ 655 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 656 { 657 PetscErrorCode ierr; 658 659 PetscFunctionBegin; 660 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 661 PetscFunctionReturn(0); 662 } 663 664 EXTERN_C_BEGIN 665 #undef __FUNCT__ 666 #define __FUNCT__ "MatNestGetSubMats_Nest" 667 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 668 { 669 Mat_Nest *bA = (Mat_Nest*)A->data; 670 PetscFunctionBegin; 671 if (M) { *M = bA->nr; } 672 if (N) { *N = bA->nc; } 673 if (mat) { *mat = bA->m; } 674 PetscFunctionReturn(0); 675 } 676 EXTERN_C_END 677 678 #undef __FUNCT__ 679 #define __FUNCT__ "MatNestGetSubMats" 680 /*@C 681 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 682 683 Not collective 684 685 Input Parameters: 686 . A - nest matri 687 688 Output Parameter: 689 . M - number of rows in the nest matrix 690 . N - number of cols in the nest matrix 691 . mat - 2d array of matrices 692 693 Notes: 694 695 The user should not free the array mat. 696 697 Level: developer 698 699 .seealso: MatNestGetSize(), MatNestGetSubMat() 700 @*/ 701 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 702 { 703 PetscErrorCode ierr; 704 705 PetscFunctionBegin; 706 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 707 PetscFunctionReturn(0); 708 } 709 710 EXTERN_C_BEGIN 711 #undef __FUNCT__ 712 #define __FUNCT__ "MatNestGetSize_Nest" 713 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 714 { 715 Mat_Nest *bA = (Mat_Nest*)A->data; 716 717 PetscFunctionBegin; 718 if (M) { *M = bA->nr; } 719 if (N) { *N = bA->nc; } 720 PetscFunctionReturn(0); 721 } 722 EXTERN_C_END 723 724 #undef __FUNCT__ 725 #define __FUNCT__ "MatNestGetSize" 726 /*@C 727 MatNestGetSize - Returns the size of the nest matrix. 728 729 Not collective 730 731 Input Parameters: 732 . A - nest matrix 733 734 Output Parameter: 735 . M - number of rows in the nested mat 736 . N - number of cols in the nested mat 737 738 Notes: 739 740 Level: developer 741 742 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 743 @*/ 744 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 745 { 746 PetscErrorCode ierr; 747 748 PetscFunctionBegin; 749 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 750 PetscFunctionReturn(0); 751 } 752 753 /* constructor */ 754 #undef __FUNCT__ 755 #define __FUNCT__ "MatNestSetOps_Private" 756 static PetscErrorCode MatNestSetOps_Private(struct _MatOps* ops) 757 { 758 PetscFunctionBegin; 759 /* 0*/ 760 ops->setvalues = 0; 761 ops->getrow = 0; 762 ops->restorerow = 0; 763 ops->mult = MatMult_Nest; 764 ops->multadd = 0; 765 /* 5*/ 766 ops->multtranspose = MatMultTranspose_Nest; 767 ops->multtransposeadd = 0; 768 ops->solve = 0; 769 ops->solveadd = 0; 770 ops->solvetranspose = 0; 771 /*10*/ 772 ops->solvetransposeadd = 0; 773 ops->lufactor = 0; 774 ops->choleskyfactor = 0; 775 ops->sor = 0; 776 ops->transpose = 0; 777 /*15*/ 778 ops->getinfo = 0; 779 ops->equal = 0; 780 ops->getdiagonal = 0; 781 ops->diagonalscale = 0; 782 ops->norm = 0; 783 /*20*/ 784 ops->assemblybegin = MatAssemblyBegin_Nest; 785 ops->assemblyend = MatAssemblyEnd_Nest; 786 ops->setoption = 0; 787 ops->zeroentries = MatZeroEntries_Nest; 788 /*24*/ 789 ops->zerorows = 0; 790 ops->lufactorsymbolic = 0; 791 ops->lufactornumeric = 0; 792 ops->choleskyfactorsymbolic = 0; 793 ops->choleskyfactornumeric = 0; 794 /*29*/ 795 ops->setuppreallocation = 0; 796 ops->ilufactorsymbolic = 0; 797 ops->iccfactorsymbolic = 0; 798 ops->getarray = 0; 799 ops->restorearray = 0; 800 /*34*/ 801 ops->duplicate = MatDuplicate_Nest; 802 ops->forwardsolve = 0; 803 ops->backwardsolve = 0; 804 ops->ilufactor = 0; 805 ops->iccfactor = 0; 806 /*39*/ 807 ops->axpy = 0; 808 ops->getsubmatrices = 0; 809 ops->increaseoverlap = 0; 810 ops->getvalues = 0; 811 ops->copy = 0; 812 /*44*/ 813 ops->getrowmax = 0; 814 ops->scale = 0; 815 ops->shift = 0; 816 ops->diagonalset = 0; 817 ops->zerorowscolumns = 0; 818 /*49*/ 819 ops->setblocksize = 0; 820 ops->getrowij = 0; 821 ops->restorerowij = 0; 822 ops->getcolumnij = 0; 823 ops->restorecolumnij = 0; 824 /*54*/ 825 ops->fdcoloringcreate = 0; 826 ops->coloringpatch = 0; 827 ops->setunfactored = 0; 828 ops->permute = 0; 829 ops->setvaluesblocked = 0; 830 /*59*/ 831 ops->getsubmatrix = MatGetSubMatrix_Nest; 832 ops->destroy = MatDestroy_Nest; 833 ops->view = MatView_Nest; 834 ops->convertfrom = 0; 835 ops->usescaledform = 0; 836 /*64*/ 837 ops->scalesystem = 0; 838 ops->unscalesystem = 0; 839 ops->setlocaltoglobalmapping = 0; 840 ops->setvalueslocal = 0; 841 ops->zerorowslocal = 0; 842 /*69*/ 843 ops->getrowmaxabs = 0; 844 ops->getrowminabs = 0; 845 ops->convert = 0;/*MatConvert_Nest;*/ 846 ops->setcoloring = 0; 847 ops->setvaluesadic = 0; 848 /* 74 */ 849 ops->setvaluesadifor = 0; 850 ops->fdcoloringapply = 0; 851 ops->setfromoptions = 0; 852 ops->multconstrained = 0; 853 ops->multtransposeconstrained = 0; 854 /*79*/ 855 ops->permutesparsify = 0; 856 ops->mults = 0; 857 ops->solves = 0; 858 ops->getinertia = 0; 859 ops->load = 0; 860 /*84*/ 861 ops->issymmetric = 0; 862 ops->ishermitian = 0; 863 ops->isstructurallysymmetric = 0; 864 ops->dummy4 = 0; 865 ops->getvecs = MatGetVecs_Nest; 866 /*89*/ 867 ops->matmult = 0;/*MatMatMult_Nest;*/ 868 ops->matmultsymbolic = 0; 869 ops->matmultnumeric = 0; 870 ops->ptap = 0; 871 ops->ptapsymbolic = 0; 872 /*94*/ 873 ops->ptapnumeric = 0; 874 ops->matmulttranspose = 0; 875 ops->matmulttransposesymbolic = 0; 876 ops->matmulttransposenumeric = 0; 877 ops->ptapsymbolic_seqaij = 0; 878 /*99*/ 879 ops->ptapnumeric_seqaij = 0; 880 ops->ptapsymbolic_mpiaij = 0; 881 ops->ptapnumeric_mpiaij = 0; 882 ops->conjugate = 0; 883 ops->setsizes = 0; 884 /*104*/ 885 ops->setvaluesrow = 0; 886 ops->realpart = 0; 887 ops->imaginarypart = 0; 888 ops->getrowuppertriangular = 0; 889 ops->restorerowuppertriangular = 0; 890 /*109*/ 891 ops->matsolve = 0; 892 ops->getredundantmatrix = 0; 893 ops->getrowmin = 0; 894 ops->getcolumnvector = 0; 895 ops->missingdiagonal = 0; 896 /* 114 */ 897 ops->getseqnonzerostructure = 0; 898 ops->create = 0; 899 ops->getghosts = 0; 900 ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 901 ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 902 /* 119 */ 903 ops->multdiagonalblock = 0; 904 ops->hermitiantranspose = 0; 905 ops->multhermitiantranspose = 0; 906 ops->multhermitiantransposeadd = 0; 907 ops->getmultiprocblock = 0; 908 /* 124 */ 909 ops->dummy1 = 0; 910 ops->dummy2 = 0; 911 ops->dummy3 = 0; 912 ops->dummy4 = 0; 913 ops->getsubmatricesparallel = 0; 914 915 PetscFunctionReturn(0); 916 } 917 918 #undef __FUNCT__ 919 #define __FUNCT__ "MatSetUp_Nest_Private" 920 static PetscErrorCode MatSetUp_Nest_Private(Mat A,PetscInt nr,PetscInt nc,const Mat *sub) 921 { 922 Mat_Nest *ctx = (Mat_Nest*)A->data; 923 PetscInt i,j; 924 PetscErrorCode ierr; 925 926 PetscFunctionBegin; 927 if(ctx->setup_called) PetscFunctionReturn(0); 928 929 ctx->nr = nr; 930 ctx->nc = nc; 931 932 /* Create space */ 933 ierr = PetscMalloc(sizeof(Mat*)*ctx->nr,&ctx->m);CHKERRQ(ierr); 934 for (i=0; i<ctx->nr; i++) { 935 ierr = PetscMalloc(sizeof(Mat)*ctx->nc,&ctx->m[i]);CHKERRQ(ierr); 936 } 937 for (i=0; i<ctx->nr; i++) { 938 for (j=0; j<ctx->nc; j++) { 939 ctx->m[i][j] = sub[i*nc+j]; 940 if (sub[i*nc+j]) { 941 ierr = PetscObjectReference((PetscObject)sub[i*nc+j]);CHKERRQ(ierr); 942 } 943 } 944 } 945 946 ierr = PetscMalloc(sizeof(IS)*ctx->nr,&ctx->isglobal.row);CHKERRQ(ierr); 947 ierr = PetscMalloc(sizeof(IS)*ctx->nc,&ctx->isglobal.col);CHKERRQ(ierr); 948 949 ierr = PetscMalloc(sizeof(PetscInt)*ctx->nr,&ctx->row_len);CHKERRQ(ierr); 950 ierr = PetscMalloc(sizeof(PetscInt)*ctx->nc,&ctx->col_len);CHKERRQ(ierr); 951 for (i=0;i<ctx->nr;i++) ctx->row_len[i]=-1; 952 for (j=0;j<ctx->nc;j++) ctx->col_len[j]=-1; 953 954 ctx->setup_called = PETSC_TRUE; 955 956 PetscFunctionReturn(0); 957 } 958 959 960 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 961 /* 962 nprocessors = NP 963 Nest x^T = ( (g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1) ) 964 proc 0: => (g_0,h_0,) 965 proc 1: => (g_1,h_1,) 966 ... 967 proc nprocs-1: => (g_NP-1,h_NP-1,) 968 969 proc 0: proc 1: proc nprocs-1: 970 is[0] = ( 0,1,2,...,nlocal(g_0)-1 ) ( 0,1,...,nlocal(g_1)-1 ) ( 0,1,...,nlocal(g_NP-1) ) 971 972 proc 0: 973 is[1] = ( nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1 ) 974 proc 1: 975 is[1] = ( nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1 ) 976 977 proc NP-1: 978 is[1] = ( nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1 ) 979 */ 980 #undef __FUNCT__ 981 #define __FUNCT__ "MatSetUp_NestIS_Private" 982 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 983 { 984 Mat_Nest *vs = (Mat_Nest*)A->data; 985 PetscInt i,j,offset,n,bs; 986 PetscErrorCode ierr; 987 Mat sub; 988 989 PetscFunctionBegin; 990 if (is_row) { /* valid IS is passed in */ 991 /* refs on is[] are incremeneted */ 992 for (i=0; i<vs->nr; i++) { 993 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 994 vs->isglobal.row[i] = is_row[i]; 995 } 996 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 997 offset = A->rmap->rstart; 998 for (i=0; i<vs->nr; i++) { 999 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1000 ierr = MatGetLocalSize(sub,&n,PETSC_NULL);CHKERRQ(ierr); 1001 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1002 ierr = ISCreateStride(((PetscObject)sub)->comm,n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1003 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1004 offset += n; 1005 } 1006 } 1007 1008 if (is_col) { /* valid IS is passed in */ 1009 /* refs on is[] are incremeneted */ 1010 for (j=0; j<vs->nc; j++) { 1011 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1012 vs->isglobal.col[j] = is_col[j]; 1013 } 1014 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1015 offset = A->cmap->rstart; 1016 for (j=0; j<vs->nc; j++) { 1017 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1018 ierr = MatGetLocalSize(sub,PETSC_NULL,&n);CHKERRQ(ierr); 1019 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1020 ierr = ISCreateStride(((PetscObject)sub)->comm,n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1021 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1022 offset += n; 1023 } 1024 } 1025 1026 /* Set up the local ISs */ 1027 ierr = PetscMalloc(vs->nr*sizeof(IS),&vs->islocal.row);CHKERRQ(ierr); 1028 ierr = PetscMalloc(vs->nc*sizeof(IS),&vs->islocal.col);CHKERRQ(ierr); 1029 for (i=0,offset=0; i<vs->nr; i++) { 1030 IS isloc; 1031 ISLocalToGlobalMapping rmap; 1032 PetscInt nlocal,bs; 1033 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1034 ierr = MatGetLocalToGlobalMapping(sub,&rmap,PETSC_NULL);CHKERRQ(ierr); 1035 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1036 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1037 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1038 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1039 vs->islocal.row[i] = isloc; 1040 offset += nlocal; 1041 } 1042 for (i=0,offset=0; i<vs->nr; i++) { 1043 IS isloc; 1044 ISLocalToGlobalMapping cmap; 1045 PetscInt nlocal,bs; 1046 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1047 ierr = MatGetLocalToGlobalMapping(sub,PETSC_NULL,&cmap);CHKERRQ(ierr); 1048 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1049 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1050 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1051 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1052 vs->islocal.col[i] = isloc; 1053 offset += nlocal; 1054 } 1055 PetscFunctionReturn(0); 1056 } 1057 1058 #undef __FUNCT__ 1059 #define __FUNCT__ "MatCreateNest" 1060 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1061 { 1062 Mat A; 1063 Mat_Nest *s; 1064 PetscInt i,m,n,M,N; 1065 PetscErrorCode ierr; 1066 1067 PetscFunctionBegin; 1068 if (nr < 0) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1069 if (nr && is_row) { 1070 PetscValidPointer(is_row,3); 1071 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1072 } 1073 if (nc < 0) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1074 if (nc && is_row) { 1075 PetscValidPointer(is_col,5); 1076 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1077 } 1078 if (nr*nc) PetscValidPointer(a,6); 1079 PetscValidPointer(B,7); 1080 1081 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1082 1083 ierr = PetscMalloc( sizeof(Mat_Nest), &s );CHKERRQ(ierr); 1084 ierr = PetscMemzero( s, sizeof(Mat_Nest) );CHKERRQ(ierr); 1085 A->data = (void*)s; 1086 s->setup_called = PETSC_FALSE; 1087 s->nr = s->nc = -1; 1088 s->m = PETSC_NULL; 1089 1090 /* define the operations */ 1091 ierr = MatNestSetOps_Private(A->ops);CHKERRQ(ierr); 1092 A->spptr = 0; 1093 A->same_nonzero = PETSC_FALSE; 1094 A->assembled = PETSC_FALSE; 1095 1096 ierr = PetscObjectChangeTypeName((PetscObject) A, MATNEST );CHKERRQ(ierr); 1097 ierr = MatSetUp_Nest_Private(A,nr,nc,a);CHKERRQ(ierr); 1098 1099 m = n = 0; 1100 M = N = 0; 1101 ierr = MatNestGetSize_Recursive(A,PETSC_TRUE,&M,&N);CHKERRQ(ierr); 1102 ierr = MatNestGetSize_Recursive(A,PETSC_FALSE,&m,&n);CHKERRQ(ierr); 1103 1104 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1105 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1106 ierr = PetscLayoutSetBlockSize(A->rmap,1);CHKERRQ(ierr); 1107 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1108 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1109 ierr = PetscLayoutSetBlockSize(A->cmap,1);CHKERRQ(ierr); 1110 1111 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1112 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1113 1114 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1115 1116 /* expose Nest api's */ 1117 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMat_C","MatNestGetSubMat_Nest",MatNestGetSubMat_Nest);CHKERRQ(ierr); 1118 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMats_C","MatNestGetSubMats_Nest",MatNestGetSubMats_Nest);CHKERRQ(ierr); 1119 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSize_C","MatNestGetSize_Nest",MatNestGetSize_Nest);CHKERRQ(ierr); 1120 1121 *B = A; 1122 PetscFunctionReturn(0); 1123 } 1124 1125