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