Actual source code: symbadbrdn.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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;
 13: 
 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:   }
 24: 
 25:   VecCheckSameSize(F, 2, dX, 3);
 26:   VecCheckMatCompatible(B, dX, 3, F, 2);
 27: 
 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:   }
 52: 
 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;
 82: 
 83: 
 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:   }
 94: 
 95:   VecCheckSameSize(X, 2, Z, 3);
 96:   VecCheckMatCompatible(B, X, 2, Z, 3);
 97: 
 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:   }
154: 
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 == SYMBRDN_SCALE_DIAG) {
186:     dbase = (Mat_LMVM*)lsb->D->data;
187:     dctx = (Mat_DiagBrdn*)dbase->ctx;
188:     dctx->forward = PETSC_FALSE;
189:   }
190:   if (Scale_Table[0][0] == '#') SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale_Table starts with illegal hashtag character!");  /* Dummy use of Scale_Table to prevent unused variable warnings in this translation unit */
191:   return(0);
192: }

194: /*------------------------------------------------------------*/

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

202:   MatCreate_LMVMSymBrdn(B);
203:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBADBRDN);
204:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBadBrdn;
205:   B->ops->solve = MatSolve_LMVMSymBadBrdn;
206: 
207:   lmvm = (Mat_LMVM*)B->data;
208:   lmvm->ops->mult = MatMult_LMVMSymBadBrdn;
209:   return(0);
210: }

212: /*------------------------------------------------------------*/

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

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

230:    Collective

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

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

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

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

253:    Level: intermediate

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

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