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