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