Actual source code: badbrdn.c

      1: #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h>

      3: /*------------------------------------------------------------*/

      5: /*
      6:   The solution method is the matrix-free implementation of the inverse Hessian in
      7:   Equation 6 on page 312 of Griewank "Broyden Updating, The Good and The Bad!"
      8:   (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).

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

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

     16:   for i=0,1,2,...,k
     17:     # Q[i] = (B_i)^{-1} * Y[i]
     18:     tau = (Y[i]^T F) / (Y[i]^T Y[i])
     19:     dX <- dX + (tau * (S[i] - Q[i]))
     20:   end
     21:  */

     23: static PetscErrorCode MatSolve_LMVMBadBrdn(Mat B, Vec F, Vec dX)
     24: {
     25:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
     26:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
     27:   PetscErrorCode    ierr;
     28:   PetscInt          i, j;
     29:   PetscScalar       yjtyi, ytf;

     32:   VecCheckSameSize(F, 2, dX, 3);
     33:   VecCheckMatCompatible(B, dX, 3, F, 2);

     35:   if (lbb->needQ) {
     36:     /* Pre-compute (Q[i] = (B_i)^{-1} * Y[i]) */
     37:     for (i = 0; i <= lmvm->k; ++i) {
     38:       MatLMVMApplyJ0Inv(B, lmvm->Y[i], lbb->Q[i]);
     39:       for (j = 0; j <= i-1; ++j) {
     40:         VecDot(lmvm->Y[j], lmvm->Y[i], &yjtyi);
     41:         VecAXPBYPCZ(lbb->Q[i], PetscRealPart(yjtyi)/lbb->yty[j], -PetscRealPart(yjtyi)/lbb->yty[j], 1.0, lmvm->S[j], lbb->Q[j]);
     42:       }
     43:     }
     44:     lbb->needQ = PETSC_FALSE;
     45:   }

     47:   MatLMVMApplyJ0Inv(B, F, dX);
     48:   for (i = 0; i <= lmvm->k; ++i) {
     49:     VecDot(lmvm->Y[i], F, &ytf);
     50:     VecAXPBYPCZ(dX, PetscRealPart(ytf)/lbb->yty[i], -PetscRealPart(ytf)/lbb->yty[i], 1.0, lmvm->S[i], lbb->Q[i]);
     51:   }
     52:   return(0);
     53: }

     55: /*------------------------------------------------------------*/

     57: /*
     58:   The forward product is the matrix-free implementation of the direct update in
     59:   Equation 6 on page 302 of Griewank "Broyden Updating, The Good and The Bad!"
     60:   (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).

     62:   P[i] = (B_i)*S[i] terms are computed ahead of time whenever
     63:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
     64:   repeated calls of MatMult inside KSP solvers without unnecessarily
     65:   recomputing P[i] terms in expensive nested-loops.

     67:   Z <- J0 * X

     69:   for i=0,1,2,...,k
     70:     # P[i] = B_i * S[i]
     71:     tau = (Y[i]^T X) / (Y[i]^T S[i])
     72:     dX <- dX + (tau * (Y[i] - P[i]))
     73:   end
     74:  */

     76: static PetscErrorCode MatMult_LMVMBadBrdn(Mat B, Vec X, Vec Z)
     77: {
     78:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
     79:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
     80:   PetscErrorCode    ierr;
     81:   PetscInt          i, j;
     82:   PetscScalar       yjtsi, ytx;

     85:   VecCheckSameSize(X, 2, Z, 3);
     86:   VecCheckMatCompatible(B, X, 2, Z, 3);

     88:   if (lbb->needP) {
     89:     /* Pre-compute (P[i] = (B_i) * S[i]) */
     90:     for (i = 0; i <= lmvm->k; ++i) {
     91:       MatLMVMApplyJ0Fwd(B, lmvm->S[i], lbb->P[i]);
     92:       for (j = 0; j <= i-1; ++j) {
     93:         VecDot(lmvm->Y[j], lmvm->S[i], &yjtsi);
     94:         VecAXPBYPCZ(lbb->P[i], PetscRealPart(yjtsi)/lbb->yts[j], -PetscRealPart(yjtsi)/lbb->yts[j], 1.0, lmvm->Y[j], lbb->P[j]);
     95:       }
     96:     }
     97:     lbb->needP = PETSC_FALSE;
     98:   }

    100:   MatLMVMApplyJ0Fwd(B, X, Z);
    101:   for (i = 0; i <= lmvm->k; ++i) {
    102:     VecDot(lmvm->Y[i], X, &ytx);
    103:     VecAXPBYPCZ(Z, PetscRealPart(ytx)/lbb->yts[i], -PetscRealPart(ytx)/lbb->yts[i], 1.0, lmvm->Y[i], lbb->P[i]);
    104:   }
    105:   return(0);
    106: }

    108: /*------------------------------------------------------------*/

    110: static PetscErrorCode MatUpdate_LMVMBadBrdn(Mat B, Vec X, Vec F)
    111: {
    112:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    113:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
    114:   PetscErrorCode    ierr;
    115:   PetscInt          old_k, i;
    116:   PetscScalar       yty, yts;

    119:   if (!lmvm->m) return(0);
    120:   if (lmvm->prev_set) {
    121:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
    122:     VecAYPX(lmvm->Xprev, -1.0, X);
    123:     VecAYPX(lmvm->Fprev, -1.0, F);
    124:     /* Accept the update */
    125:     lbb->needP = lbb->needQ = PETSC_TRUE;
    126:     old_k = lmvm->k;
    127:     MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
    128:     /* If we hit the memory limit, shift the yty and yts arrays */
    129:     if (old_k == lmvm->k) {
    130:       for (i = 0; i <= lmvm->k-1; ++i) {
    131:         lbb->yty[i] = lbb->yty[i+1];
    132:         lbb->yts[i] = lbb->yts[i+1];
    133:       }
    134:     }
    135:     /* Accumulate the latest yTy and yTs dot products */
    136:     VecDotBegin(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
    137:     VecDotBegin(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts);
    138:     VecDotEnd(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
    139:     VecDotEnd(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts);
    140:     lbb->yty[lmvm->k] = PetscRealPart(yty);
    141:     lbb->yts[lmvm->k] = PetscRealPart(yts);
    142:   }
    143:   /* Save the solution and function to be used in the next update */
    144:   VecCopy(X, lmvm->Xprev);
    145:   VecCopy(F, lmvm->Fprev);
    146:   lmvm->prev_set = PETSC_TRUE;
    147:   return(0);
    148: }

    150: /*------------------------------------------------------------*/

    152: static PetscErrorCode MatCopy_LMVMBadBrdn(Mat B, Mat M, MatStructure str)
    153: {
    154:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
    155:   Mat_Brdn          *bctx = (Mat_Brdn*)bdata->ctx;
    156:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
    157:   Mat_Brdn          *mctx = (Mat_Brdn*)mdata->ctx;
    158:   PetscErrorCode    ierr;
    159:   PetscInt          i;

    162:   mctx->needP = bctx->needP;
    163:   mctx->needQ = bctx->needQ;
    164:   for (i=0; i<=bdata->k; ++i) {
    165:     mctx->yty[i] = bctx->yty[i];
    166:     mctx->yts[i] = bctx->yts[i];
    167:     VecCopy(bctx->P[i], mctx->P[i]);
    168:     VecCopy(bctx->Q[i], mctx->Q[i]);
    169:   }
    170:   return(0);
    171: }

    173: /*------------------------------------------------------------*/

    175: static PetscErrorCode MatReset_LMVMBadBrdn(Mat B, PetscBool destructive)
    176: {
    177:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    178:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
    179:   PetscErrorCode    ierr;

    182:   lbb->needP = lbb->needQ = PETSC_TRUE;
    183:   if (destructive && lbb->allocated) {
    184:     PetscFree2(lbb->yty, lbb->yts);
    185:     VecDestroyVecs(lmvm->m, &lbb->P);
    186:     VecDestroyVecs(lmvm->m, &lbb->Q);
    187:     lbb->allocated = PETSC_FALSE;
    188:   }
    189:   MatReset_LMVM(B, destructive);
    190:   return(0);
    191: }

    193: /*------------------------------------------------------------*/

    195: static PetscErrorCode MatAllocate_LMVMBadBrdn(Mat B, Vec X, Vec F)
    196: {
    197:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    198:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
    199:   PetscErrorCode    ierr;

    202:   MatAllocate_LMVM(B, X, F);
    203:   if (!lbb->allocated) {
    204:     PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts);
    205:     if (lmvm->m > 0) {
    206:       VecDuplicateVecs(X, lmvm->m, &lbb->P);
    207:       VecDuplicateVecs(X, lmvm->m, &lbb->Q);
    208:     }
    209:     lbb->allocated = PETSC_TRUE;
    210:   }
    211:   return(0);
    212: }

    214: /*------------------------------------------------------------*/

    216: static PetscErrorCode MatDestroy_LMVMBadBrdn(Mat B)
    217: {
    218:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    219:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
    220:   PetscErrorCode    ierr;

    223:   if (lbb->allocated) {
    224:     PetscFree2(lbb->yty, lbb->yts);
    225:     VecDestroyVecs(lmvm->m, &lbb->P);
    226:     VecDestroyVecs(lmvm->m, &lbb->Q);
    227:     lbb->allocated = PETSC_FALSE;
    228:   }
    229:   PetscFree(lmvm->ctx);
    230:   MatDestroy_LMVM(B);
    231:   return(0);
    232: }

    234: /*------------------------------------------------------------*/

    236: static PetscErrorCode MatSetUp_LMVMBadBrdn(Mat B)
    237: {
    238:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    239:   Mat_Brdn          *lbb = (Mat_Brdn*)lmvm->ctx;
    240:   PetscErrorCode    ierr;

    243:   MatSetUp_LMVM(B);
    244:   if (!lbb->allocated) {
    245:     PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts);
    246:     if (lmvm->m > 0) {
    247:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->P);
    248:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->Q);
    249:     }
    250:     lbb->allocated = PETSC_TRUE;
    251:   }
    252:   return(0);
    253: }

    255: /*------------------------------------------------------------*/

    257: PetscErrorCode MatCreate_LMVMBadBrdn(Mat B)
    258: {
    259:   Mat_LMVM          *lmvm;
    260:   Mat_Brdn          *lbb;
    261:   PetscErrorCode    ierr;

    264:   MatCreate_LMVM(B);
    265:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMBADBROYDEN);
    266:   B->ops->setup = MatSetUp_LMVMBadBrdn;
    267:   B->ops->destroy = MatDestroy_LMVMBadBrdn;
    268:   B->ops->solve = MatSolve_LMVMBadBrdn;

    270:   lmvm = (Mat_LMVM*)B->data;
    271:   lmvm->square = PETSC_TRUE;
    272:   lmvm->ops->allocate = MatAllocate_LMVMBadBrdn;
    273:   lmvm->ops->reset = MatReset_LMVMBadBrdn;
    274:   lmvm->ops->mult = MatMult_LMVMBadBrdn;
    275:   lmvm->ops->update = MatUpdate_LMVMBadBrdn;
    276:   lmvm->ops->copy = MatCopy_LMVMBadBrdn;

    278:   PetscNewLog(B, &lbb);
    279:   lmvm->ctx = (void*)lbb;
    280:   lbb->allocated = PETSC_FALSE;
    281:   lbb->needP = lbb->needQ = PETSC_TRUE;
    282:   return(0);
    283: }

    285: /*------------------------------------------------------------*/

    287: /*@
    288:    MatCreateLMVMBadBroyden - Creates a limited-memory modified (aka "bad") Broyden-type
    289:    approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be
    290:    symmetric or positive-definite.

    292:    The provided local and global sizes must match the solution and function vectors
    293:    used with MatLMVMUpdate() and MatSolve(). The resulting L-BadBrdn matrix will have
    294:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
    295:    parallel. To use the L-BadBrdn matrix with other vector types, the matrix must be
    296:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
    297:    This ensures that the internal storage and work vectors are duplicated from the
    298:    correct type of vector.

    300:    Collective

    302:    Input Parameters:
    303: +  comm - MPI communicator, set to PETSC_COMM_SELF
    304: .  n - number of local rows for storage vectors
    305: -  N - global size of the storage vectors

    307:    Output Parameter:
    308: .  B - the matrix

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

    313:    Options Database Keys:
    314: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
    315: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
    316: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    317: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
    318: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    319: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
    320: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

    322:    Level: intermediate

    324: .seealso: MatCreate(), MATLMVM, MATLMVMBADBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
    325:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
    326: @*/
    327: PetscErrorCode MatCreateLMVMBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    328: {
    329:   PetscErrorCode    ierr;

    332:   MatCreate(comm, B);
    333:   MatSetSizes(*B, n, n, N, N);
    334:   MatSetType(*B, MATLMVMBADBROYDEN);
    335:   MatSetUp(*B);
    336:   return(0);
    337: }