Actual source code: ex40.c
petsc-3.3-p7 2013-05-11
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 <petscdm.h>
11: #include <petscdmadda.h>
13: PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig);
17: int main(int Argc,char **Args)
18: {
19: PetscBool flg;
20: PetscInt n = 6,i;
21: PetscScalar rho = 1.0;
22: PetscReal h;
23: PetscReal beta = 1.0;
24: DM adda;
25: PetscInt nodes[2];
26: PetscBool periodic[2];
27: PetscInt refine[2];
28: PetscRandom rctx;
29: PetscMPIInt comm_size;
30: Mat H;
31: PetscInt *lcs, *lce;
32: PetscInt x, y;
33: PetscReal r1, r2;
34: PetscScalar uxy1, uxy2;
35: ADDAIdx sxy, sxy_m;
36: PetscScalar val, valconj;
37: Mat HtH;
38: Vec b, Htb;
39: Vec xvec;
40: KSP kspmg;
41: PC pcmg;
42: PetscErrorCode ierr;
44: PetscInitialize(&Argc,&Args,(char *)0,help);
45: PetscOptionsGetInt(PETSC_NULL,"-size",&n,&flg);
46: PetscOptionsGetReal(PETSC_NULL,"-beta",&beta,&flg);
47: PetscOptionsGetScalar(PETSC_NULL,"-rho",&rho,&flg);
49: /* Set the fudge parameters, we scale the whole thing by 1/(2*h) later */
50: h = 1.;
51: rho *= 1./(2.*h);
52:
53: /* Geometry info */
54: for(i=0; i<2; i++) {
55: nodes[i] = n;
56: periodic[i] = PETSC_TRUE;
57: refine[i] = 3;
58: }
59: DMADDACreate(PETSC_COMM_WORLD, 2, nodes, PETSC_NULL, 2,periodic, &adda);
60: DMADDASetRefinement(adda, refine, 2);
61:
62: /* Random numbers */
63: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
64: PetscRandomSetFromOptions(rctx);
66: /* Single or multi processor ? */
67: MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);
69: /* construct matrix */
70: if( comm_size == 1 ) {
71: DMCreateMatrix(adda, MATSEQAIJ, &H);
72: } else {
73: DMCreateMatrix(adda, MATMPIAIJ, &H);
74: }
76: /* get local corners for this processor, user is responsible for freeing lcs,lce */
77: DMADDAGetCorners(adda, &lcs, &lce);
79: /* Allocate space for the indices that we use to construct the matrix */
80: PetscMalloc(2*sizeof(PetscInt), &(sxy.x));
81: PetscMalloc(2*sizeof(PetscInt), &(sxy_m.x));
83: /* Assemble the matrix */
84: for( x=lcs[0]; x<lce[0]; x++ ) {
85: for( y=lcs[1]; y<lce[1]; y++ ) {
86: /* each lattice point sets only the *forward* pointing parameters (right, down),
87: i.e. Nabla_1^+ and Nabla_2^+.
88: In this way we can use only local random number creation. That means
89: we also have to set the corresponding backward pointing entries. */
90: /* Compute some normally distributed random numbers via Box-Muller */
91: PetscRandomGetValueReal(rctx, &r1);
92: r1 = 1.-r1; /* to change from [0,1) to (0,1], which we need for the log */
93: PetscRandomGetValueReal(rctx, &r2);
94: PetscReal R = sqrt(-2.*log(r1));
95: PetscReal c = cos(2.*PETSC_PI*r2);
96: PetscReal s = sin(2.*PETSC_PI*r2);
98: /* use those to set the field */
99: uxy1 = PetscExpScalar( ((PetscScalar) (R*c/beta))*PETSC_i);
100: uxy2 = PetscExpScalar( ((PetscScalar) (R*s/beta))*PETSC_i);
101:
102: sxy.x[0] = x; sxy.x[1] = y; /* the point where we are */
104: /* center action */
105: sxy.d = 0; /* spin 0, 0 */
106: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy, &rho, ADD_VALUES);
107: sxy.d = 1; /* spin 1, 1 */
108: val = -rho;
109: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy, &val, ADD_VALUES);
110:
111: sxy_m.x[0] = x+1; sxy_m.x[1] = y; /* right action */
112: sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
113: val = -uxy1; valconj = PetscConj(val);
114: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
115: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
116: sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
117: val = -uxy1; valconj = PetscConj(val);
118: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
119: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
120: sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
121: val = uxy1; valconj = PetscConj(val);
122: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
123: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
124: sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
125: val = uxy1; valconj = PetscConj(val);
126: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
127: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
129: sxy_m.x[0] = x; sxy_m.x[1] = y+1; /* down action */
130: sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
131: val = -uxy2; valconj = PetscConj(val);
132: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
133: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
134: sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
135: val = -PETSC_i*uxy2; valconj = PetscConj(val);
136: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
137: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
138: sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
139: val = -PETSC_i*uxy2; valconj = PetscConj(val);
140: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
141: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
142: sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
143: val = PetscConj(uxy2); valconj = PetscConj(val);
144: DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
145: DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
146: }
147: }
148:
149: PetscFree(sxy.x);
150: PetscFree(sxy_m.x);
152: PetscFree(lcs);
153: PetscFree(lce);
155: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
158: /* scale H */
159: MatScale(H, 1./(2.*h));
161: /* construct normal equations */
162: MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &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: /* ADDACreatematrix(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: DMCreateGlobalVector(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);
206: /* construct solver */
207: KSPCreate(PETSC_COMM_WORLD,&kspmg);
208: KSPSetType(kspmg, KSPCG);
210: KSPGetPC(kspmg,&pcmg);
211: PCSetType(pcmg,PCASA);
213: /* maybe user wants to override some of the choices */
214: KSPSetFromOptions(kspmg);
216: KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);
218: PCASASetDM(pcmg, (DM) adda);
219: DMDestroy(&adda);
221: PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);
223: VecDuplicate(b, &xvec);
224: VecSet(xvec, 0.0);
226: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227: Solve the linear system
228: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230: KSPSolve(kspmg, Htb, xvec);
232: /* VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */
234: KSPDestroy(&kspmg);
236: VecDestroy(&xvec);
237: VecDestroy(&b);
238: VecDestroy(&Htb);
239: MatDestroy(&H);
240: MatDestroy(&HtH);
242: PetscRandomDestroy(&rctx);
243: PetscFinalize();
244: return 0;
245: }
247: /* --------------------------------------------------------------------- */
250: PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig) {
251: PetscErrorCode ierr;
252: PetscRandom rctx; /* random number generator context */
253: Vec x0, x, x_1, tmp;
254: PetscScalar lambda_its, lambda_its_1;
255: PetscReal norm;
256: Mat G;
257: PetscInt i;
258:
260: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
261: PetscRandomSetFromOptions(rctx);
263: /* compute G = I-1/norm(A)*A */
264: MatNorm(A, NORM_1, &norm);
265: MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &G);
266: MatShift(G, -norm);
267: MatScale(G, -1./norm);
269: MatGetVecs(G, &x_1, &x);
270: VecSetRandom(x, rctx);
271: VecDuplicate(x, &x0);
272: VecCopy(x, x0);
274: MatMult(G, x, x_1);
275: for(i=0; i<its; i++) {
276: tmp = x; x = x_1; x_1 = tmp;
277: MatMult(G, x, x_1);
278: }
279: VecDot(x0, x, &lambda_its);
280: VecDot(x0, x_1, &lambda_its_1);
282: *eig = norm*(1.-lambda_its_1/lambda_its);
284: VecDestroy(&x0);
285: VecDestroy(&x);
286: VecDestroy(&x_1);
287: PetscRandomDestroy(&rctx);
288: MatDestroy(&G);
290: return(0);
291: }