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