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 MatGetVecs_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 = MPI_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,1.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__ "MatGetVecs_Nest" 530 static PetscErrorCode MatGetVecs_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 = PetscMalloc(sizeof(Vec) * 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 = MatGetVecs(bA->m[i][j],&R[j],NULL);CHKERRQ(ierr); 548 break; 549 } 550 } 551 if (i==bA->nr) { 552 /* have an empty column */ 553 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 554 } 555 } 556 ierr = VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);CHKERRQ(ierr); 557 /* hand back control to the nest vector */ 558 for (j=0; j<bA->nc; j++) { 559 ierr = VecDestroy(&R[j]);CHKERRQ(ierr); 560 } 561 ierr = PetscFree(R);CHKERRQ(ierr); 562 } 563 564 if (left) { 565 /* allocate L */ 566 ierr = PetscMalloc(sizeof(Vec) * bA->nr, &L);CHKERRQ(ierr); 567 /* Create the left vectors */ 568 for (i=0; i<bA->nr; i++) { 569 for (j=0; j<bA->nc; j++) { 570 if (bA->m[i][j]) { 571 ierr = MatGetVecs(bA->m[i][j],NULL,&L[i]);CHKERRQ(ierr); 572 break; 573 } 574 } 575 if (j==bA->nc) { 576 /* have an empty row */ 577 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 578 } 579 } 580 581 ierr = VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);CHKERRQ(ierr); 582 for (i=0; i<bA->nr; i++) { 583 ierr = VecDestroy(&L[i]);CHKERRQ(ierr); 584 } 585 586 ierr = PetscFree(L);CHKERRQ(ierr); 587 } 588 PetscFunctionReturn(0); 589 } 590 591 #undef __FUNCT__ 592 #define __FUNCT__ "MatView_Nest" 593 static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer) 594 { 595 Mat_Nest *bA = (Mat_Nest*)A->data; 596 PetscBool isascii; 597 PetscInt i,j; 598 PetscErrorCode ierr; 599 600 PetscFunctionBegin; 601 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 602 if (isascii) { 603 604 PetscViewerASCIIPrintf(viewer,"Matrix object: \n"); 605 PetscViewerASCIIPushTab(viewer); /* push0 */ 606 PetscViewerASCIIPrintf(viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc); 607 608 PetscViewerASCIIPrintf(viewer,"MatNest structure: \n"); 609 for (i=0; i<bA->nr; i++) { 610 for (j=0; j<bA->nc; j++) { 611 MatType type; 612 char name[256] = "",prefix[256] = ""; 613 PetscInt NR,NC; 614 PetscBool isNest = PETSC_FALSE; 615 616 if (!bA->m[i][j]) { 617 PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j); 618 continue; 619 } 620 ierr = MatGetSize(bA->m[i][j],&NR,&NC);CHKERRQ(ierr); 621 ierr = MatGetType(bA->m[i][j], &type);CHKERRQ(ierr); 622 if (((PetscObject)bA->m[i][j])->name) {ierr = PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);CHKERRQ(ierr);} 623 if (((PetscObject)bA->m[i][j])->prefix) {ierr = PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);CHKERRQ(ierr);} 624 ierr = PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);CHKERRQ(ierr); 625 626 ierr = PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);CHKERRQ(ierr); 627 628 if (isNest) { 629 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); /* push1 */ 630 ierr = MatView(bA->m[i][j],viewer);CHKERRQ(ierr); 631 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); /* pop1 */ 632 } 633 } 634 } 635 PetscViewerASCIIPopTab(viewer); /* pop0 */ 636 } 637 PetscFunctionReturn(0); 638 } 639 640 #undef __FUNCT__ 641 #define __FUNCT__ "MatZeroEntries_Nest" 642 static PetscErrorCode MatZeroEntries_Nest(Mat A) 643 { 644 Mat_Nest *bA = (Mat_Nest*)A->data; 645 PetscInt i,j; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 for (i=0; i<bA->nr; i++) { 650 for (j=0; j<bA->nc; j++) { 651 if (!bA->m[i][j]) continue; 652 ierr = MatZeroEntries(bA->m[i][j]);CHKERRQ(ierr); 653 } 654 } 655 PetscFunctionReturn(0); 656 } 657 658 #undef __FUNCT__ 659 #define __FUNCT__ "MatCopy_Nest" 660 static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str) 661 { 662 Mat_Nest *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data; 663 PetscInt i,j,nr = bA->nr,nc = bA->nc; 664 PetscErrorCode ierr; 665 666 PetscFunctionBegin; 667 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); 668 for (i=0; i<nr; i++) { 669 for (j=0; j<nc; j++) { 670 if (bA->m[i][j] && bB->m[i][j]) { 671 ierr = MatCopy(bA->m[i][j],bB->m[i][j],str);CHKERRQ(ierr); 672 } 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); 673 } 674 } 675 PetscFunctionReturn(0); 676 } 677 678 #undef __FUNCT__ 679 #define __FUNCT__ "MatDuplicate_Nest" 680 static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B) 681 { 682 Mat_Nest *bA = (Mat_Nest*)A->data; 683 Mat *b; 684 PetscInt i,j,nr = bA->nr,nc = bA->nc; 685 PetscErrorCode ierr; 686 687 PetscFunctionBegin; 688 ierr = PetscMalloc1(nr*nc,&b);CHKERRQ(ierr); 689 for (i=0; i<nr; i++) { 690 for (j=0; j<nc; j++) { 691 if (bA->m[i][j]) { 692 ierr = MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);CHKERRQ(ierr); 693 } else { 694 b[i*nc+j] = NULL; 695 } 696 } 697 } 698 ierr = MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);CHKERRQ(ierr); 699 /* Give the new MatNest exclusive ownership */ 700 for (i=0; i<nr*nc; i++) { 701 ierr = MatDestroy(&b[i]);CHKERRQ(ierr); 702 } 703 ierr = PetscFree(b);CHKERRQ(ierr); 704 705 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 706 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 707 PetscFunctionReturn(0); 708 } 709 710 /* nest api */ 711 #undef __FUNCT__ 712 #define __FUNCT__ "MatNestGetSubMat_Nest" 713 PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat) 714 { 715 Mat_Nest *bA = (Mat_Nest*)A->data; 716 717 PetscFunctionBegin; 718 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 719 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 720 *mat = bA->m[idxm][jdxm]; 721 PetscFunctionReturn(0); 722 } 723 724 #undef __FUNCT__ 725 #define __FUNCT__ "MatNestGetSubMat" 726 /*@ 727 MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix. 728 729 Not collective 730 731 Input Parameters: 732 + A - nest matrix 733 . idxm - index of the matrix within the nest matrix 734 - jdxm - index of the matrix within the nest matrix 735 736 Output Parameter: 737 . sub - matrix at index idxm,jdxm within the nest matrix 738 739 Level: developer 740 741 .seealso: MatNestGetSize(), MatNestGetSubMats() 742 @*/ 743 PetscErrorCode MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub) 744 { 745 PetscErrorCode ierr; 746 747 PetscFunctionBegin; 748 ierr = PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));CHKERRQ(ierr); 749 PetscFunctionReturn(0); 750 } 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "MatNestSetSubMat_Nest" 754 PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat) 755 { 756 Mat_Nest *bA = (Mat_Nest*)A->data; 757 PetscInt m,n,M,N,mi,ni,Mi,Ni; 758 PetscErrorCode ierr; 759 760 PetscFunctionBegin; 761 if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1); 762 if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1); 763 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 764 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 765 ierr = ISGetLocalSize(bA->isglobal.row[idxm],&mi);CHKERRQ(ierr); 766 ierr = ISGetSize(bA->isglobal.row[idxm],&Mi);CHKERRQ(ierr); 767 ierr = ISGetLocalSize(bA->isglobal.col[jdxm],&ni);CHKERRQ(ierr); 768 ierr = ISGetSize(bA->isglobal.col[jdxm],&Ni);CHKERRQ(ierr); 769 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); 770 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); 771 772 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 773 ierr = MatDestroy(&bA->m[idxm][jdxm]);CHKERRQ(ierr); 774 bA->m[idxm][jdxm] = mat; 775 PetscFunctionReturn(0); 776 } 777 778 #undef __FUNCT__ 779 #define __FUNCT__ "MatNestSetSubMat" 780 /*@ 781 MatNestSetSubMat - Set a single submatrix in the nest matrix. 782 783 Logically collective on the submatrix communicator 784 785 Input Parameters: 786 + A - nest matrix 787 . idxm - index of the matrix within the nest matrix 788 . jdxm - index of the matrix within the nest matrix 789 - sub - matrix at index idxm,jdxm within the nest matrix 790 791 Notes: 792 The new submatrix must have the same size and communicator as that block of the nest. 793 794 This increments the reference count of the submatrix. 795 796 Level: developer 797 798 .seealso: MatNestSetSubMats(), MatNestGetSubMat() 799 @*/ 800 PetscErrorCode MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub) 801 { 802 PetscErrorCode ierr; 803 804 PetscFunctionBegin; 805 ierr = PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));CHKERRQ(ierr); 806 PetscFunctionReturn(0); 807 } 808 809 #undef __FUNCT__ 810 #define __FUNCT__ "MatNestGetSubMats_Nest" 811 PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 812 { 813 Mat_Nest *bA = (Mat_Nest*)A->data; 814 815 PetscFunctionBegin; 816 if (M) *M = bA->nr; 817 if (N) *N = bA->nc; 818 if (mat) *mat = bA->m; 819 PetscFunctionReturn(0); 820 } 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "MatNestGetSubMats" 824 /*@C 825 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix. 826 827 Not collective 828 829 Input Parameters: 830 . A - nest matrix 831 832 Output Parameter: 833 + M - number of rows in the nest matrix 834 . N - number of cols in the nest matrix 835 - mat - 2d array of matrices 836 837 Notes: 838 839 The user should not free the array mat. 840 841 Level: developer 842 843 .seealso: MatNestGetSize(), MatNestGetSubMat() 844 @*/ 845 PetscErrorCode MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat) 846 { 847 PetscErrorCode ierr; 848 849 PetscFunctionBegin; 850 ierr = PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));CHKERRQ(ierr); 851 PetscFunctionReturn(0); 852 } 853 854 #undef __FUNCT__ 855 #define __FUNCT__ "MatNestGetSize_Nest" 856 PetscErrorCode MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N) 857 { 858 Mat_Nest *bA = (Mat_Nest*)A->data; 859 860 PetscFunctionBegin; 861 if (M) *M = bA->nr; 862 if (N) *N = bA->nc; 863 PetscFunctionReturn(0); 864 } 865 866 #undef __FUNCT__ 867 #define __FUNCT__ "MatNestGetSize" 868 /*@ 869 MatNestGetSize - Returns the size of the nest matrix. 870 871 Not collective 872 873 Input Parameters: 874 . A - nest matrix 875 876 Output Parameter: 877 + M - number of rows in the nested mat 878 - N - number of cols in the nested mat 879 880 Notes: 881 882 Level: developer 883 884 .seealso: MatNestGetSubMat(), MatNestGetSubMats() 885 @*/ 886 PetscErrorCode MatNestGetSize(Mat A,PetscInt *M,PetscInt *N) 887 { 888 PetscErrorCode ierr; 889 890 PetscFunctionBegin; 891 ierr = PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));CHKERRQ(ierr); 892 PetscFunctionReturn(0); 893 } 894 895 #undef __FUNCT__ 896 #define __FUNCT__ "MatNestGetISs_Nest" 897 static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[]) 898 { 899 Mat_Nest *vs = (Mat_Nest*)A->data; 900 PetscInt i; 901 902 PetscFunctionBegin; 903 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i]; 904 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i]; 905 PetscFunctionReturn(0); 906 } 907 908 #undef __FUNCT__ 909 #define __FUNCT__ "MatNestGetISs" 910 /*@C 911 MatNestGetISs - Returns the index sets partitioning the row and column spaces 912 913 Not collective 914 915 Input Parameters: 916 . A - nest matrix 917 918 Output Parameter: 919 + rows - array of row index sets 920 - cols - array of column index sets 921 922 Level: advanced 923 924 Notes: 925 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 926 927 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs() 928 @*/ 929 PetscErrorCode MatNestGetISs(Mat A,IS rows[],IS cols[]) 930 { 931 PetscErrorCode ierr; 932 933 PetscFunctionBegin; 934 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 935 ierr = PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 936 PetscFunctionReturn(0); 937 } 938 939 #undef __FUNCT__ 940 #define __FUNCT__ "MatNestGetLocalISs_Nest" 941 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[]) 942 { 943 Mat_Nest *vs = (Mat_Nest*)A->data; 944 PetscInt i; 945 946 PetscFunctionBegin; 947 if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i]; 948 if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i]; 949 PetscFunctionReturn(0); 950 } 951 952 #undef __FUNCT__ 953 #define __FUNCT__ "MatNestGetLocalISs" 954 /*@C 955 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces 956 957 Not collective 958 959 Input Parameters: 960 . A - nest matrix 961 962 Output Parameter: 963 + rows - array of row index sets (or NULL to ignore) 964 - cols - array of column index sets (or NULL to ignore) 965 966 Level: advanced 967 968 Notes: 969 The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs. 970 971 .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs() 972 @*/ 973 PetscErrorCode MatNestGetLocalISs(Mat A,IS rows[],IS cols[]) 974 { 975 PetscErrorCode ierr; 976 977 PetscFunctionBegin; 978 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 979 ierr = PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));CHKERRQ(ierr); 980 PetscFunctionReturn(0); 981 } 982 983 #undef __FUNCT__ 984 #define __FUNCT__ "MatNestSetVecType_Nest" 985 PetscErrorCode MatNestSetVecType_Nest(Mat A,VecType vtype) 986 { 987 PetscErrorCode ierr; 988 PetscBool flg; 989 990 PetscFunctionBegin; 991 ierr = PetscStrcmp(vtype,VECNEST,&flg);CHKERRQ(ierr); 992 /* In reality, this only distinguishes VECNEST and "other" */ 993 if (flg) A->ops->getvecs = MatGetVecs_Nest; 994 else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0; 995 PetscFunctionReturn(0); 996 } 997 998 #undef __FUNCT__ 999 #define __FUNCT__ "MatNestSetVecType" 1000 /*@C 1001 MatNestSetVecType - Sets the type of Vec returned by MatGetVecs() 1002 1003 Not collective 1004 1005 Input Parameters: 1006 + A - nest matrix 1007 - vtype - type to use for creating vectors 1008 1009 Notes: 1010 1011 Level: developer 1012 1013 .seealso: MatGetVecs() 1014 @*/ 1015 PetscErrorCode MatNestSetVecType(Mat A,VecType vtype) 1016 { 1017 PetscErrorCode ierr; 1018 1019 PetscFunctionBegin; 1020 ierr = PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));CHKERRQ(ierr); 1021 PetscFunctionReturn(0); 1022 } 1023 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "MatNestSetSubMats_Nest" 1026 PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1027 { 1028 Mat_Nest *s = (Mat_Nest*)A->data; 1029 PetscInt i,j,m,n,M,N; 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 s->nr = nr; 1034 s->nc = nc; 1035 1036 /* Create space for submatrices */ 1037 ierr = PetscMalloc(sizeof(Mat*)*nr,&s->m);CHKERRQ(ierr); 1038 for (i=0; i<nr; i++) { 1039 ierr = PetscMalloc(sizeof(Mat)*nc,&s->m[i]);CHKERRQ(ierr); 1040 } 1041 for (i=0; i<nr; i++) { 1042 for (j=0; j<nc; j++) { 1043 s->m[i][j] = a[i*nc+j]; 1044 if (a[i*nc+j]) { 1045 ierr = PetscObjectReference((PetscObject)a[i*nc+j]);CHKERRQ(ierr); 1046 } 1047 } 1048 } 1049 1050 ierr = MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);CHKERRQ(ierr); 1051 1052 ierr = PetscMalloc(sizeof(PetscInt)*nr,&s->row_len);CHKERRQ(ierr); 1053 ierr = PetscMalloc(sizeof(PetscInt)*nc,&s->col_len);CHKERRQ(ierr); 1054 for (i=0; i<nr; i++) s->row_len[i]=-1; 1055 for (j=0; j<nc; j++) s->col_len[j]=-1; 1056 1057 ierr = MatNestGetSizes_Private(A,&m,&n,&M,&N);CHKERRQ(ierr); 1058 1059 ierr = PetscLayoutSetSize(A->rmap,M);CHKERRQ(ierr); 1060 ierr = PetscLayoutSetLocalSize(A->rmap,m);CHKERRQ(ierr); 1061 ierr = PetscLayoutSetSize(A->cmap,N);CHKERRQ(ierr); 1062 ierr = PetscLayoutSetLocalSize(A->cmap,n);CHKERRQ(ierr); 1063 1064 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1065 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1066 1067 ierr = PetscCalloc2(nr,&s->left,nc,&s->right);CHKERRQ(ierr); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "MatNestSetSubMats" 1073 /*@ 1074 MatNestSetSubMats - Sets the nested submatrices 1075 1076 Collective on Mat 1077 1078 Input Parameter: 1079 + N - nested matrix 1080 . nr - number of nested row blocks 1081 . is_row - index sets for each nested row block, or NULL to make contiguous 1082 . nc - number of nested column blocks 1083 . is_col - index sets for each nested column block, or NULL to make contiguous 1084 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1085 1086 Level: advanced 1087 1088 .seealso: MatCreateNest(), MATNEST 1089 @*/ 1090 PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[]) 1091 { 1092 PetscErrorCode ierr; 1093 PetscInt i; 1094 1095 PetscFunctionBegin; 1096 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1097 if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative"); 1098 if (nr && is_row) { 1099 PetscValidPointer(is_row,3); 1100 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_row[i],IS_CLASSID,3); 1101 } 1102 if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative"); 1103 if (nc && is_col) { 1104 PetscValidPointer(is_col,5); 1105 for (i=0; i<nr; i++) PetscValidHeaderSpecific(is_col[i],IS_CLASSID,5); 1106 } 1107 if (nr*nc) PetscValidPointer(a,6); 1108 ierr = PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));CHKERRQ(ierr); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 #undef __FUNCT__ 1113 #define __FUNCT__ "MatNestCreateAggregateL2G_Private" 1114 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog,ISLocalToGlobalMapping *ltogb) 1115 { 1116 PetscErrorCode ierr; 1117 PetscBool flg; 1118 PetscInt i,j,m,mi,*ix; 1119 1120 PetscFunctionBegin; 1121 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1122 if (islocal[i]) { 1123 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1124 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1125 } else { 1126 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1127 } 1128 m += mi; 1129 } 1130 if (flg) { 1131 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1132 for (i=0,n=0; i<n; i++) { 1133 ISLocalToGlobalMapping smap = NULL; 1134 VecScatter scat; 1135 IS isreq; 1136 Vec lvec,gvec; 1137 union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x; 1138 Mat sub; 1139 1140 if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers"); 1141 if (colflg) { 1142 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1143 } else { 1144 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1145 } 1146 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr);} 1147 if (islocal[i]) { 1148 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1149 } else { 1150 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1151 } 1152 for (j=0; j<mi; j++) ix[m+j] = j; 1153 if (smap) {ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);} 1154 /* 1155 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1156 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1157 The approach here is ugly because it uses VecScatter to move indices. 1158 */ 1159 ierr = VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);CHKERRQ(ierr); 1160 ierr = VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);CHKERRQ(ierr); 1161 ierr = ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);CHKERRQ(ierr); 1162 ierr = VecScatterCreate(gvec,isreq,lvec,NULL,&scat);CHKERRQ(ierr); 1163 ierr = VecGetArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1164 for (j=0; j<mi; j++) x[j].integer = ix[m+j]; 1165 ierr = VecRestoreArray(gvec,(PetscScalar**)&x);CHKERRQ(ierr); 1166 ierr = VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1167 ierr = VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1168 ierr = VecGetArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1169 for (j=0; j<mi; j++) ix[m+j] = x[j].integer; 1170 ierr = VecRestoreArray(lvec,(PetscScalar**)&x);CHKERRQ(ierr); 1171 ierr = VecDestroy(&lvec);CHKERRQ(ierr); 1172 ierr = VecDestroy(&gvec);CHKERRQ(ierr); 1173 ierr = ISDestroy(&isreq);CHKERRQ(ierr); 1174 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 1175 m += mi; 1176 } 1177 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1178 *ltogb = NULL; 1179 } else { 1180 *ltog = NULL; 1181 *ltogb = NULL; 1182 } 1183 PetscFunctionReturn(0); 1184 } 1185 1186 1187 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1188 /* 1189 nprocessors = NP 1190 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1191 proc 0: => (g_0,h_0,) 1192 proc 1: => (g_1,h_1,) 1193 ... 1194 proc nprocs-1: => (g_NP-1,h_NP-1,) 1195 1196 proc 0: proc 1: proc nprocs-1: 1197 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1198 1199 proc 0: 1200 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1201 proc 1: 1202 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1203 1204 proc NP-1: 1205 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1206 */ 1207 #undef __FUNCT__ 1208 #define __FUNCT__ "MatSetUp_NestIS_Private" 1209 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1210 { 1211 Mat_Nest *vs = (Mat_Nest*)A->data; 1212 PetscInt i,j,offset,n,nsum,bs; 1213 PetscErrorCode ierr; 1214 Mat sub = NULL; 1215 1216 PetscFunctionBegin; 1217 ierr = PetscMalloc(sizeof(IS)*nr,&vs->isglobal.row);CHKERRQ(ierr); 1218 ierr = PetscMalloc(sizeof(IS)*nc,&vs->isglobal.col);CHKERRQ(ierr); 1219 if (is_row) { /* valid IS is passed in */ 1220 /* refs on is[] are incremeneted */ 1221 for (i=0; i<vs->nr; i++) { 1222 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1223 1224 vs->isglobal.row[i] = is_row[i]; 1225 } 1226 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1227 nsum = 0; 1228 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1229 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1230 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1231 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1232 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1233 nsum += n; 1234 } 1235 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1236 offset -= nsum; 1237 for (i=0; i<vs->nr; i++) { 1238 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1239 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1240 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1241 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1242 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1243 offset += n; 1244 } 1245 } 1246 1247 if (is_col) { /* valid IS is passed in */ 1248 /* refs on is[] are incremeneted */ 1249 for (j=0; j<vs->nc; j++) { 1250 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1251 1252 vs->isglobal.col[j] = is_col[j]; 1253 } 1254 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1255 offset = A->cmap->rstart; 1256 nsum = 0; 1257 for (j=0; j<vs->nc; j++) { 1258 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1259 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1260 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1261 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1262 nsum += n; 1263 } 1264 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1265 offset -= nsum; 1266 for (j=0; j<vs->nc; j++) { 1267 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1268 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1269 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1270 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1271 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1272 offset += n; 1273 } 1274 } 1275 1276 /* Set up the local ISs */ 1277 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1278 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1279 for (i=0,offset=0; i<vs->nr; i++) { 1280 IS isloc; 1281 ISLocalToGlobalMapping rmap = NULL; 1282 PetscInt nlocal,bs; 1283 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1284 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1285 if (rmap) { 1286 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1287 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1288 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1289 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1290 } else { 1291 nlocal = 0; 1292 isloc = NULL; 1293 } 1294 vs->islocal.row[i] = isloc; 1295 offset += nlocal; 1296 } 1297 for (i=0,offset=0; i<vs->nc; i++) { 1298 IS isloc; 1299 ISLocalToGlobalMapping cmap = NULL; 1300 PetscInt nlocal,bs; 1301 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1302 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1303 if (cmap) { 1304 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1305 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1306 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1307 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1308 } else { 1309 nlocal = 0; 1310 isloc = NULL; 1311 } 1312 vs->islocal.col[i] = isloc; 1313 offset += nlocal; 1314 } 1315 1316 /* Set up the aggregate ISLocalToGlobalMapping */ 1317 { 1318 ISLocalToGlobalMapping rmap,rmapb,cmap,cmapb; 1319 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap,&rmapb);CHKERRQ(ierr); 1320 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap,&cmapb);CHKERRQ(ierr); 1321 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1322 if (rmapb && cmapb) {ierr = MatSetLocalToGlobalMappingBlock(A,rmapb,cmapb);CHKERRQ(ierr);} 1323 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1324 ierr = ISLocalToGlobalMappingDestroy(&rmapb);CHKERRQ(ierr); 1325 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1326 ierr = ISLocalToGlobalMappingDestroy(&cmapb);CHKERRQ(ierr); 1327 } 1328 1329 #if defined(PETSC_USE_DEBUG) 1330 for (i=0; i<vs->nr; i++) { 1331 for (j=0; j<vs->nc; j++) { 1332 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1333 Mat B = vs->m[i][j]; 1334 if (!B) continue; 1335 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1336 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1337 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1338 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1339 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1340 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1341 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); 1342 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); 1343 } 1344 } 1345 #endif 1346 1347 /* Set A->assembled if all non-null blocks are currently assembled */ 1348 for (i=0; i<vs->nr; i++) { 1349 for (j=0; j<vs->nc; j++) { 1350 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1351 } 1352 } 1353 A->assembled = PETSC_TRUE; 1354 PetscFunctionReturn(0); 1355 } 1356 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "MatCreateNest" 1359 /*@C 1360 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1361 1362 Collective on Mat 1363 1364 Input Parameter: 1365 + comm - Communicator for the new Mat 1366 . nr - number of nested row blocks 1367 . is_row - index sets for each nested row block, or NULL to make contiguous 1368 . nc - number of nested column blocks 1369 . is_col - index sets for each nested column block, or NULL to make contiguous 1370 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1371 1372 Output Parameter: 1373 . B - new matrix 1374 1375 Level: advanced 1376 1377 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST 1378 @*/ 1379 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1380 { 1381 Mat A; 1382 PetscErrorCode ierr; 1383 1384 PetscFunctionBegin; 1385 *B = 0; 1386 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1387 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1388 ierr = MatSetUp(A);CHKERRQ(ierr); 1389 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1390 *B = A; 1391 PetscFunctionReturn(0); 1392 } 1393 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "MatConvert_Nest_AIJ" 1396 PETSC_EXTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1397 { 1398 PetscErrorCode ierr; 1399 Mat_Nest *nest = (Mat_Nest*)A->data; 1400 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz; 1401 Mat C; 1402 1403 PetscFunctionBegin; 1404 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1405 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1406 switch (reuse) { 1407 case MAT_INITIAL_MATRIX: 1408 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1409 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1410 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1411 *newmat = C; 1412 break; 1413 case MAT_REUSE_MATRIX: 1414 C = *newmat; 1415 break; 1416 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1417 } 1418 1419 /* Preallocation */ 1420 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1421 onnz = dnnz + m; 1422 for (k=0; k<m; k++) { 1423 dnnz[k] = 0; 1424 onnz[k] = 0; 1425 } 1426 for (j=0; j<nest->nc; ++j) { 1427 IS bNis; 1428 PetscInt bN; 1429 const PetscInt *bNindices; 1430 /* Using global column indices and ISAllGather() is not scalable. */ 1431 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1432 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1433 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1434 for (i=0; i<nest->nr; ++i) { 1435 PetscSF bmsf; 1436 PetscSFNode *bmedges; 1437 Mat B; 1438 PetscInt bm, *bmdnnz, br; 1439 const PetscInt *bmindices; 1440 B = nest->m[i][j]; 1441 if (!B) continue; 1442 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1443 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1444 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1445 ierr = PetscMalloc1(2*bm,&bmedges);CHKERRQ(ierr); 1446 ierr = PetscMalloc1(2*bm,&bmdnnz);CHKERRQ(ierr); 1447 for (k = 0; k < 2*bm; ++k) bmdnnz[k] = 0; 1448 /* 1449 Locate the owners for all of the locally-owned global row indices for this row block. 1450 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1451 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1452 */ 1453 for (br = 0; br < bm; ++br) { 1454 PetscInt row = bmindices[br], rowowner = 0, brncols, col, colowner = 0; 1455 const PetscInt *brcols; 1456 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1457 PetscInt rowownerm; /* local row size on row's owning rank. */ 1458 1459 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1460 rowownerm = A->rmap->range[rowowner+1]-A->rmap->range[rowowner]; 1461 1462 bmedges[br].rank = rowowner; bmedges[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1463 bmedges[br].rank = rowowner; bmedges[br].index = rowrel+rowownerm; /* edge from bmonnz to onnz */ 1464 /* Now actually compute the data -- bmdnnz and bmonnz by looking at the global columns in the br row of this block. */ 1465 /* Note that this is not a pessimistic bound only because we assume the index sets embedding the blocks do not overlap. */ 1466 ierr = MatGetRow(B,br,&brncols,&brcols,NULL);CHKERRQ(ierr); 1467 for (k=0; k<brncols; k++) { 1468 col = bNindices[brcols[k]]; 1469 ierr = PetscLayoutFindOwnerIndex(A->cmap,col,&colowner,NULL);CHKERRQ(ierr); 1470 if (colowner == rowowner) bmdnnz[br]++; 1471 else onnz[br]++; 1472 } 1473 ierr = MatRestoreRow(B,br,&brncols,&brcols,NULL);CHKERRQ(ierr); 1474 } 1475 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1476 /* bsf will have to take care of disposing of bedges. */ 1477 ierr = PetscSFSetGraph(bmsf,m,2*bm,NULL,PETSC_COPY_VALUES,bmedges,PETSC_OWN_POINTER);CHKERRQ(ierr); 1478 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,bmdnnz,dnnz,MPIU_SUM);CHKERRQ(ierr); 1479 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,bmdnnz,dnnz,MPIU_SUM);CHKERRQ(ierr); 1480 ierr = PetscFree(bmdnnz);CHKERRQ(ierr); 1481 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1482 } 1483 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1484 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1485 } 1486 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1487 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1488 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1489 1490 /* Fill by row */ 1491 for (j=0; j<nest->nc; ++j) { 1492 /* Using global column indices and ISAllGather() is not scalable. */ 1493 IS bNis; 1494 PetscInt bN; 1495 const PetscInt *bNindices; 1496 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1497 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1498 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1499 for (i=0; i<nest->nr; ++i) { 1500 Mat B; 1501 PetscInt bm, br; 1502 const PetscInt *bmindices; 1503 B = nest->m[i][j]; 1504 if (!B) continue; 1505 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1506 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1507 for (br = 0; br < bm; ++br) { 1508 PetscInt row = bmindices[br], brncols, *cols; 1509 const PetscInt *brcols; 1510 const PetscScalar *brcoldata; 1511 ierr = MatGetRow(B,br,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1512 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1513 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1514 /* 1515 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1516 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1517 */ 1518 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1519 ierr = MatRestoreRow(B,br,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1520 ierr = PetscFree(cols);CHKERRQ(ierr); 1521 } 1522 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1523 } 1524 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1525 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1526 } 1527 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1528 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1529 PetscFunctionReturn(0); 1530 } 1531 1532 /*MC 1533 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1534 1535 Level: intermediate 1536 1537 Notes: 1538 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1539 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1540 It is usually used with DMComposite and DMCreateMatrix() 1541 1542 .seealso: MatCreate(), MatType, MatCreateNest() 1543 M*/ 1544 #undef __FUNCT__ 1545 #define __FUNCT__ "MatCreate_Nest" 1546 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1547 { 1548 Mat_Nest *s; 1549 PetscErrorCode ierr; 1550 1551 PetscFunctionBegin; 1552 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1553 A->data = (void*)s; 1554 1555 s->nr = -1; 1556 s->nc = -1; 1557 s->m = NULL; 1558 s->splitassembly = PETSC_FALSE; 1559 1560 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1561 1562 A->ops->mult = MatMult_Nest; 1563 A->ops->multadd = MatMultAdd_Nest; 1564 A->ops->multtranspose = MatMultTranspose_Nest; 1565 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1566 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1567 A->ops->assemblyend = MatAssemblyEnd_Nest; 1568 A->ops->zeroentries = MatZeroEntries_Nest; 1569 A->ops->copy = MatCopy_Nest; 1570 A->ops->duplicate = MatDuplicate_Nest; 1571 A->ops->getsubmatrix = MatGetSubMatrix_Nest; 1572 A->ops->destroy = MatDestroy_Nest; 1573 A->ops->view = MatView_Nest; 1574 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1575 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1576 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1577 A->ops->getdiagonal = MatGetDiagonal_Nest; 1578 A->ops->diagonalscale = MatDiagonalScale_Nest; 1579 A->ops->scale = MatScale_Nest; 1580 A->ops->shift = MatShift_Nest; 1581 1582 A->spptr = 0; 1583 A->assembled = PETSC_FALSE; 1584 1585 /* expose Nest api's */ 1586 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1587 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1588 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1589 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1590 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1591 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1592 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1593 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1594 1595 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1596 PetscFunctionReturn(0); 1597 } 1598