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