Actual source code: ex129.c
petsc-3.8.4 2018-03-24
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: {
25: PetscMPIInt size;
26: Vec x,b,y,b1;
27: DM da;
28: Mat A,F,RHS,X,C1;
29: MatFactorInfo info;
30: IS perm,iperm;
31: PetscInt dof =1,M=8,m,n,nrhs;
32: PetscScalar one = 1.0;
33: PetscReal norm,tol = 1000*PETSC_MACHINE_EPSILON;
34: PetscBool InplaceLU=PETSC_FALSE;
36: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
37: MPI_Comm_size(PETSC_COMM_WORLD,&size);
38: if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only\n");
39: PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
40: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
42: DMDACreate(PETSC_COMM_WORLD,&da);
43: DMSetDimension(da,3);
44: DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
45: DMDASetStencilType(da,DMDA_STENCIL_STAR);
46: DMDASetSizes(da,M,M,M);
47: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
48: DMDASetDof(da,dof);
49: DMDASetStencilWidth(da,1);
50: DMDASetOwnershipRanges(da,NULL,NULL,NULL);
51: DMSetMatType(da,MATBAIJ);
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,&A);
61: ComputeMatrix(da,A);
62: MatGetSize(A,&m,&n);
63: nrhs = 2;
64: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
65: ComputeRHSMatrix(m,nrhs,&RHS);
66: MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X);
68: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
70: PetscOptionsGetBool(NULL,NULL,"-inplacelu",&InplaceLU,NULL);
71: MatFactorInfoInitialize(&info);
72: if (!InplaceLU) {
73: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
74: info.fill = 5.0;
75: MatLUFactorSymbolic(F,A,perm,iperm,&info);
76: MatLUFactorNumeric(F,A,&info);
77: } else { /* Test inplace factorization */
78: MatDuplicate(A,MAT_COPY_VALUES,&F);
79: MatLUFactor(F,perm,iperm,&info);
80: }
82: VecDuplicate(y,&b1);
84: /* MatSolve */
85: MatSolve(F,b,x);
86: MatMult(A,x,b1);
87: VecAXPY(b1,-1.0,b);
88: VecNorm(b1,NORM_2,&norm);
89: if (norm > tol) {
90: PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %g\n",(double)norm);
91: }
93: /* MatSolveTranspose */
94: MatSolveTranspose(F,b,x);
95: MatMultTranspose(A,x,b1);
96: VecAXPY(b1,-1.0,b);
97: VecNorm(b1,NORM_2,&norm);
98: if (norm > tol) {
99: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %g\n",(double)norm);
100: }
102: /* MatSolveAdd */
103: MatSolveAdd(F,b,y,x);
104: MatMult(A,y,b1);
105: VecScale(b1,-1.0);
106: MatMultAdd(A,x,b1,b1);
107: VecAXPY(b1,-1.0,b);
108: VecNorm(b1,NORM_2,&norm);
109: if (norm > tol) {
110: PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %g\n",(double)norm);
111: }
113: /* MatSolveTransposeAdd */
114: MatSolveTransposeAdd(F,b,y,x);
115: MatMultTranspose(A,y,b1);
116: VecScale(b1,-1.0);
117: MatMultTransposeAdd(A,x,b1,b1);
118: VecAXPY(b1,-1.0,b);
119: VecNorm(b1,NORM_2,&norm);
120: if (norm > tol) {
121: PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %g\n",(double)norm);
122: }
124: /* MatMatSolve */
125: MatMatSolve(F,RHS,X);
126: MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
127: MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
128: MatNorm(C1,NORM_FROBENIUS,&norm);
129: if (norm > tol) {
130: PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %g\n",(double)norm);
131: }
133: VecDestroy(&x);
134: VecDestroy(&b);
135: VecDestroy(&b1);
136: VecDestroy(&y);
137: MatDestroy(&A);
138: MatDestroy(&F);
139: MatDestroy(&RHS);
140: MatDestroy(&C1);
141: MatDestroy(&X);
142: ISDestroy(&perm);
143: ISDestroy(&iperm);
144: DMDestroy(&da);
145: PetscFinalize();
146: return ierr;
147: }
149: PetscErrorCode ComputeRHS(DM da,Vec b)
150: {
152: PetscInt mx,my,mz;
153: PetscScalar h;
156: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
157: h = 1.0/((mx-1)*(my-1)*(mz-1));
158: VecSet(b,h);
159: return(0);
160: }
162: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat *C)
163: {
165: PetscRandom rand;
166: Mat RHS;
167: PetscScalar *array,rval;
168: PetscInt i,k;
171: MatCreate(PETSC_COMM_WORLD,&RHS);
172: MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
173: MatSetType(RHS,MATSEQDENSE);
174: MatSetUp(RHS);
176: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
177: PetscRandomSetFromOptions(rand);
178: MatDenseGetArray(RHS,&array);
179: for (i=0; i<m; i++) {
180: PetscRandomGetValue(rand,&rval);
181: array[i] = rval;
182: }
183: if (nrhs > 1) {
184: for (k=1; k<nrhs; k++) {
185: for (i=0; i<m; i++) {
186: array[m*k+i] = array[i];
187: }
188: }
189: }
190: MatDenseRestoreArray(RHS,&array);
191: MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
192: MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
193: *C = RHS;
194: PetscRandomDestroy(&rand);
195: return(0);
196: }
199: PetscErrorCode ComputeMatrix(DM da,Mat B)
200: {
202: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
203: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
204: MatStencil row,col;
205: PetscRandom rand;
208: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
209: PetscRandomSetSeed(rand,1);
210: PetscRandomSetInterval(rand,-.001,.001);
211: PetscRandomSetFromOptions(rand);
213: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
214: /* For simplicity, this example only works on mx=my=mz */
215: 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);
217: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
218: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
220: PetscMalloc1(2*dof*dof+1,&v);
221: v_neighbor = v + dof*dof;
222: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
223: k3 = 0;
224: for (k1=0; k1<dof; k1++) {
225: for (k2=0; k2<dof; k2++) {
226: if (k1 == k2) {
227: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
228: v_neighbor[k3] = -HxHydHz;
229: } else {
230: PetscRandomGetValue(rand,&r1);
231: PetscRandomGetValue(rand,&r2);
233: v[k3] = r1;
234: v_neighbor[k3] = r2;
235: }
236: k3++;
237: }
238: }
239: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
241: for (k=zs; k<zs+zm; k++) {
242: for (j=ys; j<ys+ym; j++) {
243: for (i=xs; i<xs+xm; i++) {
244: row.i = i; row.j = j; row.k = k;
245: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boudary points */
246: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
247: } else { /* interior points */
248: /* center */
249: col.i = i; col.j = j; col.k = k;
250: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
252: /* x neighbors */
253: col.i = i-1; col.j = j; col.k = k;
254: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
255: col.i = i+1; col.j = j; col.k = k;
256: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
258: /* y neighbors */
259: col.i = i; col.j = j-1; col.k = k;
260: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
261: col.i = i; col.j = j+1; col.k = k;
262: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
264: /* z neighbors */
265: col.i = i; col.j = j; col.k = k-1;
266: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
267: col.i = i; col.j = j; col.k = k+1;
268: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
269: }
270: }
271: }
272: }
273: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
274: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
275: PetscFree(v);
276: PetscRandomDestroy(&rand);
277: return(0);
278: }