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