Actual source code: ex40.c

petsc-3.4.4 2014-03-13
  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;

 44:   PetscInitialize(&Argc,&Args,(char*)0,help);
 45:   PetscOptionsGetInt(NULL,"-size",&n,&flg);
 46:   PetscOptionsGetReal(NULL,"-beta",&beta,&flg);
 47:   PetscOptionsGetScalar(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);

 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, NULL, 2,periodic, &adda);
 60:   DMADDASetRefinement(adda, refine, 2);

 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 = PetscSqrtReal(-2.*log(r1));
 95:       PetscReal c = PetscCosReal(2.*PETSC_PI*r2);
 96:       PetscReal s = PetscSinReal(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);

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);

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);

115:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
116:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);

118:       sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
119:       val   = -uxy1; valconj = PetscConj(val);

121:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
122:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);

124:       sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
125:       val   = uxy1; valconj = PetscConj(val);

127:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
128:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);

130:       sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
131:       val   = uxy1; valconj = PetscConj(val);

133:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
134:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);

136:       sxy_m.x[0] = x; sxy_m.x[1] = y+1; /* down action */
137:       sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
138:       val   = -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);

143:       sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
144:       val   = -PETSC_i*uxy2; valconj = PetscConj(val);

146:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
147:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);

149:       sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
150:       val   = -PETSC_i*uxy2; valconj = PetscConj(val);

152:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
153:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);

155:       sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
156:       val   = PetscConj(uxy2); valconj = PetscConj(val);

158:       DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);
159:       DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);
160:     }
161:   }

163:   PetscFree(sxy.x);
164:   PetscFree(sxy_m.x);

166:   PetscFree(lcs);
167:   PetscFree(lce);

169:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
170:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);

172:   /* scale H */
173:   MatScale(H, 1./(2.*h));

175:   /* construct normal equations */
176:   MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);

178:   PetscScalar mineval;
179:   computeMinEigVal(HtH, 1000, &mineval);
180:   PetscPrintf(PETSC_COMM_WORLD, "Minimum eigenvalue of H^{dag} H is %f\n", PetscAbsScalar(mineval));

182:   /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
183: /*   Mat perm; */
184: /*   ADDACreatematrix(adda, MATSEQAIJ, &perm); */
185: /*   PetscInt row, col; */
186: /*   PetscScalar one = 1.0; */
187: /*   for (PetscInt i=0; i<n; i++) { */
188: /*     for (PetscInt j=0; j<n; j++) { */
189: /*       row = (i*n+j)*2; col = i*n+j; */
190: /*       MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
191: /*       row = (i*n+j)*2+1; col = i*n+j + n*n; */
192: /*       MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES); */
193: /*     } */
194: /*   } */
195: /*   MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY); */
196: /*   MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY); */

198: /*   Mat Hperm; */
199: /*   MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm); */
200: /*   PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n"); */
201: /*   MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */

203: /*   Mat HtHperm; */
204: /*   MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm); */
205: /*   PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n"); */
206: /*   MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */

208:   /* right hand side */
209:   DMCreateGlobalVector(adda, &b);
210:   VecSet(b,0.0);
211:   PetscInt    ix[1]   = {0};
212:   PetscScalar vals[1] = {1.0};
213:   VecSetValues(b, 1, ix, vals, INSERT_VALUES);
214:   VecAssemblyBegin(b);
215:   VecAssemblyEnd(b);
216: /*   VecSetRandom(b, rctx); */
217:   VecDuplicate(b, &Htb);
218:   MatMultTranspose(H, b, Htb);

220:   /* construct solver */
221:   KSPCreate(PETSC_COMM_WORLD,&kspmg);
222:   KSPSetType(kspmg, KSPCG);

224:   KSPGetPC(kspmg,&pcmg);
225:   PCSetType(pcmg,PCASA);

227:   /* maybe user wants to override some of the choices */
228:   KSPSetFromOptions(kspmg);

230:   KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);

232:   PCSetDM(pcmg,adda);
233:   DMDestroy(&adda);

235:   PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);

237:   VecDuplicate(b, &xvec);
238:   VecSet(xvec, 0.0);

240:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241:                       Solve the linear system
242:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

244:   KSPSolve(kspmg, Htb, xvec);

246: /*   VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD)); */

248:   KSPDestroy(&kspmg);

250:   VecDestroy(&xvec);
251:   VecDestroy(&b);
252:   VecDestroy(&Htb);
253:   MatDestroy(&H);
254:   MatDestroy(&HtH);

256:   PetscRandomDestroy(&rctx);
257:   PetscFinalize();
258:   return 0;
259: }

261: /* --------------------------------------------------------------------- */
264: PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig)
265: {
267:   PetscRandom    rctx;      /* random number generator context */
268:   Vec            x0, x, x_1, tmp;
269:   PetscScalar    lambda_its, lambda_its_1;
270:   PetscReal      norm;
271:   Mat            G;
272:   PetscInt       i;

275:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
276:   PetscRandomSetFromOptions(rctx);

278:   /* compute G = I-1/norm(A)*A */
279:   MatNorm(A, NORM_1, &norm);
280:   MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &G);
281:   MatShift(G, -norm);
282:   MatScale(G, -1./norm);

284:   MatGetVecs(G, &x_1, &x);
285:   VecSetRandom(x, rctx);
286:   VecDuplicate(x, &x0);
287:   VecCopy(x, x0);

289:   MatMult(G, x, x_1);
290:   for (i=0; i<its; i++) {
291:     tmp  = x; x = x_1; x_1 = tmp;
292:     MatMult(G, x, x_1);
293:   }
294:   VecDot(x0, x, &lambda_its);
295:   VecDot(x0, x_1, &lambda_its_1);

297:   *eig = norm*(1.-lambda_its_1/lambda_its);

299:   VecDestroy(&x0);
300:   VecDestroy(&x);
301:   VecDestroy(&x_1);
302:   PetscRandomDestroy(&rctx);
303:   MatDestroy(&G);
304:   return(0);
305: }