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 petscda.h
16: #include petscmg.h
24: int main(int argc,char **args)
25: {
26: PetscErrorCode ierr;
27: PetscMPIInt size;
28: Vec x,b,y,b1;
29: DA da;
30: Mat A,F,C,X,C1;
31: MatFactorInfo info;
32: IS perm,iperm;
33: PetscInt dof=1,M=-8,m,n,nrhs;
34: PetscScalar one = 1.0;
35: PetscReal norm;
37: PetscInitialize(&argc,&args,(char *)0,help);
38: MPI_Comm_size(PETSC_COMM_WORLD,&size);
39: if(size != 1) SETERRQ(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: DACreate(PETSC_COMM_WORLD,&da);
44: DASetDim(da,3);
45: DASetPeriodicity(da,DA_NONPERIODIC);
46: DASetStencilType(da,DA_STENCIL_STAR);
47: DASetSizes(da,M,M,M);
48: DASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
49: DASetDof(da,dof);
50: DASetStencilWidth(da,1);
51: DASetVertexDivision(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
52: DASetFromOptions(da);
54: DACreateGlobalVector(da,&x);
55: DACreateGlobalVector(da,&b);
56: VecDuplicate(b,&y);
57: ComputeRHS(da,b);
58: VecSet(y,one);
59: DAGetMatrix(da,MATBAIJ,&A);
60: ComputeMatrix(da,A);
61: MatGetSize(A,&m,&n);
62: nrhs = 2;
63: PetscOptionsGetInt(PETSC_NULL,"-nrhs",&nrhs,PETSC_NULL);
64: ComputeRHSMatrix(m,nrhs,&C);
65: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
66:
68: MatGetOrdering(A,MATORDERING_ND,&perm,&iperm);
69: MatGetFactor(A,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&F);
70: MatFactorInfoInitialize(&info);
71: info.fill = 5.0;
72: MatLUFactorSymbolic(F,A,perm,iperm,&info);
73: MatLUFactorNumeric(F,A,&info);
74: VecDuplicate(y,&b1);
75:
76: /* MatSolve */
77: MatSolve(F,b,x);
78: MatMult(A,x,b1);
79: VecAXPY(b1,-1.0,b);
80: VecNorm(b1,NORM_2,&norm);
81: PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %A\n",norm);
82:
83: /* MatSolveTranspose */
84: MatSolveTranspose(F,b,x);
85: MatMultTranspose(A,x,b1);
86: VecAXPY(b1,-1.0,b);
87: VecNorm(b1,NORM_2,&norm);
88: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %A\n",norm);
89:
90: /* MatSolveAdd */
91: MatSolveAdd(F,b,y,x);
92: MatMult(A,y,b1);
93: VecScale(b1,-1.0);
94: MatMultAdd(A,x,b1,b1);
95: VecAXPY(b1,-1.0,b);
96: VecNorm(b1,NORM_2,&norm);
97: PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %A\n",norm);
98:
99: /* MatSolveTransposeAdd */
100: MatSolveTransposeAdd(F,b,y,x);
101: MatMultTranspose(A,y,b1);
102: VecScale(b1,-1.0);
103: MatMultTransposeAdd(A,x,b1,b1);
104: VecAXPY(b1,-1.0,b);
105: VecNorm(b1,NORM_2,&norm);
106: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %A\n",norm);
107:
108: /* MatMatSolve */
109: MatMatSolve(F,C,X);
110: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
111: MatAXPY(C1,-1.0,C,SAME_NONZERO_PATTERN);
112: MatNorm(C1,NORM_FROBENIUS,&norm);
113: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %A\n",norm);
115: VecDestroy(x);
116: VecDestroy(b);
117: VecDestroy(b1);
118: VecDestroy(y);
119: MatDestroy(A);
120: MatDestroy(F);
121: MatDestroy(C);
122: MatDestroy(C1);
123: MatDestroy(X);
124: ISDestroy(perm);
125: ISDestroy(iperm);
126: DADestroy(da);
127: PetscFinalize();
128: return 0;
129: }
133: PetscErrorCode ComputeRHS(DA da,Vec b)
134: {
136: PetscInt mx,my,mz;
137: PetscScalar h;
140: DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
141: h = 1.0/((mx-1)*(my-1)*(mz-1));
142: VecSet(b,h);
143: return(0);
144: }
148: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat* C)
149: {
151: PetscRandom rand;
152: Mat RHS;
153: PetscScalar *array,rval;
154: PetscInt i,k;
157: MatCreate(PETSC_COMM_WORLD,&RHS);
158: MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
159: MatSetType(RHS,MATSEQDENSE);
160:
161: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
162: PetscRandomSetFromOptions(rand);
163: MatGetArray(RHS,&array);
164: for (i=0; i<m; i++){
165: PetscRandomGetValue(rand,&rval);
166: array[i] = rval;
167: }
168: if (nrhs > 1){
169: for (k=1; k<nrhs; k++){
170: for (i=0; i<m; i++){
171: array[m*k+i] = array[i];
172: }
173: }
174: }
175: MatRestoreArray(RHS,&array);
176: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
177: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
178: *C = RHS;
179: PetscRandomDestroy(rand);
180: return(0);
181: }
183:
186: PetscErrorCode ComputeMatrix(DA da,Mat B)
187: {
189: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
190: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
191: MatStencil row,col;
192: PetscRandom rand;
195: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
196: PetscRandomSetType(rand,PETSCRAND);
197: PetscRandomSetSeed(rand,1);
198: PetscRandomSetInterval(rand,-.001,.001);
199: PetscRandomSetFromOptions(rand);
201: DAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);
202: /* For simplicity, this example only works on mx=my=mz */
203: if ( mx != my || mx != mz) SETERRQ3(1,"This example only works with mx %d = my %d = mz %d\n",mx,my,mz);
205: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
206: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
208: PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
209: v_neighbor = v + dof*dof;
210: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
211: k3 = 0;
212: for (k1=0; k1<dof; k1++){
213: for (k2=0; k2<dof; k2++){
214: if (k1 == k2){
215: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
216: v_neighbor[k3] = -HxHydHz;
217: } else {
218: PetscRandomGetValue(rand,&r1);
219: PetscRandomGetValue(rand,&r2);
220: v[k3] = r1;
221: v_neighbor[k3] = r2;
222: }
223: k3++;
224: }
225: }
226: DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
227:
228: for (k=zs; k<zs+zm; k++){
229: for (j=ys; j<ys+ym; j++){
230: for(i=xs; i<xs+xm; i++){
231: row.i = i; row.j = j; row.k = k;
232: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
233: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
234: } else { /* interior points */
235: /* center */
236: col.i = i; col.j = j; col.k = k;
237: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
238:
239: /* x neighbors */
240: col.i = i-1; col.j = j; col.k = k;
241: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
242: col.i = i+1; col.j = j; col.k = k;
243: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
244:
245: /* y neighbors */
246: col.i = i; col.j = j-1; col.k = k;
247: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
248: col.i = i; col.j = j+1; col.k = k;
249: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
250:
251: /* z neighbors */
252: col.i = i; col.j = j; col.k = k-1;
253: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
254: col.i = i; col.j = j; col.k = k+1;
255: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
256: }
257: }
258: }
259: }
260: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
261: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
262: PetscFree(v);
263: PetscRandomDestroy(rand);
264: return(0);
265: }