Actual source code: ex129.c
petsc-3.6.4 2016-04-12
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*);
24: int main(int argc,char **args)
25: {
27: PetscMPIInt size;
28: Vec x,b,y,b1;
29: DM da;
30: Mat A,F,RHS,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,tol=1.e-13;
36: PetscBool InplaceLU=PETSC_FALSE;
38: PetscInitialize(&argc,&args,(char*)0,help);
39: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only\n");
41: PetscOptionsGetInt(NULL,"-dof",&dof,NULL);
42: PetscOptionsGetInt(NULL,"-M",&M,NULL);
44: DMDACreate(PETSC_COMM_WORLD,&da);
45: DMSetDimension(da,3);
46: DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
47: DMDASetStencilType(da,DMDA_STENCIL_STAR);
48: DMDASetSizes(da,M,M,M);
49: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
50: DMDASetDof(da,dof);
51: DMDASetStencilWidth(da,1);
52: DMDASetOwnershipRanges(da,NULL,NULL,NULL);
53: DMSetFromOptions(da);
54: DMSetUp(da);
56: DMCreateGlobalVector(da,&x);
57: DMCreateGlobalVector(da,&b);
58: VecDuplicate(b,&y);
59: ComputeRHS(da,b);
60: VecSet(y,one);
61: DMSetMatType(da,MATBAIJ);
62: DMCreateMatrix(da,&A);
63: ComputeMatrix(da,A);
64: MatGetSize(A,&m,&n);
65: nrhs = 2;
66: PetscOptionsGetInt(NULL,"-nrhs",&nrhs,NULL);
67: ComputeRHSMatrix(m,nrhs,&RHS);
68: MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X);
70: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
73: PetscOptionsGetBool(NULL,"-inplacelu",&InplaceLU,NULL);
74: MatFactorInfoInitialize(&info);
75: if (!InplaceLU) {
76: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
77: info.fill = 5.0;
78: MatLUFactorSymbolic(F,A,perm,iperm,&info);
79: MatLUFactorNumeric(F,A,&info);
80: } else { /* Test inplace factorization */
81: MatDuplicate(A,MAT_COPY_VALUES,&F);
82: /* or create F without DMDA
83: MatType type;
84: PetscInt i,ncols;
85: const PetscInt *cols;
86: const PetscScalar *vals;
87: MatGetSize(A,&m,&n);
88: MatGetType(A,&type);
89: MatCreate(PetscObjectComm((PetscObject)A),&F);
90: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,m,n);
91: MatSetType(F,type);
92: MatSetFromOptions(F);
93: for (i=0; i<m; i++) {
94: MatGetRow(A,i,&ncols,&cols,&vals);
95: MatSetValues(F,1,&i,ncols,cols,vals,INSERT_VALUES);
96: }
97: MatAssemblyBegin(F,MAT_FINAL_ASSEMBLY);
98: MatAssemblyEnd(F,MAT_FINAL_ASSEMBLY);
99: */
100: MatLUFactor(F,perm,iperm,&info);
101: }
103: VecDuplicate(y,&b1);
105: /* MatSolve */
106: MatSolve(F,b,x);
107: MatMult(A,x,b1);
108: VecAXPY(b1,-1.0,b);
109: VecNorm(b1,NORM_2,&norm);
110: if (norm > tol) {
111: PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %g\n",(double)norm);
112: }
114: /* MatSolveTranspose */
115: MatSolveTranspose(F,b,x);
116: MatMultTranspose(A,x,b1);
117: VecAXPY(b1,-1.0,b);
118: VecNorm(b1,NORM_2,&norm);
119: if (norm > tol) {
120: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %g\n",(double)norm);
121: }
123: /* MatSolveAdd */
124: MatSolveAdd(F,b,y,x);
125: MatMult(A,y,b1);
126: VecScale(b1,-1.0);
127: MatMultAdd(A,x,b1,b1);
128: VecAXPY(b1,-1.0,b);
129: VecNorm(b1,NORM_2,&norm);
130: if (norm > tol) {
131: PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %g\n",(double)norm);
132: }
134: /* MatSolveTransposeAdd */
135: MatSolveTransposeAdd(F,b,y,x);
136: MatMultTranspose(A,y,b1);
137: VecScale(b1,-1.0);
138: MatMultTransposeAdd(A,x,b1,b1);
139: VecAXPY(b1,-1.0,b);
140: VecNorm(b1,NORM_2,&norm);
141: if (norm > tol) {
142: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %g\n",(double)norm);
143: }
145: /* MatMatSolve */
146: MatMatSolve(F,RHS,X);
147: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
148: MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
149: MatNorm(C1,NORM_FROBENIUS,&norm);
150: if (norm > tol) {
151: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %g\n",(double)norm);
152: }
154: VecDestroy(&x);
155: VecDestroy(&b);
156: VecDestroy(&b1);
157: VecDestroy(&y);
158: MatDestroy(&A);
159: MatDestroy(&F);
160: MatDestroy(&RHS);
161: MatDestroy(&C1);
162: MatDestroy(&X);
163: ISDestroy(&perm);
164: ISDestroy(&iperm);
165: DMDestroy(&da);
166: PetscFinalize();
167: return 0;
168: }
172: PetscErrorCode ComputeRHS(DM da,Vec b)
173: {
175: PetscInt mx,my,mz;
176: PetscScalar h;
179: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
180: h = 1.0/((mx-1)*(my-1)*(mz-1));
181: VecSet(b,h);
182: return(0);
183: }
187: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat *C)
188: {
190: PetscRandom rand;
191: Mat RHS;
192: PetscScalar *array,rval;
193: PetscInt i,k;
196: MatCreate(PETSC_COMM_WORLD,&RHS);
197: MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
198: MatSetType(RHS,MATSEQDENSE);
199: MatSetUp(RHS);
201: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
202: PetscRandomSetFromOptions(rand);
203: MatDenseGetArray(RHS,&array);
204: for (i=0; i<m; i++) {
205: PetscRandomGetValue(rand,&rval);
206: array[i] = rval;
207: }
208: if (nrhs > 1) {
209: for (k=1; k<nrhs; k++) {
210: for (i=0; i<m; i++) {
211: array[m*k+i] = array[i];
212: }
213: }
214: }
215: MatDenseRestoreArray(RHS,&array);
216: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
217: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
218: *C = RHS;
219: PetscRandomDestroy(&rand);
220: return(0);
221: }
226: PetscErrorCode ComputeMatrix(DM da,Mat B)
227: {
229: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
230: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
231: MatStencil row,col;
232: PetscRandom rand;
235: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
236: PetscRandomSetType(rand,PETSCRAND);
237: PetscRandomSetSeed(rand,1);
238: PetscRandomSetInterval(rand,-.001,.001);
239: PetscRandomSetFromOptions(rand);
241: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
242: /* For simplicity, this example only works on mx=my=mz */
243: 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);
245: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
246: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
248: PetscMalloc1(2*dof*dof+1,&v);
249: v_neighbor = v + dof*dof;
250: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
251: k3 = 0;
252: for (k1=0; k1<dof; k1++) {
253: for (k2=0; k2<dof; k2++) {
254: if (k1 == k2) {
255: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
256: v_neighbor[k3] = -HxHydHz;
257: } else {
258: PetscRandomGetValue(rand,&r1);
259: PetscRandomGetValue(rand,&r2);
261: v[k3] = r1;
262: v_neighbor[k3] = r2;
263: }
264: k3++;
265: }
266: }
267: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
269: for (k=zs; k<zs+zm; k++) {
270: for (j=ys; j<ys+ym; j++) {
271: for (i=xs; i<xs+xm; i++) {
272: row.i = i; row.j = j; row.k = k;
273: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boudary points */
274: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
275: } else { /* interior points */
276: /* center */
277: col.i = i; col.j = j; col.k = k;
278: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
280: /* x neighbors */
281: col.i = i-1; col.j = j; col.k = k;
282: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
283: col.i = i+1; col.j = j; col.k = k;
284: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
286: /* y neighbors */
287: col.i = i; col.j = j-1; col.k = k;
288: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
289: col.i = i; col.j = j+1; col.k = k;
290: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
292: /* z neighbors */
293: col.i = i; col.j = j; col.k = k-1;
294: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
295: col.i = i; col.j = j; col.k = k+1;
296: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
297: }
298: }
299: }
300: }
301: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
302: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
303: PetscFree(v);
304: PetscRandomDestroy(&rand);
305: return(0);
306: }