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