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; 526 PetscErrorCode ierr; 527 528 // TODO 529 PetscFunctionBegin; 530 531 ierr = PetscMalloc(sizeof(Mat*)*bA->nr,&b);CHKERRQ(ierr); 532 for (i=0; i<bA->nr; i++) { 533 ierr = PetscMalloc(sizeof(Mat)*bA->nc,&b[i]);CHKERRQ(ierr); 534 } 535 for (i=0; i<bA->nr; i++) { 536 for (j=0; j<bA->nc; j++) { 537 if (!bA->m[i][j]) continue; 538 ierr = MatDuplicate(bA->m[i][j],op,&b[i][j]);CHKERRQ(ierr); 539 } 540 } 541 ierr = MatCreateNest(((PetscObject)A)->comm,bA->nr,bA->nc,bA->is_row,bA->is_col,bA->m,B);CHKERRQ(ierr); 542 /* hand back control to the nest */ 543 for (i=0; i<bA->nr; i++) { 544 for (j=0; j<bA->nc; j++) { 545 if (!bA->m[i][j]) continue; 546 ierr = MatDestroy(b[i][j]);CHKERRQ(ierr); 547 } 548 } 549 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 550 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 551 552 for (i=0; i<bA->nr; i++) { 553 ierr = PetscFree(b[i]);CHKERRQ(ierr); 554 } 555 ierr = PetscFree(b);CHKERRQ(ierr); 556 557 PetscFunctionReturn(0); 558 } 559 560 /* nest api */ 561 EXTERN_C_BEGIN 562 #undef __FUNCT__ 563 #define __FUNCT__ "MatNestGetSubMat_Nest" 564 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 565 { 566 Mat_Nest *bA = (Mat_Nest*)A->data; 567 PetscFunctionBegin; 568 if (idxm >= bA->nr) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 569 if (jdxm >= bA->nc) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 570 *mat = bA->m[idxm][jdxm]; 571 PetscFunctionReturn(0); 572 } 573 EXTERN_C_END 574 575 #undef __FUNCT__ 576 #define __FUNCT__ "MatNestGetSubMat" 577 /*@C 578 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 579 580 Not collective 581 582 Input Parameters: 583 . A - nest matrix 584 . idxm - index of the matrix within the nest matrix 585 . jdxm - index of the matrix within the nest matrix 586 587 Output Parameter: 588 . sub - matrix at index idxm,jdxm within the nest matrix 589 590 Notes: 591 592 Level: developer 593 594 .seealso: MatNestGetSize(), MatNestGetSubMats() 595 @*/ 596 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 597 { 598 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,Mat*); 599 600 PetscFunctionBegin; 601 ierr = PetscObjectQueryFunction((PetscObject)A,"MatNestGetSubMat_C",(void (**)(void))&f);CHKERRQ(ierr); 602 if (f) { 603 ierr = (*f)(A,idxm,jdxm,sub);CHKERRQ(ierr); 604 } 605 PetscFunctionReturn(0); 606 } 607 608 EXTERN_C_BEGIN 609 #undef __FUNCT__ 610 #define __FUNCT__ "MatNestGetSubMats_Nest" 611 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 612 { 613 Mat_Nest *bA = (Mat_Nest*)A->data; 614 PetscFunctionBegin; 615 if (M) { *M = bA->nr; } 616 if (N) { *N = bA->nc; } 617 if (mat) { *mat = bA->m; } 618 PetscFunctionReturn(0); 619 } 620 EXTERN_C_END 621 622 #undef __FUNCT__ 623 #define __FUNCT__ "MatNestGetSubMats" 624 /*@C 625 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 626 627 Not collective 628 629 Input Parameters: 630 . A - nest matri 631 632 Output Parameter: 633 . M - number of rows in the nest matrix 634 . N - number of cols in the nest matrix 635 . mat - 2d array of matrices 636 637 Notes: 638 639 The user should not free the array mat. 640 641 Level: developer 642 643 .seealso: MatNestGetSize(), MatNestGetSubMat() 644 @*/ 645 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 646 { 647 PetscErrorCode ierr,(*f)(Mat,PetscInt*,PetscInt*,Mat***); 648 649 PetscFunctionBegin; 650 ierr = PetscObjectQueryFunction((PetscObject)A,"MatNestGetSubMats_C",(void (**)(void))&f);CHKERRQ(ierr); 651 if (f) { 652 ierr = (*f)(A,M,N,mat);CHKERRQ(ierr); 653 } 654 PetscFunctionReturn(0); 655 } 656 657 EXTERN_C_BEGIN 658 #undef __FUNCT__ 659 #define __FUNCT__ "MatNestGetSize_Nest" 660 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 661 { 662 Mat_Nest *bA = (Mat_Nest*)A->data; 663 664 PetscFunctionBegin; 665 if (M) { *M = bA->nr; } 666 if (N) { *N = bA->nc; } 667 PetscFunctionReturn(0); 668 } 669 EXTERN_C_END 670 671 #undef __FUNCT__ 672 #define __FUNCT__ "MatNestGetSize" 673 /*@C 674 MatNestGetSize - Returns the size of the nest matrix. 675 676 Not collective 677 678 Input Parameters: 679 . A - nest matrix 680 681 Output Parameter: 682 . M - number of rows in the nested mat 683 . N - number of cols in the nested mat 684 685 Notes: 686 687 Level: developer 688 689 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 690 @*/ 691 PetscErrorCode PETSCMAT_DLLEXPORT MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 692 { 693 PetscErrorCode ierr,(*f)(Mat,PetscInt*,PetscInt*); 694 695 PetscFunctionBegin; 696 ierr = PetscObjectQueryFunction((PetscObject)A,"MatNestGetSize_C",(void (**)(void))&f);CHKERRQ(ierr); 697 if (f) { 698 ierr = (*f)(A,M,N);CHKERRQ(ierr); 699 } 700 PetscFunctionReturn(0); 701 } 702 703 /* constructor */ 704 #undef __FUNCT__ 705 #define __FUNCT__ "MatNestSetOps_Private" 706 static PetscErrorCode MatNestSetOps_Private(struct _MatOps* ops) 707 { 708 PetscFunctionBegin; 709 /* 0*/ 710 ops->setvalues = 0; 711 ops->getrow = 0; 712 ops->restorerow = 0; 713 ops->mult = MatMult_Nest; 714 ops->multadd = 0; 715 /* 5*/ 716 ops->multtranspose = MatMultTranspose_Nest; 717 ops->multtransposeadd = 0; 718 ops->solve = 0; 719 ops->solveadd = 0; 720 ops->solvetranspose = 0; 721 /*10*/ 722 ops->solvetransposeadd = 0; 723 ops->lufactor = 0; 724 ops->choleskyfactor = 0; 725 ops->sor = 0; 726 ops->transpose = 0; 727 /*15*/ 728 ops->getinfo = 0; 729 ops->equal = 0; 730 ops->getdiagonal = 0; 731 ops->diagonalscale = 0; 732 ops->norm = 0; 733 /*20*/ 734 ops->assemblybegin = MatAssemblyBegin_Nest; 735 ops->assemblyend = MatAssemblyEnd_Nest; 736 ops->setoption = 0; 737 ops->zeroentries = MatZeroEntries_Nest; 738 /*24*/ 739 ops->zerorows = 0; 740 ops->lufactorsymbolic = 0; 741 ops->lufactornumeric = 0; 742 ops->choleskyfactorsymbolic = 0; 743 ops->choleskyfactornumeric = 0; 744 /*29*/ 745 ops->setuppreallocation = 0; 746 ops->ilufactorsymbolic = 0; 747 ops->iccfactorsymbolic = 0; 748 ops->getarray = 0; 749 ops->restorearray = 0; 750 /*34*/ 751 ops->duplicate = MatDuplicate_Nest; 752 ops->forwardsolve = 0; 753 ops->backwardsolve = 0; 754 ops->ilufactor = 0; 755 ops->iccfactor = 0; 756 /*39*/ 757 ops->axpy = 0; 758 ops->getsubmatrices = 0; 759 ops->increaseoverlap = 0; 760 ops->getvalues = 0; 761 ops->copy = 0; 762 /*44*/ 763 ops->getrowmax = 0; 764 ops->scale = 0; 765 ops->shift = 0; 766 ops->diagonalset = 0; 767 ops->zerorowscolumns = 0; 768 /*49*/ 769 ops->setblocksize = 0; 770 ops->getrowij = 0; 771 ops->restorerowij = 0; 772 ops->getcolumnij = 0; 773 ops->restorecolumnij = 0; 774 /*54*/ 775 ops->fdcoloringcreate = 0; 776 ops->coloringpatch = 0; 777 ops->setunfactored = 0; 778 ops->permute = 0; 779 ops->setvaluesblocked = 0; 780 /*59*/ 781 ops->getsubmatrix = 0; 782 ops->destroy = MatDestroy_Nest; 783 ops->view = MatView_Nest; 784 ops->convertfrom = 0; 785 ops->usescaledform = 0; 786 /*64*/ 787 ops->scalesystem = 0; 788 ops->unscalesystem = 0; 789 ops->setlocaltoglobalmapping = 0; 790 ops->setvalueslocal = 0; 791 ops->zerorowslocal = 0; 792 /*69*/ 793 ops->getrowmaxabs = 0; 794 ops->getrowminabs = 0; 795 ops->convert = 0;/*MatConvert_Nest;*/ 796 ops->setcoloring = 0; 797 ops->setvaluesadic = 0; 798 /* 74 */ 799 ops->setvaluesadifor = 0; 800 ops->fdcoloringapply = 0; 801 ops->setfromoptions = 0; 802 ops->multconstrained = 0; 803 ops->multtransposeconstrained = 0; 804 /*79*/ 805 ops->permutesparsify = 0; 806 ops->mults = 0; 807 ops->solves = 0; 808 ops->getinertia = 0; 809 ops->load = 0; 810 /*84*/ 811 ops->issymmetric = 0; 812 ops->ishermitian = 0; 813 ops->isstructurallysymmetric = 0; 814 ops->dummy4 = 0; 815 ops->getvecs = MatGetVecs_Nest; 816 /*89*/ 817 ops->matmult = 0;/*MatMatMult_Nest;*/ 818 ops->matmultsymbolic = 0; 819 ops->matmultnumeric = 0; 820 ops->ptap = 0; 821 ops->ptapsymbolic = 0; 822 /*94*/ 823 ops->ptapnumeric = 0; 824 ops->matmulttranspose = 0; 825 ops->matmulttransposesymbolic = 0; 826 ops->matmulttransposenumeric = 0; 827 ops->ptapsymbolic_seqaij = 0; 828 /*99*/ 829 ops->ptapnumeric_seqaij = 0; 830 ops->ptapsymbolic_mpiaij = 0; 831 ops->ptapnumeric_mpiaij = 0; 832 ops->conjugate = 0; 833 ops->setsizes = 0; 834 /*104*/ 835 ops->setvaluesrow = 0; 836 ops->realpart = 0; 837 ops->imaginarypart = 0; 838 ops->getrowuppertriangular = 0; 839 ops->restorerowuppertriangular = 0; 840 /*109*/ 841 ops->matsolve = 0; 842 ops->getredundantmatrix = 0; 843 ops->getrowmin = 0; 844 ops->getcolumnvector = 0; 845 ops->missingdiagonal = 0; 846 /* 114 */ 847 ops->getseqnonzerostructure = 0; 848 ops->create = 0; 849 ops->getghosts = 0; 850 ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 851 ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 852 /* 119 */ 853 ops->multdiagonalblock = 0; 854 ops->hermitiantranspose = 0; 855 ops->multhermitiantranspose = 0; 856 ops->multhermitiantransposeadd = 0; 857 ops->getmultiprocblock = 0; 858 /* 124 */ 859 ops->dummy1 = 0; 860 ops->dummy2 = 0; 861 ops->dummy3 = 0; 862 ops->dummy4 = 0; 863 ops->getsubmatricesparallel = 0; 864 865 PetscFunctionReturn(0); 866 } 867 868 #undef __FUNCT__ 869 #define __FUNCT__ "MatSetUp_Nest_Private" 870 static PetscErrorCode MatSetUp_Nest_Private(Mat A,PetscInt nr,PetscInt nc,Mat **sub) 871 { 872 Mat_Nest *ctx = (Mat_Nest*)A->data; 873 PetscInt i,j; 874 PetscErrorCode ierr; 875 876 PetscFunctionBegin; 877 if(ctx->setup_called) PetscFunctionReturn(0); 878 879 ctx->nr = nr; 880 ctx->nc = nc; 881 if (ctx->nr < 0) { 882 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Cannot create MAT_NEST with < 0 row Nests." ); 883 } 884 if (ctx->nc < 0) { 885 SETERRQ( ((PetscObject)A)->comm, PETSC_ERR_ARG_WRONG, "Cannot create MAT_NEST with < 0 col Nests." ); 886 } 887 888 /* Create space */ 889 ierr = PetscMalloc(sizeof(Mat*)*ctx->nr,&ctx->m);CHKERRQ(ierr); 890 for (i=0; i<ctx->nr; i++) { 891 ierr = PetscMalloc(sizeof(Mat)*ctx->nc,&ctx->m[i]);CHKERRQ(ierr); 892 } 893 for (i=0; i<ctx->nr; i++) { 894 for (j=0; j<ctx->nc; j++) { 895 ctx->m[i][j] = sub[i][j]; 896 if (sub[i][j]) { 897 ierr = PetscObjectReference((PetscObject)sub[i][j]);CHKERRQ(ierr); 898 } 899 } 900 } 901 902 ierr = PetscMalloc(sizeof(IS)*ctx->nr,&ctx->is_row);CHKERRQ(ierr); 903 ierr = PetscMalloc(sizeof(IS)*ctx->nc,&ctx->is_col);CHKERRQ(ierr); 904 905 ierr = PetscMalloc(sizeof(PetscInt)*ctx->nr,&ctx->row_len);CHKERRQ(ierr); 906 ierr = PetscMalloc(sizeof(PetscInt)*ctx->nc,&ctx->col_len);CHKERRQ(ierr); 907 for (i=0;i<ctx->nr;i++) ctx->row_len[i]=-1; 908 for (j=0;j<ctx->nc;j++) ctx->col_len[j]=-1; 909 910 ctx->setup_called = PETSC_TRUE; 911 912 PetscFunctionReturn(0); 913 } 914 915 916 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 917 /* 918 nprocessors = NP 919 Nest x^T = ( (g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1) ) 920 proc 0: => (g_0,h_0,) 921 proc 1: => (g_1,h_1,) 922 ... 923 proc nprocs-1: => (g_NP-1,h_NP-1,) 924 925 proc 0: proc 1: proc nprocs-1: 926 is[0] = ( 0,1,2,...,nlocal(g_0)-1 ) ( 0,1,...,nlocal(g_1)-1 ) ( 0,1,...,nlocal(g_NP-1) ) 927 928 proc 0: 929 is[1] = ( nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1 ) 930 proc 1: 931 is[1] = ( nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1 ) 932 933 proc NP-1: 934 is[1] = ( nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1 ) 935 */ 936 #undef __FUNCT__ 937 #define __FUNCT__ "MatSetUp_NestIS_Private" 938 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,PetscInt nc,IS is_row[],IS is_col[]) 939 { 940 Mat_Nest *ctx = (Mat_Nest*)A->data; 941 PetscInt i,j,offset,n,start,index; 942 PetscErrorCode ierr; 943 944 PetscFunctionBegin; 945 if (is_row) { /* valid IS is passed in */ 946 /* refs on is[] are incremeneted */ 947 for (i=0; i<ctx->nr; i++) { 948 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 949 ctx->is_row[i] = is_row[i]; 950 } 951 } else { 952 for (j=0; j<ctx->nc; j++) { 953 if (ctx->m[0][j]) { index = j; break; } 954 } 955 ierr = MatGetOwnershipRange(ctx->m[0][index],&start,PETSC_NULL);CHKERRQ(ierr); 956 offset = start; 957 for (i=0; i<ctx->nr; i++) { 958 959 ierr = MatGetLocalSize(ctx->m[i][index],&n,PETSC_NULL);CHKERRQ(ierr); 960 ierr = ISCreateStride(((PetscObject)ctx->m[i][index])->comm,n,offset,1,&ctx->is_row[i]);CHKERRQ(ierr); 961 offset = offset + n; 962 } 963 } 964 965 if (is_col) { /* valid IS is passed in */ 966 /* refs on is[] are incremeneted */ 967 for (j=0; j<ctx->nc; j++) { 968 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 969 ctx->is_col[j] = is_col[j]; 970 } 971 } else { 972 for (i=0; i<ctx->nr; i++) { 973 if (ctx->m[i][0]) { index = i; break; } 974 } 975 ierr = MatGetOwnershipRange(ctx->m[index][0],&start,PETSC_NULL);CHKERRQ(ierr); 976 offset = start; 977 for (j=0; j<ctx->nc; j++) { 978 for (i=0; i<ctx->nr; i++) { 979 if (ctx->m[i][j]) { index = i; break; } 980 } 981 ierr = MatGetLocalSize(ctx->m[index][j],PETSC_NULL,&n);CHKERRQ(ierr); 982 ierr = ISCreateStride(((PetscObject)ctx->m[index][j])->comm,n,offset,1,&ctx->is_col[j]);CHKERRQ(ierr); 983 offset = offset + n; 984 } 985 } 986 PetscFunctionReturn(0); 987 } 988 989 #undef __FUNCT__ 990 #define __FUNCT__ "MatCreateNest" 991 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNest(MPI_Comm comm,PetscInt nr,PetscInt nc,IS is_row[],IS is_col[],Mat **a,Mat *B) 992 { 993 Mat A; 994 Mat_Nest *s; 995 PetscInt m,n,M,N; 996 PetscErrorCode ierr; 997 998 PetscFunctionBegin; 999 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1000 1001 ierr = PetscMalloc( sizeof(Mat_Nest), &s );CHKERRQ(ierr); 1002 ierr = PetscMemzero( s, sizeof(Mat_Nest) );CHKERRQ(ierr); 1003 A->data = (void*)s; 1004 s->setup_called = PETSC_FALSE; 1005 s->nr = s->nc = -1; 1006 s->m = PETSC_NULL; 1007 1008 /* define the operations */ 1009 ierr = MatNestSetOps_Private(A->ops);CHKERRQ(ierr); 1010 A->spptr = 0; 1011 A->same_nonzero = PETSC_FALSE; 1012 A->assembled = PETSC_FALSE; 1013 1014 ierr = PetscObjectChangeTypeName((PetscObject) A, MATNEST );CHKERRQ(ierr); 1015 ierr = MatSetUp_Nest_Private(A,nr,nc,a);CHKERRQ(ierr); 1016 1017 m = n = 0; 1018 M = N = 0; 1019 ierr = MatNestGetSize_Recursive(A,PETSC_TRUE,&M,&N);CHKERRQ(ierr); 1020 ierr = MatNestGetSize_Recursive(A,PETSC_FALSE,&m,&n);CHKERRQ(ierr); 1021 1022 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1023 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1024 ierr = PetscLayoutSetBlockSize(A->rmap,1);CHKERRQ(ierr); 1025 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1026 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1027 ierr = PetscLayoutSetBlockSize(A->cmap,1);CHKERRQ(ierr); 1028 1029 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1030 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1031 1032 ierr = MatSetUp_NestIS_Private(A,nr,nc,is_row,is_col);CHKERRQ(ierr); 1033 1034 /* expose Nest api's */ 1035 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMat_C","MatNestGetSubMat_Nest",MatNestGetSubMat_Nest);CHKERRQ(ierr); 1036 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSubMats_C","MatNestGetSubMats_Nest",MatNestGetSubMats_Nest);CHKERRQ(ierr); 1037 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatNestGetSize_C","MatNestGetSize_Nest",MatNestGetSize_Nest);CHKERRQ(ierr); 1038 1039 *B = A; 1040 PetscFunctionReturn(0); 1041 } 1042 1043