1 2 #include <../src/mat/impls/nest/matnestimpl.h> /*I "petscmat.h" I*/ 3 #include <petscsf.h> 4 5 static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]); 6 static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left); 7 8 /* private functions */ 9 #undef __FUNCT__ 10 #define __FUNCT__ "MatNestGetSizes_Private" 11 static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N) 12 { 13 Mat_Nest *bA = (Mat_Nest*)A->data; 14 PetscInt i,j; 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 *m = *n = *M = *N = 0; 19 for (i=0; i<bA->nr; i++) { /* rows */ 20 PetscInt sm,sM; 21 ierr = ISGetLocalSize(bA->isglobal.row[i],&sm);CHKERRQ(ierr); 22 ierr = ISGetSize(bA->isglobal.row[i],&sM);CHKERRQ(ierr); 23 *m += sm; 24 *M += sM; 25 } 26 for (j=0; j<bA->nc; j++) { /* cols */ 27 PetscInt sn,sN; 28 ierr = ISGetLocalSize(bA->isglobal.col[j],&sn);CHKERRQ(ierr); 29 ierr = ISGetSize(bA->isglobal.col[j],&sN);CHKERRQ(ierr); 30 *n += sn; 31 *N += sN; 32 } 33 PetscFunctionReturn(0); 34 } 35 36 /* operations */ 37 #undef __FUNCT__ 38 #define __FUNCT__ "MatMult_Nest" 39 static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y) 40 { 41 Mat_Nest *bA = (Mat_Nest*)A->data; 42 Vec *bx = bA->right,*by = bA->left; 43 PetscInt i,j,nr = bA->nr,nc = bA->nc; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 for (i=0; i<nr; i++) {ierr = VecGetSubVector(y,bA->isglobal.row[i],&by[i]);CHKERRQ(ierr);} 48 for (i=0; i<nc; i++) {ierr = VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 49 for (i=0; i<nr; i++) { 50 ierr = VecZeroEntries(by[i]);CHKERRQ(ierr); 51 for (j=0; j<nc; j++) { 52 if (!bA->m[i][j]) continue; 53 /* y[i] <- y[i] + A[i][j] * x[j] */ 54 ierr = MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);CHKERRQ(ierr); 55 } 56 } 57 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);CHKERRQ(ierr);} 58 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 59 PetscFunctionReturn(0); 60 } 61 62 #undef __FUNCT__ 63 #define __FUNCT__ "MatMultAdd_Nest" 64 static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z) 65 { 66 Mat_Nest *bA = (Mat_Nest*)A->data; 67 Vec *bx = bA->right,*bz = bA->left; 68 PetscInt i,j,nr = bA->nr,nc = bA->nc; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 for (i=0; i<nr; i++) {ierr = VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);CHKERRQ(ierr);} 73 for (i=0; i<nc; i++) {ierr = VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 74 for (i=0; i<nr; i++) { 75 if (y != z) { 76 Vec by; 77 ierr = VecGetSubVector(y,bA->isglobal.row[i],&by);CHKERRQ(ierr); 78 ierr = VecCopy(by,bz[i]);CHKERRQ(ierr); 79 ierr = VecRestoreSubVector(y,bA->isglobal.row[i],&by);CHKERRQ(ierr); 80 } 81 for (j=0; j<nc; j++) { 82 if (!bA->m[i][j]) continue; 83 /* y[i] <- y[i] + A[i][j] * x[j] */ 84 ierr = MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);CHKERRQ(ierr); 85 } 86 } 87 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);CHKERRQ(ierr);} 88 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);CHKERRQ(ierr);} 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "MatMultTranspose_Nest" 94 static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y) 95 { 96 Mat_Nest *bA = (Mat_Nest*)A->data; 97 Vec *bx = bA->left,*by = bA->right; 98 PetscInt i,j,nr = bA->nr,nc = bA->nc; 99 PetscErrorCode ierr; 100 101 PetscFunctionBegin; 102 for (i=0; i<nr; i++) {ierr = VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 103 for (i=0; i<nc; i++) {ierr = VecGetSubVector(y,bA->isglobal.col[i],&by[i]);CHKERRQ(ierr);} 104 for (j=0; j<nc; j++) { 105 ierr = VecZeroEntries(by[j]);CHKERRQ(ierr); 106 for (i=0; i<nr; i++) { 107 if (!bA->m[i][j]) continue; 108 /* y[j] <- y[j] + (A[i][j])^T * x[i] */ 109 ierr = MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);CHKERRQ(ierr); 110 } 111 } 112 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 113 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);CHKERRQ(ierr);} 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "MatMultTransposeAdd_Nest" 119 static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z) 120 { 121 Mat_Nest *bA = (Mat_Nest*)A->data; 122 Vec *bx = bA->left,*bz = bA->right; 123 PetscInt i,j,nr = bA->nr,nc = bA->nc; 124 PetscErrorCode ierr; 125 126 PetscFunctionBegin; 127 for (i=0; i<nr; i++) {ierr = VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 128 for (i=0; i<nc; i++) {ierr = VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);CHKERRQ(ierr);} 129 for (j=0; j<nc; j++) { 130 if (y != z) { 131 Vec by; 132 ierr = VecGetSubVector(y,bA->isglobal.col[j],&by);CHKERRQ(ierr); 133 ierr = VecCopy(by,bz[j]);CHKERRQ(ierr); 134 ierr = VecRestoreSubVector(y,bA->isglobal.col[j],&by);CHKERRQ(ierr); 135 } 136 for (i=0; i<nr; i++) { 137 if (!bA->m[i][j]) continue; 138 /* z[j] <- y[j] + (A[i][j])^T * x[i] */ 139 ierr = MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);CHKERRQ(ierr); 140 } 141 } 142 for (i=0; i<nr; i++) {ierr = VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);CHKERRQ(ierr);} 143 for (i=0; i<nc; i++) {ierr = VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);CHKERRQ(ierr);} 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "MatNestDestroyISList" 149 static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list) 150 { 151 PetscErrorCode ierr; 152 IS *lst = *list; 153 PetscInt i; 154 155 PetscFunctionBegin; 156 if (!lst) PetscFunctionReturn(0); 157 for (i=0; i<n; i++) if (lst[i]) {ierr = ISDestroy(&lst[i]);CHKERRQ(ierr);} 158 ierr = PetscFree(lst);CHKERRQ(ierr); 159 *list = NULL; 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatDestroy_Nest" 165 static PetscErrorCode MatDestroy_Nest(Mat A) 166 { 167 Mat_Nest *vs = (Mat_Nest*)A->data; 168 PetscInt i,j; 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 /* release the matrices and the place holders */ 173 ierr = MatNestDestroyISList(vs->nr,&vs->isglobal.row);CHKERRQ(ierr); 174 ierr = MatNestDestroyISList(vs->nc,&vs->isglobal.col);CHKERRQ(ierr); 175 ierr = MatNestDestroyISList(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 176 ierr = MatNestDestroyISList(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 177 178 ierr = PetscFree(vs->row_len);CHKERRQ(ierr); 179 ierr = PetscFree(vs->col_len);CHKERRQ(ierr); 180 181 ierr = PetscFree2(vs->left,vs->right);CHKERRQ(ierr); 182 183 /* release the matrices and the place holders */ 184 if (vs->m) { 185 for (i=0; i<vs->nr; i++) { 186 for (j=0; j<vs->nc; j++) { 187 ierr = MatDestroy(&vs->m[i][j]);CHKERRQ(ierr); 188 } 189 ierr = PetscFree(vs->m[i]);CHKERRQ(ierr); 190 } 191 ierr = PetscFree(vs->m);CHKERRQ(ierr); 192 } 193 ierr = PetscFree(A->data);CHKERRQ(ierr); 194 195 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);CHKERRQ(ierr); 196 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);CHKERRQ(ierr); 197 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);CHKERRQ(ierr); 198 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);CHKERRQ(ierr); 199 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);CHKERRQ(ierr); 200 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);CHKERRQ(ierr); 201 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);CHKERRQ(ierr); 202 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatAssemblyBegin_Nest" 208 static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type) 209 { 210 Mat_Nest *vs = (Mat_Nest*)A->data; 211 PetscInt i,j; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 for (i=0; i<vs->nr; i++) { 216 for (j=0; j<vs->nc; j++) { 217 if (vs->m[i][j]) { 218 ierr = MatAssemblyBegin(vs->m[i][j],type);CHKERRQ(ierr); 219 if (!vs->splitassembly) { 220 /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested 221 * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was 222 * already performing an assembly, but the result would by more complicated and appears to offer less 223 * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an 224 * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives". 225 */ 226 ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr); 227 } 228 } 229 } 230 } 231 PetscFunctionReturn(0); 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "MatAssemblyEnd_Nest" 236 static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type) 237 { 238 Mat_Nest *vs = (Mat_Nest*)A->data; 239 PetscInt i,j; 240 PetscErrorCode ierr; 241 242 PetscFunctionBegin; 243 for (i=0; i<vs->nr; i++) { 244 for (j=0; j<vs->nc; j++) { 245 if (vs->m[i][j]) { 246 if (vs->splitassembly) { 247 ierr = MatAssemblyEnd(vs->m[i][j],type);CHKERRQ(ierr); 248 } 249 } 250 } 251 } 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatNestFindNonzeroSubMatRow" 257 static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B) 258 { 259 PetscErrorCode ierr; 260 Mat_Nest *vs = (Mat_Nest*)A->data; 261 PetscInt j; 262 Mat sub; 263 264 PetscFunctionBegin; 265 sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */ 266 for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j]; 267 if (sub) {ierr = MatSetUp(sub);CHKERRQ(ierr);} /* Ensure that the sizes are available */ 268 *B = sub; 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "MatNestFindNonzeroSubMatCol" 274 static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B) 275 { 276 PetscErrorCode ierr; 277 Mat_Nest *vs = (Mat_Nest*)A->data; 278 PetscInt i; 279 Mat sub; 280 281 PetscFunctionBegin; 282 sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */ 283 for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col]; 284 if (sub) {ierr = MatSetUp(sub);CHKERRQ(ierr);} /* Ensure that the sizes are available */ 285 *B = sub; 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatNestFindIS" 291 static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found) 292 { 293 PetscErrorCode ierr; 294 PetscInt i; 295 PetscBool flg; 296 297 PetscFunctionBegin; 298 PetscValidPointer(list,3); 299 PetscValidHeaderSpecific(is,IS_CLASSID,4); 300 PetscValidIntPointer(found,5); 301 *found = -1; 302 for (i=0; i<n; i++) { 303 if (!list[i]) continue; 304 ierr = ISEqual(list[i],is,&flg);CHKERRQ(ierr); 305 if (flg) { 306 *found = i; 307 PetscFunctionReturn(0); 308 } 309 } 310 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set"); 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "MatNestGetRow" 316 /* Get a block row as a new MatNest */ 317 static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B) 318 { 319 Mat_Nest *vs = (Mat_Nest*)A->data; 320 char keyname[256]; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 *B = NULL; 325 ierr = PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);CHKERRQ(ierr); 326 ierr = PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);CHKERRQ(ierr); 327 if (*B) PetscFunctionReturn(0); 328 329 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);CHKERRQ(ierr); 330 331 (*B)->assembled = A->assembled; 332 333 ierr = PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);CHKERRQ(ierr); 334 ierr = PetscObjectDereference((PetscObject)*B);CHKERRQ(ierr); /* Leave the only remaining reference in the composition */ 335 PetscFunctionReturn(0); 336 } 337 338 #undef __FUNCT__ 339 #define __FUNCT__ "MatNestFindSubMat" 340 static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B) 341 { 342 Mat_Nest *vs = (Mat_Nest*)A->data; 343 PetscErrorCode ierr; 344 PetscInt row,col; 345 PetscBool same,isFullCol,isFullColGlobal; 346 347 PetscFunctionBegin; 348 /* Check if full column space. This is a hack */ 349 isFullCol = PETSC_FALSE; 350 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);CHKERRQ(ierr); 351 if (same) { 352 PetscInt n,first,step,i,an,am,afirst,astep; 353 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 354 ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr); 355 isFullCol = PETSC_TRUE; 356 for (i=0,an=A->cmap->rstart; i<vs->nc; i++) { 357 ierr = ISStrideGetInfo(is->col[i],&afirst,&astep);CHKERRQ(ierr); 358 ierr = ISGetLocalSize(is->col[i],&am);CHKERRQ(ierr); 359 if (afirst != an || astep != step) isFullCol = PETSC_FALSE; 360 an += am; 361 } 362 if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE; 363 } 364 ierr = MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));CHKERRQ(ierr); 365 366 if (isFullColGlobal) { 367 PetscInt row; 368 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr); 369 ierr = MatNestGetRow(A,row,B);CHKERRQ(ierr); 370 } else { 371 ierr = MatNestFindIS(A,vs->nr,is->row,isrow,&row);CHKERRQ(ierr); 372 ierr = MatNestFindIS(A,vs->nc,is->col,iscol,&col);CHKERRQ(ierr); 373 *B = vs->m[row][col]; 374 } 375 PetscFunctionReturn(0); 376 } 377 378 #undef __FUNCT__ 379 #define __FUNCT__ "MatGetSubMatrix_Nest" 380 static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B) 381 { 382 PetscErrorCode ierr; 383 Mat_Nest *vs = (Mat_Nest*)A->data; 384 Mat sub; 385 386 PetscFunctionBegin; 387 ierr = MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);CHKERRQ(ierr); 388 switch (reuse) { 389 case MAT_INITIAL_MATRIX: 390 if (sub) { ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr); } 391 *B = sub; 392 break; 393 case MAT_REUSE_MATRIX: 394 if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call"); 395 break; 396 case MAT_IGNORE_MATRIX: /* Nothing to do */ 397 break; 398 } 399 PetscFunctionReturn(0); 400 } 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "MatGetLocalSubMatrix_Nest" 404 PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 405 { 406 PetscErrorCode ierr; 407 Mat_Nest *vs = (Mat_Nest*)A->data; 408 Mat sub; 409 410 PetscFunctionBegin; 411 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 412 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */ 413 if (sub) {ierr = PetscObjectReference((PetscObject)sub);CHKERRQ(ierr);} 414 *B = sub; 415 PetscFunctionReturn(0); 416 } 417 418 #undef __FUNCT__ 419 #define __FUNCT__ "MatRestoreLocalSubMatrix_Nest" 420 static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B) 421 { 422 PetscErrorCode ierr; 423 Mat_Nest *vs = (Mat_Nest*)A->data; 424 Mat sub; 425 426 PetscFunctionBegin; 427 ierr = MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);CHKERRQ(ierr); 428 if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten"); 429 if (sub) { 430 if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times"); 431 ierr = MatDestroy(B);CHKERRQ(ierr); 432 } 433 PetscFunctionReturn(0); 434 } 435 436 #undef __FUNCT__ 437 #define __FUNCT__ "MatGetDiagonal_Nest" 438 static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v) 439 { 440 Mat_Nest *bA = (Mat_Nest*)A->data; 441 PetscInt i; 442 PetscErrorCode ierr; 443 444 PetscFunctionBegin; 445 for (i=0; i<bA->nr; i++) { 446 Vec bv; 447 ierr = VecGetSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 448 if (bA->m[i][i]) { 449 ierr = MatGetDiagonal(bA->m[i][i],bv);CHKERRQ(ierr); 450 } else { 451 ierr = VecSet(bv,0.0);CHKERRQ(ierr); 452 } 453 ierr = VecRestoreSubVector(v,bA->isglobal.row[i],&bv);CHKERRQ(ierr); 454 } 455 PetscFunctionReturn(0); 456 } 457 458 #undef __FUNCT__ 459 #define __FUNCT__ "MatDiagonalScale_Nest" 460 static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r) 461 { 462 Mat_Nest *bA = (Mat_Nest*)A->data; 463 Vec bl,*br; 464 PetscInt i,j; 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 ierr = PetscCalloc1(bA->nc,&br);CHKERRQ(ierr); 469 if (r) { 470 for (j=0; j<bA->nc; j++) {ierr = VecGetSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);} 471 } 472 bl = NULL; 473 for (i=0; i<bA->nr; i++) { 474 if (l) { 475 ierr = VecGetSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr); 476 } 477 for (j=0; j<bA->nc; j++) { 478 if (bA->m[i][j]) { 479 ierr = MatDiagonalScale(bA->m[i][j],bl,br[j]);CHKERRQ(ierr); 480 } 481 } 482 if (l) { 483 ierr = VecRestoreSubVector(l,bA->isglobal.row[i],&bl);CHKERRQ(ierr); 484 } 485 } 486 if (r) { 487 for (j=0; j<bA->nc; j++) {ierr = VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);CHKERRQ(ierr);} 488 } 489 ierr = PetscFree(br);CHKERRQ(ierr); 490 PetscFunctionReturn(0); 491 } 492 493 #undef __FUNCT__ 494 #define __FUNCT__ "MatScale_Nest" 495 static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a) 496 { 497 Mat_Nest *bA = (Mat_Nest*)A->data; 498 PetscInt i,j; 499 PetscErrorCode ierr; 500 501 PetscFunctionBegin; 502 for (i=0; i<bA->nr; i++) { 503 for (j=0; j<bA->nc; j++) { 504 if (bA->m[i][j]) { 505 ierr = MatScale(bA->m[i][j],a);CHKERRQ(ierr); 506 } 507 } 508 } 509 PetscFunctionReturn(0); 510 } 511 512 #undef __FUNCT__ 513 #define __FUNCT__ "MatShift_Nest" 514 static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a) 515 { 516 Mat_Nest *bA = (Mat_Nest*)A->data; 517 PetscInt i; 518 PetscErrorCode ierr; 519 520 PetscFunctionBegin; 521 for (i=0; i<bA->nr; i++) { 522 if (!bA->m[i][i]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for shifting an empty diagonal block, insert a matrix in block (%D,%D)",i,i); 523 ierr = MatShift(bA->m[i][i],a);CHKERRQ(ierr); 524 } 525 PetscFunctionReturn(0); 526 } 527 528 #undef __FUNCT__ 529 #define __FUNCT__ "MatCreateVecs_Nest" 530 static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left) 531 { 532 Mat_Nest *bA = (Mat_Nest*)A->data; 533 Vec *L,*R; 534 MPI_Comm comm; 535 PetscInt i,j; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 540 if (right) { 541 /* allocate R */ 542 ierr = PetscMalloc1(bA->nc, &R);CHKERRQ(ierr); 543 /* Create the right vectors */ 544 for (j=0; j<bA->nc; j++) { 545 for (i=0; i<bA->nr; i++) { 546 if (bA->m[i][j]) { 547 ierr = MatCreateVecs(bA->m[i][j],&R[j],NULL);CHKERRQ(ierr); 548 break; 549 } 550 } 551 if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 552 } 553 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr); 554 /* hand back control to the nest vector */ 555 for (j=0; j<bA->nc; j++) { 556 ierr = VecDestroy(&R[j]);CHKERRQ(ierr); 557 } 558 ierr = PetscFree(R);CHKERRQ(ierr); 559 } 560 561 if (left) { 562 /* allocate L */ 563 ierr = PetscMalloc1(bA->nr, &L);CHKERRQ(ierr); 564 /* Create the left vectors */ 565 for (i=0; i<bA->nr; i++) { 566 for (j=0; j<bA->nc; j++) { 567 if (bA->m[i][j]) { 568 ierr = MatCreateVecs(bA->m[i][j],NULL,&L[i]);CHKERRQ(ierr); 569 break; 570 } 571 } 572 if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 573 } 574 575 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr); 576 for (i=0; i<bA->nr; i++) { 577 ierr = VecDestroy(&L[i]);CHKERRQ(ierr); 578 } 579 580 ierr = PetscFree(L);CHKERRQ(ierr); 581 } 582 PetscFunctionReturn(0); 583 } 584 585 #undef __FUNCT__ 586 #define __FUNCT__ "MatView_Nest" 587 static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer) 588 { 589 Mat_Nest *bA = (Mat_Nest*)A->data; 590 PetscBool isascii; 591 PetscInt i,j; 592 PetscErrorCode ierr; 593 594 PetscFunctionBegin; 595 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 596 if (isascii) { 597 598 PetscViewerASCIIPrintf(viewer,"Matrix object: \n"); 599 PetscViewerASCIIPushTab(viewer); /* push0 */ 600 PetscViewerASCIIPrintf(viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc); 601 602 PetscViewerASCIIPrintf(viewer,"MatNest structure: \n"); 603 for (i=0; i<bA->nr; i++) { 604 for (j=0; j<bA->nc; j++) { 605 MatType type; 606 char name[256] = "",prefix[256] = ""; 607 PetscInt NR,NC; 608 PetscBool isNest = PETSC_FALSE; 609 610 if (!bA->m[i][j]) { 611 PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j); 612 continue; 613 } 614 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr); 615 ierr = MatGetType(bA->m[i][j], &type);CHKERRQ(ierr); 616 if (((PetscObject)bA->m[i][j])->name) {ierr = PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);CHKERRQ(ierr);} 617 if (((PetscObject)bA->m[i][j])->prefix) {ierr = PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);CHKERRQ(ierr);} 618 ierr = PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr); 619 620 ierr = PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);CHKERRQ(ierr); 621 622 if (isNest) { 623 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); /* push1 */ 624 ierr = MatView(bA->m[i][j],viewer);CHKERRQ(ierr); 625 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop1 */ 626 } 627 } 628 } 629 PetscViewerASCIIPopTab(viewer); /* pop0 */ 630 } 631 PetscFunctionReturn(0); 632 } 633 634 #undef __FUNCT__ 635 #define __FUNCT__ "MatZeroEntries_Nest" 636 static PetscErrorCode MatZeroEntries_Nest(Mat A) 637 { 638 Mat_Nest *bA = (Mat_Nest*)A->data; 639 PetscInt i,j; 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 for (i=0; i<bA->nr; i++) { 644 for (j=0; j<bA->nc; j++) { 645 if (!bA->m[i][j]) continue; 646 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr); 647 } 648 } 649 PetscFunctionReturn(0); 650 } 651 652 #undef __FUNCT__ 653 #define __FUNCT__ "MatCopy_Nest" 654 static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str) 655 { 656 Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data; 657 PetscInt i,j,nr = bA->nr,nc = bA->nc; 658 PetscErrorCode ierr; 659 660 PetscFunctionBegin; 661 if (nr != bB->nr || nc != bB->nc) SETERRQ4(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot copy a Mat_Nest of block size (%D,%D) to a Mat_Nest of block size (%D,%D)",bB->nr,bB->nc,nr,nc); 662 for (i=0; i<nr; i++) { 663 for (j=0; j<nc; j++) { 664 if (bA->m[i][j] && bB->m[i][j]) { 665 ierr = MatCopy(bA->m[i][j],bB->m[i][j],str);CHKERRQ(ierr); 666 } else if (bA->m[i][j] || bB->m[i][j]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j); 667 } 668 } 669 PetscFunctionReturn(0); 670 } 671 672 #undef __FUNCT__ 673 #define __FUNCT__ "MatDuplicate_Nest" 674 static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 675 { 676 Mat_Nest *bA = (Mat_Nest*)A->data; 677 Mat *b; 678 PetscInt i,j,nr = bA->nr,nc = bA->nc; 679 PetscErrorCode ierr; 680 681 PetscFunctionBegin; 682 ierr = PetscMalloc1(nr*nc,&b);CHKERRQ(ierr); 683 for (i=0; i<nr; i++) { 684 for (j=0; j<nc; j++) { 685 if (bA->m[i][j]) { 686 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 687 } else { 688 b[i*nc+j] = NULL; 689 } 690 } 691 } 692 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 693 /* Give the new MatNest exclusive ownership */ 694 for (i=0; i<nr*nc; i++) { 695 ierr = MatDestroy(&b[i]);CHKERRQ(ierr); 696 } 697 ierr = PetscFree(b);CHKERRQ(ierr); 698 699 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 700 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 701 PetscFunctionReturn(0); 702 } 703 704 /* nest api */ 705 #undef __FUNCT__ 706 #define __FUNCT__ "MatNestGetSubMat_Nest" 707 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 708 { 709 Mat_Nest *bA = (Mat_Nest*)A->data; 710 711 PetscFunctionBegin; 712 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 713 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 714 *mat = bA->m[idxm][jdxm]; 715 PetscFunctionReturn(0); 716 } 717 718 #undef __FUNCT__ 719 #define __FUNCT__ "MatNestGetSubMat" 720 /*@ 721 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 722 723 Not collective 724 725 Input Parameters: 726 + A - nest matrix 727 . idxm - index of the matrix within the nest matrix 728 - jdxm - index of the matrix within the nest matrix 729 730 Output Parameter: 731 . sub - matrix at index idxm,jdxm within the nest matrix 732 733 Level: developer 734 735 .seealso: MatNestGetSize(), MatNestGetSubMats() 736 @*/ 737 PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 738 { 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 743 PetscFunctionReturn(0); 744 } 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "MatNestSetSubMat_Nest" 748 PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat) 749 { 750 Mat_Nest *bA = (Mat_Nest*)A->data; 751 PetscInt m,n,M,N,mi,ni,Mi,Ni; 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 756 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 757 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 758 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 759 ierr = ISGetLocalSize(bA->isglobal.row[idxm],&mi);CHKERRQ(ierr); 760 ierr = ISGetSize(bA->isglobal.row[idxm],&Mi);CHKERRQ(ierr); 761 ierr = ISGetLocalSize(bA->isglobal.col[jdxm],&ni);CHKERRQ(ierr); 762 ierr = ISGetSize(bA->isglobal.col[jdxm],&Ni);CHKERRQ(ierr); 763 if (M != Mi || N != Ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix dimension (%D,%D) incompatible with nest block (%D,%D)",M,N,Mi,Ni); 764 if (m != mi || n != ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix local dimension (%D,%D) incompatible with nest block (%D,%D)",m,n,mi,ni); 765 766 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 767 ierr = MatDestroy(&bA->m[idxm][jdxm]);CHKERRQ(ierr); 768 bA->m[idxm][jdxm] = mat; 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNCT__ 773 #define __FUNCT__ "MatNestSetSubMat" 774 /*@ 775 MatNestSetSubMat - Set a single submatrix in the nest matrix. 776 777 Logically collective on the submatrix communicator 778 779 Input Parameters: 780 + A - nest matrix 781 . idxm - index of the matrix within the nest matrix 782 . jdxm - index of the matrix within the nest matrix 783 - sub - matrix at index idxm,jdxm within the nest matrix 784 785 Notes: 786 The new submatrix must have the same size and communicator as that block of the nest. 787 788 This increments the reference count of the submatrix. 789 790 Level: developer 791 792 .seealso: MatNestSetSubMats(), MatNestGetSubMat() 793 @*/ 794 PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub) 795 { 796 PetscErrorCode ierr; 797 798 PetscFunctionBegin; 799 ierr = PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));CHKERRQ(ierr); 800 PetscFunctionReturn(0); 801 } 802 803 #undef __FUNCT__ 804 #define __FUNCT__ "MatNestGetSubMats_Nest" 805 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 806 { 807 Mat_Nest *bA = (Mat_Nest*)A->data; 808 809 PetscFunctionBegin; 810 if (M) *M = bA->nr; 811 if (N) *N = bA->nc; 812 if (mat) *mat = bA->m; 813 PetscFunctionReturn(0); 814 } 815 816 #undef __FUNCT__ 817 #define __FUNCT__ "MatNestGetSubMats" 818 /*@C 819 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 820 821 Not collective 822 823 Input Parameters: 824 . A - nest matrix 825 826 Output Parameter: 827 + M - number of rows in the nest matrix 828 . N - number of cols in the nest matrix 829 - mat - 2d array of matrices 830 831 Notes: 832 833 The user should not free the array mat. 834 835 Level: developer 836 837 .seealso: MatNestGetSize(), MatNestGetSubMat() 838 @*/ 839 PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 840 { 841 PetscErrorCode ierr; 842 843 PetscFunctionBegin; 844 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "MatNestGetSize_Nest" 850 PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 851 { 852 Mat_Nest *bA = (Mat_Nest*)A->data; 853 854 PetscFunctionBegin; 855 if (M) *M = bA->nr; 856 if (N) *N = bA->nc; 857 PetscFunctionReturn(0); 858 } 859 860 #undef __FUNCT__ 861 #define __FUNCT__ "MatNestGetSize" 862 /*@ 863 MatNestGetSize - Returns the size of the nest matrix. 864 865 Not collective 866 867 Input Parameters: 868 . A - nest matrix 869 870 Output Parameter: 871 + M - number of rows in the nested mat 872 - N - number of cols in the nested mat 873 874 Notes: 875 876 Level: developer 877 878 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 879 @*/ 880 PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 881 { 882 PetscErrorCode ierr; 883 884 PetscFunctionBegin; 885 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 #undef __FUNCT__ 890 #define __FUNCT__ "MatNestGetISs_Nest" 891 static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[]) 892 { 893 Mat_Nest *vs = (Mat_Nest*)A->data; 894 PetscInt i; 895 896 PetscFunctionBegin; 897 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i]; 898 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i]; 899 PetscFunctionReturn(0); 900 } 901 902 #undef __FUNCT__ 903 #define __FUNCT__ "MatNestGetISs" 904 /*@C 905 MatNestGetISs - Returns the index sets partitioning the row and column spaces 906 907 Not collective 908 909 Input Parameters: 910 . A - nest matrix 911 912 Output Parameter: 913 + rows - array of row index sets 914 - cols - array of column index sets 915 916 Level: advanced 917 918 Notes: 919 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 920 921 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs() 922 @*/ 923 PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[]) 924 { 925 PetscErrorCode ierr; 926 927 PetscFunctionBegin; 928 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 929 ierr = PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 930 PetscFunctionReturn(0); 931 } 932 933 #undef __FUNCT__ 934 #define __FUNCT__ "MatNestGetLocalISs_Nest" 935 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[]) 936 { 937 Mat_Nest *vs = (Mat_Nest*)A->data; 938 PetscInt i; 939 940 PetscFunctionBegin; 941 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i]; 942 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i]; 943 PetscFunctionReturn(0); 944 } 945 946 #undef __FUNCT__ 947 #define __FUNCT__ "MatNestGetLocalISs" 948 /*@C 949 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces 950 951 Not collective 952 953 Input Parameters: 954 . A - nest matrix 955 956 Output Parameter: 957 + rows - array of row index sets (or NULL to ignore) 958 - cols - array of column index sets (or NULL to ignore) 959 960 Level: advanced 961 962 Notes: 963 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 964 965 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs() 966 @*/ 967 PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[]) 968 { 969 PetscErrorCode ierr; 970 971 PetscFunctionBegin; 972 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 973 ierr = PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 974 PetscFunctionReturn(0); 975 } 976 977 #undef __FUNCT__ 978 #define __FUNCT__ "MatNestSetVecType_Nest" 979 PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype) 980 { 981 PetscErrorCode ierr; 982 PetscBool flg; 983 984 PetscFunctionBegin; 985 ierr = PetscStrcmp(vtype,VECNEST,&flg);CHKERRQ(ierr); 986 /* In reality, this only distinguishes VECNEST and "other" */ 987 if (flg) A->ops->getvecs = MatCreateVecs_Nest; 988 else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0; 989 PetscFunctionReturn(0); 990 } 991 992 #undef __FUNCT__ 993 #define __FUNCT__ "MatNestSetVecType" 994 /*@C 995 MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs() 996 997 Not collective 998 999 Input Parameters: 1000 + A - nest matrix 1001 - vtype - type to use for creating vectors 1002 1003 Notes: 1004 1005 Level: developer 1006 1007 .seealso: MatCreateVecs() 1008 @*/ 1009 PetscErrorCode MatNestSetVecType(Mat A,VecType vtype) 1010 { 1011 PetscErrorCode ierr; 1012 1013 PetscFunctionBegin; 1014 ierr = PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));CHKERRQ(ierr); 1015 PetscFunctionReturn(0); 1016 } 1017 1018 #undef __FUNCT__ 1019 #define __FUNCT__ "MatNestSetSubMats_Nest" 1020 PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1021 { 1022 Mat_Nest *s = (Mat_Nest*)A->data; 1023 PetscInt i,j,m,n,M,N; 1024 PetscErrorCode ierr; 1025 1026 PetscFunctionBegin; 1027 s->nr = nr; 1028 s->nc = nc; 1029 1030 /* Create space for submatrices */ 1031 ierr = PetscMalloc1(nr,&s->m);CHKERRQ(ierr); 1032 for (i=0; i<nr; i++) { 1033 ierr = PetscMalloc1(nc,&s->m[i]);CHKERRQ(ierr); 1034 } 1035 for (i=0; i<nr; i++) { 1036 for (j=0; j<nc; j++) { 1037 s->m[i][j] = a[i*nc+j]; 1038 if (a[i*nc+j]) { 1039 ierr = PetscObjectReference((PetscObject)a[i*nc+j]);CHKERRQ(ierr); 1040 } 1041 } 1042 } 1043 1044 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1045 1046 ierr = PetscMalloc1(nr,&s->row_len);CHKERRQ(ierr); 1047 ierr = PetscMalloc1(nc,&s->col_len);CHKERRQ(ierr); 1048 for (i=0; i<nr; i++) s->row_len[i]=-1; 1049 for (j=0; j<nc; j++) s->col_len[j]=-1; 1050 1051 ierr = MatNestGetSizes_Private(A,&m,&n,&M,&N);CHKERRQ(ierr); 1052 1053 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1054 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1055 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1056 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1057 1058 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1059 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1060 1061 ierr = PetscCalloc2(nr,&s->left,nc,&s->right);CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 #undef __FUNCT__ 1066 #define __FUNCT__ "MatNestSetSubMats" 1067 /*@ 1068 MatNestSetSubMats - Sets the nested submatrices 1069 1070 Collective on Mat 1071 1072 Input Parameter: 1073 + N - nested matrix 1074 . nr - number of nested row blocks 1075 . is_row - index sets for each nested row block, or NULL to make contiguous 1076 . nc - number of nested column blocks 1077 . is_col - index sets for each nested column block, or NULL to make contiguous 1078 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1079 1080 Level: advanced 1081 1082 .seealso: MatCreateNest(), MATNEST 1083 @*/ 1084 PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1085 { 1086 PetscErrorCode ierr; 1087 PetscInt i; 1088 1089 PetscFunctionBegin; 1090 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1091 if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1092 if (nr && is_row) { 1093 PetscValidPointer(is_row,3); 1094 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1095 } 1096 if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1097 if (nc && is_col) { 1098 PetscValidPointer(is_col,5); 1099 for (i=0; i<nc; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1100 } 1101 if (nr*nc) PetscValidPointer(a,6); 1102 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr); 1103 PetscFunctionReturn(0); 1104 } 1105 1106 #undef __FUNCT__ 1107 #define __FUNCT__ "MatNestCreateAggregateL2G_Private" 1108 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog) 1109 { 1110 PetscErrorCode ierr; 1111 PetscBool flg; 1112 PetscInt i,j,m,mi,*ix; 1113 1114 PetscFunctionBegin; 1115 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1116 if (islocal[i]) { 1117 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1118 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1119 } else { 1120 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1121 } 1122 m += mi; 1123 } 1124 if (flg) { 1125 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1126 for (i=0,n=0; i<n; i++) { 1127 ISLocalToGlobalMapping smap = NULL; 1128 VecScatter scat; 1129 IS isreq; 1130 Vec lvec,gvec; 1131 union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x; 1132 Mat sub; 1133 1134 if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers"); 1135 if (colflg) { 1136 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1137 } else { 1138 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1139 } 1140 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr);} 1141 if (islocal[i]) { 1142 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1143 } else { 1144 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1145 } 1146 for (j=0; j<mi; j++) ix[m+j] = j; 1147 if (smap) {ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);} 1148 /* 1149 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1150 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1151 The approach here is ugly because it uses VecScatter to move indices. 1152 */ 1153 ierr = VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);CHKERRQ(ierr); 1154 ierr = VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);CHKERRQ(ierr); 1155 ierr = ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);CHKERRQ(ierr); 1156 ierr = VecScatterCreate(gvec,isreq,lvec,NULL,&scat);CHKERRQ(ierr); 1157 ierr = VecGetArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1158 for (j=0; j<mi; j++) x[j].integer = ix[m+j]; 1159 ierr = VecRestoreArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1160 ierr = VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1161 ierr = VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1162 ierr = VecGetArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1163 for (j=0; j<mi; j++) ix[m+j] = x[j].integer; 1164 ierr = VecRestoreArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1165 ierr = VecDestroy(&lvec);CHKERRQ(ierr); 1166 ierr = VecDestroy(&gvec);CHKERRQ(ierr); 1167 ierr = ISDestroy(&isreq);CHKERRQ(ierr); 1168 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 1169 m += mi; 1170 } 1171 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1172 } else { 1173 *ltog = NULL; 1174 } 1175 PetscFunctionReturn(0); 1176 } 1177 1178 1179 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1180 /* 1181 nprocessors = NP 1182 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1183 proc 0: => (g_0,h_0,) 1184 proc 1: => (g_1,h_1,) 1185 ... 1186 proc nprocs-1: => (g_NP-1,h_NP-1,) 1187 1188 proc 0: proc 1: proc nprocs-1: 1189 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1190 1191 proc 0: 1192 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1193 proc 1: 1194 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1195 1196 proc NP-1: 1197 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1198 */ 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatSetUp_NestIS_Private" 1201 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1202 { 1203 Mat_Nest *vs = (Mat_Nest*)A->data; 1204 PetscInt i,j,offset,n,nsum,bs; 1205 PetscErrorCode ierr; 1206 Mat sub = NULL; 1207 1208 PetscFunctionBegin; 1209 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1210 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1211 if (is_row) { /* valid IS is passed in */ 1212 /* refs on is[] are incremeneted */ 1213 for (i=0; i<vs->nr; i++) { 1214 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1215 1216 vs->isglobal.row[i] = is_row[i]; 1217 } 1218 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1219 nsum = 0; 1220 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1221 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1222 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1223 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1224 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1225 nsum += n; 1226 } 1227 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1228 offset -= nsum; 1229 for (i=0; i<vs->nr; i++) { 1230 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1231 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1232 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1233 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1234 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1235 offset += n; 1236 } 1237 } 1238 1239 if (is_col) { /* valid IS is passed in */ 1240 /* refs on is[] are incremeneted */ 1241 for (j=0; j<vs->nc; j++) { 1242 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1243 1244 vs->isglobal.col[j] = is_col[j]; 1245 } 1246 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1247 offset = A->cmap->rstart; 1248 nsum = 0; 1249 for (j=0; j<vs->nc; j++) { 1250 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1251 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1252 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1253 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1254 nsum += n; 1255 } 1256 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1257 offset -= nsum; 1258 for (j=0; j<vs->nc; j++) { 1259 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1260 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1261 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1262 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1263 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1264 offset += n; 1265 } 1266 } 1267 1268 /* Set up the local ISs */ 1269 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1270 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1271 for (i=0,offset=0; i<vs->nr; i++) { 1272 IS isloc; 1273 ISLocalToGlobalMapping rmap = NULL; 1274 PetscInt nlocal,bs; 1275 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1276 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1277 if (rmap) { 1278 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1279 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1280 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1281 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1282 } else { 1283 nlocal = 0; 1284 isloc = NULL; 1285 } 1286 vs->islocal.row[i] = isloc; 1287 offset += nlocal; 1288 } 1289 for (i=0,offset=0; i<vs->nc; i++) { 1290 IS isloc; 1291 ISLocalToGlobalMapping cmap = NULL; 1292 PetscInt nlocal,bs; 1293 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1294 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1295 if (cmap) { 1296 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1297 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1298 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1299 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1300 } else { 1301 nlocal = 0; 1302 isloc = NULL; 1303 } 1304 vs->islocal.col[i] = isloc; 1305 offset += nlocal; 1306 } 1307 1308 /* Set up the aggregate ISLocalToGlobalMapping */ 1309 { 1310 ISLocalToGlobalMapping rmap,cmap; 1311 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1312 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1313 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1314 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1315 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1316 } 1317 1318 #if defined(PETSC_USE_DEBUG) 1319 for (i=0; i<vs->nr; i++) { 1320 for (j=0; j<vs->nc; j++) { 1321 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1322 Mat B = vs->m[i][j]; 1323 if (!B) continue; 1324 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1325 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1326 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1327 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1328 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1329 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1330 if (M != Mi || N != Ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Global sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",M,N,i,j,Mi,Ni); 1331 if (m != mi || n != ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Local sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",m,n,i,j,mi,ni); 1332 } 1333 } 1334 #endif 1335 1336 /* Set A->assembled if all non-null blocks are currently assembled */ 1337 for (i=0; i<vs->nr; i++) { 1338 for (j=0; j<vs->nc; j++) { 1339 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1340 } 1341 } 1342 A->assembled = PETSC_TRUE; 1343 PetscFunctionReturn(0); 1344 } 1345 1346 #undef __FUNCT__ 1347 #define __FUNCT__ "MatCreateNest" 1348 /*@C 1349 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1350 1351 Collective on Mat 1352 1353 Input Parameter: 1354 + comm - Communicator for the new Mat 1355 . nr - number of nested row blocks 1356 . is_row - index sets for each nested row block, or NULL to make contiguous 1357 . nc - number of nested column blocks 1358 . is_col - index sets for each nested column block, or NULL to make contiguous 1359 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1360 1361 Output Parameter: 1362 . B - new matrix 1363 1364 Level: advanced 1365 1366 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST 1367 @*/ 1368 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1369 { 1370 Mat A; 1371 PetscErrorCode ierr; 1372 1373 PetscFunctionBegin; 1374 *B = 0; 1375 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1376 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1377 ierr = MatSetUp(A);CHKERRQ(ierr); 1378 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1379 *B = A; 1380 PetscFunctionReturn(0); 1381 } 1382 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "MatConvert_Nest_AIJ" 1385 PETSC_EXTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1386 { 1387 PetscErrorCode ierr; 1388 Mat_Nest *nest = (Mat_Nest*)A->data; 1389 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1390 PetscInt cstart,cend; 1391 Mat C; 1392 1393 PetscFunctionBegin; 1394 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1395 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1396 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1397 switch (reuse) { 1398 case MAT_INITIAL_MATRIX: 1399 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1400 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1401 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1402 *newmat = C; 1403 break; 1404 case MAT_REUSE_MATRIX: 1405 C = *newmat; 1406 break; 1407 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1408 } 1409 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1410 onnz = dnnz + m; 1411 for (k=0; k<m; k++) { 1412 dnnz[k] = 0; 1413 onnz[k] = 0; 1414 } 1415 for (j=0; j<nest->nc; ++j) { 1416 IS bNis; 1417 PetscInt bN; 1418 const PetscInt *bNindices; 1419 /* Using global column indices and ISAllGather() is not scalable. */ 1420 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1421 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1422 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1423 for (i=0; i<nest->nr; ++i) { 1424 PetscSF bmsf; 1425 PetscSFNode *iremote; 1426 Mat B; 1427 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1428 const PetscInt *bmindices; 1429 B = nest->m[i][j]; 1430 if (!B) continue; 1431 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1432 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1433 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1434 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1435 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1436 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1437 for (k = 0; k < bm; ++k){ 1438 sub_dnnz[k] = 0; 1439 sub_onnz[k] = 0; 1440 } 1441 /* 1442 Locate the owners for all of the locally-owned global row indices for this row block. 1443 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1444 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1445 */ 1446 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1447 for (br = 0; br < bm; ++br) { 1448 PetscInt row = bmindices[br], rowowner = 0, brncols, col; 1449 const PetscInt *brcols; 1450 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1451 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1452 /* how many roots */ 1453 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1454 /* get nonzero pattern */ 1455 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1456 for (k=0; k<brncols; k++) { 1457 col = bNindices[brcols[k]]; 1458 if(col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]){ 1459 sub_dnnz[br]++; 1460 }else{ 1461 sub_onnz[br]++; 1462 } 1463 } 1464 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1465 } 1466 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1467 /* bsf will have to take care of disposing of bedges. */ 1468 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1469 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1470 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1471 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1472 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1473 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1474 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1475 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1476 } 1477 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1478 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1479 } 1480 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1481 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1482 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1483 1484 /* Fill by row */ 1485 for (j=0; j<nest->nc; ++j) { 1486 /* Using global column indices and ISAllGather() is not scalable. */ 1487 IS bNis; 1488 PetscInt bN; 1489 const PetscInt *bNindices; 1490 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1491 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1492 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1493 for (i=0; i<nest->nr; ++i) { 1494 Mat B; 1495 PetscInt bm, br; 1496 const PetscInt *bmindices; 1497 B = nest->m[i][j]; 1498 if (!B) continue; 1499 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1500 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1501 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1502 for (br = 0; br < bm; ++br) { 1503 PetscInt row = bmindices[br], brncols, *cols; 1504 const PetscInt *brcols; 1505 const PetscScalar *brcoldata; 1506 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1507 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1508 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1509 /* 1510 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1511 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1512 */ 1513 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1514 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1515 ierr = PetscFree(cols);CHKERRQ(ierr); 1516 } 1517 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1518 } 1519 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1520 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1521 } 1522 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1523 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1524 PetscFunctionReturn(0); 1525 } 1526 1527 /*MC 1528 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1529 1530 Level: intermediate 1531 1532 Notes: 1533 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1534 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1535 It is usually used with DMComposite and DMCreateMatrix() 1536 1537 .seealso: MatCreate(), MatType, MatCreateNest() 1538 M*/ 1539 #undef __FUNCT__ 1540 #define __FUNCT__ "MatCreate_Nest" 1541 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1542 { 1543 Mat_Nest *s; 1544 PetscErrorCode ierr; 1545 1546 PetscFunctionBegin; 1547 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1548 A->data = (void*)s; 1549 1550 s->nr = -1; 1551 s->nc = -1; 1552 s->m = NULL; 1553 s->splitassembly = PETSC_FALSE; 1554 1555 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1556 1557 A->ops->mult = MatMult_Nest; 1558 A->ops->multadd = MatMultAdd_Nest; 1559 A->ops->multtranspose = MatMultTranspose_Nest; 1560 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1561 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1562 A->ops->assemblyend = MatAssemblyEnd_Nest; 1563 A->ops->zeroentries = MatZeroEntries_Nest; 1564 A->ops->copy = MatCopy_Nest; 1565 A->ops->duplicate = MatDuplicate_Nest; 1566 A->ops->getsubmatrix = MatGetSubMatrix_Nest; 1567 A->ops->destroy = MatDestroy_Nest; 1568 A->ops->view = MatView_Nest; 1569 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1570 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1571 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1572 A->ops->getdiagonal = MatGetDiagonal_Nest; 1573 A->ops->diagonalscale = MatDiagonalScale_Nest; 1574 A->ops->scale = MatScale_Nest; 1575 A->ops->shift = MatShift_Nest; 1576 1577 A->spptr = 0; 1578 A->assembled = PETSC_FALSE; 1579 1580 /* expose Nest api's */ 1581 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1582 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1583 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1584 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1585 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1586 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1587 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1588 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1589 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1590 1591 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1592 PetscFunctionReturn(0); 1593 } 1594