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