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*/