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