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