Actual source code: symbadbrdn.c

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

      4: /*------------------------------------------------------------*/

      6: static PetscErrorCode MatSolve_LMVMSymBadBrdn(Mat B, Vec F, Vec dX)
      7: {
      8:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
      9:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
     10:   PetscErrorCode    ierr;
     11:   PetscInt          i, j;
     12:   PetscScalar       yjtqi, sjtyi, wtyi, ytx, stf, wtf, ytq;

     15:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
     16:   if (lsb->phi == 0.0) {
     17:     MatSolve_LMVMBFGS(B, F, dX);
     18:     return(0);
     19:   }
     20:   if (lsb->phi == 1.0) {
     21:     MatSolve_LMVMDFP(B, F, dX);
     22:     return(0);
     23:   }

     25:   VecCheckSameSize(F, 2, dX, 3);
     26:   VecCheckMatCompatible(B, dX, 3, F, 2);

     28:   if (lsb->needQ) {
     29:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
     30:     for (i = 0; i <= lmvm->k; ++i) {
     31:       MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);
     32:       for (j = 0; j <= i-1; ++j) {
     33:         /* Compute the necessary dot products */
     34:         VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);
     35:         VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
     36:         VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);
     37:         VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
     38:         /* Compute the pure DFP component of the inverse application*/
     39:         VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi)/lsb->ytq[j], PetscRealPart(sjtyi)/lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]);
     40:         /* Tack on the convexly scaled extras to the inverse application*/
     41:         if (lsb->psi[j] > 0.0) {
     42:           VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);
     43:           VecDot(lsb->work, lmvm->Y[i], &wtyi);
     44:           VecAXPY(lsb->Q[i], lsb->phi*lsb->ytq[j]*PetscRealPart(wtyi), lsb->work);
     45:         }
     46:       }
     47:       VecDot(lmvm->Y[i], lsb->Q[i], &ytq);
     48:       lsb->ytq[i] = PetscRealPart(ytq);
     49:     }
     50:     lsb->needQ = PETSC_FALSE;
     51:   }

     53:   /* Start the outer iterations for ((B^{-1}) * dX) */
     54:   MatSymBrdnApplyJ0Inv(B, F, dX);
     55:   for (i = 0; i <= lmvm->k; ++i) {
     56:     /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
     57:     VecDotBegin(lmvm->Y[i], dX, &ytx);
     58:     VecDotBegin(lmvm->S[i], F, &stf);
     59:     VecDotEnd(lmvm->Y[i], dX, &ytx);
     60:     VecDotEnd(lmvm->S[i], F, &stf);
     61:     /* Compute the pure DFP component */
     62:     VecAXPBYPCZ(dX, -PetscRealPart(ytx)/lsb->ytq[i], PetscRealPart(stf)/lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]);
     63:     /* Tack on the convexly scaled extras */
     64:     VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]);
     65:     VecDot(lsb->work, F, &wtf);
     66:     VecAXPY(dX, lsb->phi*lsb->ytq[i]*PetscRealPart(wtf), lsb->work);
     67:   }

     69:   return(0);
     70: }

     72: /*------------------------------------------------------------*/

     74: static PetscErrorCode MatMult_LMVMSymBadBrdn(Mat B, Vec X, Vec Z)
     75: {
     76:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
     77:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
     78:   PetscErrorCode    ierr;
     79:   PetscInt          i, j;
     80:   PetscReal         numer;
     81:   PetscScalar       sjtpi, sjtyi, yjtsi, yjtqi, wtsi, wtyi, stz, ytx, ytq, wtx, stp;

     85:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
     86:   if (lsb->phi == 0.0) {
     87:     MatMult_LMVMBFGS(B, X, Z);
     88:     return(0);
     89:   }
     90:   if (lsb->phi == 1.0) {
     91:     MatMult_LMVMDFP(B, X, Z);
     92:     return(0);
     93:   }

     95:   VecCheckSameSize(X, 2, Z, 3);
     96:   VecCheckMatCompatible(B, X, 2, Z, 3);

     98:   if (lsb->needQ) {
     99:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
    100:     for (i = 0; i <= lmvm->k; ++i) {
    101:       MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);
    102:       for (j = 0; j <= i-1; ++j) {
    103:         /* Compute the necessary dot products */
    104:         VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);
    105:         VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
    106:         VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);
    107:         VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
    108:         /* Compute the pure DFP component of the inverse application*/
    109:         VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi)/lsb->ytq[j], PetscRealPart(sjtyi)/lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]);
    110:         /* Tack on the convexly scaled extras to the inverse application*/
    111:         if (lsb->psi[j] > 0.0) {
    112:           VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);
    113:           VecDot(lsb->work, lmvm->Y[i], &wtyi);
    114:           VecAXPY(lsb->Q[i], lsb->phi*lsb->ytq[j]*PetscRealPart(wtyi), lsb->work);
    115:         }
    116:       }
    117:       VecDot(lmvm->Y[i], lsb->Q[i], &ytq);
    118:       lsb->ytq[i] = PetscRealPart(ytq);
    119:     }
    120:     lsb->needQ = PETSC_FALSE;
    121:   }
    122:   if (lsb->needP) {
    123:     /* Start the loop for (P[k] = (B_k) * S[k]) */
    124:     for (i = 0; i <= lmvm->k; ++i) {
    125:       MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
    126:       for (j = 0; j <= i-1; ++j) {
    127:         /* Compute the necessary dot products */
    128:         VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
    129:         VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
    130:         VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
    131:         VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
    132:         /* Compute the pure BFGS component of the forward product */
    133:         VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
    134:         /* Tack on the convexly scaled extras to the forward product */
    135:         if (lsb->phi > 0.0) {
    136:           VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
    137:           VecDot(lsb->work, lmvm->S[i], &wtsi);
    138:           VecAXPY(lsb->P[i], lsb->psi[j]*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);
    139:         }
    140:       }
    141:       VecDot(lmvm->S[i], lsb->P[i], &stp);
    142:       lsb->stp[i] = PetscRealPart(stp);
    143:       if (lsb->phi == 1.0) {
    144:         lsb->psi[i] = 0.0;
    145:       } else if (lsb->phi == 0.0) {
    146:         lsb->psi[i] = 1.0;
    147:       } else {
    148:         numer = (1.0 - lsb->phi)*lsb->yts[i]*lsb->yts[i];
    149:         lsb->psi[i] = numer / (numer + (lsb->phi*lsb->ytq[i]*lsb->stp[i]));
    150:       }
    151:     }
    152:     lsb->needP = PETSC_FALSE;
    153:   }

    155:   /* Start the outer iterations for (B * X) */
    156:   MatSymBrdnApplyJ0Fwd(B, X, Z);
    157:   for (i = 0; i <= lmvm->k; ++i) {
    158:     /* Compute the necessary dot products */
    159:     VecDotBegin(lmvm->S[i], Z, &stz);
    160:     VecDotBegin(lmvm->Y[i], X, &ytx);
    161:     VecDotEnd(lmvm->S[i], Z, &stz);
    162:     VecDotEnd(lmvm->Y[i], X, &ytx);
    163:     /* Compute the pure BFGS component */
    164:     VecAXPBYPCZ(Z, -PetscRealPart(stz)/lsb->stp[i], PetscRealPart(ytx)/lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
    165:     /* Tack on the convexly scaled extras */
    166:     VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
    167:     VecDot(lsb->work, X, &wtx);
    168:     VecAXPY(Z, lsb->psi[i]*lsb->stp[i]*PetscRealPart(wtx), lsb->work);
    169:   }
    170:   return(0);
    171: }

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

    175: static PetscErrorCode MatSetFromOptions_LMVMSymBadBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
    176: {
    177:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    178:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
    179:   Mat_LMVM          *dbase;
    180:   Mat_DiagBrdn      *dctx;
    181:   PetscErrorCode    ierr;

    184:   MatSetFromOptions_LMVMSymBrdn(PetscOptionsObject, B);
    185:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
    186:     dbase = (Mat_LMVM*)lsb->D->data;
    187:     dctx = (Mat_DiagBrdn*)dbase->ctx;
    188:     dctx->forward = PETSC_FALSE;
    189:   }
    190:   return(0);
    191: }

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

    195: PetscErrorCode MatCreate_LMVMSymBadBrdn(Mat B)
    196: {
    197:   Mat_LMVM          *lmvm;
    198:   PetscErrorCode    ierr;

    201:   MatCreate_LMVMSymBrdn(B);
    202:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBADBROYDEN);
    203:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBadBrdn;
    204:   B->ops->solve = MatSolve_LMVMSymBadBrdn;

    206:   lmvm = (Mat_LMVM*)B->data;
    207:   lmvm->ops->mult = MatMult_LMVMSymBadBrdn;
    208:   return(0);
    209: }

    211: /*------------------------------------------------------------*/

    213: /*@
    214:    MatCreateLMVMSymBadBroyden - Creates a limited-memory Symmetric "Bad" Broyden-type matrix used
    215:    for approximating Jacobians. L-SymBadBrdn is a convex combination of L-DFP and
    216:    L-BFGS such that `^{-1} = (1 - phi)*BFGS^{-1} + phi*DFP^{-1}. The combination factor
    217:    phi is restricted to the range [0, 1], where the L-SymBadBrdn matrix is guaranteed
    218:    to be symmetric positive-definite. Note that this combination is on the inverses and not
    219:    on the forwards. For forward convex combinations, use the L-SymBrdn matrix.

    221:    The provided local and global sizes must match the solution and function vectors
    222:    used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
    223:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
    224:    parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be
    225:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
    226:    This ensures that the internal storage and work vectors are duplicated from the
    227:    correct type of vector.

    229:    Collective

    231:    Input Parameters:
    232: +  comm - MPI communicator, set to PETSC_COMM_SELF
    233: .  n - number of local rows for storage vectors
    234: -  N - global size of the storage vectors

    236:    Output Parameter:
    237: .  B - the matrix

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

    242:    Options Database Keys:
    243: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
    244: .   -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
    245: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
    246: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    247: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
    248: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    249: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
    250: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

    252:    Level: intermediate

    254: .seealso: MatCreate(), MATLMVM, MATLMVMSYMBROYDEN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
    255:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
    256: @*/
    257: PetscErrorCode MatCreateLMVMSymBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    258: {
    259:   PetscErrorCode    ierr;

    262:   MatCreate(comm, B);
    263:   MatSetSizes(*B, n, n, N, N);
    264:   MatSetType(*B, MATLMVMSYMBADBROYDEN);
    265:   MatSetUp(*B);
    266:   return(0);
    267: }