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 <petscdmda.h>
23: int main(int argc,char **args)
24: {
25: PetscErrorCode ierr;
26: PetscMPIInt size;
27: Vec x,b,y,b1;
28: DM da;
29: Mat A,F,RHS,X,C1;
30: MatFactorInfo info;
31: IS perm,iperm;
32: PetscInt dof=1,M=-8,m,n,nrhs;
33: PetscScalar one = 1.0;
34: PetscReal norm;
35: PetscBool InplaceLU=PETSC_FALSE;
37: PetscInitialize(&argc,&args,(char *)0,help);
38: MPI_Comm_size(PETSC_COMM_WORLD,&size);
39: if(size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only\n");
40: PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
41: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
43: DMDACreate(PETSC_COMM_WORLD,&da);
44: DMDASetDim(da,3);
45: DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);
46: DMDASetStencilType(da,DMDA_STENCIL_STAR);
47: DMDASetSizes(da,M,M,M);
48: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
49: DMDASetDof(da,dof);
50: DMDASetStencilWidth(da,1);
51: DMDASetOwnershipRanges(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
52: DMSetFromOptions(da);
53: DMSetUp(da);
55: DMCreateGlobalVector(da,&x);
56: DMCreateGlobalVector(da,&b);
57: VecDuplicate(b,&y);
58: ComputeRHS(da,b);
59: VecSet(y,one);
60: DMGetMatrix(da,MATBAIJ,&A);
61: ComputeMatrix(da,A);
62: MatGetSize(A,&m,&n);
63: nrhs = 2;
64: PetscOptionsGetInt(PETSC_NULL,"-nrhs",&nrhs,PETSC_NULL);
65: ComputeRHSMatrix(m,nrhs,&RHS);
66: MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X);
68: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
69:
70:
71: PetscOptionsGetBool(PETSC_NULL,"-inplacelu",&InplaceLU,PETSC_NULL);
72: MatFactorInfoInitialize(&info);
73: if (!InplaceLU){
74: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
75: info.fill = 5.0;
76: MatLUFactorSymbolic(F,A,perm,iperm,&info);
77: MatLUFactorNumeric(F,A,&info);
78: } else { /* Test inplace factorization */
79: MatDuplicate(A,MAT_COPY_VALUES,&F);
80: /* or create F without DMDA
81: const MatType type;
82: PetscInt i,ncols;
83: const PetscInt *cols;
84: const PetscScalar *vals;
85: MatGetSize(A,&m,&n);
86: MatGetType(A,&type);
87: MatCreate(((PetscObject)A)->comm,&F);
88: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,m,n);
89: MatSetType(F,type);
90: MatSetFromOptions(F);
91: for (i=0; i<m; i++) {
92: MatGetRow(A,i,&ncols,&cols,&vals);
93: MatSetValues(F,1,&i,ncols,cols,vals,INSERT_VALUES);
94: }
95: MatAssemblyBegin(F,MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(F,MAT_FINAL_ASSEMBLY);
97: */
98: MatLUFactor(F,perm,iperm,&info);
99: }
101: VecDuplicate(y,&b1);
102:
103: /* MatSolve */
104: MatSolve(F,b,x);
105: MatMult(A,x,b1);
106: VecAXPY(b1,-1.0,b);
107: VecNorm(b1,NORM_2,&norm);
108: PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %A\n",norm);
109:
110: /* MatSolveTranspose */
111: MatSolveTranspose(F,b,x);
112: MatMultTranspose(A,x,b1);
113: VecAXPY(b1,-1.0,b);
114: VecNorm(b1,NORM_2,&norm);
115: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %A\n",norm);
116:
117: /* MatSolveAdd */
118: MatSolveAdd(F,b,y,x);
119: MatMult(A,y,b1);
120: VecScale(b1,-1.0);
121: MatMultAdd(A,x,b1,b1);
122: VecAXPY(b1,-1.0,b);
123: VecNorm(b1,NORM_2,&norm);
124: PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %A\n",norm);
125:
126: /* MatSolveTransposeAdd */
127: MatSolveTransposeAdd(F,b,y,x);
128: MatMultTranspose(A,y,b1);
129: VecScale(b1,-1.0);
130: MatMultTransposeAdd(A,x,b1,b1);
131: VecAXPY(b1,-1.0,b);
132: VecNorm(b1,NORM_2,&norm);
133: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %A\n",norm);
134:
135: /* MatMatSolve */
136: MatMatSolve(F,RHS,X);
137: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
138: MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
139: MatNorm(C1,NORM_FROBENIUS,&norm);
140: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %A\n",norm);
143: VecDestroy(&x);
144: VecDestroy(&b);
145: VecDestroy(&b1);
146: VecDestroy(&y);
147: MatDestroy(&A);
148: MatDestroy(&F);
149: MatDestroy(&RHS);
150: MatDestroy(&C1);
151: MatDestroy(&X);
152: ISDestroy(&perm);
153: ISDestroy(&iperm);
154: DMDestroy(&da);
155: PetscFinalize();
156: return 0;
157: }
161: PetscErrorCode ComputeRHS(DM da,Vec b)
162: {
164: PetscInt mx,my,mz;
165: PetscScalar h;
168: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
169: h = 1.0/((mx-1)*(my-1)*(mz-1));
170: VecSet(b,h);
171: return(0);
172: }
176: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat* C)
177: {
179: PetscRandom rand;
180: Mat RHS;
181: PetscScalar *array,rval;
182: PetscInt i,k;
185: MatCreate(PETSC_COMM_WORLD,&RHS);
186: MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
187: MatSetType(RHS,MATSEQDENSE);
188:
189: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
190: PetscRandomSetFromOptions(rand);
191: MatGetArray(RHS,&array);
192: for (i=0; i<m; i++){
193: PetscRandomGetValue(rand,&rval);
194: array[i] = rval;
195: }
196: if (nrhs > 1){
197: for (k=1; k<nrhs; k++){
198: for (i=0; i<m; i++){
199: array[m*k+i] = array[i];
200: }
201: }
202: }
203: MatRestoreArray(RHS,&array);
204: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
205: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
206: *C = RHS;
207: PetscRandomDestroy(&rand);
208: return(0);
209: }
211:
214: PetscErrorCode ComputeMatrix(DM da,Mat B)
215: {
217: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
218: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
219: MatStencil row,col;
220: PetscRandom rand;
223: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
224: PetscRandomSetType(rand,PETSCRAND);
225: PetscRandomSetSeed(rand,1);
226: PetscRandomSetInterval(rand,-.001,.001);
227: PetscRandomSetFromOptions(rand);
229: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
230: /* For simplicity, this example only works on mx=my=mz */
231: if ( mx != my || mx != mz) SETERRQ3(PETSC_COMM_SELF,1,"This example only works with mx %d = my %d = mz %d\n",mx,my,mz);
233: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
234: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
236: PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
237: v_neighbor = v + dof*dof;
238: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
239: k3 = 0;
240: for (k1=0; k1<dof; k1++){
241: for (k2=0; k2<dof; k2++){
242: if (k1 == k2){
243: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
244: v_neighbor[k3] = -HxHydHz;
245: } else {
246: PetscRandomGetValue(rand,&r1);
247: PetscRandomGetValue(rand,&r2);
248: v[k3] = r1;
249: v_neighbor[k3] = r2;
250: }
251: k3++;
252: }
253: }
254: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
255:
256: for (k=zs; k<zs+zm; k++){
257: for (j=ys; j<ys+ym; j++){
258: for(i=xs; i<xs+xm; i++){
259: row.i = i; row.j = j; row.k = k;
260: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
261: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
262: } else { /* interior points */
263: /* center */
264: col.i = i; col.j = j; col.k = k;
265: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
266:
267: /* x neighbors */
268: col.i = i-1; col.j = j; col.k = k;
269: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
270: col.i = i+1; col.j = j; col.k = k;
271: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
272:
273: /* y neighbors */
274: col.i = i; col.j = j-1; col.k = k;
275: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
276: col.i = i; col.j = j+1; col.k = k;
277: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
278:
279: /* z neighbors */
280: col.i = i; col.j = j; col.k = k-1;
281: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
282: col.i = i; col.j = j; col.k = k+1;
283: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
284: }
285: }
286: }
287: }
288: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
289: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
290: PetscFree(v);
291: PetscRandomDestroy(&rand);
292: return(0);
293: }