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__ "MatDestroy_Nest" 240 PetscErrorCode MatDestroy_Nest(Mat A) 241 { 242 Mat_Nest *vs = (Mat_Nest*)A->data; 243 PetscInt i,j; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 /* release the matrices and the place holders */ 248 if (vs->isglobal.row) { 249 for (i=0; i<vs->nr; i++) { 250 if(vs->isglobal.row[i]) ierr = ISDestroy(vs->isglobal.row[i]);CHKERRQ(ierr); 251 } 252 ierr = PetscFree(vs->isglobal.row);CHKERRQ(ierr); 253 } 254 if (vs->isglobal.col) { 255 for (j=0; j<vs->nc; j++) { 256 if(vs->isglobal.col[j]) ierr = ISDestroy(vs->isglobal.col[j]);CHKERRQ(ierr); 257 } 258 ierr = PetscFree(vs->isglobal.col);CHKERRQ(ierr); 259 } 260 261 ierr = PetscFree(vs->row_len);CHKERRQ(ierr); 262 ierr = PetscFree(vs->col_len);CHKERRQ(ierr); 263 264 /* release the matrices and the place holders */ 265 if (vs->m) { 266 for (i=0; i<vs->nr; i++) { 267 for (j=0; j<vs->nc; j++) { 268 269 if (vs->m[i][j]) { 270 ierr = MatDestroy(vs->m[i][j]);CHKERRQ(ierr); 271 vs->m[i][j] = PETSC_NULL; 272 } 273 } 274 ierr = PetscFree( vs->m[i] );CHKERRQ(ierr); 275 vs->m[i] = PETSC_NULL; 276 } 277 ierr = PetscFree(vs->m);CHKERRQ(ierr); 278 vs->m = PETSC_NULL; 279 } 280 ierr = PetscFree(vs);CHKERRQ(ierr); 281 282 PetscFunctionReturn(0); 283 } 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "MatAssemblyBegin_Nest" 287 PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type) 288 { 289 Mat_Nest *vs = (Mat_Nest*)A->data; 290 PetscInt i,j; 291 PetscErrorCode ierr; 292 293 PetscFunctionBegin; 294 for (i=0; i<vs->nr; i++) { 295 for (j=0; j<vs->nc; j++) { 296 if (vs->m[i][j]) { ierr = MatAssemblyBegin(vs->m[i][j],type);CHKERRQ(ierr); } 297 } 298 } 299 PetscFunctionReturn(0); 300 } 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "MatAssemblyEnd_Nest" 304 PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type) 305 { 306 Mat_Nest *vs = (Mat_Nest*)A->data; 307 PetscInt i,j; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 for (i=0; i<vs->nr; i++) { 312 for (j=0; j<vs->nc; j++) { 313 if (vs->m[i][j]) { ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr); } 314 } 315 } 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "MatNestFindNonzeroSubMatRow" 321 static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B) 322 { 323 Mat_Nest *vs = (Mat_Nest*)A->data; 324 PetscInt j; 325 Mat sub; 326 327 PetscFunctionBegin; 328 sub = vs->m[row][row]; /* Prefer to find on the diagonal */ 329 for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j]; 330 if (!sub) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",row); 331 *B = sub; 332 PetscFunctionReturn(0); 333 } 334 335 #undef __FUNCT__ 336 #define __FUNCT__ "MatNestFindNonzeroSubMatCol" 337 static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B) 338 { 339 Mat_Nest *vs = (Mat_Nest*)A->data; 340 PetscInt i; 341 Mat sub; 342 343 PetscFunctionBegin; 344 sub = vs->m[col][col]; /* Prefer to find on the diagonal */ 345 for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col]; 346 if (!sub) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",col); 347 *B = sub; 348 PetscFunctionReturn(0); 349 } 350 351 #undef __FUNCT__ 352 #define __FUNCT__ "MatNestFindIS" 353 static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found) 354 { 355 PetscErrorCode ierr; 356 PetscInt i; 357 PetscBool flg; 358 359 PetscFunctionBegin; 360 PetscValidPointer(list,3); 361 PetscValidHeaderSpecific(is,IS_CLASSID,4); 362 PetscValidIntPointer(found,5); 363 *found = -1; 364 for (i=0; i<n; i++) { 365 ierr = ISEqual(list[i],is,&flg);CHKERRQ(ierr); 366 if (flg) { 367 *found = i; 368 PetscFunctionReturn(0); 369 } 370 } 371 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_INCOMP,"Could not find index set"); 372 PetscFunctionReturn(0); 373 } 374 375 #undef __FUNCT__ 376 #define __FUNCT__ "MatNestFindSubMat" 377 static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B) 378 { 379 PetscErrorCode ierr; 380 Mat_Nest *vs = (Mat_Nest*)A->data; 381 PetscInt row,col; 382 383 PetscFunctionBegin; 384 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr); 385 ierr = MatNestFindIS(A,vs->nc,is->col,iscol,&col);CHKERRQ(ierr); 386 *B = vs->m[row][col]; 387 PetscFunctionReturn(0); 388 } 389 390 #undef __FUNCT__ 391 #define __FUNCT__ "MatGetSubMatrix_Nest" 392 static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B) 393 { 394 PetscErrorCode ierr; 395 Mat_Nest *vs = (Mat_Nest*)A->data; 396 Mat sub; 397 398 PetscFunctionBegin; 399 ierr = MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);CHKERRQ(ierr); 400 switch (reuse) { 401 case MAT_INITIAL_MATRIX: 402 ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr); 403 *B = sub; 404 break; 405 case MAT_REUSE_MATRIX: 406 if (sub != *B) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call"); 407 break; 408 case MAT_IGNORE_MATRIX: /* Nothing to do */ 409 break; 410 } 411 PetscFunctionReturn(0); 412 } 413 414 #undef __FUNCT__ 415 #define __FUNCT__ "MatGetLocalSubMatrix_Nest" 416 PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 417 { 418 PetscErrorCode ierr; 419 Mat_Nest *vs = (Mat_Nest*)A->data; 420 Mat sub; 421 422 PetscFunctionBegin; 423 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 424 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */ 425 if (sub) {ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr);} 426 *B = sub; 427 PetscFunctionReturn(0); 428 } 429 430 #undef __FUNCT__ 431 #define __FUNCT__ "MatRestoreLocalSubMatrix_Nest" 432 PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 433 { 434 PetscErrorCode ierr; 435 Mat_Nest *vs = (Mat_Nest*)A->data; 436 Mat sub; 437 438 PetscFunctionBegin; 439 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 440 if (*B != sub) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten"); 441 if (sub) { 442 if (((PetscObject)sub)->refct <= 1) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times"); 443 ierr = MatDestroy(*B);CHKERRQ(ierr); 444 } 445 *B = PETSC_NULL; 446 PetscFunctionReturn(0); 447 } 448 449 #undef __FUNCT__ 450 #define __FUNCT__ "MatGetVecs_Nest" 451 PetscErrorCode MatGetVecs_Nest(Mat A,Vec *right,Vec *left) 452 { 453 Mat_Nest *bA = (Mat_Nest*)A->data; 454 Vec *L,*R; 455 MPI_Comm comm; 456 PetscInt i,j; 457 PetscErrorCode ierr; 458 459 PetscFunctionBegin; 460 comm = ((PetscObject)A)->comm; 461 if (right) { 462 /* allocate R */ 463 ierr = PetscMalloc( sizeof(Vec) * bA->nc, &R );CHKERRQ(ierr); 464 /* Create the right vectors */ 465 for (j=0; j<bA->nc; j++) { 466 for (i=0; i<bA->nr; i++) { 467 if (bA->m[i][j]) { 468 ierr = MatGetVecs(bA->m[i][j],&R[j],PETSC_NULL);CHKERRQ(ierr); 469 break; 470 } 471 } 472 if (i==bA->nr) { 473 /* have an empty column */ 474 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 475 } 476 } 477 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr); 478 /* hand back control to the nest vector */ 479 for (j=0; j<bA->nc; j++) { 480 ierr = VecDestroy(R[j]);CHKERRQ(ierr); 481 } 482 ierr = PetscFree(R);CHKERRQ(ierr); 483 } 484 485 if (left) { 486 /* allocate L */ 487 ierr = PetscMalloc( sizeof(Vec) * bA->nr, &L );CHKERRQ(ierr); 488 /* Create the left vectors */ 489 for (i=0; i<bA->nr; i++) { 490 for (j=0; j<bA->nc; j++) { 491 if (bA->m[i][j]) { 492 ierr = MatGetVecs(bA->m[i][j],PETSC_NULL,&L[i]);CHKERRQ(ierr); 493 break; 494 } 495 } 496 if (j==bA->nc) { 497 /* have an empty row */ 498 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 499 } 500 } 501 502 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr); 503 for (i=0; i<bA->nr; i++) { 504 ierr = VecDestroy(L[i]);CHKERRQ(ierr); 505 } 506 507 ierr = PetscFree(L);CHKERRQ(ierr); 508 } 509 PetscFunctionReturn(0); 510 } 511 512 #undef __FUNCT__ 513 #define __FUNCT__ "MatView_Nest" 514 PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer) 515 { 516 Mat_Nest *bA = (Mat_Nest*)A->data; 517 PetscBool isascii; 518 PetscInt i,j; 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 523 if (isascii) { 524 525 PetscViewerASCIIPrintf(viewer,"Matrix object: \n" ); 526 PetscViewerASCIIPushTab( viewer ); /* push0 */ 527 PetscViewerASCIIPrintf( viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc); 528 529 PetscViewerASCIIPrintf(viewer,"MatNest structure: \n" ); 530 for (i=0; i<bA->nr; i++) { 531 for (j=0; j<bA->nc; j++) { 532 const MatType type; 533 const char *name; 534 PetscInt NR,NC; 535 PetscBool isNest = PETSC_FALSE; 536 537 if (!bA->m[i][j]) { 538 PetscViewerASCIIPrintf( viewer, "(%D,%D) : PETSC_NULL \n",i,j); 539 continue; 540 } 541 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr); 542 ierr = MatGetType( bA->m[i][j], &type );CHKERRQ(ierr); 543 name = ((PetscObject)bA->m[i][j])->prefix; 544 ierr = PetscTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr); 545 546 PetscViewerASCIIPrintf( viewer, "(%D,%D) : name=\"%s\", type=%s, rows=%D, cols=%D \n",i,j,name,type,NR,NC); 547 548 if (isNest) { 549 PetscViewerASCIIPushTab( viewer ); /* push1 */ 550 ierr = MatView( bA->m[i][j], viewer );CHKERRQ(ierr); 551 PetscViewerASCIIPopTab(viewer); /* pop1 */ 552 } 553 } 554 } 555 PetscViewerASCIIPopTab(viewer); /* pop0 */ 556 } 557 PetscFunctionReturn(0); 558 } 559 560 #undef __FUNCT__ 561 #define __FUNCT__ "MatZeroEntries_Nest" 562 PetscErrorCode MatZeroEntries_Nest(Mat A) 563 { 564 Mat_Nest *bA = (Mat_Nest*)A->data; 565 PetscInt i,j; 566 PetscErrorCode ierr; 567 568 PetscFunctionBegin; 569 for (i=0; i<bA->nr; i++) { 570 for (j=0; j<bA->nc; j++) { 571 if (!bA->m[i][j]) continue; 572 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr); 573 } 574 } 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "MatDuplicate_Nest" 580 PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 581 { 582 Mat_Nest *bA = (Mat_Nest*)A->data; 583 Mat *b; 584 PetscInt i,j,nr = bA->nr,nc = bA->nc; 585 PetscErrorCode ierr; 586 587 PetscFunctionBegin; 588 ierr = PetscMalloc(nr*nc*sizeof(Mat),&b);CHKERRQ(ierr); 589 for (i=0; i<nr; i++) { 590 for (j=0; j<nc; j++) { 591 if (bA->m[i][j]) { 592 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 593 } else { 594 b[i*nc+j] = PETSC_NULL; 595 } 596 } 597 } 598 ierr = MatCreateNest(((PetscObject)A)->comm,nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 599 /* Give the new MatNest exclusive ownership */ 600 for (i=0; i<nr*nc; i++) { 601 if (b[i]) {ierr = MatDestroy(b[i]);CHKERRQ(ierr);} 602 } 603 ierr = PetscFree(b);CHKERRQ(ierr); 604 605 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 606 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 607 PetscFunctionReturn(0); 608 } 609 610 /* nest api */ 611 EXTERN_C_BEGIN 612 #undef __FUNCT__ 613 #define __FUNCT__ "MatNestGetSubMat_Nest" 614 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 615 { 616 Mat_Nest *bA = (Mat_Nest*)A->data; 617 PetscFunctionBegin; 618 if (idxm >= bA->nr) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 619 if (jdxm >= bA->nc) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 620 *mat = bA->m[idxm][jdxm]; 621 PetscFunctionReturn(0); 622 } 623 EXTERN_C_END 624 625 #undef __FUNCT__ 626 #define __FUNCT__ "MatNestGetSubMat" 627 /*@C 628 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 629 630 Not collective 631 632 Input Parameters: 633 . A - nest matrix 634 . idxm - index of the matrix within the nest matrix 635 . jdxm - index of the matrix within the nest matrix 636 637 Output Parameter: 638 . sub - matrix at index idxm,jdxm within the nest matrix 639 640 Notes: 641 642 Level: developer 643 644 .seealso: MatNestGetSize(), MatNestGetSubMats() 645 @*/ 646 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 647 { 648 PetscErrorCode ierr; 649 650 PetscFunctionBegin; 651 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 EXTERN_C_BEGIN 656 #undef __FUNCT__ 657 #define __FUNCT__ "MatNestGetSubMats_Nest" 658 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 659 { 660 Mat_Nest *bA = (Mat_Nest*)A->data; 661 PetscFunctionBegin; 662 if (M) { *M = bA->nr; } 663 if (N) { *N = bA->nc; } 664 if (mat) { *mat = bA->m; } 665 PetscFunctionReturn(0); 666 } 667 EXTERN_C_END 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "MatNestGetSubMats" 671 /*@C 672 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 673 674 Not collective 675 676 Input Parameters: 677 . A - nest matri 678 679 Output Parameter: 680 . M - number of rows in the nest matrix 681 . N - number of cols in the nest matrix 682 . mat - 2d array of matrices 683 684 Notes: 685 686 The user should not free the array mat. 687 688 Level: developer 689 690 .seealso: MatNestGetSize(), MatNestGetSubMat() 691 @*/ 692 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 693 { 694 PetscErrorCode ierr; 695 696 PetscFunctionBegin; 697 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 698 PetscFunctionReturn(0); 699 } 700 701 EXTERN_C_BEGIN 702 #undef __FUNCT__ 703 #define __FUNCT__ "MatNestGetSize_Nest" 704 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 705 { 706 Mat_Nest *bA = (Mat_Nest*)A->data; 707 708 PetscFunctionBegin; 709 if (M) { *M = bA->nr; } 710 if (N) { *N = bA->nc; } 711 PetscFunctionReturn(0); 712 } 713 EXTERN_C_END 714 715 #undef __FUNCT__ 716 #define __FUNCT__ "MatNestGetSize" 717 /*@C 718 MatNestGetSize - Returns the size of the nest matrix. 719 720 Not collective 721 722 Input Parameters: 723 . A - nest matrix 724 725 Output Parameter: 726 . M - number of rows in the nested mat 727 . N - number of cols in the nested mat 728 729 Notes: 730 731 Level: developer 732 733 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 734 @*/ 735 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 736 { 737 PetscErrorCode ierr; 738 739 PetscFunctionBegin; 740 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 741 PetscFunctionReturn(0); 742 } 743 744 /* constructor */ 745 #undef __FUNCT__ 746 #define __FUNCT__ "MatNestSetOps_Private" 747 static PetscErrorCode MatNestSetOps_Private(struct _MatOps* ops) 748 { 749 PetscFunctionBegin; 750 /* 0*/ 751 ops->setvalues = 0; 752 ops->getrow = 0; 753 ops->restorerow = 0; 754 ops->mult = MatMult_Nest; 755 ops->multadd = 0; 756 /* 5*/ 757 ops->multtranspose = MatMultTranspose_Nest; 758 ops->multtransposeadd = 0; 759 ops->solve = 0; 760 ops->solveadd = 0; 761 ops->solvetranspose = 0; 762 /*10*/ 763 ops->solvetransposeadd = 0; 764 ops->lufactor = 0; 765 ops->choleskyfactor = 0; 766 ops->sor = 0; 767 ops->transpose = 0; 768 /*15*/ 769 ops->getinfo = 0; 770 ops->equal = 0; 771 ops->getdiagonal = 0; 772 ops->diagonalscale = 0; 773 ops->norm = 0; 774 /*20*/ 775 ops->assemblybegin = MatAssemblyBegin_Nest; 776 ops->assemblyend = MatAssemblyEnd_Nest; 777 ops->setoption = 0; 778 ops->zeroentries = MatZeroEntries_Nest; 779 /*24*/ 780 ops->zerorows = 0; 781 ops->lufactorsymbolic = 0; 782 ops->lufactornumeric = 0; 783 ops->choleskyfactorsymbolic = 0; 784 ops->choleskyfactornumeric = 0; 785 /*29*/ 786 ops->setuppreallocation = 0; 787 ops->ilufactorsymbolic = 0; 788 ops->iccfactorsymbolic = 0; 789 ops->getarray = 0; 790 ops->restorearray = 0; 791 /*34*/ 792 ops->duplicate = MatDuplicate_Nest; 793 ops->forwardsolve = 0; 794 ops->backwardsolve = 0; 795 ops->ilufactor = 0; 796 ops->iccfactor = 0; 797 /*39*/ 798 ops->axpy = 0; 799 ops->getsubmatrices = 0; 800 ops->increaseoverlap = 0; 801 ops->getvalues = 0; 802 ops->copy = 0; 803 /*44*/ 804 ops->getrowmax = 0; 805 ops->scale = 0; 806 ops->shift = 0; 807 ops->diagonalset = 0; 808 ops->zerorowscolumns = 0; 809 /*49*/ 810 ops->setblocksize = 0; 811 ops->getrowij = 0; 812 ops->restorerowij = 0; 813 ops->getcolumnij = 0; 814 ops->restorecolumnij = 0; 815 /*54*/ 816 ops->fdcoloringcreate = 0; 817 ops->coloringpatch = 0; 818 ops->setunfactored = 0; 819 ops->permute = 0; 820 ops->setvaluesblocked = 0; 821 /*59*/ 822 ops->getsubmatrix = MatGetSubMatrix_Nest; 823 ops->destroy = MatDestroy_Nest; 824 ops->view = MatView_Nest; 825 ops->convertfrom = 0; 826 ops->usescaledform = 0; 827 /*64*/ 828 ops->scalesystem = 0; 829 ops->unscalesystem = 0; 830 ops->setlocaltoglobalmapping = 0; 831 ops->setvalueslocal = 0; 832 ops->zerorowslocal = 0; 833 /*69*/ 834 ops->getrowmaxabs = 0; 835 ops->getrowminabs = 0; 836 ops->convert = 0;/*MatConvert_Nest;*/ 837 ops->setcoloring = 0; 838 ops->setvaluesadic = 0; 839 /* 74 */ 840 ops->setvaluesadifor = 0; 841 ops->fdcoloringapply = 0; 842 ops->setfromoptions = 0; 843 ops->multconstrained = 0; 844 ops->multtransposeconstrained = 0; 845 /*79*/ 846 ops->permutesparsify = 0; 847 ops->mults = 0; 848 ops->solves = 0; 849 ops->getinertia = 0; 850 ops->load = 0; 851 /*84*/ 852 ops->issymmetric = 0; 853 ops->ishermitian = 0; 854 ops->isstructurallysymmetric = 0; 855 ops->dummy4 = 0; 856 ops->getvecs = MatGetVecs_Nest; 857 /*89*/ 858 ops->matmult = 0;/*MatMatMult_Nest;*/ 859 ops->matmultsymbolic = 0; 860 ops->matmultnumeric = 0; 861 ops->ptap = 0; 862 ops->ptapsymbolic = 0; 863 /*94*/ 864 ops->ptapnumeric = 0; 865 ops->matmulttranspose = 0; 866 ops->matmulttransposesymbolic = 0; 867 ops->matmulttransposenumeric = 0; 868 ops->ptapsymbolic_seqaij = 0; 869 /*99*/ 870 ops->ptapnumeric_seqaij = 0; 871 ops->ptapsymbolic_mpiaij = 0; 872 ops->ptapnumeric_mpiaij = 0; 873 ops->conjugate = 0; 874 ops->setsizes = 0; 875 /*104*/ 876 ops->setvaluesrow = 0; 877 ops->realpart = 0; 878 ops->imaginarypart = 0; 879 ops->getrowuppertriangular = 0; 880 ops->restorerowuppertriangular = 0; 881 /*109*/ 882 ops->matsolve = 0; 883 ops->getredundantmatrix = 0; 884 ops->getrowmin = 0; 885 ops->getcolumnvector = 0; 886 ops->missingdiagonal = 0; 887 /* 114 */ 888 ops->getseqnonzerostructure = 0; 889 ops->create = 0; 890 ops->getghosts = 0; 891 ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 892 ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 893 /* 119 */ 894 ops->multdiagonalblock = 0; 895 ops->hermitiantranspose = 0; 896 ops->multhermitiantranspose = 0; 897 ops->multhermitiantransposeadd = 0; 898 ops->getmultiprocblock = 0; 899 /* 124 */ 900 ops->dummy1 = 0; 901 ops->dummy2 = 0; 902 ops->dummy3 = 0; 903 ops->dummy4 = 0; 904 ops->getsubmatricesparallel = 0; 905 906 PetscFunctionReturn(0); 907 } 908 909 #undef __FUNCT__ 910 #define __FUNCT__ "MatSetUp_Nest_Private" 911 static PetscErrorCode MatSetUp_Nest_Private(Mat A,PetscInt nr,PetscInt nc,const Mat *sub) 912 { 913 Mat_Nest *ctx = (Mat_Nest*)A->data; 914 PetscInt i,j; 915 PetscErrorCode ierr; 916 917 PetscFunctionBegin; 918 if(ctx->setup_called) PetscFunctionReturn(0); 919 920 ctx->nr = nr; 921 ctx->nc = nc; 922 923 /* Create space */ 924 ierr = PetscMalloc(sizeof(Mat*)*ctx->nr,&ctx->m);CHKERRQ(ierr); 925 for (i=0; i<ctx->nr; i++) { 926 ierr = PetscMalloc(sizeof(Mat)*ctx->nc,&ctx->m[i]);CHKERRQ(ierr); 927 } 928 for (i=0; i<ctx->nr; i++) { 929 for (j=0; j<ctx->nc; j++) { 930 ctx->m[i][j] = sub[i*nc+j]; 931 if (sub[i*nc+j]) { 932 ierr = PetscObjectReference((PetscObject)sub[i*nc+j]);CHKERRQ(ierr); 933 } 934 } 935 } 936 937 ierr = PetscMalloc(sizeof(IS)*ctx->nr,&ctx->isglobal.row);CHKERRQ(ierr); 938 ierr = PetscMalloc(sizeof(IS)*ctx->nc,&ctx->isglobal.col);CHKERRQ(ierr); 939 940 ierr = PetscMalloc(sizeof(PetscInt)*ctx->nr,&ctx->row_len);CHKERRQ(ierr); 941 ierr = PetscMalloc(sizeof(PetscInt)*ctx->nc,&ctx->col_len);CHKERRQ(ierr); 942 for (i=0;i<ctx->nr;i++) ctx->row_len[i]=-1; 943 for (j=0;j<ctx->nc;j++) ctx->col_len[j]=-1; 944 945 ctx->setup_called = PETSC_TRUE; 946 947 PetscFunctionReturn(0); 948 } 949 950 951 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 952 /* 953 nprocessors = NP 954 Nest x^T = ( (g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1) ) 955 proc 0: => (g_0,h_0,) 956 proc 1: => (g_1,h_1,) 957 ... 958 proc nprocs-1: => (g_NP-1,h_NP-1,) 959 960 proc 0: proc 1: proc nprocs-1: 961 is[0] = ( 0,1,2,...,nlocal(g_0)-1 ) ( 0,1,...,nlocal(g_1)-1 ) ( 0,1,...,nlocal(g_NP-1) ) 962 963 proc 0: 964 is[1] = ( nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1 ) 965 proc 1: 966 is[1] = ( nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1 ) 967 968 proc NP-1: 969 is[1] = ( nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1 ) 970 */ 971 #undef __FUNCT__ 972 #define __FUNCT__ "MatSetUp_NestIS_Private" 973 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 974 { 975 Mat_Nest *ctx = (Mat_Nest*)A->data; 976 PetscInt i,j,offset,n,bs; 977 PetscErrorCode ierr; 978 Mat sub; 979 980 PetscFunctionBegin; 981 if (is_row) { /* valid IS is passed in */ 982 /* refs on is[] are incremeneted */ 983 for (i=0; i<ctx->nr; i++) { 984 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 985 ctx->isglobal.row[i] = is_row[i]; 986 } 987 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 988 offset = A->rmap->rstart; 989 for (i=0; i<ctx->nr; i++) { 990 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 991 ierr = MatGetLocalSize(sub,&n,PETSC_NULL);CHKERRQ(ierr); 992 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 993 ierr = ISCreateStride(((PetscObject)sub)->comm,n,offset,1,&ctx->isglobal.row[i]);CHKERRQ(ierr); 994 ierr = ISSetBlockSize(ctx->isglobal.row[i],bs);CHKERRQ(ierr); 995 offset += n; 996 } 997 } 998 999 if (is_col) { /* valid IS is passed in */ 1000 /* refs on is[] are incremeneted */ 1001 for (j=0; j<ctx->nc; j++) { 1002 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1003 ctx->isglobal.col[j] = is_col[j]; 1004 } 1005 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1006 offset = A->cmap->rstart; 1007 for (j=0; j<ctx->nc; j++) { 1008 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1009 ierr = MatGetLocalSize(sub,PETSC_NULL,&n);CHKERRQ(ierr); 1010 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1011 ierr = ISCreateStride(((PetscObject)sub)->comm,n,offset,1,&ctx->isglobal.col[j]);CHKERRQ(ierr); 1012 ierr = ISSetBlockSize(ctx->isglobal.col[j],bs);CHKERRQ(ierr); 1013 offset += n; 1014 } 1015 } 1016 PetscFunctionReturn(0); 1017 } 1018 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "MatCreateNest" 1021 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1022 { 1023 Mat A; 1024 Mat_Nest *s; 1025 PetscInt i,m,n,M,N; 1026 PetscErrorCode ierr; 1027 1028 PetscFunctionBegin; 1029 if (nr < 0) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1030 if (nr && is_row) { 1031 PetscValidPointer(is_row,3); 1032 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1033 } 1034 if (nc < 0) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1035 if (nc && is_row) { 1036 PetscValidPointer(is_col,5); 1037 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1038 } 1039 if (nr*nc) PetscValidPointer(a,6); 1040 PetscValidPointer(B,7); 1041 1042 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1043 1044 ierr = PetscMalloc( sizeof(Mat_Nest), &s );CHKERRQ(ierr); 1045 ierr = PetscMemzero( s, sizeof(Mat_Nest) );CHKERRQ(ierr); 1046 A->data = (void*)s; 1047 s->setup_called = PETSC_FALSE; 1048 s->nr = s->nc = -1; 1049 s->m = PETSC_NULL; 1050 1051 /* define the operations */ 1052 ierr = MatNestSetOps_Private(A->ops);CHKERRQ(ierr); 1053 A->spptr = 0; 1054 A->same_nonzero = PETSC_FALSE; 1055 A->assembled = PETSC_FALSE; 1056 1057 ierr = PetscObjectChangeTypeName((PetscObject) A, MATNEST );CHKERRQ(ierr); 1058 ierr = MatSetUp_Nest_Private(A,nr,nc,a);CHKERRQ(ierr); 1059 1060 m = n = 0; 1061 M = N = 0; 1062 ierr = MatNestGetSize_Recursive(A,PETSC_TRUE,&M,&N);CHKERRQ(ierr); 1063 ierr = MatNestGetSize_Recursive(A,PETSC_FALSE,&m,&n);CHKERRQ(ierr); 1064 1065 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1066 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1067 ierr = PetscLayoutSetBlockSize(A->rmap,1);CHKERRQ(ierr); 1068 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1069 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1070 ierr = PetscLayoutSetBlockSize(A->cmap,1);CHKERRQ(ierr); 1071 1072 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1073 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1074 1075 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1076 1077 /* expose Nest api's */ 1078 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMat_C","MatNestGetSubMat_Nest",MatNestGetSubMat_Nest);CHKERRQ(ierr); 1079 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMats_C","MatNestGetSubMats_Nest",MatNestGetSubMats_Nest);CHKERRQ(ierr); 1080 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSize_C","MatNestGetSize_Nest",MatNestGetSize_Nest);CHKERRQ(ierr); 1081 1082 *B = A; 1083 PetscFunctionReturn(0); 1084 } 1085 1086