Actual source code: brdn.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
  7:   representation in 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 = (S[i]^T dX) / (S[i]^T Q[i])
 19:     dX <- dX + (tau * (S[i] - Q[i]))
 20:   end
 21:  */

 23: static PetscErrorCode MatSolve_LMVMBrdn(Mat B, Vec F, Vec dX)
 24: {
 25:   Mat_LMVM   *lmvm  = (Mat_LMVM *)B->data;
 26:   Mat_Brdn   *lbrdn = (Mat_Brdn *)lmvm->ctx;
 27:   PetscInt    i, j;
 28:   PetscScalar sjtqi, stx, stq;

 30:   VecCheckSameSize(F, 2, dX, 3);
 31:   VecCheckMatCompatible(B, dX, 3, F, 2);

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

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

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

 57: /*
 58:   The forward product is the matrix-free implementation of Equation 2 in
 59:   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 = (S[i]^T X) / (S[i]^T S[i])
 72:     dX <- dX + (tau * (Y[i] - P[i]))
 73:   end
 74:  */

 76: static PetscErrorCode MatMult_LMVMBrdn(Mat B, Vec X, Vec Z)
 77: {
 78:   Mat_LMVM   *lmvm  = (Mat_LMVM *)B->data;
 79:   Mat_Brdn   *lbrdn = (Mat_Brdn *)lmvm->ctx;
 80:   PetscInt    i, j;
 81:   PetscScalar sjtsi, stx;

 83:   VecCheckSameSize(X, 2, Z, 3);
 84:   VecCheckMatCompatible(B, X, 2, Z, 3);

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

 98:   MatLMVMApplyJ0Fwd(B, X, Z);
 99:   for (i = 0; i <= lmvm->k; ++i) {
100:     VecDot(lmvm->S[i], X, &stx);
101:     VecAXPBYPCZ(Z, PetscRealPart(stx) / lbrdn->sts[i], -PetscRealPart(stx) / lbrdn->sts[i], 1.0, lmvm->Y[i], lbrdn->P[i]);
102:   }
103:   return 0;
104: }

106: /*------------------------------------------------------------*/

108: static PetscErrorCode MatUpdate_LMVMBrdn(Mat B, Vec X, Vec F)
109: {
110:   Mat_LMVM   *lmvm  = (Mat_LMVM *)B->data;
111:   Mat_Brdn   *lbrdn = (Mat_Brdn *)lmvm->ctx;
112:   PetscInt    old_k, i;
113:   PetscScalar sts;

115:   if (!lmvm->m) return 0;
116:   if (lmvm->prev_set) {
117:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
118:     VecAYPX(lmvm->Xprev, -1.0, X);
119:     VecAYPX(lmvm->Fprev, -1.0, F);
120:     /* Accept the update */
121:     lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
122:     old_k                       = lmvm->k;
123:     MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
124:     /* If we hit the memory limit, shift the sts array */
125:     if (old_k == lmvm->k) {
126:       for (i = 0; i <= lmvm->k - 1; ++i) lbrdn->sts[i] = lbrdn->sts[i + 1];
127:     }
128:     VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &sts);
129:     lbrdn->sts[lmvm->k] = PetscRealPart(sts);
130:   }
131:   /* Save the solution and function to be used in the next update */
132:   VecCopy(X, lmvm->Xprev);
133:   VecCopy(F, lmvm->Fprev);
134:   lmvm->prev_set = PETSC_TRUE;
135:   return 0;
136: }

138: /*------------------------------------------------------------*/

140: static PetscErrorCode MatCopy_LMVMBrdn(Mat B, Mat M, MatStructure str)
141: {
142:   Mat_LMVM *bdata = (Mat_LMVM *)B->data;
143:   Mat_Brdn *bctx  = (Mat_Brdn *)bdata->ctx;
144:   Mat_LMVM *mdata = (Mat_LMVM *)M->data;
145:   Mat_Brdn *mctx  = (Mat_Brdn *)mdata->ctx;
146:   PetscInt  i;

148:   mctx->needP = bctx->needP;
149:   mctx->needQ = bctx->needQ;
150:   for (i = 0; i <= bdata->k; ++i) {
151:     mctx->sts[i] = bctx->sts[i];
152:     mctx->stq[i] = bctx->stq[i];
153:     VecCopy(bctx->P[i], mctx->P[i]);
154:     VecCopy(bctx->Q[i], mctx->Q[i]);
155:   }
156:   return 0;
157: }

159: /*------------------------------------------------------------*/

161: static PetscErrorCode MatReset_LMVMBrdn(Mat B, PetscBool destructive)
162: {
163:   Mat_LMVM *lmvm  = (Mat_LMVM *)B->data;
164:   Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;

166:   lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
167:   if (destructive && lbrdn->allocated) {
168:     PetscFree2(lbrdn->sts, lbrdn->stq);
169:     VecDestroyVecs(lmvm->m, &lbrdn->P);
170:     VecDestroyVecs(lmvm->m, &lbrdn->Q);
171:     lbrdn->allocated = PETSC_FALSE;
172:   }
173:   MatReset_LMVM(B, destructive);
174:   return 0;
175: }

177: /*------------------------------------------------------------*/

179: static PetscErrorCode MatAllocate_LMVMBrdn(Mat B, Vec X, Vec F)
180: {
181:   Mat_LMVM *lmvm  = (Mat_LMVM *)B->data;
182:   Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;

184:   MatAllocate_LMVM(B, X, F);
185:   if (!lbrdn->allocated) {
186:     PetscMalloc2(lmvm->m, &lbrdn->sts, lmvm->m, &lbrdn->stq);
187:     if (lmvm->m > 0) {
188:       VecDuplicateVecs(X, lmvm->m, &lbrdn->P);
189:       VecDuplicateVecs(X, lmvm->m, &lbrdn->Q);
190:     }
191:     lbrdn->allocated = PETSC_TRUE;
192:   }
193:   return 0;
194: }

196: /*------------------------------------------------------------*/

198: static PetscErrorCode MatDestroy_LMVMBrdn(Mat B)
199: {
200:   Mat_LMVM *lmvm  = (Mat_LMVM *)B->data;
201:   Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;

203:   if (lbrdn->allocated) {
204:     PetscFree2(lbrdn->sts, lbrdn->stq);
205:     VecDestroyVecs(lmvm->m, &lbrdn->P);
206:     VecDestroyVecs(lmvm->m, &lbrdn->Q);
207:     lbrdn->allocated = PETSC_FALSE;
208:   }
209:   PetscFree(lmvm->ctx);
210:   MatDestroy_LMVM(B);
211:   return 0;
212: }

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

216: static PetscErrorCode MatSetUp_LMVMBrdn(Mat B)
217: {
218:   Mat_LMVM *lmvm  = (Mat_LMVM *)B->data;
219:   Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;

221:   MatSetUp_LMVM(B);
222:   if (!lbrdn->allocated) {
223:     PetscMalloc2(lmvm->m, &lbrdn->sts, lmvm->m, &lbrdn->stq);
224:     if (lmvm->m > 0) {
225:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbrdn->P);
226:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbrdn->Q);
227:     }
228:     lbrdn->allocated = PETSC_TRUE;
229:   }
230:   return 0;
231: }

233: /*------------------------------------------------------------*/

235: PetscErrorCode MatCreate_LMVMBrdn(Mat B)
236: {
237:   Mat_LMVM *lmvm;
238:   Mat_Brdn *lbrdn;

240:   MatCreate_LMVM(B);
241:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMBROYDEN);
242:   B->ops->setup   = MatSetUp_LMVMBrdn;
243:   B->ops->destroy = MatDestroy_LMVMBrdn;
244:   B->ops->solve   = MatSolve_LMVMBrdn;

246:   lmvm                = (Mat_LMVM *)B->data;
247:   lmvm->square        = PETSC_TRUE;
248:   lmvm->ops->allocate = MatAllocate_LMVMBrdn;
249:   lmvm->ops->reset    = MatReset_LMVMBrdn;
250:   lmvm->ops->mult     = MatMult_LMVMBrdn;
251:   lmvm->ops->update   = MatUpdate_LMVMBrdn;
252:   lmvm->ops->copy     = MatCopy_LMVMBrdn;

254:   PetscNew(&lbrdn);
255:   lmvm->ctx        = (void *)lbrdn;
256:   lbrdn->allocated = PETSC_FALSE;
257:   lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
258:   return 0;
259: }

261: /*------------------------------------------------------------*/

263: /*@
264:    MatCreateLMVMBroyden - Creates a limited-memory "good" Broyden-type approximation
265:    matrix used for a Jacobian. L-Brdn is not guaranteed to be symmetric or
266:    positive-definite.

268:    The provided local and global sizes must match the solution and function vectors
269:    used with `MatLMVMUpdate()` and `MatSolve()`. The resulting L-Brdn matrix will have
270:    storage vectors allocated with `VecCreateSeq()` in serial and `VecCreateMPI()` in
271:    parallel. To use the L-Brdn matrix with other vector types, the matrix must be
272:    created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
273:    This ensures that the internal storage and work vectors are duplicated from the
274:    correct type of vector.

276:    Collective

278:    Input Parameters:
279: +  comm - MPI communicator
280: .  n - number of local rows for storage vectors
281: -  N - global size of the storage vectors

283:    Output Parameter:
284: .  B - the matrix

286:    Level: intermediate

288:    Note:
289:    It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
290:    paradigm instead of this routine directly.

292: .seealso: [](chapter_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
293:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBadBrdn()`, `MatCreateLMVMSymBrdn()`
294: @*/
295: PetscErrorCode MatCreateLMVMBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
296: {
297:   MatCreate(comm, B);
298:   MatSetSizes(*B, n, n, N, N);
299:   MatSetType(*B, MATLMVMBROYDEN);
300:   MatSetUp(*B);
301:   return 0;
302: }