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