Actual source code: ex129.c


  2: /*
  3:   Laplacian in 3D. Use for testing MatSolve routines.
  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[] = "This example is for testing different MatSolve routines :MatSolve(), MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), and MatMatSolve().\n\
 13: Example usage: ./ex129 -mat_type aij -dof 2\n\n";

 15: #include <petscdm.h>
 16: #include <petscdmda.h>

 18: extern PetscErrorCode ComputeMatrix(DM,Mat);
 19: extern PetscErrorCode ComputeRHS(DM,Vec);
 20: extern PetscErrorCode ComputeRHSMatrix(PetscInt,PetscInt,Mat*);

 22: int main(int argc,char **args)
 23: {
 24:   PetscMPIInt    size;
 25:   Vec            x,b,y,b1;
 26:   DM             da;
 27:   Mat            A,F,RHS,X,C1;
 28:   MatFactorInfo  info;
 29:   IS             perm,iperm;
 30:   PetscInt       dof =1,M=8,m,n,nrhs;
 31:   PetscScalar    one = 1.0;
 32:   PetscReal      norm,tol = 1000*PETSC_MACHINE_EPSILON;
 33:   PetscBool      InplaceLU=PETSC_FALSE;

 35:   PetscInitialize(&argc,&args,(char*)0,help);
 36:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 38:   PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
 39:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);

 41:   DMDACreate(PETSC_COMM_WORLD,&da);
 42:   DMSetDimension(da,3);
 43:   DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
 44:   DMDASetStencilType(da,DMDA_STENCIL_STAR);
 45:   DMDASetSizes(da,M,M,M);
 46:   DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 47:   DMDASetDof(da,dof);
 48:   DMDASetStencilWidth(da,1);
 49:   DMDASetOwnershipRanges(da,NULL,NULL,NULL);
 50:   DMSetMatType(da,MATBAIJ);
 51:   DMSetFromOptions(da);
 52:   DMSetUp(da);

 54:   DMCreateGlobalVector(da,&x);
 55:   DMCreateGlobalVector(da,&b);
 56:   VecDuplicate(b,&y);
 57:   ComputeRHS(da,b);
 58:   VecSet(y,one);
 59:   DMCreateMatrix(da,&A);
 60:   ComputeMatrix(da,A);
 61:   MatGetSize(A,&m,&n);
 62:   nrhs = 2;
 63:   PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
 64:   ComputeRHSMatrix(m,nrhs,&RHS);
 65:   MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X);

 67:   MatGetOrdering(A,MATORDERINGND,&perm,&iperm);

 69:   PetscOptionsGetBool(NULL,NULL,"-inplacelu",&InplaceLU,NULL);
 70:   MatFactorInfoInitialize(&info);
 71:   if (!InplaceLU) {
 72:     MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
 73:     info.fill = 5.0;
 74:     MatLUFactorSymbolic(F,A,perm,iperm,&info);
 75:     MatLUFactorNumeric(F,A,&info);
 76:   } else { /* Test inplace factorization */
 77:     MatDuplicate(A,MAT_COPY_VALUES,&F);
 78:     MatLUFactor(F,perm,iperm,&info);
 79:   }

 81:   VecDuplicate(y,&b1);

 83:   /* MatSolve */
 84:   MatSolve(F,b,x);
 85:   MatMult(A,x,b1);
 86:   VecAXPY(b1,-1.0,b);
 87:   VecNorm(b1,NORM_2,&norm);
 88:   if (norm > tol) {
 89:     PetscPrintf(PETSC_COMM_WORLD,"MatSolve              : Error of norm %g\n",(double)norm);
 90:   }

 92:   /* MatSolveTranspose */
 93:   MatSolveTranspose(F,b,x);
 94:   MatMultTranspose(A,x,b1);
 95:   VecAXPY(b1,-1.0,b);
 96:   VecNorm(b1,NORM_2,&norm);
 97:   if (norm > tol) {
 98:     PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose     : Error of norm %g\n",(double)norm);
 99:   }

101:   /* MatSolveAdd */
102:   MatSolveAdd(F,b,y,x);
103:   MatMult(A,y,b1);
104:   VecScale(b1,-1.0);
105:   MatMultAdd(A,x,b1,b1);
106:   VecAXPY(b1,-1.0,b);
107:   VecNorm(b1,NORM_2,&norm);
108:   if (norm > tol) {
109:     PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd           : Error of norm %g\n",(double)norm);
110:   }

112:   /* MatSolveTransposeAdd */
113:   MatSolveTransposeAdd(F,b,y,x);
114:   MatMultTranspose(A,y,b1);
115:   VecScale(b1,-1.0);
116:   MatMultTransposeAdd(A,x,b1,b1);
117:   VecAXPY(b1,-1.0,b);
118:   VecNorm(b1,NORM_2,&norm);
119:   if (norm > tol) {
120:     PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd  : Error of norm %g\n",(double)norm);
121:   }

123:   /* MatMatSolve */
124:   MatMatSolve(F,RHS,X);
125:   MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
126:   MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
127:   MatNorm(C1,NORM_FROBENIUS,&norm);
128:   if (norm > tol) {
129:     PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve           : Error of norm %g\n",(double)norm);
130:   }

132:   VecDestroy(&x);
133:   VecDestroy(&b);
134:   VecDestroy(&b1);
135:   VecDestroy(&y);
136:   MatDestroy(&A);
137:   MatDestroy(&F);
138:   MatDestroy(&RHS);
139:   MatDestroy(&C1);
140:   MatDestroy(&X);
141:   ISDestroy(&perm);
142:   ISDestroy(&iperm);
143:   DMDestroy(&da);
144:   PetscFinalize();
145:   return 0;
146: }

148: PetscErrorCode ComputeRHS(DM da,Vec b)
149: {
150:   PetscInt       mx,my,mz;
151:   PetscScalar    h;

153:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
154:   h    = 1.0/((mx-1)*(my-1)*(mz-1));
155:   VecSet(b,h);
156:   return 0;
157: }

159: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat *C)
160: {
161:   PetscRandom    rand;
162:   Mat            RHS;
163:   PetscScalar    *array,rval;
164:   PetscInt       i,k;

166:   MatCreate(PETSC_COMM_WORLD,&RHS);
167:   MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
168:   MatSetType(RHS,MATSEQDENSE);
169:   MatSetUp(RHS);

171:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
172:   PetscRandomSetFromOptions(rand);
173:   MatDenseGetArray(RHS,&array);
174:   for (i=0; i<m; i++) {
175:     PetscRandomGetValue(rand,&rval);
176:     array[i] = rval;
177:   }
178:   if (nrhs > 1) {
179:     for (k=1; k<nrhs; k++) {
180:       for (i=0; i<m; i++) {
181:         array[m*k+i] = array[i];
182:       }
183:     }
184:   }
185:   MatDenseRestoreArray(RHS,&array);
186:   MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
187:   MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
188:   *C   = RHS;
189:   PetscRandomDestroy(&rand);
190:   return 0;
191: }

193: PetscErrorCode ComputeMatrix(DM da,Mat B)
194: {
195:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
196:   PetscScalar    *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
197:   MatStencil     row,col;
198:   PetscRandom    rand;

200:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
201:   PetscRandomSetSeed(rand,1);
202:   PetscRandomSetInterval(rand,-.001,.001);
203:   PetscRandomSetFromOptions(rand);

205:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
206:   /* For simplicity, this example only works on mx=my=mz */

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

212:   PetscMalloc1(2*dof*dof+1,&v);
213:   v_neighbor = v + dof*dof;
214:   PetscArrayzero(v,2*dof*dof+1);
215:   k3         = 0;
216:   for (k1=0; k1<dof; k1++) {
217:     for (k2=0; k2<dof; k2++) {
218:       if (k1 == k2) {
219:         v[k3]          = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
220:         v_neighbor[k3] = -HxHydHz;
221:       } else {
222:         PetscRandomGetValue(rand,&r1);
223:         PetscRandomGetValue(rand,&r2);

225:         v[k3]          = r1;
226:         v_neighbor[k3] = r2;
227:       }
228:       k3++;
229:     }
230:   }
231:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

233:   for (k=zs; k<zs+zm; k++) {
234:     for (j=ys; j<ys+ym; j++) {
235:       for (i=xs; i<xs+xm; i++) {
236:         row.i = i; row.j = j; row.k = k;
237:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boundary points */
238:           MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
239:         } else { /* interior points */
240:           /* center */
241:           col.i = i; col.j = j; col.k = k;
242:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);

244:           /* x neighbors */
245:           col.i = i-1; col.j = j; col.k = k;
246:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
247:           col.i = i+1; col.j = j; col.k = k;
248:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);

250:           /* y neighbors */
251:           col.i = i; col.j = j-1; col.k = k;
252:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
253:           col.i = i; col.j = j+1; col.k = k;
254:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);

256:           /* z neighbors */
257:           col.i = i; col.j = j; col.k = k-1;
258:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
259:           col.i = i; col.j = j; col.k = k+1;
260:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
261:         }
262:       }
263:     }
264:   }
265:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
266:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
267:   PetscFree(v);
268:   PetscRandomDestroy(&rand);
269:   return 0;
270: }

272: /*TEST

274:    test:
275:       args: -dm_mat_type aij -dof 1
276:       output_file: output/ex129.out

278:    test:
279:       suffix: 2
280:       args: -dm_mat_type aij -dof 1 -inplacelu
281:       output_file: output/ex129.out

283: TEST*/