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