1 #define PETSCMAT_DLL 2 3 /* 4 This provides a matrix that consists of Mats 5 */ 6 7 #include "src/mat/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 #define CHUNKSIZE 15 12 13 typedef struct { 14 SEQAIJHEADER(Mat); 15 SEQBAIJHEADER; 16 Mat *diags; 17 18 Vec left,right,middle,workb; /* dummy vectors to perform local parts of product */ 19 } Mat_BlockMat; 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "MatRelax_BlockMat_Symmetric" 23 PetscErrorCode MatRelax_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_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 40 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat"); 41 if (omega != 1.0) SETERRQ(PETSC_ERR_SUP,"No support yet for omega not equal to 1.0"); 42 if (fshift) SETERRQ(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_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]], MATORDERING_ND,&row,&col);CHKERRQ(ierr); 51 ierr = MatCholeskyFactorSymbolic(a->a[a->diag[i]],row,&info,a->diags+i);CHKERRQ(ierr); 52 ierr = MatCholeskyFactorNumeric(a->a[a->diag[i]],&info,a->diags+i);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 = VecGetArray(a->workb,(PetscScalar**)&b);CHKERRQ(ierr); 65 66 /* need to add code for when initial guess is zero, see MatRelax_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 MatRelax_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 = VecRestoreArray(a->workb,(PetscScalar**)&b);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 #undef __FUNCT__ 129 #define __FUNCT__ "MatRelax_BlockMat" 130 PetscErrorCode MatRelax_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_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 147 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat"); 148 if (omega != 1.0) SETERRQ(PETSC_ERR_SUP,"No support yet for omega not equal to 1.0"); 149 if (fshift) SETERRQ(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]], MATORDERING_ND,&row,&col);CHKERRQ(ierr); 156 ierr = MatLUFactorSymbolic(a->a[a->diag[i]],row,col,&info,a->diags+i);CHKERRQ(ierr); 157 ierr = MatLUFactorNumeric(a->a[a->diag[i]],&info,a->diags+i);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 = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 167 168 /* need to add code for when initial guess is zero, see MatRelax_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 = VecRestoreArray(bb,(PetscScalar**)&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 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_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_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) { 262 value = v[l + k*n]; 263 } else { 264 value = v[k + l*m]; 265 } 266 if (col <= lastcol) low = 0; else high = nrow; 267 lastcol = col; 268 while (high-low > 7) { 269 t = (low+high)/2; 270 if (rp[t] > bcol) high = t; 271 else low = t; 272 } 273 for (i=low; i<high; i++) { 274 if (rp[i] > bcol) break; 275 if (rp[i] == bcol) { 276 /* printf("row %d col %d found i %d\n",brow,bcol,i);*/ 277 goto noinsert1; 278 } 279 } 280 if (nonew == 1) goto noinsert1; 281 if (nonew == -1) SETERRQ2(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 /* printf("N %d i %d\n",N,i);*/ 286 for (ii=N; ii>=i; ii--) { 287 rp[ii+1] = rp[ii]; 288 ap[ii+1] = ap[ii]; 289 } 290 if (N>=i) ap[i] = 0; 291 rp[i] = bcol; 292 a->nz++; 293 noinsert1:; 294 if (!*(ap+i)) { 295 /* printf("create matrix at i %d rw %d col %d\n",i,brow,bcol);*/ 296 if (A->symmetric && brow == bcol) { 297 /* don't use SBAIJ since want to reorder in sparse factorization */ 298 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,0,ap+i);CHKERRQ(ierr); 299 } else { 300 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,0,ap+i);CHKERRQ(ierr); 301 } 302 } 303 /* printf("numerical value at i %d row %d col %d cidx %d ridx %d value %g\n",i,brow,bcol,cidx,ridx,value);*/ 304 ierr = MatSetValues(ap[i],1,&ridx,1,&cidx,&value,is);CHKERRQ(ierr); 305 low = i; 306 } 307 /* printf("nrow for row %d %d\n",nrow,brow);*/ 308 ailen[brow] = nrow; 309 } 310 A->same_nonzero = PETSC_FALSE; 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "MatLoad_BlockMat" 316 PetscErrorCode MatLoad_BlockMat(PetscViewer viewer, MatType type,Mat *A) 317 { 318 PetscErrorCode ierr; 319 Mat tmpA; 320 PetscInt i,m,n,bs = 1,ncols; 321 const PetscInt *cols; 322 const PetscScalar *values; 323 PetscTruth flg; 324 325 PetscFunctionBegin; 326 ierr = MatLoad_SeqAIJ(viewer,MATSEQAIJ,&tmpA);CHKERRQ(ierr); 327 328 ierr = MatGetLocalSize(tmpA,&m,&n);CHKERRQ(ierr); 329 ierr = PetscOptionsBegin(PETSC_COMM_SELF,PETSC_NULL,"Options for loading BlockMat matrix","Mat");CHKERRQ(ierr); 330 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,PETSC_NULL);CHKERRQ(ierr); 331 ierr = PetscOptionsEnd();CHKERRQ(ierr); 332 ierr = MatCreateBlockMat(PETSC_COMM_SELF,m,n,bs,A);CHKERRQ(ierr); 333 ierr = PetscOptionsBegin(PETSC_COMM_SELF,PETSC_NULL,"Options for loading BlockMat matrix","Mat");CHKERRQ(ierr); 334 ierr = PetscOptionsName("-matload_symmetric","Store the matrix as symmetric","MatLoad",&flg);CHKERRQ(ierr); 335 if (flg) { 336 ierr = MatSetOption(*A,MAT_SYMMETRIC);CHKERRQ(ierr); 337 } 338 ierr = PetscOptionsEnd();CHKERRQ(ierr); 339 for (i=0; i<m; i++) { 340 ierr = MatGetRow(tmpA,i,&ncols,&cols,&values);CHKERRQ(ierr); 341 ierr = MatSetValues(*A,1,&i,ncols,cols,values,INSERT_VALUES);CHKERRQ(ierr); 342 ierr = MatRestoreRow(tmpA,i,&ncols,&cols,&values);CHKERRQ(ierr); 343 } 344 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 345 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 #undef __FUNCT__ 350 #define __FUNCT__ "MatView_BlockMat" 351 PetscErrorCode MatView_BlockMat(Mat A,PetscViewer viewer) 352 { 353 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 354 PetscErrorCode ierr; 355 const char *name; 356 PetscViewerFormat format; 357 358 PetscFunctionBegin; 359 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 360 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 361 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 362 ierr = PetscViewerASCIIPrintf(viewer,"Nonzero block matrices = %D \n",a->nz);CHKERRQ(ierr); 363 } 364 PetscFunctionReturn(0); 365 } 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "MatDestroy_BlockMat" 369 PetscErrorCode MatDestroy_BlockMat(Mat mat) 370 { 371 PetscErrorCode ierr; 372 Mat_BlockMat *bmat = (Mat_BlockMat*)mat->data; 373 PetscInt i; 374 375 PetscFunctionBegin; 376 if (bmat->right) { 377 ierr = VecDestroy(bmat->right);CHKERRQ(ierr); 378 } 379 if (bmat->left) { 380 ierr = VecDestroy(bmat->left);CHKERRQ(ierr); 381 } 382 if (bmat->middle) { 383 ierr = VecDestroy(bmat->middle);CHKERRQ(ierr); 384 } 385 if (bmat->workb) { 386 ierr = VecDestroy(bmat->workb);CHKERRQ(ierr); 387 } 388 if (bmat->diags) { 389 for (i=0; i<mat->rmap.n/mat->rmap.bs; i++) { 390 if (bmat->diags[i]) {ierr = MatDestroy(bmat->diags[i]);CHKERRQ(ierr);} 391 } 392 } 393 if (bmat->a) { 394 for (i=0; i<bmat->nz; i++) { 395 if (bmat->a[i]) {ierr = MatDestroy(bmat->a[i]);CHKERRQ(ierr);} 396 } 397 } 398 ierr = MatSeqXAIJFreeAIJ(mat,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);CHKERRQ(ierr); 399 ierr = PetscFree(bmat);CHKERRQ(ierr); 400 PetscFunctionReturn(0); 401 } 402 403 #undef __FUNCT__ 404 #define __FUNCT__ "MatMult_BlockMat" 405 PetscErrorCode MatMult_BlockMat(Mat A,Vec x,Vec y) 406 { 407 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 408 PetscErrorCode ierr; 409 PetscScalar *xx,*yy; 410 PetscInt *aj,i,*ii,jrow,m = A->rmap.n/A->rmap.bs,bs = A->rmap.bs,n,j; 411 Mat *aa; 412 413 PetscFunctionBegin; 414 CHKMEMQ; 415 /* 416 Standard CSR multiply except each entry is a Mat 417 */ 418 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 419 420 ierr = VecSet(y,0.0);CHKERRQ(ierr); 421 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 422 aj = bmat->j; 423 aa = bmat->a; 424 ii = bmat->i; 425 for (i=0; i<m; i++) { 426 jrow = ii[i]; 427 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 428 n = ii[i+1] - jrow; 429 for (j=0; j<n; j++) { 430 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 431 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 432 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 433 jrow++; 434 } 435 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 436 } 437 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 438 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 439 CHKMEMQ; 440 PetscFunctionReturn(0); 441 } 442 443 #undef __FUNCT__ 444 #define __FUNCT__ "MatMult_BlockMat_Symmetric" 445 PetscErrorCode MatMult_BlockMat_Symmetric(Mat A,Vec x,Vec y) 446 { 447 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 448 PetscErrorCode ierr; 449 PetscScalar *xx,*yy; 450 PetscInt *aj,i,*ii,jrow,m = A->rmap.n/A->rmap.bs,bs = A->rmap.bs,n,j; 451 Mat *aa; 452 453 PetscFunctionBegin; 454 CHKMEMQ; 455 /* 456 Standard CSR multiply except each entry is a Mat 457 */ 458 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 459 460 ierr = VecSet(y,0.0);CHKERRQ(ierr); 461 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 462 aj = bmat->j; 463 aa = bmat->a; 464 ii = bmat->i; 465 for (i=0; i<m; i++) { 466 jrow = ii[i]; 467 n = ii[i+1] - jrow; 468 ierr = VecPlaceArray(bmat->left,yy + bs*i);CHKERRQ(ierr); 469 ierr = VecPlaceArray(bmat->middle,xx + bs*i);CHKERRQ(ierr); 470 /* if we ALWAYS required a diagonal entry then could remove this if test */ 471 if (aj[jrow] == i) { 472 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); 473 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 474 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 475 jrow++; 476 n--; 477 } 478 for (j=0; j<n; j++) { 479 ierr = VecPlaceArray(bmat->right,xx + bs*aj[jrow]);CHKERRQ(ierr); /* upper triangular part */ 480 ierr = MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);CHKERRQ(ierr); 481 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 482 483 ierr = VecPlaceArray(bmat->right,yy + bs*aj[jrow]);CHKERRQ(ierr); /* lower triangular part */ 484 ierr = MatMultTransposeAdd(aa[jrow],bmat->middle,bmat->right,bmat->right);CHKERRQ(ierr); 485 ierr = VecResetArray(bmat->right);CHKERRQ(ierr); 486 jrow++; 487 } 488 ierr = VecResetArray(bmat->left);CHKERRQ(ierr); 489 ierr = VecResetArray(bmat->middle);CHKERRQ(ierr); 490 } 491 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 492 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 493 CHKMEMQ; 494 PetscFunctionReturn(0); 495 } 496 497 #undef __FUNCT__ 498 #define __FUNCT__ "MatMultAdd_BlockMat" 499 PetscErrorCode MatMultAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 500 { 501 PetscFunctionBegin; 502 PetscFunctionReturn(0); 503 } 504 505 #undef __FUNCT__ 506 #define __FUNCT__ "MatMultTranspose_BlockMat" 507 PetscErrorCode MatMultTranspose_BlockMat(Mat A,Vec x,Vec y) 508 { 509 PetscFunctionBegin; 510 PetscFunctionReturn(0); 511 } 512 513 #undef __FUNCT__ 514 #define __FUNCT__ "MatMultTransposeAdd_BlockMat" 515 PetscErrorCode MatMultTransposeAdd_BlockMat(Mat A,Vec x,Vec y,Vec z) 516 { 517 PetscFunctionBegin; 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNCT__ 522 #define __FUNCT__ "MatSetBlockSize_BlockMat" 523 PetscErrorCode MatSetBlockSize_BlockMat(Mat A,PetscInt bs) 524 { 525 Mat_BlockMat *bmat = (Mat_BlockMat*)A->data; 526 PetscErrorCode ierr; 527 PetscInt nz = 10,i; 528 529 PetscFunctionBegin; 530 if (A->rmap.n % bs) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Blocksize %D does not divide number of rows %D",bs,A->rmap.n); 531 if (A->cmap.n % bs) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Blocksize %D does not divide number of columns %D",bs,A->cmap.n); 532 A->rmap.bs = A->cmap.bs = bs; 533 534 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,PETSC_NULL,&bmat->right);CHKERRQ(ierr); 535 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,PETSC_NULL,&bmat->middle);CHKERRQ(ierr); 536 ierr = VecCreateSeq(PETSC_COMM_SELF,bs,&bmat->left);CHKERRQ(ierr); 537 538 ierr = PetscMalloc2(A->rmap.n,PetscInt,&bmat->imax,A->rmap.n,PetscInt,&bmat->ilen);CHKERRQ(ierr); 539 for (i=0; i<A->rmap.n; i++) bmat->imax[i] = nz; 540 nz = nz*A->rmap.n; 541 542 bmat->mbs = A->rmap.n/A->rmap.bs; 543 544 /* bmat->ilen will count nonzeros in each row so far. */ 545 for (i=0; i<bmat->mbs; i++) { bmat->ilen[i] = 0;} 546 547 /* allocate the matrix space */ 548 ierr = PetscMalloc3(nz,Mat,&bmat->a,nz,PetscInt,&bmat->j,A->rmap.n+1,PetscInt,&bmat->i);CHKERRQ(ierr); 549 bmat->i[0] = 0; 550 for (i=1; i<bmat->mbs+1; i++) { 551 bmat->i[i] = bmat->i[i-1] + bmat->imax[i-1]; 552 } 553 bmat->singlemalloc = PETSC_TRUE; 554 bmat->free_a = PETSC_TRUE; 555 bmat->free_ij = PETSC_TRUE; 556 557 bmat->nz = 0; 558 bmat->maxnz = nz; 559 A->info.nz_unneeded = (double)bmat->maxnz; 560 561 PetscFunctionReturn(0); 562 } 563 564 /* 565 Adds diagonal pointers to sparse matrix structure. 566 */ 567 #undef __FUNCT__ 568 #define __FUNCT__ "MatMarkDiagonal_BlockMat" 569 PetscErrorCode MatMarkDiagonal_BlockMat(Mat A) 570 { 571 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 572 PetscErrorCode ierr; 573 PetscInt i,j,mbs = A->rmap.n/A->rmap.bs; 574 575 PetscFunctionBegin; 576 if (!a->diag) { 577 ierr = PetscMalloc(mbs*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 578 } 579 for (i=0; i<mbs; i++) { 580 a->diag[i] = a->i[i+1]; 581 for (j=a->i[i]; j<a->i[i+1]; j++) { 582 if (a->j[j] == i) { 583 a->diag[i] = j; 584 break; 585 } 586 } 587 } 588 PetscFunctionReturn(0); 589 } 590 591 #undef __FUNCT__ 592 #define __FUNCT__ "MatGetSubMatrix_BlockMat" 593 PetscErrorCode MatGetSubMatrix_BlockMat(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 594 { 595 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 596 Mat_SeqAIJ *c; 597 PetscErrorCode ierr; 598 PetscInt i,k,first,step,lensi,nrows,ncols; 599 PetscInt *j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 600 PetscScalar *a_new,value; 601 Mat C,*aa = a->a; 602 PetscTruth stride,equal; 603 604 PetscFunctionBegin; 605 ierr = ISEqual(isrow,iscol,&equal);CHKERRQ(ierr); 606 if (!equal) SETERRQ(PETSC_ERR_SUP,"Only for idential column and row indices"); 607 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 608 if (!stride) SETERRQ(PETSC_ERR_SUP,"Only for stride indices"); 609 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 610 if (step != A->rmap.bs) SETERRQ(PETSC_ERR_SUP,"Can only select one entry from each block"); 611 612 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 613 ncols = nrows; 614 615 /* create submatrix */ 616 if (scall == MAT_REUSE_MATRIX) { 617 PetscInt n_cols,n_rows; 618 C = *B; 619 ierr = MatGetSize(C,&n_rows,&n_cols);CHKERRQ(ierr); 620 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 621 ierr = MatZeroEntries(C);CHKERRQ(ierr); 622 } else { 623 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 624 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 625 if (A->symmetric) { 626 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 627 } else { 628 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 629 } 630 ierr = MatSeqAIJSetPreallocation(C,0,ailen);CHKERRQ(ierr); 631 ierr = MatSeqSBAIJSetPreallocation(C,1,0,ailen);CHKERRQ(ierr); 632 } 633 c = (Mat_SeqAIJ*)C->data; 634 635 /* loop over rows inserting into submatrix */ 636 a_new = c->a; 637 j_new = c->j; 638 i_new = c->i; 639 640 for (i=0; i<nrows; i++) { 641 ii = ai[i]; 642 lensi = ailen[i]; 643 for (k=0; k<lensi; k++) { 644 *j_new++ = *aj++; 645 ierr = MatGetValue(*aa++,first,first,value);CHKERRQ(ierr); 646 *a_new++ = value; 647 } 648 i_new[i+1] = i_new[i] + lensi; 649 c->ilen[i] = lensi; 650 } 651 652 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 653 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 654 *B = C; 655 PetscFunctionReturn(0); 656 } 657 658 #undef __FUNCT__ 659 #define __FUNCT__ "MatAssemblyEnd_BlockMat" 660 PetscErrorCode MatAssemblyEnd_BlockMat(Mat A,MatAssemblyType mode) 661 { 662 Mat_BlockMat *a = (Mat_BlockMat*)A->data; 663 PetscErrorCode ierr; 664 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 665 PetscInt m = a->mbs,*ip,N,*ailen = a->ilen,rmax = 0; 666 Mat *aa = a->a,*ap; 667 668 PetscFunctionBegin; 669 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 670 671 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 672 for (i=1; i<m; i++) { 673 /* move each row back by the amount of empty slots (fshift) before it*/ 674 fshift += imax[i-1] - ailen[i-1]; 675 rmax = PetscMax(rmax,ailen[i]); 676 if (fshift) { 677 ip = aj + ai[i] ; 678 ap = aa + ai[i] ; 679 N = ailen[i]; 680 for (j=0; j<N; j++) { 681 ip[j-fshift] = ip[j]; 682 ap[j-fshift] = ap[j]; 683 } 684 } 685 ai[i] = ai[i-1] + ailen[i-1]; 686 } 687 if (m) { 688 fshift += imax[m-1] - ailen[m-1]; 689 ai[m] = ai[m-1] + ailen[m-1]; 690 } 691 /* reset ilen and imax for each row */ 692 for (i=0; i<m; i++) { 693 ailen[i] = imax[i] = ai[i+1] - ai[i]; 694 } 695 a->nz = ai[m]; 696 for (i=0; i<a->nz; i++) { 697 #if defined(PETSC_USE_DEBUG) 698 if (!aa[i]) SETERRQ3(PETSC_ERR_PLIB,"Null matrix at location %D column %D nz %D",i,aj[i],a->nz); 699 #endif 700 ierr = MatAssemblyBegin(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 701 ierr = MatAssemblyEnd(aa[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 702 } 703 CHKMEMQ; 704 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap.n,fshift,a->nz);CHKERRQ(ierr); 705 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 706 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 707 a->reallocs = 0; 708 A->info.nz_unneeded = (double)fshift; 709 a->rmax = rmax; 710 711 A->same_nonzero = PETSC_TRUE; 712 ierr = MatMarkDiagonal_BlockMat(A);CHKERRQ(ierr); 713 PetscFunctionReturn(0); 714 } 715 716 #undef __FUNCT__ 717 #define __FUNCT__ "MatSetOption_BlockMat" 718 PetscErrorCode MatSetOption_BlockMat(Mat A,MatOption opt) 719 { 720 PetscFunctionBegin; 721 if (opt == MAT_SYMMETRIC) { 722 A->ops->relax = MatRelax_BlockMat_Symmetric; 723 A->ops->mult = MatMult_BlockMat_Symmetric; 724 } else { 725 PetscInfo1(A,"Unused matrix option %s\n",MatOptions[opt]); 726 } 727 PetscFunctionReturn(0); 728 } 729 730 731 static struct _MatOps MatOps_Values = {MatSetValues_BlockMat, 732 0, 733 0, 734 MatMult_BlockMat, 735 /* 4*/ MatMultAdd_BlockMat, 736 MatMultTranspose_BlockMat, 737 MatMultTransposeAdd_BlockMat, 738 0, 739 0, 740 0, 741 /*10*/ 0, 742 0, 743 0, 744 MatRelax_BlockMat, 745 0, 746 /*15*/ 0, 747 0, 748 0, 749 0, 750 0, 751 /*20*/ 0, 752 MatAssemblyEnd_BlockMat, 753 0, 754 MatSetOption_BlockMat, 755 0, 756 /*25*/ 0, 757 0, 758 0, 759 0, 760 0, 761 /*30*/ 0, 762 0, 763 0, 764 0, 765 0, 766 /*35*/ 0, 767 0, 768 0, 769 0, 770 0, 771 /*40*/ 0, 772 0, 773 0, 774 0, 775 0, 776 /*45*/ 0, 777 0, 778 0, 779 0, 780 0, 781 /*50*/ MatSetBlockSize_BlockMat, 782 0, 783 0, 784 0, 785 0, 786 /*55*/ 0, 787 0, 788 0, 789 0, 790 0, 791 /*60*/ MatGetSubMatrix_BlockMat, 792 MatDestroy_BlockMat, 793 MatView_BlockMat, 794 0, 795 0, 796 /*65*/ 0, 797 0, 798 0, 799 0, 800 0, 801 /*70*/ 0, 802 0, 803 0, 804 0, 805 0, 806 /*75*/ 0, 807 0, 808 0, 809 0, 810 0, 811 /*80*/ 0, 812 0, 813 0, 814 0, 815 MatLoad_BlockMat, 816 /*85*/ 0, 817 0, 818 0, 819 0, 820 0, 821 /*90*/ 0, 822 0, 823 0, 824 0, 825 0, 826 /*95*/ 0, 827 0, 828 0, 829 0}; 830 831 /*MC 832 MATBLOCKMAT - A matrix that is defined by a set of Mat's that represents a sparse block matrix 833 consisting of (usually) sparse blocks. 834 835 Level: advanced 836 837 .seealso: MatCreateBlockMat() 838 839 M*/ 840 841 EXTERN_C_BEGIN 842 #undef __FUNCT__ 843 #define __FUNCT__ "MatCreate_BlockMat" 844 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_BlockMat(Mat A) 845 { 846 Mat_BlockMat *b; 847 PetscErrorCode ierr; 848 849 PetscFunctionBegin; 850 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 851 ierr = PetscNew(Mat_BlockMat,&b);CHKERRQ(ierr); 852 853 A->data = (void*)b; 854 855 ierr = PetscMapInitialize(A->comm,&A->rmap);CHKERRQ(ierr); 856 ierr = PetscMapInitialize(A->comm,&A->cmap);CHKERRQ(ierr); 857 858 A->assembled = PETSC_TRUE; 859 A->preallocated = PETSC_FALSE; 860 ierr = PetscObjectChangeTypeName((PetscObject)A,MATBLOCKMAT);CHKERRQ(ierr); 861 862 ierr = PetscOptionsBegin(A->comm,A->prefix,"Matrix Option","Mat");CHKERRQ(ierr); 863 ierr = PetscOptionsEnd(); 864 865 PetscFunctionReturn(0); 866 } 867 EXTERN_C_END 868 869 #undef __FUNCT__ 870 #define __FUNCT__ "MatCreateBlockMat" 871 /*@C 872 MatCreateBlockMat - Creates a new matrix based sparse Mat storage 873 874 Collective on MPI_Comm 875 876 Input Parameters: 877 + comm - MPI communicator 878 . m - number of rows 879 . n - number of columns 880 - bs - size of each submatrix 881 882 883 Output Parameter: 884 . A - the matrix 885 886 Level: intermediate 887 888 PETSc requires that matrices and vectors being used for certain 889 operations are partitioned accordingly. For example, when 890 creating a bmat matrix, A, that supports parallel matrix-vector 891 products using MatMult(A,x,y) the user should set the number 892 of local matrix rows to be the number of local elements of the 893 corresponding result vector, y. Note that this is information is 894 required for use of the matrix interface routines, even though 895 the bmat matrix may not actually be physically partitioned. 896 For example, 897 898 .keywords: matrix, bmat, create 899 900 .seealso: MATBLOCKMAT 901 @*/ 902 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateBlockMat(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt bs,Mat *A) 903 { 904 PetscErrorCode ierr; 905 906 PetscFunctionBegin; 907 ierr = MatCreate(comm,A);CHKERRQ(ierr); 908 ierr = MatSetSizes(*A,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 909 ierr = MatSetType(*A,MATBLOCKMAT);CHKERRQ(ierr); 910 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 911 PetscFunctionReturn(0); 912 } 913 914 915 916