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