Actual source code: ex39.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 <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: }