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 CHKMEMQ; 411 412 ierr = PetscFree3(lens,ii,ilens);CHKERRQ(ierr); 413 ierr = PetscFree(llens);CHKERRQ(ierr); 414 415 /* copy over the matrix, one row at a time */ 416 for (i=0; i<m; i++) { 417 ierr = MatGetRow(tmpA,i,&ncols,&cols,&values);CHKERRQ(ierr); 418 ierr = MatSetValues(newmat,1,&i,ncols,cols,values,INSERT_VALUES);CHKERRQ(ierr); 419 ierr = MatRestoreRow(tmpA,i,&ncols,&cols,&values);CHKERRQ(ierr); 420 } 421 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 422 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 423 PetscFunctionReturn(0); 424 } 425 426 #undef __FUNCT__ 427 #define __FUNCT__ "MatView_BlockMat" 428 PetscErrorCode MatView_BlockMat(Mat A,PetscViewer viewer) 429 { 430 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 431 PetscErrorCode ierr; 432 const char *name; 433 PetscViewerFormat format; 434 435 PetscFunctionBegin; 436 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 437 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 438 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 439 ierr = PetscViewerASCIIPrintf(viewer,"Nonzero block matrices = %D \n",a->nz);CHKERRQ(ierr); 440 if (A->symmetric) { 441 ierr = PetscViewerASCIIPrintf(viewer,"Only upper triangular part of symmetric matrix is stored\n");CHKERRQ(ierr); 442 } 443 } 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNCT__ 448 #define __FUNCT__ "MatDestroy_BlockMat" 449 PetscErrorCode MatDestroy_BlockMat(Mat mat) 450 { 451 PetscErrorCode ierr; 452 Mat_BlockMat *bmat = (Mat_BlockMat*)mat->data; 453 PetscInt i; 454 455 PetscFunctionBegin; 456 ierr = VecDestroy(&bmat->right);CHKERRQ(ierr); 457 ierr = VecDestroy(&bmat->left);CHKERRQ(ierr); 458 ierr = VecDestroy(&bmat->middle);CHKERRQ(ierr); 459 ierr = VecDestroy(&bmat->workb);CHKERRQ(ierr); 460 if (bmat->diags) { 461 for (i=0; i<mat->rmap->n/mat->rmap->bs; i++) { 462 ierr = MatDestroy(&bmat->diags[i]);CHKERRQ(ierr); 463 } 464 } 465 if (bmat->a) { 466 for (i=0; i<bmat->nz; i++) { 467 ierr = MatDestroy(&bmat->a[i]);CHKERRQ(ierr); 468 } 469 } 470 ierr = MatSeqXAIJFreeAIJ(mat,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);CHKERRQ(ierr); 471 ierr = PetscFree(mat->data);CHKERRQ(ierr); 472 PetscFunctionReturn(0); 473 } 474 475 #undef __FUNCT__ 476 #define __FUNCT__ "MatMult_BlockMat" 477 PetscErrorCode MatMult_BlockMat(Mat A,Vec x,Vec y) 478 { 479 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 480 PetscErrorCode ierr; 481 PetscScalar *xx,*yy; 482 PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j; 483 Mat *aa; 484 485 PetscFunctionBegin; 486 CHKMEMQ; 487 /* 488 Standard CSR multiply except each entry is a Mat 489 */ 490 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 491 492 ierr = VecSet(y,0.0);CHKERRQ(ierr); 493 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 494 aj = bmat->j; 495 aa = bmat->a; 496 ii = bmat->i; 497 for (i=0; i<m; i++) { 498 jrow = ii[i]; 499 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 500 n = ii[i+1] - jrow; 501 for (j=0; j<n; j++) { 502 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 503 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 504 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 505 jrow++; 506 } 507 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 508 } 509 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 510 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 511 CHKMEMQ; 512 PetscFunctionReturn(0); 513 } 514 515 #undef __FUNCT__ 516 #define __FUNCT__ "MatMult_BlockMat_Symmetric" 517 PetscErrorCode MatMult_BlockMat_Symmetric(Mat A,Vec x,Vec y) 518 { 519 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 520 PetscErrorCode ierr; 521 PetscScalar *xx,*yy; 522 PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j; 523 Mat *aa; 524 525 PetscFunctionBegin; 526 CHKMEMQ; 527 /* 528 Standard CSR multiply except each entry is a Mat 529 */ 530 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 531 532 ierr = VecSet(y,0.0);CHKERRQ(ierr); 533 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 534 aj = bmat->j; 535 aa = bmat->a; 536 ii = bmat->i; 537 for (i=0; i<m; i++) { 538 jrow = ii[i]; 539 n = ii[i+1] - jrow; 540 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 541 ierr = VecPlaceArray(bmat->middle,xx + bs*i);CHKERRQ(ierr); 542 /* if we ALWAYS required a diagonal entry then could remove this if test */ 543 if (aj[jrow] == i) { 544 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 545 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 546 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 547 jrow++; 548 n--; 549 } 550 for (j=0; j<n; j++) { 551 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); /* upper triangular part */ 552 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 553 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 554 555 ierr = VecPlaceArray(bmat->right,yy + bs*aj[jrow]);CHKERRQ(ierr); /* lower triangular part */ 556 ierr = MatMultTransposeAdd(aa[jrow],bmat->middle,bmat->right,bmat->right);CHKERRQ(ierr); 557 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 558 jrow++; 559 } 560 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 561 ierr = VecResetArray(bmat->middle);CHKERRQ(ierr); 562 } 563 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 564 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 565 CHKMEMQ; 566 PetscFunctionReturn(0); 567 } 568 569 #undef __FUNCT__ 570 #define __FUNCT__ "MatMultAdd_BlockMat" 571 PetscErrorCode MatMultAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 572 { 573 PetscFunctionBegin; 574 PetscFunctionReturn(0); 575 } 576 577 #undef __FUNCT__ 578 #define __FUNCT__ "MatMultTranspose_BlockMat" 579 PetscErrorCode MatMultTranspose_BlockMat(Mat A,Vec x,Vec y) 580 { 581 PetscFunctionBegin; 582 PetscFunctionReturn(0); 583 } 584 585 #undef __FUNCT__ 586 #define __FUNCT__ "MatMultTransposeAdd_BlockMat" 587 PetscErrorCode MatMultTransposeAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 588 { 589 PetscFunctionBegin; 590 PetscFunctionReturn(0); 591 } 592 593 /* 594 Adds diagonal pointers to sparse matrix structure. 595 */ 596 #undef __FUNCT__ 597 #define __FUNCT__ "MatMarkDiagonal_BlockMat" 598 PetscErrorCode MatMarkDiagonal_BlockMat(Mat A) 599 { 600 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 601 PetscErrorCode ierr; 602 PetscInt i,j,mbs = A->rmap->n/A->rmap->bs; 603 604 PetscFunctionBegin; 605 if (!a->diag) { 606 ierr = PetscMalloc(mbs*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 607 } 608 for (i=0; i<mbs; i++) { 609 a->diag[i] = a->i[i+1]; 610 for (j=a->i[i]; j<a->i[i+1]; j++) { 611 if (a->j[j] == i) { 612 a->diag[i] = j; 613 break; 614 } 615 } 616 } 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "MatGetSubMatrix_BlockMat" 622 PetscErrorCode MatGetSubMatrix_BlockMat(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 623 { 624 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 625 Mat_SeqAIJ *c; 626 PetscErrorCode ierr; 627 PetscInt i,k,first,step,lensi,nrows,ncols; 628 PetscInt *j_new,*i_new,*aj = a->j,*ailen = a->ilen; 629 PetscScalar *a_new; 630 Mat C,*aa = a->a; 631 PetscBool stride,equal; 632 633 PetscFunctionBegin; 634 ierr = ISEqual(isrow,iscol,&equal);CHKERRQ(ierr); 635 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices"); 636 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 637 if (!stride) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for stride indices"); 638 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 639 if (step != A->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only select one entry from each block"); 640 641 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 642 ncols = nrows; 643 644 /* create submatrix */ 645 if (scall == MAT_REUSE_MATRIX) { 646 PetscInt n_cols,n_rows; 647 C = *B; 648 ierr = MatGetSize(C,&n_rows,&n_cols);CHKERRQ(ierr); 649 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 650 ierr = MatZeroEntries(C);CHKERRQ(ierr); 651 } else { 652 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 653 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 654 if (A->symmetric) { 655 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 656 } else { 657 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 658 } 659 ierr = MatSeqAIJSetPreallocation(C,0,ailen);CHKERRQ(ierr); 660 ierr = MatSeqSBAIJSetPreallocation(C,1,0,ailen);CHKERRQ(ierr); 661 } 662 c = (Mat_SeqAIJ*)C->data; 663 664 /* loop over rows inserting into submatrix */ 665 a_new = c->a; 666 j_new = c->j; 667 i_new = c->i; 668 669 for (i=0; i<nrows; i++) { 670 lensi = ailen[i]; 671 for (k=0; k<lensi; k++) { 672 *j_new++ = *aj++; 673 ierr = MatGetValue(*aa++,first,first,a_new++);CHKERRQ(ierr); 674 } 675 i_new[i+1] = i_new[i] + lensi; 676 c->ilen[i] = lensi; 677 } 678 679 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 680 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 681 *B = C; 682 PetscFunctionReturn(0); 683 } 684 685 #undef __FUNCT__ 686 #define __FUNCT__ "MatAssemblyEnd_BlockMat" 687 PetscErrorCode MatAssemblyEnd_BlockMat(Mat A,MatAssemblyType mode) 688 { 689 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 690 PetscErrorCode ierr; 691 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 692 PetscInt m = a->mbs,*ip,N,*ailen = a->ilen,rmax = 0; 693 Mat *aa = a->a,*ap; 694 695 PetscFunctionBegin; 696 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 697 698 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 699 for (i=1; i<m; i++) { 700 /* move each row back by the amount of empty slots (fshift) before it*/ 701 fshift += imax[i-1] - ailen[i-1]; 702 rmax = PetscMax(rmax,ailen[i]); 703 if (fshift) { 704 ip = aj + ai[i]; 705 ap = aa + ai[i]; 706 N = ailen[i]; 707 for (j=0; j<N; j++) { 708 ip[j-fshift] = ip[j]; 709 ap[j-fshift] = ap[j]; 710 } 711 } 712 ai[i] = ai[i-1] + ailen[i-1]; 713 } 714 if (m) { 715 fshift += imax[m-1] - ailen[m-1]; 716 ai[m] = ai[m-1] + ailen[m-1]; 717 } 718 /* reset ilen and imax for each row */ 719 for (i=0; i<m; i++) { 720 ailen[i] = imax[i] = ai[i+1] - ai[i]; 721 } 722 a->nz = ai[m]; 723 for (i=0; i<a->nz; i++) { 724 #if defined(PETSC_USE_DEBUG) 725 if (!aa[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Null matrix at location %D column %D nz %D",i,aj[i],a->nz); 726 #endif 727 ierr = MatAssemblyBegin(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 728 ierr = MatAssemblyEnd(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 729 } 730 CHKMEMQ; 731 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); 732 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 733 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 734 735 A->info.mallocs += a->reallocs; 736 a->reallocs = 0; 737 A->info.nz_unneeded = (double)fshift; 738 a->rmax = rmax; 739 740 A->same_nonzero = PETSC_TRUE; 741 ierr = MatMarkDiagonal_BlockMat(A);CHKERRQ(ierr); 742 PetscFunctionReturn(0); 743 } 744 745 #undef __FUNCT__ 746 #define __FUNCT__ "MatSetOption_BlockMat" 747 PetscErrorCode MatSetOption_BlockMat(Mat A,MatOption opt,PetscBool flg) 748 { 749 PetscFunctionBegin; 750 if (opt == MAT_SYMMETRIC && flg) { 751 A->ops->sor = MatSOR_BlockMat_Symmetric; 752 A->ops->mult = MatMult_BlockMat_Symmetric; 753 } else { 754 PetscInfo1(A,"Unused matrix option %s\n",MatOptions[opt]); 755 } 756 PetscFunctionReturn(0); 757 } 758 759 760 static struct _MatOps MatOps_Values = { 761 MatSetValues_BlockMat, 762 0, 763 0, 764 MatMult_BlockMat, 765 /* 4*/ MatMultAdd_BlockMat, 766 MatMultTranspose_BlockMat, 767 MatMultTransposeAdd_BlockMat, 768 0, 769 0, 770 0, 771 /* 10*/ 0, 772 0, 773 0, 774 MatSOR_BlockMat, 775 0, 776 /* 15*/ 0, 777 0, 778 0, 779 0, 780 0, 781 /* 20*/ 0, 782 MatAssemblyEnd_BlockMat, 783 MatSetOption_BlockMat, 784 0, 785 /* 24*/ 0, 786 0, 787 0, 788 0, 789 0, 790 /* 29*/ 0, 791 0, 792 0, 793 0, 794 0, 795 /* 34*/ 0, 796 0, 797 0, 798 0, 799 0, 800 /* 39*/ 0, 801 0, 802 0, 803 0, 804 0, 805 /* 44*/ 0, 806 0, 807 0, 808 0, 809 0, 810 /* 49*/ 0, 811 0, 812 0, 813 0, 814 0, 815 /* 54*/ 0, 816 0, 817 0, 818 0, 819 0, 820 /* 59*/ MatGetSubMatrix_BlockMat, 821 MatDestroy_BlockMat, 822 MatView_BlockMat, 823 0, 824 0, 825 /* 64*/ 0, 826 0, 827 0, 828 0, 829 0, 830 /* 69*/ 0, 831 0, 832 0, 833 0, 834 0, 835 /* 74*/ 0, 836 0, 837 0, 838 0, 839 0, 840 /* 79*/ 0, 841 0, 842 0, 843 0, 844 MatLoad_BlockMat, 845 /* 84*/ 0, 846 0, 847 0, 848 0, 849 0, 850 /* 89*/ 0, 851 0, 852 0, 853 0, 854 0, 855 /* 94*/ 0, 856 0, 857 0, 858 0, 859 0, 860 /* 99*/ 0, 861 0, 862 0, 863 0, 864 0, 865 /*104*/ 0, 866 0, 867 0, 868 0, 869 0, 870 /*109*/ 0, 871 0, 872 0, 873 0, 874 0, 875 /*114*/ 0, 876 0, 877 0, 878 0, 879 0, 880 /*119*/ 0, 881 0, 882 0, 883 0 884 }; 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "MatBlockMatSetPreallocation" 888 /*@C 889 MatBlockMatSetPreallocation - For good matrix assembly performance 890 the user should preallocate the matrix storage by setting the parameter nz 891 (or the array nnz). By setting these parameters accurately, performance 892 during matrix assembly can be increased by more than a factor of 50. 893 894 Collective on MPI_Comm 895 896 Input Parameters: 897 + B - The matrix 898 . bs - size of each block in matrix 899 . nz - number of nonzeros per block row (same for all rows) 900 - nnz - array containing the number of nonzeros in the various block rows 901 (possibly different for each row) or NULL 902 903 Notes: 904 If nnz is given then nz is ignored 905 906 Specify the preallocated storage with either nz or nnz (not both). 907 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 908 allocation. For large problems you MUST preallocate memory or you 909 will get TERRIBLE performance, see the users' manual chapter on matrices. 910 911 Level: intermediate 912 913 .seealso: MatCreate(), MatCreateBlockMat(), MatSetValues() 914 915 @*/ 916 PetscErrorCode MatBlockMatSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 917 { 918 PetscErrorCode ierr; 919 920 PetscFunctionBegin; 921 ierr = PetscTryMethod(B,"MatBlockMatSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 922 PetscFunctionReturn(0); 923 } 924 925 EXTERN_C_BEGIN 926 #undef __FUNCT__ 927 #define __FUNCT__ "MatBlockMatSetPreallocation_BlockMat" 928 PetscErrorCode MatBlockMatSetPreallocation_BlockMat(Mat A,PetscInt bs,PetscInt nz,PetscInt *nnz) 929 { 930 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 931 PetscErrorCode ierr; 932 PetscInt i; 933 934 PetscFunctionBegin; 935 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 936 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 937 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 938 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 939 ierr = PetscLayoutGetBlockSize(A->rmap,&bs);CHKERRQ(ierr); 940 941 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 942 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 943 if (nnz) { 944 for (i=0; i<A->rmap->n/bs; i++) { 945 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]); 946 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); 947 } 948 } 949 bmat->mbs = A->rmap->n/bs; 950 951 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->right);CHKERRQ(ierr); 952 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->middle);CHKERRQ(ierr); 953 ierr = VecCreateSeq(PETSC_COMM_SELF,bs,&bmat->left);CHKERRQ(ierr); 954 955 if (!bmat->imax) { 956 ierr = PetscMalloc2(A->rmap->n,PetscInt,&bmat->imax,A->rmap->n,PetscInt,&bmat->ilen);CHKERRQ(ierr); 957 ierr = PetscLogObjectMemory(A,2*A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 958 } 959 if (nnz) { 960 nz = 0; 961 for (i=0; i<A->rmap->n/A->rmap->bs; i++) { 962 bmat->imax[i] = nnz[i]; 963 nz += nnz[i]; 964 } 965 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently requires block row by row preallocation"); 966 967 /* bmat->ilen will count nonzeros in each row so far. */ 968 for (i=0; i<bmat->mbs; i++) bmat->ilen[i] = 0; 969 970 /* allocate the matrix space */ 971 ierr = MatSeqXAIJFreeAIJ(A,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);CHKERRQ(ierr); 972 ierr = PetscMalloc3(nz,Mat,&bmat->a,nz,PetscInt,&bmat->j,A->rmap->n+1,PetscInt,&bmat->i);CHKERRQ(ierr); 973 ierr = PetscLogObjectMemory(A,(A->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 974 bmat->i[0] = 0; 975 for (i=1; i<bmat->mbs+1; i++) { 976 bmat->i[i] = bmat->i[i-1] + bmat->imax[i-1]; 977 } 978 bmat->singlemalloc = PETSC_TRUE; 979 bmat->free_a = PETSC_TRUE; 980 bmat->free_ij = PETSC_TRUE; 981 982 bmat->nz = 0; 983 bmat->maxnz = nz; 984 A->info.nz_unneeded = (double)bmat->maxnz; 985 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 986 PetscFunctionReturn(0); 987 } 988 EXTERN_C_END 989 990 /*MC 991 MATBLOCKMAT - A matrix that is defined by a set of Mat's that represents a sparse block matrix 992 consisting of (usually) sparse blocks. 993 994 Level: advanced 995 996 .seealso: MatCreateBlockMat() 997 998 M*/ 999 1000 EXTERN_C_BEGIN 1001 #undef __FUNCT__ 1002 #define __FUNCT__ "MatCreate_BlockMat" 1003 PetscErrorCode MatCreate_BlockMat(Mat A) 1004 { 1005 Mat_BlockMat *b; 1006 PetscErrorCode ierr; 1007 1008 PetscFunctionBegin; 1009 ierr = PetscNewLog(A,Mat_BlockMat,&b);CHKERRQ(ierr); 1010 A->data = (void*)b; 1011 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1012 1013 A->assembled = PETSC_TRUE; 1014 A->preallocated = PETSC_FALSE; 1015 ierr = PetscObjectChangeTypeName((PetscObject)A,MATBLOCKMAT);CHKERRQ(ierr); 1016 1017 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatBlockMatSetPreallocation_C", 1018 "MatBlockMatSetPreallocation_BlockMat", 1019 MatBlockMatSetPreallocation_BlockMat);CHKERRQ(ierr); 1020 PetscFunctionReturn(0); 1021 } 1022 EXTERN_C_END 1023 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "MatCreateBlockMat" 1026 /*@C 1027 MatCreateBlockMat - Creates a new matrix based sparse Mat storage 1028 1029 Collective on MPI_Comm 1030 1031 Input Parameters: 1032 + comm - MPI communicator 1033 . m - number of rows 1034 . n - number of columns 1035 . bs - size of each submatrix 1036 . nz - expected maximum number of nonzero blocks in row (use PETSC_DEFAULT if not known) 1037 - nnz - expected number of nonzers per block row if known (use NULL otherwise) 1038 1039 1040 Output Parameter: 1041 . A - the matrix 1042 1043 Level: intermediate 1044 1045 PETSc requires that matrices and vectors being used for certain 1046 operations are partitioned accordingly. For example, when 1047 creating a bmat matrix, A, that supports parallel matrix-vector 1048 products using MatMult(A,x,y) the user should set the number 1049 of local matrix rows to be the number of local elements of the 1050 corresponding result vector, y. Note that this is information is 1051 required for use of the matrix interface routines, even though 1052 the bmat matrix may not actually be physically partitioned. 1053 For example, 1054 1055 .keywords: matrix, bmat, create 1056 1057 .seealso: MATBLOCKMAT 1058 @*/ 1059 PetscErrorCode MatCreateBlockMat(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt bs,PetscInt nz,PetscInt *nnz, Mat *A) 1060 { 1061 PetscErrorCode ierr; 1062 1063 PetscFunctionBegin; 1064 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1065 ierr = MatSetSizes(*A,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1066 ierr = MatSetType(*A,MATBLOCKMAT);CHKERRQ(ierr); 1067 ierr = MatBlockMatSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 1072 1073