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