Actual source code: ex39.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Lattice Gauge 2D model.\n"
3: "Parameters:\n"
4: "-size n to use a grid size of n, i.e n space and n time steps\n"
5: "-beta b controls the randomness of the gauge field\n"
6: "-rho r the quark mass (?)";
8: #include <petscksp.h>
9: #include <petscpcasa.h>
10: #include <petscdmda.h>
12: PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig);
16: int main(int Argc,char **Args)
17: {
18: PetscBool flg;
19: PetscInt n = -6;
20: PetscScalar rho = 1.0;
21: PetscReal h;
22: PetscReal beta = 1.0;
23: DM da;
24: PetscRandom rctx;
25: PetscMPIInt comm_size;
26: Mat H,HtH;
27: PetscInt x, y, xs, ys, xm, ym;
28: PetscReal r1, r2;
29: PetscScalar uxy1, uxy2;
30: MatStencil sxy, sxy_m;
31: PetscScalar val, valconj;
32: Vec b, Htb,xvec;
33: KSP kspmg;
34: PC pcmg;
36: PetscInt ix[1] = {0};
37: PetscScalar vals[1] = {1.0};
39: PetscInitialize(&Argc,&Args,(char*)0,help);
40: PetscOptionsGetInt(NULL,"-size",&n,&flg);
41: PetscOptionsGetReal(NULL,"-beta",&beta,&flg);
42: PetscOptionsGetScalar(NULL,"-rho",&rho,&flg);
44: /* Set the fudge parameters, we scale the whole thing by 1/(2*h) later */
45: h = 1.;
46: rho *= 1./(2.*h);
48: /* Geometry info */
49: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, n, n,
50: PETSC_DECIDE, PETSC_DECIDE, 2 /* this is the # of dof's */,
51: 1, NULL, NULL, &da);
53: /* Random numbers */
54: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
55: PetscRandomSetFromOptions(rctx);
57: /* Single or multi processor ? */
58: MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);
60: /* construct matrix */
61: if (comm_size == 1) {
62: DMCreateMatrix(da, MATSEQAIJ, &H);
63: } else {
64: DMCreateMatrix(da, MATMPIAIJ, &H);
65: }
67: /* get local corners for this processor */
68: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
70: /* Assemble the matrix */
71: for (x=xs; x<xs+xm; x++) {
72: for (y=ys; y<ys+ym; y++) {
73: /* each lattice point sets only the *forward* pointing parameters (right, down),
74: i.e. Nabla_1^+ and Nabla_2^+.
75: In this way we can use only local random number creation. That means
76: we also have to set the corresponding backward pointing entries. */
77: /* Compute some normally distributed random numbers via Box-Muller */
78: PetscRandomGetValueReal(rctx, &r1);
79: r1 = 1.-r1; /* to change from [0,1) to (0,1], which we need for the log */
80: PetscRandomGetValueReal(rctx, &r2);
81: PetscReal R = PetscSqrtReal(-2.*PetscLogReal(r1));
82: PetscReal c = PetscCosReal(2.*PETSC_PI*r2);
83: PetscReal s = PetscSinReal(2.*PETSC_PI*r2);
85: /* use those to set the field */
86: uxy1 = PetscExpScalar(((PetscScalar) (R*c/beta))*PETSC_i);
87: uxy2 = PetscExpScalar(((PetscScalar) (R*s/beta))*PETSC_i);
89: sxy.i = x; sxy.j = y; /* the point where we are */
91: /* center action */
92: sxy.c = 0; /* spin 0, 0 */
93: MatSetValuesStencil(H, 1, &sxy, 1, &sxy, &rho, ADD_VALUES);
94: sxy.c = 1; /* spin 1, 1 */
95: val = -rho;
96: MatSetValuesStencil(H, 1, &sxy, 1, &sxy, &val, ADD_VALUES);
98: sxy_m.i = x+1; sxy_m.j = y; /* right action */
99: sxy.c = 0; sxy_m.c = 0; /* spin 0, 0 */
100: val = -uxy1; valconj = PetscConj(val);
101: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
102: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
103: sxy.c = 0; sxy_m.c = 1; /* spin 0, 1 */
104: val = -uxy1; valconj = PetscConj(val);
105: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
106: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
107: sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
108: val = uxy1; valconj = PetscConj(val);
109: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
110: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
111: sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
112: val = uxy1; valconj = PetscConj(val);
113: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
114: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
116: sxy_m.i = x; sxy_m.j = y+1; /* down action */
117: sxy.c = 0; sxy_m.c = 0; /* spin 0, 0 */
118: val = -uxy2; valconj = PetscConj(val);
119: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
120: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
121: sxy.c = 0; sxy_m.c = 1; /* spin 0, 1 */
122: val = -PETSC_i*uxy2; valconj = PetscConj(val);
123: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
124: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
125: sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
126: val = -PETSC_i*uxy2; valconj = PetscConj(val);
127: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
128: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
129: sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
130: val = PetscConj(uxy2); valconj = PetscConj(val);
131: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
132: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
133: }
134: }
136: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
137: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
139: /* scale H */
140: MatScale(H, 1./(2.*h));
142: /* it looks like H is Hermetian */
143: /* construct normal equations */
144: MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);
146: /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
147: /* Mat perm; */
148: /* DMCreateMatrix(da, MATSEQAIJ, &perm); */
149: /* PetscInt row, col; */
150: /* PetscScalar one = 1.0; */
151: /* for (PetscInt i=0; i<n; i++) { */
152: /* for (PetscInt j=0; j<n; j++) { */
153: /* row = (i*n+j)*2; col = i*n+j; */
154: /* MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
155: /* row = (i*n+j)*2+1; col = i*n+j + n*n; */
156: /* MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
157: /* } */
158: /* } */
159: /* MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY); */
160: /* MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY); */
162: /* Mat Hperm; */
163: /* MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm); */
164: /* PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n"); */
165: /* MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
167: /* Mat HtHperm; */
168: /* MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm); */
169: /* PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n"); */
170: /* MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
172: /* right hand side */
173: DMCreateGlobalVector(da, &b);
174: VecSet(b,0.0);
175: VecSetValues(b, 1, ix, vals, INSERT_VALUES);
176: VecAssemblyBegin(b);
177: VecAssemblyEnd(b);
178: /* VecSetRandom(b, rctx); */
179: VecDuplicate(b, &Htb);
180: MatMultTranspose(H, b, Htb);
182: /* construct solver */
183: KSPCreate(PETSC_COMM_WORLD,&kspmg);
184: KSPSetType(kspmg, KSPCG);
186: KSPGetPC(kspmg,&pcmg);
187: PCSetType(pcmg,PCASA);
189: /* maybe user wants to override some of the choices */
190: KSPSetFromOptions(kspmg);
192: KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);
194: DMDASetRefinementFactor(da, 3, 3, 3);
195: PCSetDM(pcmg,da);
197: PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);
199: VecDuplicate(b, &xvec);
200: VecSet(xvec, 0.0);
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Solve the linear system
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: KSPSolve(kspmg, Htb, xvec);
208: /* VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
210: KSPDestroy(&kspmg);
211: VecDestroy(&xvec);
213: /* seems to be destroyed by KSPDestroy */
214: VecDestroy(&b);
215: VecDestroy(&Htb);
216: MatDestroy(&HtH);
217: MatDestroy(&H);
219: DMDestroy(&da);
220: PetscRandomDestroy(&rctx);
221: PetscFinalize();
222: return 0;
223: }
225: /* --------------------------------------------------------------------- */
228: PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig)
229: {
231: PetscRandom rctx; /* random number generator context */
232: Vec x0, x, x_1, tmp;
233: PetscScalar lambda_its, lambda_its_1;
234: PetscInt i;
237: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
238: PetscRandomSetFromOptions(rctx);
239: MatGetVecs(A, &x_1, &x);
240: VecSetRandom(x, rctx);
241: VecDuplicate(x, &x0);
242: VecCopy(x, x0);
244: MatMult(A, x, x_1);
245: for (i=0; i<its; i++) {
246: tmp = x; x = x_1; x_1 = tmp;
247: MatMult(A, x, x_1);
248: }
249: VecDot(x0, x, &lambda_its);
250: VecDot(x0, x_1, &lambda_its_1);
252: *eig = lambda_its_1/lambda_its;
254: VecDestroy(&x0);
255: VecDestroy(&x);
256: VecDestroy(&x_1);
257: return(0);
258: }