Actual source code: dfp.c

      1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
      2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>

      4: /*
      5:   Limited-memory Davidon-Fletcher-Powell method for approximating both
      6:   the forward product and inverse application of a Jacobian.
      7:  */

      9: /*------------------------------------------------------------*/

     11: /*
     12:   The solution method (approximate inverse Jacobian application) is
     13:   matrix-vector product version of the recursive formula given in
     14:   Equation (6.15) of Nocedal and Wright "Numerical Optimization" 2nd
     15:   edition, pg 139.

     17:   Note: Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
     18:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
     19:   repeated calls of MatSolve without incurring redundant computation.

     21:   dX <- J0^{-1} * F

     23:   for i = 0,1,2,...,k
     24:     # Q[i] = (B_i)^{-1} * Y[i]
     25:     gamma = (S[i]^T F) / (Y[i]^T S[i])
     26:     zeta = (Y[i]^T dX) / (Y[i]^T Q[i])
     27:     dX <- dX + (gamma * S[i]) - (zeta * Q[i])
     28:   end
     29: */
     30: PetscErrorCode MatSolve_LMVMDFP(Mat B, Vec F, Vec dX)
     31: {
     32:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
     33:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
     34:   PetscErrorCode    ierr;
     35:   PetscInt          i, j;
     36:   PetscScalar       yjtqi, sjtyi, ytx, stf, ytq;

     41:   VecCheckSameSize(F, 2, dX, 3);
     42:   VecCheckMatCompatible(B, dX, 3, F, 2);

     44:   if (ldfp->needQ) {
     45:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
     46:     for (i = 0; i <= lmvm->k; ++i) {
     47:       MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], ldfp->Q[i]);
     48:       for (j = 0; j <= i-1; ++j) {
     49:         /* Compute the necessary dot products */
     50:         VecDotBegin(lmvm->Y[j], ldfp->Q[i], &yjtqi);
     51:         VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
     52:         VecDotEnd(lmvm->Y[j], ldfp->Q[i], &yjtqi);
     53:         VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
     54:         /* Compute the pure DFP component of the inverse application*/
     55:         VecAXPBYPCZ(ldfp->Q[i], -PetscRealPart(yjtqi)/ldfp->ytq[j], PetscRealPart(sjtyi)/ldfp->yts[j], 1.0, ldfp->Q[j], lmvm->S[j]);
     56:       }
     57:       VecDot(lmvm->Y[i], ldfp->Q[i], &ytq);
     58:       ldfp->ytq[i] = PetscRealPart(ytq);
     59:     }
     60:     ldfp->needQ = PETSC_FALSE;
     61:   }

     63:   /* Start the outer loop (i) for the recursive formula */
     64:   MatSymBrdnApplyJ0Inv(B, F, dX);
     65:   for (i = 0; i <= lmvm->k; ++i) {
     66:     /* Get all the dot products we need */
     67:     VecDotBegin(lmvm->Y[i], dX, &ytx);
     68:     VecDotBegin(lmvm->S[i], F, &stf);
     69:     VecDotEnd(lmvm->Y[i], dX, &ytx);
     70:     VecDotEnd(lmvm->S[i], F, &stf);
     71:     /* Update dX_{i+1} = (B^{-1})_{i+1} * f */
     72:     VecAXPBYPCZ(dX, -PetscRealPart(ytx)/ldfp->ytq[i], PetscRealPart(stf)/ldfp->yts[i], 1.0, ldfp->Q[i], lmvm->S[i]);
     73:   }
     74:   return(0);
     75: }

     77: /*------------------------------------------------------------*/

     79: /*
     80:   The forward product for the approximate Jacobian is the matrix-free
     81:   implementation of the recursive formula given in Equation 6.13 of
     82:   Nocedal and Wright "Numerical Optimization" 2nd edition, pg 139.

     84:   This forward product has a two-loop form similar to the BFGS two-loop
     85:   formulation for the inverse Jacobian application. However, the S and
     86:   Y vectors have interchanged roles.

     88:   work <- X

     90:   for i = k,k-1,k-2,...,0
     91:     rho[i] = 1 / (Y[i]^T S[i])
     92:     alpha[i] = rho[i] * (Y[i]^T work)
     93:     work <- work - (alpha[i] * S[i])
     94:   end

     96:   Z <- J0 * work

     98:   for i = 0,1,2,...,k
     99:     beta = rho[i] * (S[i]^T Y)
    100:     Z <- Z + ((alpha[i] - beta) * Y[i])
    101:   end
    102: */
    103: PetscErrorCode MatMult_LMVMDFP(Mat B, Vec X, Vec Z)
    104: {
    105:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    106:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
    107:   PetscErrorCode    ierr;
    108:   PetscInt          i;
    109:   PetscReal         *alpha, beta;
    110:   PetscScalar       ytx, stz;

    113:   /* Copy the function into the work vector for the first loop */
    114:   VecCopy(X, ldfp->work);

    116:   /* Start the first loop */
    117:   PetscMalloc1(lmvm->k+1, &alpha);
    118:   for (i = lmvm->k; i >= 0; --i) {
    119:     VecDot(lmvm->Y[i], ldfp->work, &ytx);
    120:     alpha[i] = PetscRealPart(ytx)/ldfp->yts[i];
    121:     VecAXPY(ldfp->work, -alpha[i], lmvm->S[i]);
    122:   }

    124:   /* Apply the forward product with initial Jacobian */
    125:   MatSymBrdnApplyJ0Fwd(B, ldfp->work, Z);

    127:   /* Start the second loop */
    128:   for (i = 0; i <= lmvm->k; ++i) {
    129:     VecDot(lmvm->S[i], Z, &stz);
    130:     beta = PetscRealPart(stz)/ldfp->yts[i];
    131:     VecAXPY(Z, alpha[i]-beta, lmvm->Y[i]);
    132:   }
    133:   PetscFree(alpha);
    134:   return(0);
    135: }

    137: /*------------------------------------------------------------*/

    139: static PetscErrorCode MatUpdate_LMVMDFP(Mat B, Vec X, Vec F)
    140: {
    141:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    142:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
    143:   Mat_LMVM          *dbase;
    144:   Mat_DiagBrdn      *dctx;
    145:   PetscErrorCode    ierr;
    146:   PetscInt          old_k, i;
    147:   PetscReal         curvtol;
    148:   PetscScalar       curvature, ytytmp, ststmp;

    151:   if (!lmvm->m) return(0);
    152:   if (lmvm->prev_set) {
    153:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
    154:     VecAYPX(lmvm->Xprev, -1.0, X);
    155:     VecAYPX(lmvm->Fprev, -1.0, F);
    156:     /* Test if the updates can be accepted */
    157:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
    158:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
    159:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
    160:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
    161:     if (PetscRealPart(ststmp) < lmvm->eps) {
    162:       curvtol = 0.0;
    163:     } else {
    164:       curvtol = lmvm->eps * PetscRealPart(ststmp);
    165:     }
    166:     if (PetscRealPart(curvature) > curvtol) {
    167:       /* Update is good, accept it */
    168:       ldfp->watchdog = 0;
    169:       ldfp->needQ = PETSC_TRUE;
    170:       old_k = lmvm->k;
    171:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
    172:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
    173:       if (old_k == lmvm->k) {
    174:         for (i = 0; i <= lmvm->k-1; ++i) {
    175:           ldfp->yts[i] = ldfp->yts[i+1];
    176:           ldfp->yty[i] = ldfp->yty[i+1];
    177:           ldfp->sts[i] = ldfp->sts[i+1];
    178:         }
    179:       }
    180:       /* Update history of useful scalars */
    181:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
    182:       ldfp->yts[lmvm->k] = PetscRealPart(curvature);
    183:       ldfp->yty[lmvm->k] = PetscRealPart(ytytmp);
    184:       ldfp->sts[lmvm->k] = PetscRealPart(ststmp);
    185:       /* Compute the scalar scale if necessary */
    186:       if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) {
    187:         MatSymBrdnComputeJ0Scalar(B);
    188:       }
    189:     } else {
    190:       /* Update is bad, skip it */
    191:       ++lmvm->nrejects;
    192:       ++ldfp->watchdog;
    193:     }
    194:   } else {
    195:     switch (ldfp->scale_type) {
    196:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    197:       dbase = (Mat_LMVM*)ldfp->D->data;
    198:       dctx = (Mat_DiagBrdn*)dbase->ctx;
    199:       VecSet(dctx->invD, ldfp->delta);
    200:       break;
    201:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    202:       ldfp->sigma = ldfp->delta;
    203:       break;
    204:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    205:       ldfp->sigma = 1.0;
    206:       break;
    207:     default:
    208:       break;
    209:     }
    210:   }

    212:   /* Update the scaling */
    213:   if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
    214:     MatLMVMUpdate(ldfp->D, X, F);
    215:   }

    217:   if (ldfp->watchdog > ldfp->max_seq_rejects) {
    218:     MatLMVMReset(B, PETSC_FALSE);
    219:     if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
    220:       MatLMVMReset(ldfp->D, PETSC_FALSE);
    221:     }
    222:   }

    224:   /* Save the solution and function to be used in the next update */
    225:   VecCopy(X, lmvm->Xprev);
    226:   VecCopy(F, lmvm->Fprev);
    227:   lmvm->prev_set = PETSC_TRUE;
    228:   return(0);
    229: }

    231: /*------------------------------------------------------------*/

    233: static PetscErrorCode MatCopy_LMVMDFP(Mat B, Mat M, MatStructure str)
    234: {
    235:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
    236:   Mat_SymBrdn       *bctx = (Mat_SymBrdn*)bdata->ctx;
    237:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
    238:   Mat_SymBrdn       *mctx = (Mat_SymBrdn*)mdata->ctx;
    239:   PetscErrorCode    ierr;
    240:   PetscInt          i;

    243:   mctx->needQ = bctx->needQ;
    244:   for (i=0; i<=bdata->k; ++i) {
    245:     mctx->ytq[i] = bctx->ytq[i];
    246:     mctx->yts[i] = bctx->yts[i];
    247:     VecCopy(bctx->Q[i], mctx->Q[i]);
    248:   }
    249:   mctx->scale_type      = bctx->scale_type;
    250:   mctx->alpha           = bctx->alpha;
    251:   mctx->beta            = bctx->beta;
    252:   mctx->rho             = bctx->rho;
    253:   mctx->sigma_hist      = bctx->sigma_hist;
    254:   mctx->watchdog        = bctx->watchdog;
    255:   mctx->max_seq_rejects = bctx->max_seq_rejects;
    256:   switch (bctx->scale_type) {
    257:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    258:     mctx->sigma = bctx->sigma;
    259:     break;
    260:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    261:     MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
    262:     break;
    263:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    264:     mctx->sigma = 1.0;
    265:     break;
    266:   default:
    267:     break;
    268:   }
    269:   return(0);
    270: }

    272: /*------------------------------------------------------------*/

    274: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
    275: {
    276:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    277:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
    278:   Mat_LMVM          *dbase;
    279:   Mat_DiagBrdn      *dctx;
    280:   PetscErrorCode    ierr;

    283:   ldfp->watchdog = 0;
    284:   ldfp->needQ = PETSC_TRUE;
    285:   if (ldfp->allocated) {
    286:     if (destructive) {
    287:       VecDestroy(&ldfp->work);
    288:       PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
    289:       VecDestroyVecs(lmvm->m, &ldfp->Q);
    290:       switch (ldfp->scale_type) {
    291:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    292:         MatLMVMReset(ldfp->D, PETSC_TRUE);
    293:         break;
    294:       default:
    295:         break;
    296:       }
    297:       ldfp->allocated = PETSC_FALSE;
    298:     } else {
    299:       switch (ldfp->scale_type) {
    300:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    301:         ldfp->sigma = ldfp->delta;
    302:         break;
    303:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    304:         MatLMVMReset(ldfp->D, PETSC_FALSE);
    305:         dbase = (Mat_LMVM*)ldfp->D->data;
    306:         dctx = (Mat_DiagBrdn*)dbase->ctx;
    307:         VecSet(dctx->invD, ldfp->delta);
    308:         break;
    309:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    310:         ldfp->sigma = 1.0;
    311:         break;
    312:       default:
    313:         break;
    314:       }
    315:     }
    316:   }
    317:   MatReset_LMVM(B, destructive);
    318:   return(0);
    319: }

    321: /*------------------------------------------------------------*/

    323: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
    324: {
    325:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    326:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
    327:   PetscErrorCode    ierr;

    330:   MatAllocate_LMVM(B, X, F);
    331:   if (!ldfp->allocated) {
    332:     VecDuplicate(X, &ldfp->work);
    333:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
    334:     if (lmvm->m > 0) {
    335:       VecDuplicateVecs(X, lmvm->m, &ldfp->Q);
    336:     }
    337:     switch (ldfp->scale_type) {
    338:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    339:       MatLMVMAllocate(ldfp->D, X, F);
    340:       break;
    341:     default:
    342:       break;
    343:     }
    344:     ldfp->allocated = PETSC_TRUE;
    345:   }
    346:   return(0);
    347: }

    349: /*------------------------------------------------------------*/

    351: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
    352: {
    353:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    354:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
    355:   PetscErrorCode    ierr;

    358:   if (ldfp->allocated) {
    359:     VecDestroy(&ldfp->work);
    360:     PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
    361:     VecDestroyVecs(lmvm->m, &ldfp->Q);
    362:     ldfp->allocated = PETSC_FALSE;
    363:   }
    364:   MatDestroy(&ldfp->D);
    365:   PetscFree(lmvm->ctx);
    366:   MatDestroy_LMVM(B);
    367:   return(0);
    368: }

    370: /*------------------------------------------------------------*/

    372: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
    373: {
    374:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    375:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
    376:   PetscErrorCode    ierr;
    377:   PetscInt          n, N;

    380:   MatSetUp_LMVM(B);
    381:   if (!ldfp->allocated) {
    382:     VecDuplicate(lmvm->Xprev, &ldfp->work);
    383:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
    384:     if (lmvm->m > 0) {
    385:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q);
    386:     }
    387:     switch (ldfp->scale_type) {
    388:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    389:       MatGetLocalSize(B, &n, &n);
    390:       MatGetSize(B, &N, &N);
    391:       MatSetSizes(ldfp->D, n, n, N, N);
    392:       MatSetUp(ldfp->D);
    393:       break;
    394:     default:
    395:       break;
    396:     }
    397:     ldfp->allocated = PETSC_TRUE;
    398:   }
    399:   return(0);
    400: }

    402: /*------------------------------------------------------------*/

    404: static PetscErrorCode MatSetFromOptions_LMVMDFP(PetscOptionItems *PetscOptionsObject, Mat B)
    405: {
    406:   PetscErrorCode               ierr;

    409:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
    410:   PetscOptionsHead(PetscOptionsObject,"DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
    411:   MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
    412:   PetscOptionsTail();
    413:   return(0);
    414: }

    416: /*------------------------------------------------------------*/

    418: PetscErrorCode MatCreate_LMVMDFP(Mat B)
    419: {
    420:   Mat_LMVM          *lmvm;
    421:   Mat_SymBrdn       *ldfp;
    422:   PetscErrorCode    ierr;

    425:   MatCreate_LMVMSymBrdn(B);
    426:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP);
    427:   B->ops->setup = MatSetUp_LMVMDFP;
    428:   B->ops->destroy = MatDestroy_LMVMDFP;
    429:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
    430:   B->ops->solve = MatSolve_LMVMDFP;

    432:   lmvm = (Mat_LMVM*)B->data;
    433:   lmvm->ops->allocate = MatAllocate_LMVMDFP;
    434:   lmvm->ops->reset = MatReset_LMVMDFP;
    435:   lmvm->ops->update = MatUpdate_LMVMDFP;
    436:   lmvm->ops->mult = MatMult_LMVMDFP;
    437:   lmvm->ops->copy = MatCopy_LMVMDFP;

    439:   ldfp = (Mat_SymBrdn*)lmvm->ctx;
    440:   ldfp->needP           = PETSC_FALSE;
    441:   ldfp->phi             = 1.0;
    442:   return(0);
    443: }

    445: /*------------------------------------------------------------*/

    447: /*@
    448:    MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix
    449:    used for approximating Jacobians. L-DFP is symmetric positive-definite by
    450:    construction, and is the dual of L-BFGS where Y and S vectors swap roles.

    452:    The provided local and global sizes must match the solution and function vectors
    453:    used with MatLMVMUpdate() and MatSolve(). The resulting L-DFP matrix will have
    454:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
    455:    parallel. To use the L-DFP matrix with other vector types, the matrix must be
    456:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
    457:    This ensures that the internal storage and work vectors are duplicated from the
    458:    correct type of vector.

    460:    Collective

    462:    Input Parameters:
    463: +  comm - MPI communicator, set to PETSC_COMM_SELF
    464: .  n - number of local rows for storage vectors
    465: -  N - global size of the storage vectors

    467:    Output Parameter:
    468: .  B - the matrix

    470:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
    471:    paradigm instead of this routine directly.

    473:    Options Database Keys:
    474: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
    475: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
    476: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    477: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
    478: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    479: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
    480: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

    482:    Level: intermediate

    484: .seealso: MatCreate(), MATLMVM, MATLMVMDFP, MatCreateLMVMBFGS(), MatCreateLMVMSR1(),
    485:            MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
    486: @*/
    487: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    488: {
    489:   PetscErrorCode    ierr;

    492:   MatCreate(comm, B);
    493:   MatSetSizes(*B, n, n, N, N);
    494:   MatSetType(*B, MATLMVMDFP);
    495:   MatSetUp(*B);
    496:   return(0);
    497: }