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