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