Actual source code: ex39.c
petsc-3.4.2 2013-07-02

      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: }