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