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