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