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