Actual source code: ex129.c
petsc-3.3-p6 2013-02-11
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>
17: extern PetscErrorCode ComputeMatrix(DM,Mat);
18: extern PetscErrorCode ComputeRHS(DM,Vec);
19: extern PetscErrorCode ComputeRHSMatrix(PetscInt,PetscInt,Mat*);
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,tol=1.e-13;
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: DMCreateMatrix(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: if (norm > tol){
109: PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %G\n",norm);
110: }
112: /* MatSolveTranspose */
113: MatSolveTranspose(F,b,x);
114: MatMultTranspose(A,x,b1);
115: VecAXPY(b1,-1.0,b);
116: VecNorm(b1,NORM_2,&norm);
117: if (norm > tol){
118: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %G\n",norm);
119: }
121: /* MatSolveAdd */
122: MatSolveAdd(F,b,y,x);
123: MatMult(A,y,b1);
124: VecScale(b1,-1.0);
125: MatMultAdd(A,x,b1,b1);
126: VecAXPY(b1,-1.0,b);
127: VecNorm(b1,NORM_2,&norm);
128: if (norm > tol){
129: PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %G\n",norm);
130: }
132: /* MatSolveTransposeAdd */
133: MatSolveTransposeAdd(F,b,y,x);
134: MatMultTranspose(A,y,b1);
135: VecScale(b1,-1.0);
136: MatMultTransposeAdd(A,x,b1,b1);
137: VecAXPY(b1,-1.0,b);
138: VecNorm(b1,NORM_2,&norm);
139: if (norm > tol){
140: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %G\n",norm);
141: }
143: /* MatMatSolve */
144: MatMatSolve(F,RHS,X);
145: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
146: MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
147: MatNorm(C1,NORM_FROBENIUS,&norm);
148: if (norm > tol){
149: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %G\n",norm);
150: }
152: VecDestroy(&x);
153: VecDestroy(&b);
154: VecDestroy(&b1);
155: VecDestroy(&y);
156: MatDestroy(&A);
157: MatDestroy(&F);
158: MatDestroy(&RHS);
159: MatDestroy(&C1);
160: MatDestroy(&X);
161: ISDestroy(&perm);
162: ISDestroy(&iperm);
163: DMDestroy(&da);
164: PetscFinalize();
165: return 0;
166: }
170: PetscErrorCode ComputeRHS(DM da,Vec b)
171: {
173: PetscInt mx,my,mz;
174: PetscScalar h;
177: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
178: h = 1.0/((mx-1)*(my-1)*(mz-1));
179: VecSet(b,h);
180: return(0);
181: }
185: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat* C)
186: {
188: PetscRandom rand;
189: Mat RHS;
190: PetscScalar *array,rval;
191: PetscInt i,k;
194: MatCreate(PETSC_COMM_WORLD,&RHS);
195: MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
196: MatSetType(RHS,MATSEQDENSE);
197: MatSetUp(RHS);
199: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
200: PetscRandomSetFromOptions(rand);
201: MatGetArray(RHS,&array);
202: for (i=0; i<m; i++){
203: PetscRandomGetValue(rand,&rval);
204: array[i] = rval;
205: }
206: if (nrhs > 1){
207: for (k=1; k<nrhs; k++){
208: for (i=0; i<m; i++){
209: array[m*k+i] = array[i];
210: }
211: }
212: }
213: MatRestoreArray(RHS,&array);
214: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
215: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
216: *C = RHS;
217: PetscRandomDestroy(&rand);
218: return(0);
219: }
221:
224: PetscErrorCode ComputeMatrix(DM da,Mat B)
225: {
227: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
228: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
229: MatStencil row,col;
230: PetscRandom rand;
233: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
234: PetscRandomSetType(rand,PETSCRAND);
235: PetscRandomSetSeed(rand,1);
236: PetscRandomSetInterval(rand,-.001,.001);
237: PetscRandomSetFromOptions(rand);
239: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
240: /* For simplicity, this example only works on mx=my=mz */
241: 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);
243: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
244: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
246: PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
247: v_neighbor = v + dof*dof;
248: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
249: k3 = 0;
250: for (k1=0; k1<dof; k1++){
251: for (k2=0; k2<dof; k2++){
252: if (k1 == k2){
253: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
254: v_neighbor[k3] = -HxHydHz;
255: } else {
256: PetscRandomGetValue(rand,&r1);
257: PetscRandomGetValue(rand,&r2);
258: v[k3] = r1;
259: v_neighbor[k3] = r2;
260: }
261: k3++;
262: }
263: }
264: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
265:
266: for (k=zs; k<zs+zm; k++){
267: for (j=ys; j<ys+ym; j++){
268: for(i=xs; i<xs+xm; i++){
269: row.i = i; row.j = j; row.k = k;
270: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
271: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
272: } else { /* interior points */
273: /* center */
274: col.i = i; col.j = j; col.k = k;
275: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
276:
277: /* x neighbors */
278: col.i = i-1; col.j = j; col.k = k;
279: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
280: col.i = i+1; col.j = j; col.k = k;
281: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
282:
283: /* y neighbors */
284: col.i = i; col.j = j-1; col.k = k;
285: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
286: col.i = i; col.j = j+1; col.k = k;
287: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
288:
289: /* z neighbors */
290: col.i = i; col.j = j; col.k = k-1;
291: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
292: col.i = i; col.j = j; col.k = k+1;
293: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
294: }
295: }
296: }
297: }
298: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
299: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
300: PetscFree(v);
301: PetscRandomDestroy(&rand);
302: return(0);
303: }