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",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 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 for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) { 1258 if (islocal[i]) { 1259 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1260 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1261 } else { 1262 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1263 } 1264 m += mi; 1265 } 1266 if (flg) { 1267 ierr = PetscMalloc1(m,&ix);CHKERRQ(ierr); 1268 for (i=0,m=0; i<n; i++) { 1269 ISLocalToGlobalMapping smap = NULL; 1270 Mat sub = NULL; 1271 PetscSF sf; 1272 PetscLayout map; 1273 PetscInt *ix2; 1274 1275 if (!colflg) { 1276 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1277 } else { 1278 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1279 } 1280 if (sub) { 1281 if (!colflg) { 1282 ierr = MatGetLocalToGlobalMapping(sub,&smap,NULL);CHKERRQ(ierr); 1283 } else { 1284 ierr = MatGetLocalToGlobalMapping(sub,NULL,&smap);CHKERRQ(ierr); 1285 } 1286 } 1287 if (islocal[i]) { 1288 ierr = ISGetSize(islocal[i],&mi);CHKERRQ(ierr); 1289 } else { 1290 ierr = ISGetSize(isglobal[i],&mi);CHKERRQ(ierr); 1291 } 1292 for (j=0; j<mi; j++) ix[m+j] = j; 1293 if (smap) {ierr = ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);CHKERRQ(ierr);} 1294 1295 /* 1296 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1297 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1298 */ 1299 ierr = PetscMalloc1(mi,&ix2);CHKERRQ(ierr); 1300 ierr = PetscSFCreate(((PetscObject)isglobal[i])->comm,&sf);CHKERRQ(ierr); 1301 ierr = PetscLayoutCreate(((PetscObject)isglobal[i])->comm,&map);CHKERRQ(ierr); 1302 ierr = PetscLayoutSetLocalSize(map,mi);CHKERRQ(ierr); 1303 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); 1304 ierr = PetscSFSetGraphLayout(sf,map,mi,NULL,PETSC_USE_POINTER,ix+m);CHKERRQ(ierr); 1305 ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr); 1306 for (j=0; j<mi; j++) ix2[j] = ix[m+j]; 1307 ierr = PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr); 1308 ierr = PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);CHKERRQ(ierr); 1309 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1310 ierr = PetscFree(ix2);CHKERRQ(ierr); 1311 m += mi; 1312 } 1313 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);CHKERRQ(ierr); 1314 } else { 1315 *ltog = NULL; 1316 } 1317 PetscFunctionReturn(0); 1318 } 1319 1320 1321 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1322 /* 1323 nprocessors = NP 1324 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1325 proc 0: => (g_0,h_0,) 1326 proc 1: => (g_1,h_1,) 1327 ... 1328 proc nprocs-1: => (g_NP-1,h_NP-1,) 1329 1330 proc 0: proc 1: proc nprocs-1: 1331 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1332 1333 proc 0: 1334 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1335 proc 1: 1336 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1337 1338 proc NP-1: 1339 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1340 */ 1341 static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[]) 1342 { 1343 Mat_Nest *vs = (Mat_Nest*)A->data; 1344 PetscInt i,j,offset,n,nsum,bs; 1345 PetscErrorCode ierr; 1346 Mat sub = NULL; 1347 1348 PetscFunctionBegin; 1349 ierr = PetscMalloc1(nr,&vs->isglobal.row);CHKERRQ(ierr); 1350 ierr = PetscMalloc1(nc,&vs->isglobal.col);CHKERRQ(ierr); 1351 if (is_row) { /* valid IS is passed in */ 1352 /* refs on is[] are incremeneted */ 1353 for (i=0; i<vs->nr; i++) { 1354 ierr = PetscObjectReference((PetscObject)is_row[i]);CHKERRQ(ierr); 1355 1356 vs->isglobal.row[i] = is_row[i]; 1357 } 1358 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1359 nsum = 0; 1360 for (i=0; i<vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1361 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1362 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i); 1363 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1364 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1365 nsum += n; 1366 } 1367 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1368 offset -= nsum; 1369 for (i=0; i<vs->nr; i++) { 1370 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1371 ierr = MatGetLocalSize(sub,&n,NULL);CHKERRQ(ierr); 1372 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1373 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);CHKERRQ(ierr); 1374 ierr = ISSetBlockSize(vs->isglobal.row[i],bs);CHKERRQ(ierr); 1375 offset += n; 1376 } 1377 } 1378 1379 if (is_col) { /* valid IS is passed in */ 1380 /* refs on is[] are incremeneted */ 1381 for (j=0; j<vs->nc; j++) { 1382 ierr = PetscObjectReference((PetscObject)is_col[j]);CHKERRQ(ierr); 1383 1384 vs->isglobal.col[j] = is_col[j]; 1385 } 1386 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1387 offset = A->cmap->rstart; 1388 nsum = 0; 1389 for (j=0; j<vs->nc; j++) { 1390 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1391 if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i); 1392 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1393 if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix"); 1394 nsum += n; 1395 } 1396 ierr = MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1397 offset -= nsum; 1398 for (j=0; j<vs->nc; j++) { 1399 ierr = MatNestFindNonzeroSubMatCol(A,j,&sub);CHKERRQ(ierr); 1400 ierr = MatGetLocalSize(sub,NULL,&n);CHKERRQ(ierr); 1401 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1402 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);CHKERRQ(ierr); 1403 ierr = ISSetBlockSize(vs->isglobal.col[j],bs);CHKERRQ(ierr); 1404 offset += n; 1405 } 1406 } 1407 1408 /* Set up the local ISs */ 1409 ierr = PetscMalloc1(vs->nr,&vs->islocal.row);CHKERRQ(ierr); 1410 ierr = PetscMalloc1(vs->nc,&vs->islocal.col);CHKERRQ(ierr); 1411 for (i=0,offset=0; i<vs->nr; i++) { 1412 IS isloc; 1413 ISLocalToGlobalMapping rmap = NULL; 1414 PetscInt nlocal,bs; 1415 ierr = MatNestFindNonzeroSubMatRow(A,i,&sub);CHKERRQ(ierr); 1416 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,&rmap,NULL);CHKERRQ(ierr);} 1417 if (rmap) { 1418 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1419 ierr = ISLocalToGlobalMappingGetSize(rmap,&nlocal);CHKERRQ(ierr); 1420 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1421 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1422 } else { 1423 nlocal = 0; 1424 isloc = NULL; 1425 } 1426 vs->islocal.row[i] = isloc; 1427 offset += nlocal; 1428 } 1429 for (i=0,offset=0; i<vs->nc; i++) { 1430 IS isloc; 1431 ISLocalToGlobalMapping cmap = NULL; 1432 PetscInt nlocal,bs; 1433 ierr = MatNestFindNonzeroSubMatCol(A,i,&sub);CHKERRQ(ierr); 1434 if (sub) {ierr = MatGetLocalToGlobalMapping(sub,NULL,&cmap);CHKERRQ(ierr);} 1435 if (cmap) { 1436 ierr = MatGetBlockSize(sub,&bs);CHKERRQ(ierr); 1437 ierr = ISLocalToGlobalMappingGetSize(cmap,&nlocal);CHKERRQ(ierr); 1438 ierr = ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);CHKERRQ(ierr); 1439 ierr = ISSetBlockSize(isloc,bs);CHKERRQ(ierr); 1440 } else { 1441 nlocal = 0; 1442 isloc = NULL; 1443 } 1444 vs->islocal.col[i] = isloc; 1445 offset += nlocal; 1446 } 1447 1448 /* Set up the aggregate ISLocalToGlobalMapping */ 1449 { 1450 ISLocalToGlobalMapping rmap,cmap; 1451 ierr = MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);CHKERRQ(ierr); 1452 ierr = MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);CHKERRQ(ierr); 1453 if (rmap && cmap) {ierr = MatSetLocalToGlobalMapping(A,rmap,cmap);CHKERRQ(ierr);} 1454 ierr = ISLocalToGlobalMappingDestroy(&rmap);CHKERRQ(ierr); 1455 ierr = ISLocalToGlobalMappingDestroy(&cmap);CHKERRQ(ierr); 1456 } 1457 1458 #if defined(PETSC_USE_DEBUG) 1459 for (i=0; i<vs->nr; i++) { 1460 for (j=0; j<vs->nc; j++) { 1461 PetscInt m,n,M,N,mi,ni,Mi,Ni; 1462 Mat B = vs->m[i][j]; 1463 if (!B) continue; 1464 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 1465 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 1466 ierr = ISGetSize(vs->isglobal.row[i],&Mi);CHKERRQ(ierr); 1467 ierr = ISGetSize(vs->isglobal.col[j],&Ni);CHKERRQ(ierr); 1468 ierr = ISGetLocalSize(vs->isglobal.row[i],&mi);CHKERRQ(ierr); 1469 ierr = ISGetLocalSize(vs->isglobal.col[j],&ni);CHKERRQ(ierr); 1470 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); 1471 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); 1472 } 1473 } 1474 #endif 1475 1476 /* Set A->assembled if all non-null blocks are currently assembled */ 1477 for (i=0; i<vs->nr; i++) { 1478 for (j=0; j<vs->nc; j++) { 1479 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(0); 1480 } 1481 } 1482 A->assembled = PETSC_TRUE; 1483 PetscFunctionReturn(0); 1484 } 1485 1486 /*@C 1487 MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately 1488 1489 Collective on Mat 1490 1491 Input Parameter: 1492 + comm - Communicator for the new Mat 1493 . nr - number of nested row blocks 1494 . is_row - index sets for each nested row block, or NULL to make contiguous 1495 . nc - number of nested column blocks 1496 . is_col - index sets for each nested column block, or NULL to make contiguous 1497 - a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL 1498 1499 Output Parameter: 1500 . B - new matrix 1501 1502 Level: advanced 1503 1504 .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST 1505 @*/ 1506 PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B) 1507 { 1508 Mat A; 1509 PetscErrorCode ierr; 1510 1511 PetscFunctionBegin; 1512 *B = 0; 1513 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 1514 ierr = MatSetType(A,MATNEST);CHKERRQ(ierr); 1515 A->preallocated = PETSC_TRUE; 1516 ierr = MatNestSetSubMats(A,nr,is_row,nc,is_col,a);CHKERRQ(ierr); 1517 *B = A; 1518 PetscFunctionReturn(0); 1519 } 1520 1521 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1522 { 1523 Mat_Nest *nest = (Mat_Nest*)A->data; 1524 Mat *trans; 1525 PetscScalar **avv; 1526 PetscScalar *vv; 1527 PetscInt **aii,**ajj; 1528 PetscInt *ii,*jj,*ci; 1529 PetscInt nr,nc,nnz,i,j; 1530 PetscBool done; 1531 PetscErrorCode ierr; 1532 1533 PetscFunctionBegin; 1534 ierr = MatGetSize(A,&nr,&nc);CHKERRQ(ierr); 1535 if (reuse == MAT_REUSE_MATRIX) { 1536 PetscInt rnr; 1537 1538 ierr = MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);CHKERRQ(ierr); 1539 if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ"); 1540 if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows"); 1541 ierr = MatSeqAIJGetArray(*newmat,&vv);CHKERRQ(ierr); 1542 } 1543 /* extract CSR for nested SeqAIJ matrices */ 1544 nnz = 0; 1545 ierr = PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);CHKERRQ(ierr); 1546 for (i=0; i<nest->nr; ++i) { 1547 for (j=0; j<nest->nc; ++j) { 1548 Mat B = nest->m[i][j]; 1549 if (B) { 1550 PetscScalar *naa; 1551 PetscInt *nii,*njj,nnr; 1552 PetscBool istrans; 1553 1554 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1555 if (istrans) { 1556 Mat Bt; 1557 1558 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1559 ierr = MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);CHKERRQ(ierr); 1560 B = trans[i*nest->nc+j]; 1561 } 1562 ierr = MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);CHKERRQ(ierr); 1563 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ"); 1564 ierr = MatSeqAIJGetArray(B,&naa);CHKERRQ(ierr); 1565 nnz += nii[nnr]; 1566 1567 aii[i*nest->nc+j] = nii; 1568 ajj[i*nest->nc+j] = njj; 1569 avv[i*nest->nc+j] = naa; 1570 } 1571 } 1572 } 1573 if (reuse != MAT_REUSE_MATRIX) { 1574 ierr = PetscMalloc1(nr+1,&ii);CHKERRQ(ierr); 1575 ierr = PetscMalloc1(nnz,&jj);CHKERRQ(ierr); 1576 ierr = PetscMalloc1(nnz,&vv);CHKERRQ(ierr); 1577 } else { 1578 if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros"); 1579 } 1580 1581 /* new row pointer */ 1582 ierr = PetscArrayzero(ii,nr+1);CHKERRQ(ierr); 1583 for (i=0; i<nest->nr; ++i) { 1584 PetscInt ncr,rst; 1585 1586 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1587 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1588 for (j=0; j<nest->nc; ++j) { 1589 if (aii[i*nest->nc+j]) { 1590 PetscInt *nii = aii[i*nest->nc+j]; 1591 PetscInt ir; 1592 1593 for (ir=rst; ir<ncr+rst; ++ir) { 1594 ii[ir+1] += nii[1]-nii[0]; 1595 nii++; 1596 } 1597 } 1598 } 1599 } 1600 for (i=0; i<nr; i++) ii[i+1] += ii[i]; 1601 1602 /* construct CSR for the new matrix */ 1603 ierr = PetscCalloc1(nr,&ci);CHKERRQ(ierr); 1604 for (i=0; i<nest->nr; ++i) { 1605 PetscInt ncr,rst; 1606 1607 ierr = ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);CHKERRQ(ierr); 1608 ierr = ISGetLocalSize(nest->isglobal.row[i],&ncr);CHKERRQ(ierr); 1609 for (j=0; j<nest->nc; ++j) { 1610 if (aii[i*nest->nc+j]) { 1611 PetscScalar *nvv = avv[i*nest->nc+j]; 1612 PetscInt *nii = aii[i*nest->nc+j]; 1613 PetscInt *njj = ajj[i*nest->nc+j]; 1614 PetscInt ir,cst; 1615 1616 ierr = ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);CHKERRQ(ierr); 1617 for (ir=rst; ir<ncr+rst; ++ir) { 1618 PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir]; 1619 1620 for (ij=0;ij<rsize;ij++) { 1621 jj[ist+ij] = *njj+cst; 1622 vv[ist+ij] = *nvv; 1623 njj++; 1624 nvv++; 1625 } 1626 ci[ir] += rsize; 1627 nii++; 1628 } 1629 } 1630 } 1631 } 1632 ierr = PetscFree(ci);CHKERRQ(ierr); 1633 1634 /* restore info */ 1635 for (i=0; i<nest->nr; ++i) { 1636 for (j=0; j<nest->nc; ++j) { 1637 Mat B = nest->m[i][j]; 1638 if (B) { 1639 PetscInt nnr = 0, k = i*nest->nc+j; 1640 1641 B = (trans[k] ? trans[k] : B); 1642 ierr = MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);CHKERRQ(ierr); 1643 if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ"); 1644 ierr = MatSeqAIJRestoreArray(B,&avv[k]);CHKERRQ(ierr); 1645 ierr = MatDestroy(&trans[k]);CHKERRQ(ierr); 1646 } 1647 } 1648 } 1649 ierr = PetscFree4(aii,ajj,avv,trans);CHKERRQ(ierr); 1650 1651 /* finalize newmat */ 1652 if (reuse == MAT_INITIAL_MATRIX) { 1653 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);CHKERRQ(ierr); 1654 } else if (reuse == MAT_INPLACE_MATRIX) { 1655 Mat B; 1656 1657 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);CHKERRQ(ierr); 1658 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1659 } 1660 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1661 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1662 { 1663 Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data); 1664 a->free_a = PETSC_TRUE; 1665 a->free_ij = PETSC_TRUE; 1666 } 1667 PetscFunctionReturn(0); 1668 } 1669 1670 PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1671 { 1672 PetscErrorCode ierr; 1673 Mat_Nest *nest = (Mat_Nest*)A->data; 1674 PetscInt m,n,M,N,i,j,k,*dnnz,*onnz,rstart; 1675 PetscInt cstart,cend; 1676 PetscMPIInt size; 1677 Mat C; 1678 1679 PetscFunctionBegin; 1680 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1681 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1682 PetscInt nf; 1683 PetscBool fast; 1684 1685 ierr = PetscStrcmp(newtype,MATAIJ,&fast);CHKERRQ(ierr); 1686 if (!fast) { 1687 ierr = PetscStrcmp(newtype,MATSEQAIJ,&fast);CHKERRQ(ierr); 1688 } 1689 for (i=0; i<nest->nr && fast; ++i) { 1690 for (j=0; j<nest->nc && fast; ++j) { 1691 Mat B = nest->m[i][j]; 1692 if (B) { 1693 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);CHKERRQ(ierr); 1694 if (!fast) { 1695 PetscBool istrans; 1696 1697 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 1698 if (istrans) { 1699 Mat Bt; 1700 1701 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 1702 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);CHKERRQ(ierr); 1703 } 1704 } 1705 } 1706 } 1707 } 1708 for (i=0, nf=0; i<nest->nr && fast; ++i) { 1709 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1710 if (fast) { 1711 PetscInt f,s; 1712 1713 ierr = ISStrideGetInfo(nest->isglobal.row[i],&f,&s);CHKERRQ(ierr); 1714 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1715 else { 1716 ierr = ISGetSize(nest->isglobal.row[i],&f);CHKERRQ(ierr); 1717 nf += f; 1718 } 1719 } 1720 } 1721 for (i=0, nf=0; i<nest->nc && fast; ++i) { 1722 ierr = PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);CHKERRQ(ierr); 1723 if (fast) { 1724 PetscInt f,s; 1725 1726 ierr = ISStrideGetInfo(nest->isglobal.col[i],&f,&s);CHKERRQ(ierr); 1727 if (f != nf || s != 1) { fast = PETSC_FALSE; } 1728 else { 1729 ierr = ISGetSize(nest->isglobal.col[i],&f);CHKERRQ(ierr); 1730 nf += f; 1731 } 1732 } 1733 } 1734 if (fast) { 1735 ierr = MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);CHKERRQ(ierr); 1736 PetscFunctionReturn(0); 1737 } 1738 } 1739 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1740 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 1741 ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr); 1742 switch (reuse) { 1743 case MAT_INITIAL_MATRIX: 1744 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1745 ierr = MatSetType(C,newtype);CHKERRQ(ierr); 1746 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 1747 *newmat = C; 1748 break; 1749 case MAT_REUSE_MATRIX: 1750 C = *newmat; 1751 break; 1752 default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse"); 1753 } 1754 ierr = PetscMalloc1(2*m,&dnnz);CHKERRQ(ierr); 1755 onnz = dnnz + m; 1756 for (k=0; k<m; k++) { 1757 dnnz[k] = 0; 1758 onnz[k] = 0; 1759 } 1760 for (j=0; j<nest->nc; ++j) { 1761 IS bNis; 1762 PetscInt bN; 1763 const PetscInt *bNindices; 1764 /* Using global column indices and ISAllGather() is not scalable. */ 1765 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1766 ierr = ISGetSize(bNis, &bN);CHKERRQ(ierr); 1767 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1768 for (i=0; i<nest->nr; ++i) { 1769 PetscSF bmsf; 1770 PetscSFNode *iremote; 1771 Mat B; 1772 PetscInt bm, *sub_dnnz,*sub_onnz, br; 1773 const PetscInt *bmindices; 1774 B = nest->m[i][j]; 1775 if (!B) continue; 1776 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1777 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1778 ierr = PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);CHKERRQ(ierr); 1779 ierr = PetscMalloc1(bm,&iremote);CHKERRQ(ierr); 1780 ierr = PetscMalloc1(bm,&sub_dnnz);CHKERRQ(ierr); 1781 ierr = PetscMalloc1(bm,&sub_onnz);CHKERRQ(ierr); 1782 for (k = 0; k < bm; ++k){ 1783 sub_dnnz[k] = 0; 1784 sub_onnz[k] = 0; 1785 } 1786 /* 1787 Locate the owners for all of the locally-owned global row indices for this row block. 1788 These determine the roots of PetscSF used to communicate preallocation data to row owners. 1789 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 1790 */ 1791 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1792 for (br = 0; br < bm; ++br) { 1793 PetscInt row = bmindices[br], rowowner = 0, brncols, col; 1794 const PetscInt *brcols; 1795 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 1796 ierr = PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);CHKERRQ(ierr); 1797 /* how many roots */ 1798 iremote[br].rank = rowowner; iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 1799 /* get nonzero pattern */ 1800 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1801 for (k=0; k<brncols; k++) { 1802 col = bNindices[brcols[k]]; 1803 if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) { 1804 sub_dnnz[br]++; 1805 } else { 1806 sub_onnz[br]++; 1807 } 1808 } 1809 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);CHKERRQ(ierr); 1810 } 1811 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1812 /* bsf will have to take care of disposing of bedges. */ 1813 ierr = PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr); 1814 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1815 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);CHKERRQ(ierr); 1816 ierr = PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1817 ierr = PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);CHKERRQ(ierr); 1818 ierr = PetscFree(sub_dnnz);CHKERRQ(ierr); 1819 ierr = PetscFree(sub_onnz);CHKERRQ(ierr); 1820 ierr = PetscSFDestroy(&bmsf);CHKERRQ(ierr); 1821 } 1822 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1823 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1824 } 1825 /* Resize preallocation if overestimated */ 1826 for (i=0;i<m;i++) { 1827 dnnz[i] = PetscMin(dnnz[i],A->cmap->n); 1828 onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n); 1829 } 1830 ierr = MatSeqAIJSetPreallocation(C,0,dnnz);CHKERRQ(ierr); 1831 ierr = MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);CHKERRQ(ierr); 1832 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1833 1834 /* Fill by row */ 1835 for (j=0; j<nest->nc; ++j) { 1836 /* Using global column indices and ISAllGather() is not scalable. */ 1837 IS bNis; 1838 PetscInt bN; 1839 const PetscInt *bNindices; 1840 ierr = ISAllGather(nest->isglobal.col[j], &bNis);CHKERRQ(ierr); 1841 ierr = ISGetSize(bNis,&bN);CHKERRQ(ierr); 1842 ierr = ISGetIndices(bNis,&bNindices);CHKERRQ(ierr); 1843 for (i=0; i<nest->nr; ++i) { 1844 Mat B; 1845 PetscInt bm, br; 1846 const PetscInt *bmindices; 1847 B = nest->m[i][j]; 1848 if (!B) continue; 1849 ierr = ISGetLocalSize(nest->isglobal.row[i],&bm);CHKERRQ(ierr); 1850 ierr = ISGetIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1851 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 1852 for (br = 0; br < bm; ++br) { 1853 PetscInt row = bmindices[br], brncols, *cols; 1854 const PetscInt *brcols; 1855 const PetscScalar *brcoldata; 1856 ierr = MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1857 ierr = PetscMalloc1(brncols,&cols);CHKERRQ(ierr); 1858 for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]]; 1859 /* 1860 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1861 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1862 */ 1863 ierr = MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);CHKERRQ(ierr); 1864 ierr = MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);CHKERRQ(ierr); 1865 ierr = PetscFree(cols);CHKERRQ(ierr); 1866 } 1867 ierr = ISRestoreIndices(nest->isglobal.row[i],&bmindices);CHKERRQ(ierr); 1868 } 1869 ierr = ISRestoreIndices(bNis,&bNindices);CHKERRQ(ierr); 1870 ierr = ISDestroy(&bNis);CHKERRQ(ierr); 1871 } 1872 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1873 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1874 PetscFunctionReturn(0); 1875 } 1876 1877 PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has) 1878 { 1879 PetscFunctionBegin; 1880 *has = PETSC_FALSE; 1881 if (op == MATOP_MULT_TRANSPOSE) { 1882 Mat_Nest *bA = (Mat_Nest*)mat->data; 1883 PetscInt i,j,nr = bA->nr,nc = bA->nc; 1884 PetscErrorCode ierr; 1885 PetscBool flg; 1886 1887 for (j=0; j<nc; j++) { 1888 for (i=0; i<nr; i++) { 1889 if (!bA->m[i][j]) continue; 1890 ierr = MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);CHKERRQ(ierr); 1891 if (!flg) PetscFunctionReturn(0); 1892 } 1893 } 1894 } 1895 if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 1896 PetscFunctionReturn(0); 1897 } 1898 1899 /*MC 1900 MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately. 1901 1902 Level: intermediate 1903 1904 Notes: 1905 This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices. 1906 It allows the use of symmetric and block formats for parts of multi-physics simulations. 1907 It is usually used with DMComposite and DMCreateMatrix() 1908 1909 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 1910 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 1911 than the nest matrix. 1912 1913 .seealso: MatCreate(), MatType, MatCreateNest() 1914 M*/ 1915 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 1916 { 1917 Mat_Nest *s; 1918 PetscErrorCode ierr; 1919 1920 PetscFunctionBegin; 1921 ierr = PetscNewLog(A,&s);CHKERRQ(ierr); 1922 A->data = (void*)s; 1923 1924 s->nr = -1; 1925 s->nc = -1; 1926 s->m = NULL; 1927 s->splitassembly = PETSC_FALSE; 1928 1929 ierr = PetscMemzero(A->ops,sizeof(*A->ops));CHKERRQ(ierr); 1930 1931 A->ops->mult = MatMult_Nest; 1932 A->ops->multadd = MatMultAdd_Nest; 1933 A->ops->multtranspose = MatMultTranspose_Nest; 1934 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 1935 A->ops->transpose = MatTranspose_Nest; 1936 A->ops->assemblybegin = MatAssemblyBegin_Nest; 1937 A->ops->assemblyend = MatAssemblyEnd_Nest; 1938 A->ops->zeroentries = MatZeroEntries_Nest; 1939 A->ops->copy = MatCopy_Nest; 1940 A->ops->axpy = MatAXPY_Nest; 1941 A->ops->duplicate = MatDuplicate_Nest; 1942 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 1943 A->ops->destroy = MatDestroy_Nest; 1944 A->ops->view = MatView_Nest; 1945 A->ops->getvecs = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 1946 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 1947 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 1948 A->ops->getdiagonal = MatGetDiagonal_Nest; 1949 A->ops->diagonalscale = MatDiagonalScale_Nest; 1950 A->ops->scale = MatScale_Nest; 1951 A->ops->shift = MatShift_Nest; 1952 A->ops->diagonalset = MatDiagonalSet_Nest; 1953 A->ops->setrandom = MatSetRandom_Nest; 1954 A->ops->hasoperation = MatHasOperation_Nest; 1955 1956 A->spptr = 0; 1957 A->assembled = PETSC_FALSE; 1958 1959 /* expose Nest api's */ 1960 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C", MatNestGetSubMat_Nest);CHKERRQ(ierr); 1961 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C", MatNestSetSubMat_Nest);CHKERRQ(ierr); 1962 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C", MatNestGetSubMats_Nest);CHKERRQ(ierr); 1963 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C", MatNestGetSize_Nest);CHKERRQ(ierr); 1964 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C", MatNestGetISs_Nest);CHKERRQ(ierr); 1965 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);CHKERRQ(ierr); 1966 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C", MatNestSetVecType_Nest);CHKERRQ(ierr); 1967 ierr = PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C", MatNestSetSubMats_Nest);CHKERRQ(ierr); 1968 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1969 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",MatConvert_Nest_AIJ);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C", MatConvert_Nest_AIJ);CHKERRQ(ierr); 1971 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);CHKERRQ(ierr); 1972 1973 ierr = PetscObjectChangeTypeName((PetscObject)A,MATNEST);CHKERRQ(ierr); 1974 PetscFunctionReturn(0); 1975 } 1976