1 2 /* 3 This provides a matrix that consists of Mats 4 */ 5 6 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/impls/baij/seq/baij.h> /* use the common AIJ data-structure */ 8 9 typedef struct { 10 SEQAIJHEADER(Mat); 11 SEQBAIJHEADER; 12 Mat *diags; 13 14 Vec left,right,middle,workb; /* dummy vectors to perform local parts of product */ 15 } Mat_BlockMat; 16 17 static PetscErrorCode MatSOR_BlockMat_Symmetric(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 18 { 19 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 20 PetscScalar *x; 21 const Mat *v; 22 const PetscScalar *b; 23 PetscErrorCode ierr; 24 PetscInt n = A->cmap->n,i,mbs = n/A->rmap->bs,j,bs = A->rmap->bs; 25 const PetscInt *idx; 26 IS row,col; 27 MatFactorInfo info; 28 Vec left = a->left,right = a->right, middle = a->middle; 29 Mat *diag; 30 31 PetscFunctionBegin; 32 its = its*lits; 33 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 34 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat"); 35 if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for omega not equal to 1.0"); 36 if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for fshift"); 37 if ((flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) && !(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)) { 38 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot do backward sweep without forward sweep"); 39 } 40 41 if (!a->diags) { 42 ierr = PetscMalloc1(mbs,&a->diags);CHKERRQ(ierr); 43 ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr); 44 for (i=0; i<mbs; i++) { 45 ierr = MatGetOrdering(a->a[a->diag[i]], MATORDERINGND,&row,&col);CHKERRQ(ierr); 46 ierr = MatCholeskyFactorSymbolic(a->diags[i],a->a[a->diag[i]],row,&info);CHKERRQ(ierr); 47 ierr = MatCholeskyFactorNumeric(a->diags[i],a->a[a->diag[i]],&info);CHKERRQ(ierr); 48 ierr = ISDestroy(&row);CHKERRQ(ierr); 49 ierr = ISDestroy(&col);CHKERRQ(ierr); 50 } 51 ierr = VecDuplicate(bb,&a->workb);CHKERRQ(ierr); 52 } 53 diag = a->diags; 54 55 ierr = VecSet(xx,0.0);CHKERRQ(ierr); 56 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 57 /* copy right hand side because it must be modified during iteration */ 58 ierr = VecCopy(bb,a->workb);CHKERRQ(ierr); 59 ierr = VecGetArrayRead(a->workb,&b);CHKERRQ(ierr); 60 61 /* need to add code for when initial guess is zero, see MatSOR_SeqAIJ */ 62 while (its--) { 63 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 64 65 for (i=0; i<mbs; i++) { 66 n = a->i[i+1] - a->i[i] - 1; 67 idx = a->j + a->i[i] + 1; 68 v = a->a + a->i[i] + 1; 69 70 ierr = VecSet(left,0.0);CHKERRQ(ierr); 71 for (j=0; j<n; j++) { 72 ierr = VecPlaceArray(right,x + idx[j]*bs);CHKERRQ(ierr); 73 ierr = MatMultAdd(v[j],right,left,left);CHKERRQ(ierr); 74 ierr = VecResetArray(right);CHKERRQ(ierr); 75 } 76 ierr = VecPlaceArray(right,b + i*bs);CHKERRQ(ierr); 77 ierr = VecAYPX(left,-1.0,right);CHKERRQ(ierr); 78 ierr = VecResetArray(right);CHKERRQ(ierr); 79 80 ierr = VecPlaceArray(right,x + i*bs);CHKERRQ(ierr); 81 ierr = MatSolve(diag[i],left,right);CHKERRQ(ierr); 82 83 /* now adjust right hand side, see MatSOR_SeqSBAIJ */ 84 for (j=0; j<n; j++) { 85 ierr = MatMultTranspose(v[j],right,left);CHKERRQ(ierr); 86 ierr = VecPlaceArray(middle,b + idx[j]*bs);CHKERRQ(ierr); 87 ierr = VecAXPY(middle,-1.0,left);CHKERRQ(ierr); 88 ierr = VecResetArray(middle);CHKERRQ(ierr); 89 } 90 ierr = VecResetArray(right);CHKERRQ(ierr); 91 92 } 93 } 94 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 95 96 for (i=mbs-1; i>=0; i--) { 97 n = a->i[i+1] - a->i[i] - 1; 98 idx = a->j + a->i[i] + 1; 99 v = a->a + a->i[i] + 1; 100 101 ierr = VecSet(left,0.0);CHKERRQ(ierr); 102 for (j=0; j<n; j++) { 103 ierr = VecPlaceArray(right,x + idx[j]*bs);CHKERRQ(ierr); 104 ierr = MatMultAdd(v[j],right,left,left);CHKERRQ(ierr); 105 ierr = VecResetArray(right);CHKERRQ(ierr); 106 } 107 ierr = VecPlaceArray(right,b + i*bs);CHKERRQ(ierr); 108 ierr = VecAYPX(left,-1.0,right);CHKERRQ(ierr); 109 ierr = VecResetArray(right);CHKERRQ(ierr); 110 111 ierr = VecPlaceArray(right,x + i*bs);CHKERRQ(ierr); 112 ierr = MatSolve(diag[i],left,right);CHKERRQ(ierr); 113 ierr = VecResetArray(right);CHKERRQ(ierr); 114 115 } 116 } 117 } 118 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 119 ierr = VecRestoreArrayRead(a->workb,&b);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 static PetscErrorCode MatSOR_BlockMat(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 124 { 125 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 126 PetscScalar *x; 127 const Mat *v; 128 const PetscScalar *b; 129 PetscErrorCode ierr; 130 PetscInt n = A->cmap->n,i,mbs = n/A->rmap->bs,j,bs = A->rmap->bs; 131 const PetscInt *idx; 132 IS row,col; 133 MatFactorInfo info; 134 Vec left = a->left,right = a->right; 135 Mat *diag; 136 137 PetscFunctionBegin; 138 its = its*lits; 139 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 140 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat"); 141 if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for omega not equal to 1.0"); 142 if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for fshift"); 143 144 if (!a->diags) { 145 ierr = PetscMalloc1(mbs,&a->diags);CHKERRQ(ierr); 146 ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr); 147 for (i=0; i<mbs; i++) { 148 ierr = MatGetOrdering(a->a[a->diag[i]], MATORDERINGND,&row,&col);CHKERRQ(ierr); 149 ierr = MatLUFactorSymbolic(a->diags[i],a->a[a->diag[i]],row,col,&info);CHKERRQ(ierr); 150 ierr = MatLUFactorNumeric(a->diags[i],a->a[a->diag[i]],&info);CHKERRQ(ierr); 151 ierr = ISDestroy(&row);CHKERRQ(ierr); 152 ierr = ISDestroy(&col);CHKERRQ(ierr); 153 } 154 } 155 diag = a->diags; 156 157 ierr = VecSet(xx,0.0);CHKERRQ(ierr); 158 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 159 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 160 161 /* need to add code for when initial guess is zero, see MatSOR_SeqAIJ */ 162 while (its--) { 163 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 164 165 for (i=0; i<mbs; i++) { 166 n = a->i[i+1] - a->i[i]; 167 idx = a->j + a->i[i]; 168 v = a->a + a->i[i]; 169 170 ierr = VecSet(left,0.0);CHKERRQ(ierr); 171 for (j=0; j<n; j++) { 172 if (idx[j] != i) { 173 ierr = VecPlaceArray(right,x + idx[j]*bs);CHKERRQ(ierr); 174 ierr = MatMultAdd(v[j],right,left,left);CHKERRQ(ierr); 175 ierr = VecResetArray(right);CHKERRQ(ierr); 176 } 177 } 178 ierr = VecPlaceArray(right,b + i*bs);CHKERRQ(ierr); 179 ierr = VecAYPX(left,-1.0,right);CHKERRQ(ierr); 180 ierr = VecResetArray(right);CHKERRQ(ierr); 181 182 ierr = VecPlaceArray(right,x + i*bs);CHKERRQ(ierr); 183 ierr = MatSolve(diag[i],left,right);CHKERRQ(ierr); 184 ierr = VecResetArray(right);CHKERRQ(ierr); 185 } 186 } 187 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 188 189 for (i=mbs-1; i>=0; i--) { 190 n = a->i[i+1] - a->i[i]; 191 idx = a->j + a->i[i]; 192 v = a->a + a->i[i]; 193 194 ierr = VecSet(left,0.0);CHKERRQ(ierr); 195 for (j=0; j<n; j++) { 196 if (idx[j] != i) { 197 ierr = VecPlaceArray(right,x + idx[j]*bs);CHKERRQ(ierr); 198 ierr = MatMultAdd(v[j],right,left,left);CHKERRQ(ierr); 199 ierr = VecResetArray(right);CHKERRQ(ierr); 200 } 201 } 202 ierr = VecPlaceArray(right,b + i*bs);CHKERRQ(ierr); 203 ierr = VecAYPX(left,-1.0,right);CHKERRQ(ierr); 204 ierr = VecResetArray(right);CHKERRQ(ierr); 205 206 ierr = VecPlaceArray(right,x + i*bs);CHKERRQ(ierr); 207 ierr = MatSolve(diag[i],left,right);CHKERRQ(ierr); 208 ierr = VecResetArray(right);CHKERRQ(ierr); 209 210 } 211 } 212 } 213 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 214 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 static PetscErrorCode MatSetValues_BlockMat(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 219 { 220 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 221 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 222 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 223 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 224 PetscErrorCode ierr; 225 PetscInt ridx,cidx; 226 PetscBool roworiented=a->roworiented; 227 MatScalar value; 228 Mat *ap,*aa = a->a; 229 230 PetscFunctionBegin; 231 for (k=0; k<m; k++) { /* loop over added rows */ 232 row = im[k]; 233 brow = row/bs; 234 if (row < 0) continue; 235 if (PetscUnlikelyDebug(row >= A->rmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 236 rp = aj + ai[brow]; 237 ap = aa + ai[brow]; 238 rmax = imax[brow]; 239 nrow = ailen[brow]; 240 low = 0; 241 high = nrow; 242 for (l=0; l<n; l++) { /* loop over added columns */ 243 if (in[l] < 0) continue; 244 if (PetscUnlikelyDebug(in[l] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 245 col = in[l]; bcol = col/bs; 246 if (A->symmetric && brow > bcol) continue; 247 ridx = row % bs; cidx = col % bs; 248 if (roworiented) value = v[l + k*n]; 249 else value = v[k + l*m]; 250 251 if (col <= lastcol) low = 0; 252 else high = nrow; 253 lastcol = col; 254 while (high-low > 7) { 255 t = (low+high)/2; 256 if (rp[t] > bcol) high = t; 257 else low = t; 258 } 259 for (i=low; i<high; i++) { 260 if (rp[i] > bcol) break; 261 if (rp[i] == bcol) goto noinsert1; 262 } 263 if (nonew == 1) goto noinsert1; 264 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 265 MatSeqXAIJReallocateAIJ(A,a->mbs,1,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,Mat); 266 N = nrow++ - 1; high++; 267 /* shift up all the later entries in this row */ 268 for (ii=N; ii>=i; ii--) { 269 rp[ii+1] = rp[ii]; 270 ap[ii+1] = ap[ii]; 271 } 272 if (N>=i) ap[i] = NULL; 273 rp[i] = bcol; 274 a->nz++; 275 A->nonzerostate++; 276 noinsert1:; 277 if (!*(ap+i)) { 278 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,NULL,ap+i);CHKERRQ(ierr); 279 } 280 ierr = MatSetValues(ap[i],1,&ridx,1,&cidx,&value,is);CHKERRQ(ierr); 281 low = i; 282 } 283 ailen[brow] = nrow; 284 } 285 PetscFunctionReturn(0); 286 } 287 288 static PetscErrorCode MatLoad_BlockMat(Mat newmat, PetscViewer viewer) 289 { 290 PetscErrorCode ierr; 291 Mat tmpA; 292 PetscInt i,j,m,n,bs = 1,ncols,*lens,currentcol,mbs,**ii,*ilens,nextcol,*llens,cnt = 0; 293 const PetscInt *cols; 294 const PetscScalar *values; 295 PetscBool flg = PETSC_FALSE,notdone; 296 Mat_SeqAIJ *a; 297 Mat_BlockMat *amat; 298 299 PetscFunctionBegin; 300 /* force binary viewer to load .info file if it has not yet done so */ 301 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 302 ierr = MatCreate(PETSC_COMM_SELF,&tmpA);CHKERRQ(ierr); 303 ierr = MatSetType(tmpA,MATSEQAIJ);CHKERRQ(ierr); 304 ierr = MatLoad_SeqAIJ(tmpA,viewer);CHKERRQ(ierr); 305 306 ierr = MatGetLocalSize(tmpA,&m,&n);CHKERRQ(ierr); 307 ierr = PetscOptionsBegin(PETSC_COMM_SELF,NULL,"Options for loading BlockMat matrix 1","Mat");CHKERRQ(ierr); 308 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 309 ierr = PetscOptionsBool("-matload_symmetric","Store the matrix as symmetric","MatLoad",flg,&flg,NULL);CHKERRQ(ierr); 310 ierr = PetscOptionsEnd();CHKERRQ(ierr); 311 312 /* Determine number of nonzero blocks for each block row */ 313 a = (Mat_SeqAIJ*) tmpA->data; 314 mbs = m/bs; 315 ierr = PetscMalloc3(mbs,&lens,bs,&ii,bs,&ilens);CHKERRQ(ierr); 316 ierr = PetscArrayzero(lens,mbs);CHKERRQ(ierr); 317 318 for (i=0; i<mbs; i++) { 319 for (j=0; j<bs; j++) { 320 ii[j] = a->j + a->i[i*bs + j]; 321 ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j]; 322 } 323 324 currentcol = -1; 325 while (PETSC_TRUE) { 326 notdone = PETSC_FALSE; 327 nextcol = 1000000000; 328 for (j=0; j<bs; j++) { 329 while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) { 330 ii[j]++; 331 ilens[j]--; 332 } 333 if (ilens[j] > 0) { 334 notdone = PETSC_TRUE; 335 nextcol = PetscMin(nextcol,ii[j][0]/bs); 336 } 337 } 338 if (!notdone) break; 339 if (!flg || (nextcol >= i)) lens[i]++; 340 currentcol = nextcol; 341 } 342 } 343 344 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 345 ierr = MatSetSizes(newmat,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 346 } 347 ierr = MatBlockMatSetPreallocation(newmat,bs,0,lens);CHKERRQ(ierr); 348 if (flg) { 349 ierr = MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 350 } 351 amat = (Mat_BlockMat*)(newmat)->data; 352 353 /* preallocate the submatrices */ 354 ierr = PetscMalloc1(bs,&llens);CHKERRQ(ierr); 355 for (i=0; i<mbs; i++) { /* loops for block rows */ 356 for (j=0; j<bs; j++) { 357 ii[j] = a->j + a->i[i*bs + j]; 358 ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j]; 359 } 360 361 currentcol = 1000000000; 362 for (j=0; j<bs; j++) { /* loop over rows in block finding first nonzero block */ 363 if (ilens[j] > 0) { 364 currentcol = PetscMin(currentcol,ii[j][0]/bs); 365 } 366 } 367 368 while (PETSC_TRUE) { /* loops over blocks in block row */ 369 notdone = PETSC_FALSE; 370 nextcol = 1000000000; 371 ierr = PetscArrayzero(llens,bs);CHKERRQ(ierr); 372 for (j=0; j<bs; j++) { /* loop over rows in block */ 373 while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) { /* loop over columns in row */ 374 ii[j]++; 375 ilens[j]--; 376 llens[j]++; 377 } 378 if (ilens[j] > 0) { 379 notdone = PETSC_TRUE; 380 nextcol = PetscMin(nextcol,ii[j][0]/bs); 381 } 382 } 383 if (cnt >= amat->maxnz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of blocks found greater than expected %D",cnt); 384 if (!flg || currentcol >= i) { 385 amat->j[cnt] = currentcol; 386 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,llens,amat->a+cnt++);CHKERRQ(ierr); 387 } 388 389 if (!notdone) break; 390 currentcol = nextcol; 391 } 392 amat->ilen[i] = lens[i]; 393 } 394 395 ierr = PetscFree3(lens,ii,ilens);CHKERRQ(ierr); 396 ierr = PetscFree(llens);CHKERRQ(ierr); 397 398 /* copy over the matrix, one row at a time */ 399 for (i=0; i<m; i++) { 400 ierr = MatGetRow(tmpA,i,&ncols,&cols,&values);CHKERRQ(ierr); 401 ierr = MatSetValues(newmat,1,&i,ncols,cols,values,INSERT_VALUES);CHKERRQ(ierr); 402 ierr = MatRestoreRow(tmpA,i,&ncols,&cols,&values);CHKERRQ(ierr); 403 } 404 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 405 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408 409 static PetscErrorCode MatView_BlockMat(Mat A,PetscViewer viewer) 410 { 411 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 412 PetscErrorCode ierr; 413 const char *name; 414 PetscViewerFormat format; 415 416 PetscFunctionBegin; 417 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 418 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 419 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 420 ierr = PetscViewerASCIIPrintf(viewer,"Nonzero block matrices = %D \n",a->nz);CHKERRQ(ierr); 421 if (A->symmetric) { 422 ierr = PetscViewerASCIIPrintf(viewer,"Only upper triangular part of symmetric matrix is stored\n");CHKERRQ(ierr); 423 } 424 } 425 PetscFunctionReturn(0); 426 } 427 428 static PetscErrorCode MatDestroy_BlockMat(Mat mat) 429 { 430 PetscErrorCode ierr; 431 Mat_BlockMat *bmat = (Mat_BlockMat*)mat->data; 432 PetscInt i; 433 434 PetscFunctionBegin; 435 ierr = VecDestroy(&bmat->right);CHKERRQ(ierr); 436 ierr = VecDestroy(&bmat->left);CHKERRQ(ierr); 437 ierr = VecDestroy(&bmat->middle);CHKERRQ(ierr); 438 ierr = VecDestroy(&bmat->workb);CHKERRQ(ierr); 439 if (bmat->diags) { 440 for (i=0; i<mat->rmap->n/mat->rmap->bs; i++) { 441 ierr = MatDestroy(&bmat->diags[i]);CHKERRQ(ierr); 442 } 443 } 444 if (bmat->a) { 445 for (i=0; i<bmat->nz; i++) { 446 ierr = MatDestroy(&bmat->a[i]);CHKERRQ(ierr); 447 } 448 } 449 ierr = MatSeqXAIJFreeAIJ(mat,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);CHKERRQ(ierr); 450 ierr = PetscFree(mat->data);CHKERRQ(ierr); 451 PetscFunctionReturn(0); 452 } 453 454 static PetscErrorCode MatMult_BlockMat(Mat A,Vec x,Vec y) 455 { 456 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 457 PetscErrorCode ierr; 458 PetscScalar *xx,*yy; 459 PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j; 460 Mat *aa; 461 462 PetscFunctionBegin; 463 /* 464 Standard CSR multiply except each entry is a Mat 465 */ 466 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 467 468 ierr = VecSet(y,0.0);CHKERRQ(ierr); 469 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 470 aj = bmat->j; 471 aa = bmat->a; 472 ii = bmat->i; 473 for (i=0; i<m; i++) { 474 jrow = ii[i]; 475 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 476 n = ii[i+1] - jrow; 477 for (j=0; j<n; j++) { 478 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 479 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 480 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 481 jrow++; 482 } 483 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 484 } 485 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 486 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 487 PetscFunctionReturn(0); 488 } 489 490 PetscErrorCode MatMult_BlockMat_Symmetric(Mat A,Vec x,Vec y) 491 { 492 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 493 PetscErrorCode ierr; 494 PetscScalar *xx,*yy; 495 PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j; 496 Mat *aa; 497 498 PetscFunctionBegin; 499 /* 500 Standard CSR multiply except each entry is a Mat 501 */ 502 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 503 504 ierr = VecSet(y,0.0);CHKERRQ(ierr); 505 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 506 aj = bmat->j; 507 aa = bmat->a; 508 ii = bmat->i; 509 for (i=0; i<m; i++) { 510 jrow = ii[i]; 511 n = ii[i+1] - jrow; 512 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 513 ierr = VecPlaceArray(bmat->middle,xx + bs*i);CHKERRQ(ierr); 514 /* if we ALWAYS required a diagonal entry then could remove this if test */ 515 if (aj[jrow] == i) { 516 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 517 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 518 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 519 jrow++; 520 n--; 521 } 522 for (j=0; j<n; j++) { 523 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); /* upper triangular part */ 524 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 525 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 526 527 ierr = VecPlaceArray(bmat->right,yy + bs*aj[jrow]);CHKERRQ(ierr); /* lower triangular part */ 528 ierr = MatMultTransposeAdd(aa[jrow],bmat->middle,bmat->right,bmat->right);CHKERRQ(ierr); 529 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 530 jrow++; 531 } 532 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 533 ierr = VecResetArray(bmat->middle);CHKERRQ(ierr); 534 } 535 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 536 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 537 PetscFunctionReturn(0); 538 } 539 540 static PetscErrorCode MatMultAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 541 { 542 PetscFunctionBegin; 543 PetscFunctionReturn(0); 544 } 545 546 static PetscErrorCode MatMultTranspose_BlockMat(Mat A,Vec x,Vec y) 547 { 548 PetscFunctionBegin; 549 PetscFunctionReturn(0); 550 } 551 552 static PetscErrorCode MatMultTransposeAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 553 { 554 PetscFunctionBegin; 555 PetscFunctionReturn(0); 556 } 557 558 /* 559 Adds diagonal pointers to sparse matrix structure. 560 */ 561 static PetscErrorCode MatMarkDiagonal_BlockMat(Mat A) 562 { 563 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 564 PetscErrorCode ierr; 565 PetscInt i,j,mbs = A->rmap->n/A->rmap->bs; 566 567 PetscFunctionBegin; 568 if (!a->diag) { 569 ierr = PetscMalloc1(mbs,&a->diag);CHKERRQ(ierr); 570 } 571 for (i=0; i<mbs; i++) { 572 a->diag[i] = a->i[i+1]; 573 for (j=a->i[i]; j<a->i[i+1]; j++) { 574 if (a->j[j] == i) { 575 a->diag[i] = j; 576 break; 577 } 578 } 579 } 580 PetscFunctionReturn(0); 581 } 582 583 static PetscErrorCode MatCreateSubMatrix_BlockMat(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 584 { 585 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 586 Mat_SeqAIJ *c; 587 PetscErrorCode ierr; 588 PetscInt i,k,first,step,lensi,nrows,ncols; 589 PetscInt *j_new,*i_new,*aj = a->j,*ailen = a->ilen; 590 PetscScalar *a_new; 591 Mat C,*aa = a->a; 592 PetscBool stride,equal; 593 594 PetscFunctionBegin; 595 ierr = ISEqual(isrow,iscol,&equal);CHKERRQ(ierr); 596 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for identical column and row indices"); 597 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 598 if (!stride) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for stride indices"); 599 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 600 if (step != A->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only select one entry from each block"); 601 602 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 603 ncols = nrows; 604 605 /* create submatrix */ 606 if (scall == MAT_REUSE_MATRIX) { 607 PetscInt n_cols,n_rows; 608 C = *B; 609 ierr = MatGetSize(C,&n_rows,&n_cols);CHKERRQ(ierr); 610 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 611 ierr = MatZeroEntries(C);CHKERRQ(ierr); 612 } else { 613 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 614 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 615 if (A->symmetric) { 616 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 617 } else { 618 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 619 } 620 ierr = MatSeqAIJSetPreallocation(C,0,ailen);CHKERRQ(ierr); 621 ierr = MatSeqSBAIJSetPreallocation(C,1,0,ailen);CHKERRQ(ierr); 622 } 623 c = (Mat_SeqAIJ*)C->data; 624 625 /* loop over rows inserting into submatrix */ 626 a_new = c->a; 627 j_new = c->j; 628 i_new = c->i; 629 630 for (i=0; i<nrows; i++) { 631 lensi = ailen[i]; 632 for (k=0; k<lensi; k++) { 633 *j_new++ = *aj++; 634 ierr = MatGetValue(*aa++,first,first,a_new++);CHKERRQ(ierr); 635 } 636 i_new[i+1] = i_new[i] + lensi; 637 c->ilen[i] = lensi; 638 } 639 640 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 641 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 642 *B = C; 643 PetscFunctionReturn(0); 644 } 645 646 static PetscErrorCode MatAssemblyEnd_BlockMat(Mat A,MatAssemblyType mode) 647 { 648 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 649 PetscErrorCode ierr; 650 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 651 PetscInt m = a->mbs,*ip,N,*ailen = a->ilen,rmax = 0; 652 Mat *aa = a->a,*ap; 653 654 PetscFunctionBegin; 655 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 656 657 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 658 for (i=1; i<m; i++) { 659 /* move each row back by the amount of empty slots (fshift) before it*/ 660 fshift += imax[i-1] - ailen[i-1]; 661 rmax = PetscMax(rmax,ailen[i]); 662 if (fshift) { 663 ip = aj + ai[i]; 664 ap = aa + ai[i]; 665 N = ailen[i]; 666 for (j=0; j<N; j++) { 667 ip[j-fshift] = ip[j]; 668 ap[j-fshift] = ap[j]; 669 } 670 } 671 ai[i] = ai[i-1] + ailen[i-1]; 672 } 673 if (m) { 674 fshift += imax[m-1] - ailen[m-1]; 675 ai[m] = ai[m-1] + ailen[m-1]; 676 } 677 /* reset ilen and imax for each row */ 678 for (i=0; i<m; i++) { 679 ailen[i] = imax[i] = ai[i+1] - ai[i]; 680 } 681 a->nz = ai[m]; 682 for (i=0; i<a->nz; i++) { 683 if (PetscUnlikelyDebug(!aa[i])) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Null matrix at location %D column %D nz %D",i,aj[i],a->nz); 684 ierr = MatAssemblyBegin(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 685 ierr = MatAssemblyEnd(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 686 } 687 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n/A->cmap->bs,fshift,a->nz);CHKERRQ(ierr); 688 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 689 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 690 691 A->info.mallocs += a->reallocs; 692 a->reallocs = 0; 693 A->info.nz_unneeded = (double)fshift; 694 a->rmax = rmax; 695 ierr = MatMarkDiagonal_BlockMat(A);CHKERRQ(ierr); 696 PetscFunctionReturn(0); 697 } 698 699 static PetscErrorCode MatSetOption_BlockMat(Mat A,MatOption opt,PetscBool flg) 700 { 701 PetscErrorCode ierr; 702 PetscFunctionBegin; 703 if (opt == MAT_SYMMETRIC && flg) { 704 A->ops->sor = MatSOR_BlockMat_Symmetric; 705 A->ops->mult = MatMult_BlockMat_Symmetric; 706 } else { 707 ierr = PetscInfo1(A,"Unused matrix option %s\n",MatOptions[opt]);CHKERRQ(ierr); 708 } 709 PetscFunctionReturn(0); 710 } 711 712 713 static struct _MatOps MatOps_Values = {MatSetValues_BlockMat, 714 NULL, 715 NULL, 716 MatMult_BlockMat, 717 /* 4*/ MatMultAdd_BlockMat, 718 MatMultTranspose_BlockMat, 719 MatMultTransposeAdd_BlockMat, 720 NULL, 721 NULL, 722 NULL, 723 /* 10*/ NULL, 724 NULL, 725 NULL, 726 MatSOR_BlockMat, 727 NULL, 728 /* 15*/ NULL, 729 NULL, 730 NULL, 731 NULL, 732 NULL, 733 /* 20*/ NULL, 734 MatAssemblyEnd_BlockMat, 735 MatSetOption_BlockMat, 736 NULL, 737 /* 24*/ NULL, 738 NULL, 739 NULL, 740 NULL, 741 NULL, 742 /* 29*/ NULL, 743 NULL, 744 NULL, 745 NULL, 746 NULL, 747 /* 34*/ NULL, 748 NULL, 749 NULL, 750 NULL, 751 NULL, 752 /* 39*/ NULL, 753 NULL, 754 NULL, 755 NULL, 756 NULL, 757 /* 44*/ NULL, 758 NULL, 759 MatShift_Basic, 760 NULL, 761 NULL, 762 /* 49*/ NULL, 763 NULL, 764 NULL, 765 NULL, 766 NULL, 767 /* 54*/ NULL, 768 NULL, 769 NULL, 770 NULL, 771 NULL, 772 /* 59*/ MatCreateSubMatrix_BlockMat, 773 MatDestroy_BlockMat, 774 MatView_BlockMat, 775 NULL, 776 NULL, 777 /* 64*/ NULL, 778 NULL, 779 NULL, 780 NULL, 781 NULL, 782 /* 69*/ NULL, 783 NULL, 784 NULL, 785 NULL, 786 NULL, 787 /* 74*/ NULL, 788 NULL, 789 NULL, 790 NULL, 791 NULL, 792 /* 79*/ NULL, 793 NULL, 794 NULL, 795 NULL, 796 MatLoad_BlockMat, 797 /* 84*/ NULL, 798 NULL, 799 NULL, 800 NULL, 801 NULL, 802 /* 89*/ NULL, 803 NULL, 804 NULL, 805 NULL, 806 NULL, 807 /* 94*/ NULL, 808 NULL, 809 NULL, 810 NULL, 811 NULL, 812 /* 99*/ NULL, 813 NULL, 814 NULL, 815 NULL, 816 NULL, 817 /*104*/ NULL, 818 NULL, 819 NULL, 820 NULL, 821 NULL, 822 /*109*/ NULL, 823 NULL, 824 NULL, 825 NULL, 826 NULL, 827 /*114*/ NULL, 828 NULL, 829 NULL, 830 NULL, 831 NULL, 832 /*119*/ NULL, 833 NULL, 834 NULL, 835 NULL, 836 NULL, 837 /*124*/ NULL, 838 NULL, 839 NULL, 840 NULL, 841 NULL, 842 /*129*/ NULL, 843 NULL, 844 NULL, 845 NULL, 846 NULL, 847 /*134*/ NULL, 848 NULL, 849 NULL, 850 NULL, 851 NULL, 852 /*139*/ NULL, 853 NULL, 854 NULL 855 }; 856 857 /*@C 858 MatBlockMatSetPreallocation - For good matrix assembly performance 859 the user should preallocate the matrix storage by setting the parameter nz 860 (or the array nnz). By setting these parameters accurately, performance 861 during matrix assembly can be increased by more than a factor of 50. 862 863 Collective 864 865 Input Parameters: 866 + B - The matrix 867 . bs - size of each block in matrix 868 . nz - number of nonzeros per block row (same for all rows) 869 - nnz - array containing the number of nonzeros in the various block rows 870 (possibly different for each row) or NULL 871 872 Notes: 873 If nnz is given then nz is ignored 874 875 Specify the preallocated storage with either nz or nnz (not both). 876 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 877 allocation. For large problems you MUST preallocate memory or you 878 will get TERRIBLE performance, see the users' manual chapter on matrices. 879 880 Level: intermediate 881 882 .seealso: MatCreate(), MatCreateBlockMat(), MatSetValues() 883 884 @*/ 885 PetscErrorCode MatBlockMatSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 886 { 887 PetscErrorCode ierr; 888 889 PetscFunctionBegin; 890 ierr = PetscTryMethod(B,"MatBlockMatSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 891 PetscFunctionReturn(0); 892 } 893 894 static PetscErrorCode MatBlockMatSetPreallocation_BlockMat(Mat A,PetscInt bs,PetscInt nz,PetscInt *nnz) 895 { 896 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 897 PetscErrorCode ierr; 898 PetscInt i; 899 900 PetscFunctionBegin; 901 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 902 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 903 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 904 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 905 ierr = PetscLayoutGetBlockSize(A->rmap,&bs);CHKERRQ(ierr); 906 907 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 908 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 909 if (nnz) { 910 for (i=0; i<A->rmap->n/bs; i++) { 911 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 912 if (nnz[i] > A->cmap->n/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],A->cmap->n/bs); 913 } 914 } 915 bmat->mbs = A->rmap->n/bs; 916 917 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->right);CHKERRQ(ierr); 918 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->middle);CHKERRQ(ierr); 919 ierr = VecCreateSeq(PETSC_COMM_SELF,bs,&bmat->left);CHKERRQ(ierr); 920 921 if (!bmat->imax) { 922 ierr = PetscMalloc2(A->rmap->n,&bmat->imax,A->rmap->n,&bmat->ilen);CHKERRQ(ierr); 923 ierr = PetscLogObjectMemory((PetscObject)A,2*A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 924 } 925 if (nnz) { 926 nz = 0; 927 for (i=0; i<A->rmap->n/A->rmap->bs; i++) { 928 bmat->imax[i] = nnz[i]; 929 nz += nnz[i]; 930 } 931 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently requires block row by row preallocation"); 932 933 /* bmat->ilen will count nonzeros in each row so far. */ 934 for (i=0; i<bmat->mbs; i++) bmat->ilen[i] = 0; 935 936 /* allocate the matrix space */ 937 ierr = MatSeqXAIJFreeAIJ(A,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);CHKERRQ(ierr); 938 ierr = PetscMalloc3(nz,&bmat->a,nz,&bmat->j,A->rmap->n+1,&bmat->i);CHKERRQ(ierr); 939 ierr = PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 940 bmat->i[0] = 0; 941 for (i=1; i<bmat->mbs+1; i++) { 942 bmat->i[i] = bmat->i[i-1] + bmat->imax[i-1]; 943 } 944 bmat->singlemalloc = PETSC_TRUE; 945 bmat->free_a = PETSC_TRUE; 946 bmat->free_ij = PETSC_TRUE; 947 948 bmat->nz = 0; 949 bmat->maxnz = nz; 950 A->info.nz_unneeded = (double)bmat->maxnz; 951 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 952 PetscFunctionReturn(0); 953 } 954 955 /*MC 956 MATBLOCKMAT - A matrix that is defined by a set of Mat's that represents a sparse block matrix 957 consisting of (usually) sparse blocks. 958 959 Level: advanced 960 961 .seealso: MatCreateBlockMat() 962 963 M*/ 964 965 PETSC_EXTERN PetscErrorCode MatCreate_BlockMat(Mat A) 966 { 967 Mat_BlockMat *b; 968 PetscErrorCode ierr; 969 970 PetscFunctionBegin; 971 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 972 A->data = (void*)b; 973 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 974 975 A->assembled = PETSC_TRUE; 976 A->preallocated = PETSC_FALSE; 977 ierr = PetscObjectChangeTypeName((PetscObject)A,MATBLOCKMAT);CHKERRQ(ierr); 978 979 ierr = PetscObjectComposeFunction((PetscObject)A,"MatBlockMatSetPreallocation_C",MatBlockMatSetPreallocation_BlockMat);CHKERRQ(ierr); 980 PetscFunctionReturn(0); 981 } 982 983 /*@C 984 MatCreateBlockMat - Creates a new matrix in which each block contains a uniform-size sequential Mat object 985 986 Collective 987 988 Input Parameters: 989 + comm - MPI communicator 990 . m - number of rows 991 . n - number of columns 992 . bs - size of each submatrix 993 . nz - expected maximum number of nonzero blocks in row (use PETSC_DEFAULT if not known) 994 - nnz - expected number of nonzers per block row if known (use NULL otherwise) 995 996 997 Output Parameter: 998 . A - the matrix 999 1000 Level: intermediate 1001 1002 Notes: 1003 Matrices of this type are nominally-sparse matrices in which each "entry" is a Mat object. Each Mat must 1004 have the same size and be sequential. The local and global sizes must be compatible with this decomposition. 1005 1006 For matrices containing parallel submatrices and variable block sizes, see MATNEST. 1007 1008 .seealso: MATBLOCKMAT, MatCreateNest() 1009 @*/ 1010 PetscErrorCode MatCreateBlockMat(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt bs,PetscInt nz,PetscInt *nnz, Mat *A) 1011 { 1012 PetscErrorCode ierr; 1013 1014 PetscFunctionBegin; 1015 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1016 ierr = MatSetSizes(*A,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1017 ierr = MatSetType(*A,MATBLOCKMAT);CHKERRQ(ierr); 1018 ierr = MatBlockMatSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1019 PetscFunctionReturn(0); 1020 } 1021