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