Actual source code: ex32.c

  2: /*
  3:   Laplacian in 3D. Use for testing BAIJ matrix. 
  4:   Modeled by the partial differential equation

  6:    - Laplacian u = 1,0 < x,y,z < 1,

  8:    with boundary conditions
  9:    u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
 10: */

 12: static char help[] = "Solves 3D Laplacian using wirebasket based multigrid.\n\n";

 14:  #include petscda.h
 15:  #include petscksp.h
 16:  #include petscmg.h


 23: int main(int argc,char **argv)
 24: {
 26:   KSP            ksp;
 27:   PC             pc;
 28:   Vec            x,b;
 29:   DA             da;
 30:   Mat            A,Atrans;
 31:   PetscInt       dof=1,M=-8;
 32:   PetscTruth     flg,trans=PETSC_FALSE;

 34:   PetscInitialize(&argc,&argv,(char *)0,help);
 35:   PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
 36:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 37:   PetscOptionsGetTruth(PETSC_NULL,"-trans",&trans,PETSC_NULL);

 39:   DACreate(PETSC_COMM_WORLD,&da);
 40:   DASetDim(da,3);
 41:   DASetPeriodicity(da,DA_NONPERIODIC);
 42:   DASetStencilType(da,DA_STENCIL_STAR);
 43:   DASetSizes(da,M,M,M);
 44:   DASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 45:   DASetDof(da,dof);
 46:   DASetStencilWidth(da,1);
 47:   DASetVertexDivision(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
 48:   DASetFromOptions(da);

 50:   DACreateGlobalVector(da,&x);
 51:   DACreateGlobalVector(da,&b);
 52:   ComputeRHS(da,b);
 53:   DAGetMatrix(da,MATBAIJ,&A);
 54:   ComputeMatrix(da,A);


 57:   /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
 58:   MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);
 59:   MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);
 60:   MatScale(A,0.5);
 61:   MatDestroy(Atrans);

 63:   /* Test sbaij matrix */
 64:   flg  = PETSC_FALSE;
 65:   PetscOptionsGetTruth(PETSC_NULL, "-test_sbaij1", &flg,PETSC_NULL);
 66:   if (flg){
 67:     Mat sA;
 68:     MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
 69:     MatDestroy(A);
 70:     A = sA;
 71:   }

 73:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 74:   KSPSetFromOptions(ksp);
 75:   KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
 76:   KSPGetPC(ksp,&pc);
 77:   PCSetDA(pc,da);
 78: 
 79:   if (trans) {
 80:     KSPSolveTranspose(ksp,b,x);
 81:   } else {
 82:     KSPSolve(ksp,b,x);
 83:   }

 85:   /* check final residual */
 86:   flg  = PETSC_FALSE;
 87:   PetscOptionsGetTruth(PETSC_NULL, "-check_final_residual", &flg,PETSC_NULL);
 88:   if (flg){
 89:     Vec            b1;
 90:     PetscReal      norm;
 91:     KSPGetSolution(ksp,&x);
 92:     VecDuplicate(b,&b1);
 93:     MatMult(A,x,b1);
 94:     VecAXPY(b1,-1.0,b);
 95:     VecNorm(b1,NORM_2,&norm);
 96:     PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
 97:     VecDestroy(b1);
 98:   }
 99: 
100:   KSPDestroy(ksp);
101:   VecDestroy(x);
102:   VecDestroy(b);
103:   MatDestroy(A);
104:   DADestroy(da);
105:   PetscFinalize();
106:   return 0;
107: }

111: PetscErrorCode ComputeRHS(DA da,Vec b)
112: {
114:   PetscInt       mx,my,mz;
115:   PetscScalar    h;

118:   DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
119:   h    = 1.0/((mx-1)*(my-1)*(mz-1));
120:   VecSet(b,h);
121:   return(0);
122: }
123: 
126: PetscErrorCode ComputeMatrix(DA da,Mat B)
127: {
129:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
130:   PetscScalar    *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
131:   MatStencil     row,col;
132: 
134:   DAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);
135:   /* For simplicity, this example only works on mx=my=mz */
136:   if ( mx != my || mx != mz) SETERRQ3(1,"This example only works with mx %d = my %d = mz %d\n",mx,my,mz);

138:   Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
139:   HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;

141:   PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
142:   v_neighbor = v + dof*dof;
143:   PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
144:   k3 = 0;
145:   for (k1=0; k1<dof; k1++){
146:     for (k2=0; k2<dof; k2++){
147:       if (k1 == k2){
148:         v[k3]          = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
149:         v_neighbor[k3] = -HxHydHz;
150:       } else {
151:         v[k3] = k1/(dof*dof); ;
152:         v_neighbor[k3] = k2/(dof*dof);
153:       }
154:       k3++;
155:     }
156:   }
157:   DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
158: 
159:   for (k=zs; k<zs+zm; k++){
160:     for (j=ys; j<ys+ym; j++){
161:       for(i=xs; i<xs+xm; i++){
162:         row.i = i; row.j = j; row.k = k;
163:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
164:           MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
165:         } else { /* interior points */
166:           /* center */
167:           col.i = i; col.j = j; col.k = k;
168:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
169: 
170:           /* x neighbors */
171:           col.i = i-1; col.j = j; col.k = k;
172:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
173:           col.i = i+1; col.j = j; col.k = k;
174:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
175: 
176:           /* y neighbors */
177:           col.i = i; col.j = j-1; col.k = k;
178:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
179:           col.i = i; col.j = j+1; col.k = k;
180:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
181: 
182:           /* z neighbors */
183:           col.i = i; col.j = j; col.k = k-1;
184:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
185:           col.i = i; col.j = j; col.k = k+1;
186:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
187:         }
188:       }
189:     }
190:   }
191:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
192:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
193:   PetscFree(v);
194:   return(0);
195: }