Actual source code: ex39.c
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 petscasa.h
11: PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig);
15: int main(int Argc,char **Args)
16: {
17: PetscTruth flg;
18: PetscInt n = -6;
19: PetscScalar rho = 1.0;
20: PetscReal h;
21: PetscReal beta = 1.0;
22: DA da;
23: PetscRandom rctx;
24: PetscMPIInt comm_size;
25: Mat H,HtH;
26: PetscInt x, y, xs, ys, xm, ym;
27: PetscReal r1, r2;
28: PetscScalar uxy1, uxy2;
29: MatStencil sxy, sxy_m;
30: PetscScalar val, valconj;
31: Vec b, Htb,xvec;
32: KSP kspmg;
33: PC pcmg;
34: PetscErrorCode ierr;
35: PetscInt ix[1] = {0};
36: PetscScalar vals[1] = {1.0};
38: PetscInitialize(&Argc,&Args,(char *)0,help);
39: PetscOptionsGetInt(PETSC_NULL,"-size",&n,&flg);
40: PetscOptionsGetReal(PETSC_NULL,"-beta",&beta,&flg);
41: PetscOptionsGetScalar(PETSC_NULL,"-rho",&rho,&flg);
43: /* Set the fudge parameters, we scale the whole thing by 1/(2*h) later */
44: h = 1.;
45: rho *= 1./(2.*h);
46:
47: /* Geometry info */
48: DACreate2d(PETSC_COMM_WORLD, DA_XYPERIODIC, DA_STENCIL_STAR, n, n,
49: PETSC_DECIDE, PETSC_DECIDE, 2 /* this is the # of dof's */,
50: 1, PETSC_NULL, PETSC_NULL, &da);
51:
52: /* Random numbers */
53: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
54: PetscRandomSetFromOptions(rctx);
56: /* Single or multi processor ? */
57: MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);
59: /* construct matrix */
60: if( comm_size == 1 ) {
61: DAGetMatrix(da, MATSEQAIJ, &H);
62: } else {
63: DAGetMatrix(da, MATMPIAIJ, &H);
64: }
66: /* get local corners for this processor */
67: DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
69: /* Assemble the matrix */
70: for( x=xs; x<xs+xm; x++ ) {
71: for( y=ys; y<ys+ym; y++ ) {
72: /* each lattice point sets only the *forward* pointing parameters (right, down),
73: i.e. Nabla_1^+ and Nabla_2^+.
74: In this way we can use only local random number creation. That means
75: we also have to set the corresponding backward pointing entries. */
76: /* Compute some normally distributed random numbers via Box-Muller */
77: PetscRandomGetValueReal(rctx, &r1);
78: r1 = 1.-r1; /* to change from [0,1) to (0,1], which we need for the log */
79: PetscRandomGetValueReal(rctx, &r2);
80: PetscReal R = sqrt(-2.*log(r1));
81: PetscReal c = cos(2.*PETSC_PI*r2);
82: PetscReal s = sin(2.*PETSC_PI*r2);
84: /* use those to set the field */
85: uxy1 = PetscExpScalar( ((PetscScalar) (R*c/beta))*PETSC_i);
86: uxy2 = PetscExpScalar( ((PetscScalar) (R*s/beta))*PETSC_i);
87:
88: sxy.i = x; sxy.j = y; /* the point where we are */
90: /* center action */
91: sxy.c = 0; /* spin 0, 0 */
92: MatSetValuesStencil(H, 1, &sxy, 1, &sxy, &rho, ADD_VALUES);
93: sxy.c = 1; /* spin 1, 1 */
94: val = -rho;
95: MatSetValuesStencil(H, 1, &sxy, 1, &sxy, &val, ADD_VALUES);
96:
97: sxy_m.i = x+1; sxy_m.j = y; /* right action */
98: sxy.c = 0; sxy_m.c = 0; /* spin 0, 0 */
99: val = -uxy1; valconj = PetscConj(val);
100: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
101: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
102: sxy.c = 0; sxy_m.c = 1; /* spin 0, 1 */
103: val = -uxy1; valconj = PetscConj(val);
104: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
105: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
106: sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
107: val = uxy1; valconj = PetscConj(val);
108: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
109: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
110: sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
111: val = uxy1; valconj = PetscConj(val);
112: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
113: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
115: sxy_m.i = x; sxy_m.j = y+1; /* down action */
116: sxy.c = 0; sxy_m.c = 0; /* spin 0, 0 */
117: val = -uxy2; valconj = PetscConj(val);
118: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
119: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
120: sxy.c = 0; sxy_m.c = 1; /* spin 0, 1 */
121: val = -PETSC_i*uxy2; valconj = PetscConj(val);
122: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
123: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
124: sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
125: val = -PETSC_i*uxy2; valconj = PetscConj(val);
126: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
127: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
128: sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
129: val = PetscConj(uxy2); valconj = PetscConj(val);
130: MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);
131: MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);
132: }
133: }
134:
135: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
136: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
138: /* scale H */
139: MatScale(H, 1./(2.*h));
141: /* it looks like H is Hermetian */
142: /* construct normal equations */
143: MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);
145: /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
146: /* Mat perm; */
147: /* DAGetMatrix(da, MATSEQAIJ, &perm); */
148: /* PetscInt row, col; */
149: /* PetscScalar one = 1.0; */
150: /* for(PetscInt i=0; i<n; i++) { */
151: /* for(PetscInt j=0; j<n; j++) { */
152: /* row = (i*n+j)*2; col = i*n+j; */
153: /* MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
154: /* row = (i*n+j)*2+1; col = i*n+j + n*n; */
155: /* MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
156: /* } */
157: /* } */
158: /* MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY); */
159: /* MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY); */
161: /* Mat Hperm; */
162: /* MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm); */
163: /* PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n"); */
164: /* MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
166: /* Mat HtHperm; */
167: /* MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm); */
168: /* PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n"); */
169: /* MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
171: /* right hand side */
172: DACreateGlobalVector(da, &b);
173: VecSet(b,0.0);
174: VecSetValues(b, 1, ix, vals, INSERT_VALUES);
175: VecAssemblyBegin(b);
176: VecAssemblyEnd(b);
177: /* VecSetRandom(b, rctx); */
178: VecDuplicate(b, &Htb);
179: MatMultTranspose(H, b, Htb);
181: /* construct solver */
182: KSPCreate(PETSC_COMM_WORLD,&kspmg);
183: KSPSetType(kspmg, KSPCG);
185: KSPGetPC(kspmg,&pcmg);
186: PCSetType(pcmg,PCASA);
188: /* maybe user wants to override some of the choices */
189: KSPSetFromOptions(kspmg);
191: KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);
193: DASetRefinementFactor(da, 3, 3, 3);
194: PCASASetDM(pcmg, (DM) da);
196: PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);
198: VecDuplicate(b, &xvec);
199: VecSet(xvec, 0.0);
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Solve the linear system
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205: KSPSolve(kspmg, Htb, xvec);
207: /* VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
209: KSPDestroy(kspmg);
210: VecDestroy(xvec);
212: /* seems to be destroyed by KSPDestroy */
213: VecDestroy(b);
214: VecDestroy(Htb);
215: MatDestroy(HtH);
216: MatDestroy(H);
218: DADestroy(da);
219: PetscRandomDestroy(rctx);
220: PetscFinalize();
221: return 0;
222: }
224: /* --------------------------------------------------------------------- */
227: PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig) {
228: PetscErrorCode ierr;
229: PetscRandom rctx; /* random number generator context */
230: Vec x0, x, x_1, tmp;
231: PetscScalar lambda_its, lambda_its_1;
232: PetscInt i;
233:
235: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
236: PetscRandomSetFromOptions(rctx);
237: MatGetVecs(A, &x_1, &x);
238: VecSetRandom(x, rctx);
239: VecDuplicate(x, &x0);
240: VecCopy(x, x0);
242: MatMult(A, x, x_1);
243: for(i=0; i<its; i++) {
244: tmp = x; x = x_1; x_1 = tmp;
245: MatMult(A, x, x_1);
246: }
247: VecDot(x0, x, &lambda_its);
248: VecDot(x0, x_1, &lambda_its_1);
250: *eig = lambda_its_1/lambda_its;
252: VecDestroy(x0);
253: VecDestroy(x);
254: VecDestroy(x_1);
256: return(0);
257: }