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