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