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(PetscObjectComm((PetscObject)A),&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 = {MatSetValues_BlockMat, 761 0, 762 0, 763 MatMult_BlockMat, 764 /* 4*/ MatMultAdd_BlockMat, 765 MatMultTranspose_BlockMat, 766 MatMultTransposeAdd_BlockMat, 767 0, 768 0, 769 0, 770 /* 10*/ 0, 771 0, 772 0, 773 MatSOR_BlockMat, 774 0, 775 /* 15*/ 0, 776 0, 777 0, 778 0, 779 0, 780 /* 20*/ 0, 781 MatAssemblyEnd_BlockMat, 782 MatSetOption_BlockMat, 783 0, 784 /* 24*/ 0, 785 0, 786 0, 787 0, 788 0, 789 /* 29*/ 0, 790 0, 791 0, 792 0, 793 0, 794 /* 34*/ 0, 795 0, 796 0, 797 0, 798 0, 799 /* 39*/ 0, 800 0, 801 0, 802 0, 803 0, 804 /* 44*/ 0, 805 0, 806 0, 807 0, 808 0, 809 /* 49*/ 0, 810 0, 811 0, 812 0, 813 0, 814 /* 54*/ 0, 815 0, 816 0, 817 0, 818 0, 819 /* 59*/ MatGetSubMatrix_BlockMat, 820 MatDestroy_BlockMat, 821 MatView_BlockMat, 822 0, 823 0, 824 /* 64*/ 0, 825 0, 826 0, 827 0, 828 0, 829 /* 69*/ 0, 830 0, 831 0, 832 0, 833 0, 834 /* 74*/ 0, 835 0, 836 0, 837 0, 838 0, 839 /* 79*/ 0, 840 0, 841 0, 842 0, 843 MatLoad_BlockMat, 844 /* 84*/ 0, 845 0, 846 0, 847 0, 848 0, 849 /* 89*/ 0, 850 0, 851 0, 852 0, 853 0, 854 /* 94*/ 0, 855 0, 856 0, 857 0, 858 0, 859 /* 99*/ 0, 860 0, 861 0, 862 0, 863 0, 864 /*104*/ 0, 865 0, 866 0, 867 0, 868 0, 869 /*109*/ 0, 870 0, 871 0, 872 0, 873 0, 874 /*114*/ 0, 875 0, 876 0, 877 0, 878 0, 879 /*119*/ 0, 880 0, 881 0, 882 0, 883 0, 884 /*124*/ 0, 885 0, 886 0, 887 0, 888 0, 889 /*129*/ 0, 890 0, 891 0, 892 0, 893 0, 894 /*134*/ 0, 895 0, 896 0, 897 0, 898 0, 899 /*139*/ 0, 900 0 901 }; 902 903 #undef __FUNCT__ 904 #define __FUNCT__ "MatBlockMatSetPreallocation" 905 /*@C 906 MatBlockMatSetPreallocation - For good matrix assembly performance 907 the user should preallocate the matrix storage by setting the parameter nz 908 (or the array nnz). By setting these parameters accurately, performance 909 during matrix assembly can be increased by more than a factor of 50. 910 911 Collective on MPI_Comm 912 913 Input Parameters: 914 + B - The matrix 915 . bs - size of each block in matrix 916 . nz - number of nonzeros per block row (same for all rows) 917 - nnz - array containing the number of nonzeros in the various block rows 918 (possibly different for each row) or NULL 919 920 Notes: 921 If nnz is given then nz is ignored 922 923 Specify the preallocated storage with either nz or nnz (not both). 924 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 925 allocation. For large problems you MUST preallocate memory or you 926 will get TERRIBLE performance, see the users' manual chapter on matrices. 927 928 Level: intermediate 929 930 .seealso: MatCreate(), MatCreateBlockMat(), MatSetValues() 931 932 @*/ 933 PetscErrorCode MatBlockMatSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 934 { 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 ierr = PetscTryMethod(B,"MatBlockMatSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 939 PetscFunctionReturn(0); 940 } 941 942 #undef __FUNCT__ 943 #define __FUNCT__ "MatBlockMatSetPreallocation_BlockMat" 944 PetscErrorCode MatBlockMatSetPreallocation_BlockMat(Mat A,PetscInt bs,PetscInt nz,PetscInt *nnz) 945 { 946 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 947 PetscErrorCode ierr; 948 PetscInt i; 949 950 PetscFunctionBegin; 951 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 952 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 953 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 954 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 955 ierr = PetscLayoutGetBlockSize(A->rmap,&bs);CHKERRQ(ierr); 956 957 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 958 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 959 if (nnz) { 960 for (i=0; i<A->rmap->n/bs; i++) { 961 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]); 962 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); 963 } 964 } 965 bmat->mbs = A->rmap->n/bs; 966 967 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->right);CHKERRQ(ierr); 968 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->middle);CHKERRQ(ierr); 969 ierr = VecCreateSeq(PETSC_COMM_SELF,bs,&bmat->left);CHKERRQ(ierr); 970 971 if (!bmat->imax) { 972 ierr = PetscMalloc2(A->rmap->n,PetscInt,&bmat->imax,A->rmap->n,PetscInt,&bmat->ilen);CHKERRQ(ierr); 973 ierr = PetscLogObjectMemory(A,2*A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 974 } 975 if (nnz) { 976 nz = 0; 977 for (i=0; i<A->rmap->n/A->rmap->bs; i++) { 978 bmat->imax[i] = nnz[i]; 979 nz += nnz[i]; 980 } 981 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently requires block row by row preallocation"); 982 983 /* bmat->ilen will count nonzeros in each row so far. */ 984 for (i=0; i<bmat->mbs; i++) bmat->ilen[i] = 0; 985 986 /* allocate the matrix space */ 987 ierr = MatSeqXAIJFreeAIJ(A,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);CHKERRQ(ierr); 988 ierr = PetscMalloc3(nz,Mat,&bmat->a,nz,PetscInt,&bmat->j,A->rmap->n+1,PetscInt,&bmat->i);CHKERRQ(ierr); 989 ierr = PetscLogObjectMemory(A,(A->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 990 bmat->i[0] = 0; 991 for (i=1; i<bmat->mbs+1; i++) { 992 bmat->i[i] = bmat->i[i-1] + bmat->imax[i-1]; 993 } 994 bmat->singlemalloc = PETSC_TRUE; 995 bmat->free_a = PETSC_TRUE; 996 bmat->free_ij = PETSC_TRUE; 997 998 bmat->nz = 0; 999 bmat->maxnz = nz; 1000 A->info.nz_unneeded = (double)bmat->maxnz; 1001 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1002 PetscFunctionReturn(0); 1003 } 1004 1005 /*MC 1006 MATBLOCKMAT - A matrix that is defined by a set of Mat's that represents a sparse block matrix 1007 consisting of (usually) sparse blocks. 1008 1009 Level: advanced 1010 1011 .seealso: MatCreateBlockMat() 1012 1013 M*/ 1014 1015 #undef __FUNCT__ 1016 #define __FUNCT__ "MatCreate_BlockMat" 1017 PETSC_EXTERN PetscErrorCode MatCreate_BlockMat(Mat A) 1018 { 1019 Mat_BlockMat *b; 1020 PetscErrorCode ierr; 1021 1022 PetscFunctionBegin; 1023 ierr = PetscNewLog(A,Mat_BlockMat,&b);CHKERRQ(ierr); 1024 A->data = (void*)b; 1025 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1026 1027 A->assembled = PETSC_TRUE; 1028 A->preallocated = PETSC_FALSE; 1029 ierr = PetscObjectChangeTypeName((PetscObject)A,MATBLOCKMAT);CHKERRQ(ierr); 1030 1031 ierr = PetscObjectComposeFunction((PetscObject)A,"MatBlockMatSetPreallocation_C","MatBlockMatSetPreallocation_BlockMat",MatBlockMatSetPreallocation_BlockMat);CHKERRQ(ierr); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 #undef __FUNCT__ 1036 #define __FUNCT__ "MatCreateBlockMat" 1037 /*@C 1038 MatCreateBlockMat - Creates a new matrix based sparse Mat storage 1039 1040 Collective on MPI_Comm 1041 1042 Input Parameters: 1043 + comm - MPI communicator 1044 . m - number of rows 1045 . n - number of columns 1046 . bs - size of each submatrix 1047 . nz - expected maximum number of nonzero blocks in row (use PETSC_DEFAULT if not known) 1048 - nnz - expected number of nonzers per block row if known (use NULL otherwise) 1049 1050 1051 Output Parameter: 1052 . A - the matrix 1053 1054 Level: intermediate 1055 1056 PETSc requires that matrices and vectors being used for certain 1057 operations are partitioned accordingly. For example, when 1058 creating a bmat matrix, A, that supports parallel matrix-vector 1059 products using MatMult(A,x,y) the user should set the number 1060 of local matrix rows to be the number of local elements of the 1061 corresponding result vector, y. Note that this is information is 1062 required for use of the matrix interface routines, even though 1063 the bmat matrix may not actually be physically partitioned. 1064 For example, 1065 1066 .keywords: matrix, bmat, create 1067 1068 .seealso: MATBLOCKMAT 1069 @*/ 1070 PetscErrorCode MatCreateBlockMat(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt bs,PetscInt nz,PetscInt *nnz, Mat *A) 1071 { 1072 PetscErrorCode ierr; 1073 1074 PetscFunctionBegin; 1075 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1076 ierr = MatSetSizes(*A,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1077 ierr = MatSetType(*A,MATBLOCKMAT);CHKERRQ(ierr); 1078 ierr = MatBlockMatSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1079 PetscFunctionReturn(0); 1080 } 1081 1082 1083 1084