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 = PetscMalloc1(mbs,&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 = PetscMalloc1(mbs,&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 for (k=0; k<m; k++) { /* loop over added rows */ 240 row = im[k]; 241 brow = row/bs; 242 if (row < 0) continue; 243 #if defined(PETSC_USE_DEBUG) 244 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); 245 #endif 246 rp = aj + ai[brow]; 247 ap = aa + ai[brow]; 248 rmax = imax[brow]; 249 nrow = ailen[brow]; 250 low = 0; 251 high = nrow; 252 for (l=0; l<n; l++) { /* loop over added columns */ 253 if (in[l] < 0) continue; 254 #if defined(PETSC_USE_DEBUG) 255 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); 256 #endif 257 col = in[l]; bcol = col/bs; 258 if (A->symmetric && brow > bcol) continue; 259 ridx = row % bs; cidx = col % bs; 260 if (roworiented) value = v[l + k*n]; 261 else value = v[k + l*m]; 262 263 if (col <= lastcol) low = 0; 264 else high = nrow; 265 lastcol = col; 266 while (high-low > 7) { 267 t = (low+high)/2; 268 if (rp[t] > bcol) high = t; 269 else low = t; 270 } 271 for (i=low; i<high; i++) { 272 if (rp[i] > bcol) break; 273 if (rp[i] == bcol) goto noinsert1; 274 } 275 if (nonew == 1) goto noinsert1; 276 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 277 MatSeqXAIJReallocateAIJ(A,a->mbs,1,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,Mat); 278 N = nrow++ - 1; high++; 279 /* shift up all the later entries in this row */ 280 for (ii=N; ii>=i; ii--) { 281 rp[ii+1] = rp[ii]; 282 ap[ii+1] = ap[ii]; 283 } 284 if (N>=i) ap[i] = 0; 285 rp[i] = bcol; 286 a->nz++; 287 A->nonzerostate++; 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 PetscFunctionReturn(0); 298 } 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "MatLoad_BlockMat" 302 PetscErrorCode MatLoad_BlockMat(Mat newmat, PetscViewer viewer) 303 { 304 PetscErrorCode ierr; 305 Mat tmpA; 306 PetscInt i,j,m,n,bs = 1,ncols,*lens,currentcol,mbs,**ii,*ilens,nextcol,*llens,cnt = 0; 307 const PetscInt *cols; 308 const PetscScalar *values; 309 PetscBool flg = PETSC_FALSE,notdone; 310 Mat_SeqAIJ *a; 311 Mat_BlockMat *amat; 312 313 PetscFunctionBegin; 314 /* force binary viewer to load .info file if it has not yet done so */ 315 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 316 ierr = MatCreate(PETSC_COMM_SELF,&tmpA);CHKERRQ(ierr); 317 ierr = MatSetType(tmpA,MATSEQAIJ);CHKERRQ(ierr); 318 ierr = MatLoad_SeqAIJ(tmpA,viewer);CHKERRQ(ierr); 319 320 ierr = MatGetLocalSize(tmpA,&m,&n);CHKERRQ(ierr); 321 ierr = PetscOptionsBegin(PETSC_COMM_SELF,NULL,"Options for loading BlockMat matrix 1","Mat");CHKERRQ(ierr); 322 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 323 ierr = PetscOptionsBool("-matload_symmetric","Store the matrix as symmetric","MatLoad",flg,&flg,NULL);CHKERRQ(ierr); 324 ierr = PetscOptionsEnd();CHKERRQ(ierr); 325 326 /* Determine number of nonzero blocks for each block row */ 327 a = (Mat_SeqAIJ*) tmpA->data; 328 mbs = m/bs; 329 ierr = PetscMalloc3(mbs,&lens,bs,&ii,bs,&ilens);CHKERRQ(ierr); 330 ierr = PetscMemzero(lens,mbs*sizeof(PetscInt));CHKERRQ(ierr); 331 332 for (i=0; i<mbs; i++) { 333 for (j=0; j<bs; j++) { 334 ii[j] = a->j + a->i[i*bs + j]; 335 ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j]; 336 } 337 338 currentcol = -1; 339 notdone = PETSC_TRUE; 340 while (PETSC_TRUE) { 341 notdone = PETSC_FALSE; 342 nextcol = 1000000000; 343 for (j=0; j<bs; j++) { 344 while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) { 345 ii[j]++; 346 ilens[j]--; 347 } 348 if (ilens[j] > 0) { 349 notdone = PETSC_TRUE; 350 nextcol = PetscMin(nextcol,ii[j][0]/bs); 351 } 352 } 353 if (!notdone) break; 354 if (!flg || (nextcol >= i)) lens[i]++; 355 currentcol = nextcol; 356 } 357 } 358 359 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 360 ierr = MatSetSizes(newmat,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 361 } 362 ierr = MatBlockMatSetPreallocation(newmat,bs,0,lens);CHKERRQ(ierr); 363 if (flg) { 364 ierr = MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 365 } 366 amat = (Mat_BlockMat*)(newmat)->data; 367 368 /* preallocate the submatrices */ 369 ierr = PetscMalloc1(bs,&llens);CHKERRQ(ierr); 370 for (i=0; i<mbs; i++) { /* loops for block rows */ 371 for (j=0; j<bs; j++) { 372 ii[j] = a->j + a->i[i*bs + j]; 373 ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j]; 374 } 375 376 currentcol = 1000000000; 377 for (j=0; j<bs; j++) { /* loop over rows in block finding first nonzero block */ 378 if (ilens[j] > 0) { 379 currentcol = PetscMin(currentcol,ii[j][0]/bs); 380 } 381 } 382 383 notdone = PETSC_TRUE; 384 while (PETSC_TRUE) { /* loops over blocks in block row */ 385 386 notdone = PETSC_FALSE; 387 nextcol = 1000000000; 388 ierr = PetscMemzero(llens,bs*sizeof(PetscInt));CHKERRQ(ierr); 389 for (j=0; j<bs; j++) { /* loop over rows in block */ 390 while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) { /* loop over columns in row */ 391 ii[j]++; 392 ilens[j]--; 393 llens[j]++; 394 } 395 if (ilens[j] > 0) { 396 notdone = PETSC_TRUE; 397 nextcol = PetscMin(nextcol,ii[j][0]/bs); 398 } 399 } 400 if (cnt >= amat->maxnz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of blocks found greater than expected %D",cnt); 401 if (!flg || currentcol >= i) { 402 amat->j[cnt] = currentcol; 403 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,llens,amat->a+cnt++);CHKERRQ(ierr); 404 } 405 406 if (!notdone) break; 407 currentcol = nextcol; 408 } 409 amat->ilen[i] = lens[i]; 410 } 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 /* 487 Standard CSR multiply except each entry is a Mat 488 */ 489 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 490 491 ierr = VecSet(y,0.0);CHKERRQ(ierr); 492 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 493 aj = bmat->j; 494 aa = bmat->a; 495 ii = bmat->i; 496 for (i=0; i<m; i++) { 497 jrow = ii[i]; 498 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 499 n = ii[i+1] - jrow; 500 for (j=0; j<n; j++) { 501 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 502 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 503 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 504 jrow++; 505 } 506 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 507 } 508 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 509 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 510 PetscFunctionReturn(0); 511 } 512 513 #undef __FUNCT__ 514 #define __FUNCT__ "MatMult_BlockMat_Symmetric" 515 PetscErrorCode MatMult_BlockMat_Symmetric(Mat A,Vec x,Vec y) 516 { 517 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 518 PetscErrorCode ierr; 519 PetscScalar *xx,*yy; 520 PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j; 521 Mat *aa; 522 523 PetscFunctionBegin; 524 /* 525 Standard CSR multiply except each entry is a Mat 526 */ 527 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 528 529 ierr = VecSet(y,0.0);CHKERRQ(ierr); 530 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 531 aj = bmat->j; 532 aa = bmat->a; 533 ii = bmat->i; 534 for (i=0; i<m; i++) { 535 jrow = ii[i]; 536 n = ii[i+1] - jrow; 537 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 538 ierr = VecPlaceArray(bmat->middle,xx + bs*i);CHKERRQ(ierr); 539 /* if we ALWAYS required a diagonal entry then could remove this if test */ 540 if (aj[jrow] == i) { 541 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 542 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 543 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 544 jrow++; 545 n--; 546 } 547 for (j=0; j<n; j++) { 548 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); /* upper triangular part */ 549 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 550 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 551 552 ierr = VecPlaceArray(bmat->right,yy + bs*aj[jrow]);CHKERRQ(ierr); /* lower triangular part */ 553 ierr = MatMultTransposeAdd(aa[jrow],bmat->middle,bmat->right,bmat->right);CHKERRQ(ierr); 554 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 555 jrow++; 556 } 557 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 558 ierr = VecResetArray(bmat->middle);CHKERRQ(ierr); 559 } 560 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 561 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 #undef __FUNCT__ 566 #define __FUNCT__ "MatMultAdd_BlockMat" 567 PetscErrorCode MatMultAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 568 { 569 PetscFunctionBegin; 570 PetscFunctionReturn(0); 571 } 572 573 #undef __FUNCT__ 574 #define __FUNCT__ "MatMultTranspose_BlockMat" 575 PetscErrorCode MatMultTranspose_BlockMat(Mat A,Vec x,Vec y) 576 { 577 PetscFunctionBegin; 578 PetscFunctionReturn(0); 579 } 580 581 #undef __FUNCT__ 582 #define __FUNCT__ "MatMultTransposeAdd_BlockMat" 583 PetscErrorCode MatMultTransposeAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 584 { 585 PetscFunctionBegin; 586 PetscFunctionReturn(0); 587 } 588 589 /* 590 Adds diagonal pointers to sparse matrix structure. 591 */ 592 #undef __FUNCT__ 593 #define __FUNCT__ "MatMarkDiagonal_BlockMat" 594 PetscErrorCode MatMarkDiagonal_BlockMat(Mat A) 595 { 596 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 597 PetscErrorCode ierr; 598 PetscInt i,j,mbs = A->rmap->n/A->rmap->bs; 599 600 PetscFunctionBegin; 601 if (!a->diag) { 602 ierr = PetscMalloc1(mbs,&a->diag);CHKERRQ(ierr); 603 } 604 for (i=0; i<mbs; i++) { 605 a->diag[i] = a->i[i+1]; 606 for (j=a->i[i]; j<a->i[i+1]; j++) { 607 if (a->j[j] == i) { 608 a->diag[i] = j; 609 break; 610 } 611 } 612 } 613 PetscFunctionReturn(0); 614 } 615 616 #undef __FUNCT__ 617 #define __FUNCT__ "MatGetSubMatrix_BlockMat" 618 PetscErrorCode MatGetSubMatrix_BlockMat(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 619 { 620 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 621 Mat_SeqAIJ *c; 622 PetscErrorCode ierr; 623 PetscInt i,k,first,step,lensi,nrows,ncols; 624 PetscInt *j_new,*i_new,*aj = a->j,*ailen = a->ilen; 625 PetscScalar *a_new; 626 Mat C,*aa = a->a; 627 PetscBool stride,equal; 628 629 PetscFunctionBegin; 630 ierr = ISEqual(isrow,iscol,&equal);CHKERRQ(ierr); 631 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices"); 632 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 633 if (!stride) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for stride indices"); 634 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 635 if (step != A->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only select one entry from each block"); 636 637 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 638 ncols = nrows; 639 640 /* create submatrix */ 641 if (scall == MAT_REUSE_MATRIX) { 642 PetscInt n_cols,n_rows; 643 C = *B; 644 ierr = MatGetSize(C,&n_rows,&n_cols);CHKERRQ(ierr); 645 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 646 ierr = MatZeroEntries(C);CHKERRQ(ierr); 647 } else { 648 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 649 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 650 if (A->symmetric) { 651 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 652 } else { 653 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 654 } 655 ierr = MatSeqAIJSetPreallocation(C,0,ailen);CHKERRQ(ierr); 656 ierr = MatSeqSBAIJSetPreallocation(C,1,0,ailen);CHKERRQ(ierr); 657 } 658 c = (Mat_SeqAIJ*)C->data; 659 660 /* loop over rows inserting into submatrix */ 661 a_new = c->a; 662 j_new = c->j; 663 i_new = c->i; 664 665 for (i=0; i<nrows; i++) { 666 lensi = ailen[i]; 667 for (k=0; k<lensi; k++) { 668 *j_new++ = *aj++; 669 ierr = MatGetValue(*aa++,first,first,a_new++);CHKERRQ(ierr); 670 } 671 i_new[i+1] = i_new[i] + lensi; 672 c->ilen[i] = lensi; 673 } 674 675 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 676 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 677 *B = C; 678 PetscFunctionReturn(0); 679 } 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "MatAssemblyEnd_BlockMat" 683 PetscErrorCode MatAssemblyEnd_BlockMat(Mat A,MatAssemblyType mode) 684 { 685 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 686 PetscErrorCode ierr; 687 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 688 PetscInt m = a->mbs,*ip,N,*ailen = a->ilen,rmax = 0; 689 Mat *aa = a->a,*ap; 690 691 PetscFunctionBegin; 692 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 693 694 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 695 for (i=1; i<m; i++) { 696 /* move each row back by the amount of empty slots (fshift) before it*/ 697 fshift += imax[i-1] - ailen[i-1]; 698 rmax = PetscMax(rmax,ailen[i]); 699 if (fshift) { 700 ip = aj + ai[i]; 701 ap = aa + ai[i]; 702 N = ailen[i]; 703 for (j=0; j<N; j++) { 704 ip[j-fshift] = ip[j]; 705 ap[j-fshift] = ap[j]; 706 } 707 } 708 ai[i] = ai[i-1] + ailen[i-1]; 709 } 710 if (m) { 711 fshift += imax[m-1] - ailen[m-1]; 712 ai[m] = ai[m-1] + ailen[m-1]; 713 } 714 /* reset ilen and imax for each row */ 715 for (i=0; i<m; i++) { 716 ailen[i] = imax[i] = ai[i+1] - ai[i]; 717 } 718 a->nz = ai[m]; 719 for (i=0; i<a->nz; i++) { 720 #if defined(PETSC_USE_DEBUG) 721 if (!aa[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Null matrix at location %D column %D nz %D",i,aj[i],a->nz); 722 #endif 723 ierr = MatAssemblyBegin(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 724 ierr = MatAssemblyEnd(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 725 } 726 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); 727 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 728 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 729 730 A->info.mallocs += a->reallocs; 731 a->reallocs = 0; 732 A->info.nz_unneeded = (double)fshift; 733 a->rmax = rmax; 734 ierr = MatMarkDiagonal_BlockMat(A);CHKERRQ(ierr); 735 PetscFunctionReturn(0); 736 } 737 738 #undef __FUNCT__ 739 #define __FUNCT__ "MatSetOption_BlockMat" 740 PetscErrorCode MatSetOption_BlockMat(Mat A,MatOption opt,PetscBool flg) 741 { 742 PetscFunctionBegin; 743 if (opt == MAT_SYMMETRIC && flg) { 744 A->ops->sor = MatSOR_BlockMat_Symmetric; 745 A->ops->mult = MatMult_BlockMat_Symmetric; 746 } else { 747 PetscInfo1(A,"Unused matrix option %s\n",MatOptions[opt]); 748 } 749 PetscFunctionReturn(0); 750 } 751 752 753 static struct _MatOps MatOps_Values = {MatSetValues_BlockMat, 754 0, 755 0, 756 MatMult_BlockMat, 757 /* 4*/ MatMultAdd_BlockMat, 758 MatMultTranspose_BlockMat, 759 MatMultTransposeAdd_BlockMat, 760 0, 761 0, 762 0, 763 /* 10*/ 0, 764 0, 765 0, 766 MatSOR_BlockMat, 767 0, 768 /* 15*/ 0, 769 0, 770 0, 771 0, 772 0, 773 /* 20*/ 0, 774 MatAssemblyEnd_BlockMat, 775 MatSetOption_BlockMat, 776 0, 777 /* 24*/ 0, 778 0, 779 0, 780 0, 781 0, 782 /* 29*/ 0, 783 0, 784 0, 785 0, 786 0, 787 /* 34*/ 0, 788 0, 789 0, 790 0, 791 0, 792 /* 39*/ 0, 793 0, 794 0, 795 0, 796 0, 797 /* 44*/ 0, 798 0, 799 MatShift_Basic, 800 0, 801 0, 802 /* 49*/ 0, 803 0, 804 0, 805 0, 806 0, 807 /* 54*/ 0, 808 0, 809 0, 810 0, 811 0, 812 /* 59*/ MatGetSubMatrix_BlockMat, 813 MatDestroy_BlockMat, 814 MatView_BlockMat, 815 0, 816 0, 817 /* 64*/ 0, 818 0, 819 0, 820 0, 821 0, 822 /* 69*/ 0, 823 0, 824 0, 825 0, 826 0, 827 /* 74*/ 0, 828 0, 829 0, 830 0, 831 0, 832 /* 79*/ 0, 833 0, 834 0, 835 0, 836 MatLoad_BlockMat, 837 /* 84*/ 0, 838 0, 839 0, 840 0, 841 0, 842 /* 89*/ 0, 843 0, 844 0, 845 0, 846 0, 847 /* 94*/ 0, 848 0, 849 0, 850 0, 851 0, 852 /* 99*/ 0, 853 0, 854 0, 855 0, 856 0, 857 /*104*/ 0, 858 0, 859 0, 860 0, 861 0, 862 /*109*/ 0, 863 0, 864 0, 865 0, 866 0, 867 /*114*/ 0, 868 0, 869 0, 870 0, 871 0, 872 /*119*/ 0, 873 0, 874 0, 875 0, 876 0, 877 /*124*/ 0, 878 0, 879 0, 880 0, 881 0, 882 /*129*/ 0, 883 0, 884 0, 885 0, 886 0, 887 /*134*/ 0, 888 0, 889 0, 890 0, 891 0, 892 /*139*/ 0, 893 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,&bmat->imax,A->rmap->n,&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,&bmat->a,nz,&bmat->j,A->rmap->n+1,&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,&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 in which each block contains a uniform-size sequential Mat object 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 Notes: Matrices of this type are nominally-sparse matrices in which each "entry" is a Mat object. Each Mat must 1051 have the same size and be sequential. The local and global sizes must be compatible with this decomposition. 1052 1053 For matrices containing parallel submatrices and variable block sizes, see MATNEST. 1054 1055 .keywords: matrix, bmat, create 1056 1057 .seealso: MATBLOCKMAT, MatCreateNest() 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