Actual source code: snesnoise.c

      1: #include <petsc/private/snesimpl.h>

      3: PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES, Vec, void **);
      4: PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES, void *, Vec, Vec, PetscReal *, PetscReal *);
      5: PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void *);

      7: /* Data used by Jorge's diff parameter computation method */
      8: typedef struct {
      9:   Vec     *workv;          /* work vectors */
     10:   FILE    *fp;             /* output file */
     11:   PetscInt function_count; /* count of function evaluations for diff param estimation */
     12:   double   fnoise_min;     /* minimum allowable noise */
     13:   double   hopt_min;       /* minimum allowable hopt */
     14:   double   h_first_try;    /* first try for h used in diff parameter estimate */
     15:   PetscInt fnoise_resets;  /* number of times we've reset the noise estimate */
     16:   PetscInt hopt_resets;    /* number of times we've reset the hopt estimate */
     17: } DIFFPAR_MORE;

     19: PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes);
     20: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);

     22: static PetscErrorCode JacMatMultCompare(SNES, Vec, Vec, double);

     24: PetscErrorCode SNESDiffParameterCreate_More(SNES snes, Vec x, void **outneP)
     25: {
     26:   DIFFPAR_MORE *neP;
     27:   Vec           w;
     28:   PetscRandom   rctx; /* random number generator context */
     29:   PetscBool     flg;
     30:   char          noise_file[PETSC_MAX_PATH_LEN];

     32:   PetscFunctionBegin;
     33:   PetscCall(PetscNew(&neP));

     35:   neP->function_count = 0;
     36:   neP->fnoise_min     = 1.0e-20;
     37:   neP->hopt_min       = 1.0e-8;
     38:   neP->h_first_try    = 1.0e-3;
     39:   neP->fnoise_resets  = 0;
     40:   neP->hopt_resets    = 0;

     42:   /* Create work vectors */
     43:   PetscCall(VecDuplicateVecs(x, 3, &neP->workv));
     44:   w = neP->workv[0];

     46:   /* Set components of vector w to random numbers */
     47:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)snes), &rctx));
     48:   PetscCall(PetscRandomSetFromOptions(rctx));
     49:   PetscCall(VecSetRandom(w, rctx));
     50:   PetscCall(PetscRandomDestroy(&rctx));

     52:   /* Open output file */
     53:   PetscCall(PetscOptionsGetString(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_noise_file", noise_file, sizeof(noise_file), &flg));
     54:   if (flg) neP->fp = fopen(noise_file, "w");
     55:   else neP->fp = fopen("noise.out", "w");
     56:   PetscCheck(neP->fp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file");
     57:   PetscCall(PetscInfo(snes, "Creating Jorge's differencing parameter context\n"));

     59:   *outneP = neP;
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: PetscErrorCode SNESDiffParameterDestroy_More(void *nePv)
     64: {
     65:   DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv;
     66:   int           err;

     68:   PetscFunctionBegin;
     69:   /* Destroy work vectors and close output file */
     70:   PetscCall(VecDestroyVecs(3, &neP->workv));
     71:   err = fclose(neP->fp);
     72:   PetscCheck(!err, PETSC_COMM_SELF, PETSC_ERR_SYS, "fclose() failed on file");
     73:   PetscCall(PetscFree(neP));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: PetscErrorCode SNESDiffParameterCompute_More(SNES snes, void *nePv, Vec x, Vec p, double *fnoise, double *hopt)
     78: {
     79:   DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv;
     80:   Vec           w, xp, fvec; /* work vectors to use in computing h */
     81:   double        zero = 0.0, hl, hu, h, fnoise_s, fder2_s;
     82:   PetscScalar   alpha;
     83:   PetscScalar   fval[7], tab[7][7], eps[7], f = -1;
     84:   double        rerrf = -1., fder2;
     85:   PetscInt      iter, k, i, j, info;
     86:   PetscInt      nf = 7; /* number of function evaluations */
     87:   PetscInt      fcount;
     88:   MPI_Comm      comm;
     89:   FILE         *fp;
     90:   PetscBool     noise_test = PETSC_FALSE;

     92:   PetscFunctionBegin;
     93:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
     94:   /* Call to SNESSetUp() just to set data structures in SNES context */
     95:   if (!snes->setupcalled) PetscCall(SNESSetUp(snes));

     97:   w    = neP->workv[0];
     98:   xp   = neP->workv[1];
     99:   fvec = neP->workv[2];
    100:   fp   = neP->fp;

    102:   /* Initialize parameters */
    103:   hl       = zero;
    104:   hu       = zero;
    105:   h        = neP->h_first_try;
    106:   fnoise_s = zero;
    107:   fder2_s  = zero;
    108:   fcount   = neP->function_count;

    110:   /* We have 5 tries to attempt to compute a good hopt value */
    111:   PetscCall(SNESGetIterationNumber(snes, &i));
    112:   PetscCall(PetscFPrintf(comm, fp, "\n ------- SNES iteration %" PetscInt_FMT " ---------\n", i));
    113:   for (iter = 0; iter < 5; iter++) {
    114:     neP->h_first_try = h;

    116:     /* Compute the nf function values needed to estimate the noise from
    117:        the difference table */
    118:     for (k = 0; k < nf; k++) {
    119:       alpha = h * (k + 1 - (nf + 1) / 2);
    120:       PetscCall(VecWAXPY(xp, alpha, p, x));
    121:       PetscCall(SNESComputeFunction(snes, xp, fvec));
    122:       neP->function_count++;
    123:       PetscCall(VecDot(fvec, w, &fval[k]));
    124:     }
    125:     f = fval[(nf + 1) / 2 - 1];

    127:     /* Construct the difference table */
    128:     for (i = 0; i < nf; i++) tab[i][0] = fval[i];

    130:     for (j = 0; j < nf - 1; j++) {
    131:       for (i = 0; i < nf - j - 1; i++) tab[i][j + 1] = tab[i + 1][j] - tab[i][j];
    132:     }

    134:     /* Print the difference table */
    135:     PetscCall(PetscFPrintf(comm, fp, "Difference Table: iter = %" PetscInt_FMT "\n", iter));
    136:     for (i = 0; i < nf; i++) {
    137:       for (j = 0; j < nf - i; j++) PetscCall(PetscFPrintf(comm, fp, " %10.2e ", tab[i][j]));
    138:       PetscCall(PetscFPrintf(comm, fp, "\n"));
    139:     }

    141:     /* Call the noise estimator */
    142:     PetscCall(SNESNoise_dnest_(&nf, fval, &h, fnoise, &fder2, hopt, &info, eps));

    144:     /* Output statements */
    145:     rerrf = *fnoise / PetscAbsScalar(f);
    146:     if (info == 1) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected"));
    147:     if (info == 2) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too small"));
    148:     if (info == 3) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too large"));
    149:     if (info == 4) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected, but unreliable hopt"));
    150:     PetscCall(PetscFPrintf(comm, fp, "Approximate epsfcn %g  %g  %g  %g  %g  %g\n", (double)eps[0], (double)eps[1], (double)eps[2], (double)eps[3], (double)eps[4], (double)eps[5]));
    151:     PetscCall(PetscFPrintf(comm, fp, "h = %g, fnoise = %g, fder2 = %g, rerrf = %g, hopt = %g\n\n", (double)h, (double)*fnoise, (double)fder2, (double)rerrf, (double)*hopt));

    153:     /* Save fnoise and fder2. */
    154:     if (*fnoise) fnoise_s = *fnoise;
    155:     if (fder2) fder2_s = fder2;

    157:     /* Check for noise detection. */
    158:     if (fnoise_s && fder2_s) {
    159:       *fnoise = fnoise_s;
    160:       fder2   = fder2_s;
    161:       *hopt   = 1.68 * sqrt(*fnoise / PetscAbsScalar(fder2));
    162:       goto theend;
    163:     } else {
    164:       /* Update hl and hu, and determine new h */
    165:       if (info == 2 || info == 4) {
    166:         hl = h;
    167:         if (hu == zero) h = 100 * h;
    168:         else h = PetscMin(100 * h, 0.1 * hu);
    169:       } else if (info == 3) {
    170:         hu = h;
    171:         h  = PetscMax(1.0e-3, sqrt(hl / hu)) * hu;
    172:       }
    173:     }
    174:   }
    175: theend:

    177:   if (*fnoise < neP->fnoise_min) {
    178:     PetscCall(PetscFPrintf(comm, fp, "Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n", (double)*fnoise, (double)neP->fnoise_min));
    179:     *fnoise = neP->fnoise_min;
    180:     neP->fnoise_resets++;
    181:   }
    182:   if (*hopt < neP->hopt_min) {
    183:     PetscCall(PetscFPrintf(comm, fp, "Resetting hopt: hopt1 = %g, hopt_min = %g\n", (double)*hopt, (double)neP->hopt_min));
    184:     *hopt = neP->hopt_min;
    185:     neP->hopt_resets++;
    186:   }

    188:   PetscCall(PetscFPrintf(comm, fp, "Errors in derivative:\n"));
    189:   PetscCall(PetscFPrintf(comm, fp, "f = %g, fnoise = %g, fder2 = %g, hopt = %g\n", (double)f, (double)*fnoise, (double)fder2, (double)*hopt));

    191:   /* For now, compute h **each** MV Mult!! */
    192:   /*
    193:   PetscCall(PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg));
    194:   if (!flg) {
    195:     Mat mat;
    196:     PetscCall(SNESGetJacobian(snes,&mat,NULL,NULL));
    197:     PetscCall(MatSNESMFMoreSetParameters(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt));
    198:   }
    199:   */
    200:   fcount = neP->function_count - fcount;
    201:   PetscCall(PetscInfo(snes, "fct_now = %" PetscInt_FMT ", fct_cum = %" PetscInt_FMT ", rerrf=%g, sqrt(noise)=%g, h_more=%g\n", fcount, neP->function_count, (double)rerrf, (double)PetscSqrtReal(*fnoise), (double)*hopt));

    203:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-noise_test", &noise_test, NULL));
    204:   if (noise_test) PetscCall(JacMatMultCompare(snes, x, p, *hopt));
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: static PetscErrorCode JacMatMultCompare(SNES snes, Vec x, Vec p, double hopt)
    209: {
    210:   Vec         yy1, yy2; /* work vectors */
    211:   PetscViewer view2;    /* viewer */
    212:   Mat         J;        /* analytic Jacobian (set as matrix from which to compute the preconditioner) */
    213:   Mat         Jmf;      /* matrix-free Jacobian (set as true system matrix) */
    214:   double      h;        /* differencing parameter */
    215:   Vec         f;
    216:   PetscScalar alpha;
    217:   PetscReal   yy1n, yy2n, enorm;
    218:   PetscInt    i;
    219:   PetscBool   printv = PETSC_FALSE;
    220:   char        filename[32];
    221:   MPI_Comm    comm;

    223:   PetscFunctionBegin;
    224:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
    225:   /* Compute function and analytic Jacobian at x */
    226:   PetscCall(SNESGetJacobian(snes, &Jmf, &J, NULL, NULL));
    227:   PetscCall(SNESComputeJacobian(snes, x, Jmf, J));
    228:   PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
    229:   PetscCall(SNESComputeFunction(snes, x, f));

    231:   /* Duplicate work vectors */
    232:   PetscCall(VecDuplicate(x, &yy2));
    233:   PetscCall(VecDuplicate(x, &yy1));

    235:   /* Compute true matrix-vector product */
    236:   PetscCall(MatMult(J, p, yy1));
    237:   PetscCall(VecNorm(yy1, NORM_2, &yy1n));

    239:   /* View product vector if desired */
    240:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_vecs", &printv, NULL));
    241:   if (printv) {
    242:     PetscCall(PetscViewerASCIIOpen(comm, "y1.out", &view2));
    243:     PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON));
    244:     PetscCall(VecView(yy1, view2));
    245:     PetscCall(PetscViewerPopFormat(view2));
    246:     PetscCall(PetscViewerDestroy(&view2));
    247:   }

    249:   /* Test Jacobian-vector product computation */
    250:   alpha = -1.0;
    251:   h     = 0.01 * hopt;
    252:   for (i = 0; i < 5; i++) {
    253:     /* Set differencing parameter for matrix-free multiplication */
    254:     PetscCall(MatSNESMFMoreSetParameters(Jmf, PETSC_DEFAULT, PETSC_DEFAULT, h));

    256:     /* Compute matrix-vector product via differencing approximation */
    257:     PetscCall(MatMult(Jmf, p, yy2));
    258:     PetscCall(VecNorm(yy2, NORM_2, &yy2n));

    260:     /* View product vector if desired */
    261:     if (printv) {
    262:       PetscCall(PetscSNPrintf(filename, PETSC_STATIC_ARRAY_LENGTH(filename), "y2.%d.out", (int)i));
    263:       PetscCall(PetscViewerASCIIOpen(comm, filename, &view2));
    264:       PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON));
    265:       PetscCall(VecView(yy2, view2));
    266:       PetscCall(PetscViewerPopFormat(view2));
    267:       PetscCall(PetscViewerDestroy(&view2));
    268:     }

    270:     /* Compute relative error */
    271:     PetscCall(VecAXPY(yy2, alpha, yy1));
    272:     PetscCall(VecNorm(yy2, NORM_2, &enorm));
    273:     enorm = enorm / yy1n;
    274:     PetscCall(PetscFPrintf(comm, stdout, "h = %g: relative error = %g\n", (double)h, (double)enorm));
    275:     h *= 10.0;
    276:   }
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }