Actual source code: ex40.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
10: #include petscda.h
12: PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig);
16: int main(int Argc,char **Args)
17: {
18: PetscTruth flg;
19: PetscInt n = 6,i;
20: PetscScalar rho = 1.0;
21: PetscReal h;
22: PetscReal beta = 1.0;
23: ADDA adda;
24: PetscInt nodes[2];
25: PetscTruth periodic[2];
26: PetscInt refine[2];
27: PetscRandom rctx;
28: PetscMPIInt comm_size;
29: Mat H;
30: PetscInt *lcs, *lce;
31: PetscInt x, y;
32: PetscReal r1, r2;
33: PetscScalar uxy1, uxy2;
34: ADDAIdx sxy, sxy_m;
35: PetscScalar val, valconj;
36: Mat HtH;
37: Vec b, Htb;
38: Vec xvec;
39: KSP kspmg;
40: PC pcmg;
41: PetscErrorCode ierr;
43: PetscInitialize(&Argc,&Args,(char *)0,help);
44: PetscOptionsGetInt(PETSC_NULL,"-size",&n,&flg);
45: PetscOptionsGetReal(PETSC_NULL,"-beta",&beta,&flg);
46: PetscOptionsGetScalar(PETSC_NULL,"-rho",&rho,&flg);
48: /* Set the fudge parameters, we scale the whole thing by 1/(2*h) later */
49: h = 1.;
50: rho *= 1./(2.*h);
51:
52: /* Geometry info */
53: for(i=0; i<2; i++) {
54: nodes[i] = n;
55: periodic[i] = PETSC_TRUE;
56: refine[i] = 3;
57: }
58: ADDACreate(PETSC_COMM_WORLD, 2, nodes, PETSC_NULL, 2 /* this is the # of dof's */,
59: periodic, &adda);
60: PetscObjectReference((PetscObject) adda);
61: ADDASetRefinement(adda, refine, 2);
62:
63: /* Random numbers */
64: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
65: PetscRandomSetFromOptions(rctx);
67: /* Single or multi processor ? */
68: MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);
70: /* construct matrix */
71: if( comm_size == 1 ) {
72: ADDAGetMatrix(adda, MATSEQAIJ, &H);
73: } else {
74: ADDAGetMatrix(adda, MATMPIAIJ, &H);
75: }
77: /* get local corners for this processor */
78: ADDAGetCorners(adda, &lcs, &lce);
80: /* Allocate space for the indices that we use to construct the matrix */
81: PetscMalloc(2*sizeof(PetscInt), &(sxy.x));
82: PetscMalloc(2*sizeof(PetscInt), &(sxy_m.x));
84: /* Assemble the matrix */
85: for( x=lcs[0]; x<lce[0]; x++ ) {
86: for( y=lcs[1]; y<lce[1]; y++ ) {
87: /* each lattice point sets only the *forward* pointing parameters (right, down),
88: i.e. Nabla_1^+ and Nabla_2^+.
89: In this way we can use only local random number creation. That means
90: we also have to set the corresponding backward pointing entries. */
91: /* Compute some normally distributed random numbers via Box-Muller */
92: PetscRandomGetValueReal(rctx, &r1);
93: r1 = 1.-r1; /* to change from [0,1) to (0,1], which we need for the log */
94: PetscRandomGetValueReal(rctx, &r2);
95: PetscReal R = sqrt(-2.*log(r1));
96: PetscReal c = cos(2.*PETSC_PI*r2);
97: PetscReal s = sin(2.*PETSC_PI*r2);
99: /* use those to set the field */
100: uxy1 = PetscExpScalar( ((PetscScalar) (R*c/beta))*PETSC_i);
101: uxy2 = PetscExpScalar( ((PetscScalar) (R*s/beta))*PETSC_i);
102:
103: sxy.x[0] = x; sxy.x[1] = y; /* the point where we are */
105: /* center action */
106: sxy.d = 0; /* spin 0, 0 */
107: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy, &rho, ADD_VALUES);
108: sxy.d = 1; /* spin 1, 1 */
109: val = -rho;
110: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy, &val, ADD_VALUES);
111:
112: sxy_m.x[0] = x+1; sxy_m.x[1] = y; /* right action */
113: sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
114: val = -uxy1; valconj = PetscConj(val);
115: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
116: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
117: sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
118: val = -uxy1; valconj = PetscConj(val);
119: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
120: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
121: sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
122: val = uxy1; valconj = PetscConj(val);
123: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
124: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
125: sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
126: val = uxy1; valconj = PetscConj(val);
127: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
128: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
130: sxy_m.x[0] = x; sxy_m.x[1] = y+1; /* down action */
131: sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
132: val = -uxy2; valconj = PetscConj(val);
133: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
134: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
135: sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
136: val = -PETSC_i*uxy2; valconj = PetscConj(val);
137: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
138: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
139: sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
140: val = -PETSC_i*uxy2; valconj = PetscConj(val);
141: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
142: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
143: sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
144: val = PetscConj(uxy2); valconj = PetscConj(val);
145: ADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
146: ADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
147: }
148: }
149:
150: PetscFree(sxy.x);
151: PetscFree(sxy_m.x);
153: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
154: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
156: /* scale H */
157: MatScale(H, 1./(2.*h));
159: /* construct normal equations */
160: MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);
161: /* reference matrix */
162: PetscObjectReference((PetscObject) HtH);
164: PetscScalar mineval;
165: computeMinEigVal(HtH, 1000, &mineval);
166: PetscPrintf(PETSC_COMM_WORLD, "Minimum eigenvalue of H^{dag} H is %f\n", PetscAbsScalar(mineval));
168: /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
169: /* Mat perm; */
170: /* ADDAGetMatrix(adda, MATSEQAIJ, &perm); */
171: /* PetscInt row, col; */
172: /* PetscScalar one = 1.0; */
173: /* for(PetscInt i=0; i<n; i++) { */
174: /* for(PetscInt j=0; j<n; j++) { */
175: /* row = (i*n+j)*2; col = i*n+j; */
176: /* MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
177: /* row = (i*n+j)*2+1; col = i*n+j + n*n; */
178: /* MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
179: /* } */
180: /* } */
181: /* MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY); */
182: /* MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY); */
184: /* Mat Hperm; */
185: /* MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm); */
186: /* PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n"); */
187: /* MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
189: /* Mat HtHperm; */
190: /* MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm); */
191: /* PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n"); */
192: /* MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
194: /* right hand side */
195: ADDACreateGlobalVector(adda, &b);
196: VecSet(b,0.0);
197: PetscInt ix[1] = {0};
198: PetscScalar vals[1] = {1.0};
199: VecSetValues(b, 1, ix, vals, INSERT_VALUES);
200: VecAssemblyBegin(b);
201: VecAssemblyEnd(b);
202: /* VecSetRandom(b, rctx); */
203: VecDuplicate(b, &Htb);
204: MatMultTranspose(H, b, Htb);
205: /* reference it */
206: PetscObjectReference((PetscObject) Htb);
208: /* construct solver */
209: KSPCreate(PETSC_COMM_WORLD,&kspmg);
210: KSPSetType(kspmg, KSPCG);
212: KSPGetPC(kspmg,&pcmg);
213: PCSetType(pcmg,PCASA);
215: /* maybe user wants to override some of the choices */
216: KSPSetFromOptions(kspmg);
218: KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);
220: PCASASetDM(pcmg, (DM) adda);
222: PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);
224: VecDuplicate(b, &xvec);
225: VecSet(xvec, 0.0);
227: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228: Solve the linear system
229: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231: KSPSolve(kspmg, Htb, xvec);
233: /* VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
235: KSPDestroy(kspmg);
237: VecDestroy(xvec);
238: VecDestroy(b);
239: VecDestroy(Htb);
240: MatDestroy(H);
241: MatDestroy(HtH);
243: ADDADestroy(adda);
244: PetscRandomDestroy(rctx);
245: PetscFinalize();
246: return 0;
247: }
249: /* --------------------------------------------------------------------- */
252: PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig) {
253: PetscErrorCode ierr;
254: PetscRandom rctx; /* random number generator context */
255: Vec x0, x, x_1, tmp;
256: PetscScalar lambda_its, lambda_its_1;
257: PetscReal norm;
258: Mat G;
259: PetscInt i;
260:
262: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
263: PetscRandomSetFromOptions(rctx);
265: /* compute G = I-1/norm(A)*A */
266: MatNorm(A, NORM_1, &norm);
267: MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &G);
268: MatShift(G, -norm);
269: MatScale(G, -1./norm);
271: MatGetVecs(G, &x_1, &x);
272: VecSetRandom(x, rctx);
273: VecDuplicate(x, &x0);
274: VecCopy(x, x0);
276: MatMult(G, x, x_1);
277: for(i=0; i<its; i++) {
278: tmp = x; x = x_1; x_1 = tmp;
279: MatMult(G, x, x_1);
280: }
281: VecDot(x0, x, &lambda_its);
282: VecDot(x0, x_1, &lambda_its_1);
284: *eig = norm*(1.-lambda_its_1/lambda_its);
286: VecDestroy(x0);
287: VecDestroy(x);
288: VecDestroy(x_1);
289: PetscRandomDestroy(rctx);
290: MatDestroy(G);
292: return(0);
293: }