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