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