Actual source code: denseqn.c
1: #include <../src/ksp/ksp/utils/lmvm/dense/denseqn.h>
2: #include <../src/ksp/ksp/utils/lmvm/blas_cyclic/blas_cyclic.h>
3: #include <petscblaslapack.h>
4: #include <petscmat.h>
5: #include <petscsys.h>
6: #include <petscsystypes.h>
7: #include <petscis.h>
8: #include <petscoptions.h>
9: #include <petscdevice.h>
10: #include <petsc/private/deviceimpl.h>
12: static PetscErrorCode MatMult_LMVMDQN(Mat, Vec, Vec);
13: static PetscErrorCode MatMult_LMVMDBFGS(Mat, Vec, Vec);
14: static PetscErrorCode MatMult_LMVMDDFP(Mat, Vec, Vec);
15: static PetscErrorCode MatSolve_LMVMDQN(Mat, Vec, Vec);
16: static PetscErrorCode MatSolve_LMVMDBFGS(Mat, Vec, Vec);
17: static PetscErrorCode MatSolve_LMVMDDFP(Mat, Vec, Vec);
19: static inline PetscInt recycle_index(PetscInt m, PetscInt idx)
20: {
21: return idx % m;
22: }
24: static inline PetscInt history_index(PetscInt m, PetscInt num_updates, PetscInt idx)
25: {
26: return (idx - num_updates) + PetscMin(m, num_updates);
27: }
29: static inline PetscInt oldest_update(PetscInt m, PetscInt idx)
30: {
31: return PetscMax(0, idx - m);
32: }
34: static PetscErrorCode MatView_LMVMDQN(Mat B, PetscViewer pv)
35: {
36: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
37: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
39: PetscBool isascii;
41: PetscFunctionBegin;
42: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
43: PetscCall(MatView_LMVM(B, pv));
44: PetscCall(SymBroydenRescaleView(lqn->rescale, pv));
45: if (isascii) PetscCall(PetscViewerASCIIPrintf(pv, "Counts: S x : %" PetscInt_FMT ", S^T x : %" PetscInt_FMT ", Y x : %" PetscInt_FMT ", Y^T x: %" PetscInt_FMT "\n", lqn->S_count, lqn->St_count, lqn->Y_count, lqn->Yt_count));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode MatLMVMDQNResetDestructive(Mat B)
50: {
51: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
52: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
54: PetscFunctionBegin;
55: PetscCall(MatDestroy(&lqn->HY));
56: PetscCall(MatDestroy(&lqn->BS));
57: PetscCall(MatDestroy(&lqn->StY_triu));
58: PetscCall(MatDestroy(&lqn->YtS_triu));
59: PetscCall(VecDestroy(&lqn->StFprev));
60: PetscCall(VecDestroy(&lqn->Fprev_ref));
61: lqn->Fprev_state = 0;
62: PetscCall(MatDestroy(&lqn->YtS_triu_strict));
63: PetscCall(MatDestroy(&lqn->StY_triu_strict));
64: PetscCall(MatDestroy(&lqn->StBS));
65: PetscCall(MatDestroy(&lqn->YtHY));
66: PetscCall(MatDestroy(&lqn->J));
67: PetscCall(MatDestroy(&lqn->temp_mat));
68: PetscCall(VecDestroy(&lqn->diag_vec));
69: PetscCall(VecDestroy(&lqn->diag_vec_recycle_order));
70: PetscCall(VecDestroy(&lqn->inv_diag_vec));
71: PetscCall(VecDestroy(&lqn->column_work));
72: PetscCall(VecDestroy(&lqn->column_work2));
73: PetscCall(VecDestroy(&lqn->rwork1));
74: PetscCall(VecDestroy(&lqn->rwork2));
75: PetscCall(VecDestroy(&lqn->rwork3));
76: PetscCall(VecDestroy(&lqn->rwork2_local));
77: PetscCall(VecDestroy(&lqn->rwork3_local));
78: PetscCall(VecDestroy(&lqn->cyclic_work_vec));
79: PetscCall(VecDestroyVecs(lmvm->m, &lqn->PQ));
80: PetscCall(PetscFree(lqn->stp));
81: PetscCall(PetscFree(lqn->yts));
82: PetscCall(PetscFree(lqn->ytq));
83: lqn->allocated = PETSC_FALSE;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode MatReset_LMVMDQN_Internal(Mat B, MatLMVMResetMode mode)
88: {
89: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
90: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
92: PetscFunctionBegin;
93: lqn->watchdog = 0;
94: lqn->needPQ = PETSC_TRUE;
95: lqn->num_updates = 0;
96: lqn->num_mult_updates = 0;
97: if (MatLMVMResetClearsBases(mode)) PetscCall(MatLMVMDQNResetDestructive(B));
98: else {
99: if (lqn->BS) PetscCall(MatZeroEntries(lqn->BS));
100: if (lqn->HY) PetscCall(MatZeroEntries(lqn->HY));
101: if (lqn->StY_triu) { /* Set to identity by default so it is invertible */
102: PetscCall(MatZeroEntries(lqn->StY_triu));
103: PetscCall(MatShift(lqn->StY_triu, 1.0));
104: }
105: if (lqn->YtS_triu) {
106: PetscCall(MatZeroEntries(lqn->YtS_triu));
107: PetscCall(MatShift(lqn->YtS_triu, 1.0));
108: }
109: if (lqn->YtS_triu_strict) PetscCall(MatZeroEntries(lqn->YtS_triu_strict));
110: if (lqn->StY_triu_strict) PetscCall(MatZeroEntries(lqn->StY_triu_strict));
111: if (lqn->StBS) {
112: PetscCall(MatZeroEntries(lqn->StBS));
113: PetscCall(MatShift(lqn->StBS, 1.0));
114: }
115: if (lqn->YtHY) {
116: PetscCall(MatZeroEntries(lqn->YtHY));
117: PetscCall(MatShift(lqn->YtHY, 1.0));
118: }
119: PetscCall(VecDestroy(&lqn->Fprev_ref));
120: lqn->Fprev_state = 0;
121: if (lqn->StFprev) PetscCall(VecZeroEntries(lqn->StFprev));
122: }
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode MatReset_LMVMDQN(Mat B, MatLMVMResetMode mode)
127: {
128: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
129: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
131: PetscFunctionBegin;
132: PetscCall(SymBroydenRescaleReset(B, lqn->rescale, mode));
133: PetscCall(MatReset_LMVMDQN_Internal(B, mode));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode MatAllocate_LMVMDQN_Internal(Mat B)
138: {
139: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
140: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
142: PetscFunctionBegin;
143: if (!lqn->allocated) {
144: if (lmvm->m > 0) {
145: PetscMPIInt rank;
146: PetscInt n, N, m, M;
147: PetscBool is_dbfgs, is_ddfp, is_dqn;
148: VecType vec_type;
149: MPI_Comm comm = PetscObjectComm((PetscObject)B);
150: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
152: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
153: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
154: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
156: PetscCallMPI(MPI_Comm_rank(comm, &rank));
157: PetscCall(MatGetSize(B, &N, NULL));
158: PetscCall(MatGetLocalSize(B, &n, NULL));
159: M = lmvm->m;
160: m = (rank == 0) ? M : 0;
162: /* For DBFGS: Create data needed for MatSolve() eagerly; data needed for MatMult() will be created on demand
163: * For DDFP : Create data needed for MatMult() eagerly; data needed for MatSolve() will be created on demand
164: * For DQN : Create all data eagerly */
165: PetscCall(VecGetType(lmvm->Xprev, &vec_type));
166: if (is_dqn) {
167: PetscCall(MatCreateDenseFromVecType(comm, vec_type, m, m, M, M, -1, NULL, &lqn->StY_triu));
168: PetscCall(MatCreateDenseFromVecType(comm, vec_type, m, m, M, M, -1, NULL, &lqn->YtS_triu));
169: PetscCall(MatCreateVecs(lqn->StY_triu, &lqn->diag_vec, &lqn->rwork1));
170: PetscCall(MatCreateVecs(lqn->StY_triu, &lqn->rwork2, &lqn->rwork3));
171: } else if (is_ddfp) {
172: PetscCall(MatCreateDenseFromVecType(comm, vec_type, m, m, M, M, -1, NULL, &lqn->YtS_triu));
173: PetscCall(MatDuplicate(Sfull, MAT_SHARE_NONZERO_PATTERN, &lqn->HY));
174: PetscCall(MatCreateVecs(lqn->YtS_triu, &lqn->diag_vec, &lqn->rwork1));
175: PetscCall(MatCreateVecs(lqn->YtS_triu, &lqn->rwork2, &lqn->rwork3));
176: } else if (is_dbfgs) {
177: PetscCall(MatCreateDenseFromVecType(comm, vec_type, m, m, M, M, -1, NULL, &lqn->StY_triu));
178: PetscCall(MatDuplicate(Sfull, MAT_SHARE_NONZERO_PATTERN, &lqn->BS));
179: PetscCall(MatCreateVecs(lqn->StY_triu, &lqn->diag_vec, &lqn->rwork1));
180: PetscCall(MatCreateVecs(lqn->StY_triu, &lqn->rwork2, &lqn->rwork3));
181: } else {
182: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "MatAllocate_LMVMDQN is only available for dense derived types. (DBFGS, DDFP, DQN");
183: }
184: /* initialize StY_triu and YtS_triu to identity, if they exist, so it is invertible */
185: if (lqn->StY_triu) {
186: PetscCall(MatZeroEntries(lqn->StY_triu));
187: PetscCall(MatShift(lqn->StY_triu, 1.0));
188: }
189: if (lqn->YtS_triu) {
190: PetscCall(MatZeroEntries(lqn->YtS_triu));
191: PetscCall(MatShift(lqn->YtS_triu, 1.0));
192: }
193: if (lqn->use_recursive && (is_dbfgs || is_ddfp)) {
194: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lqn->PQ));
195: PetscCall(VecDuplicate(lmvm->Xprev, &lqn->column_work2));
196: PetscCall(PetscMalloc1(lmvm->m, &lqn->yts));
197: if (is_dbfgs) PetscCall(PetscMalloc1(lmvm->m, &lqn->stp));
198: else if (is_ddfp) PetscCall(PetscMalloc1(lmvm->m, &lqn->ytq));
199: }
200: PetscCall(VecDuplicate(lqn->rwork2, &lqn->cyclic_work_vec));
201: PetscCall(VecZeroEntries(lqn->rwork1));
202: PetscCall(VecZeroEntries(lqn->rwork2));
203: PetscCall(VecZeroEntries(lqn->rwork3));
204: PetscCall(VecZeroEntries(lqn->diag_vec));
205: }
206: PetscCall(VecDuplicate(lmvm->Xprev, &lqn->column_work));
207: lqn->allocated = PETSC_TRUE;
208: }
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode MatSetUp_LMVMDQN(Mat B)
213: {
214: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
215: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
217: PetscFunctionBegin;
218: PetscCall(MatSetUp_LMVM(B));
219: PetscCall(SymBroydenRescaleInitializeJ0(B, lqn->rescale));
220: PetscCall(MatAllocate_LMVMDQN_Internal(B));
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: static PetscErrorCode MatSetFromOptions_LMVMDQN(Mat B, PetscOptionItems PetscOptionsObject)
225: {
226: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
227: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
228: PetscBool is_dbfgs, is_ddfp, is_dqn;
230: PetscFunctionBegin;
231: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
232: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
233: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
234: PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
235: PetscOptionsHeadBegin(PetscOptionsObject, "Dense symmetric Broyden method for approximating SPD Jacobian actions");
236: if (is_dqn) {
237: PetscCall(PetscOptionsEnum("-mat_lqn_type", "Implementation options for L-QN", "MatLMVMDenseType", MatLMVMDenseTypes, (PetscEnum)lqn->strategy, (PetscEnum *)&lqn->strategy, NULL));
238: } else if (is_dbfgs) {
239: PetscCall(PetscOptionsBool("-mat_lbfgs_recursive", "Use recursive formulation for MatMult_LMVMDBFGS, instead of Cholesky", "", lqn->use_recursive, &lqn->use_recursive, NULL));
240: PetscCall(PetscOptionsEnum("-mat_lbfgs_type", "Implementation options for L-BFGS", "MatLMVMDenseType", MatLMVMDenseTypes, (PetscEnum)lqn->strategy, (PetscEnum *)&lqn->strategy, NULL));
241: } else if (is_ddfp) {
242: PetscCall(PetscOptionsBool("-mat_ldfp_recursive", "Use recursive formulation for MatSolve_LMVMDDFP, instead of Cholesky", "", lqn->use_recursive, &lqn->use_recursive, NULL));
243: PetscCall(PetscOptionsEnum("-mat_ldfp_type", "Implementation options for L-DFP", "MatLMVMDenseType", MatLMVMDenseTypes, (PetscEnum)lqn->strategy, (PetscEnum *)&lqn->strategy, NULL));
244: } else {
245: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "MatSetFromOptions_LMVMDQN is only available for dense derived types. (DBFGS, DDFP, DQN");
246: }
247: PetscCall(SymBroydenRescaleSetFromOptions(B, lqn->rescale, PetscOptionsObject));
248: PetscOptionsHeadEnd();
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: static PetscErrorCode MatDestroy_LMVMDQN(Mat B)
253: {
254: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
255: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
257: PetscFunctionBegin;
258: PetscCall(SymBroydenRescaleDestroy(&lqn->rescale));
259: PetscCall(MatReset_LMVMDQN_Internal(B, MAT_LMVM_RESET_ALL));
260: PetscCall(PetscFree(lqn->workscalar));
261: PetscCall(PetscFree(lmvm->ctx));
262: PetscCall(MatDestroy_LMVM(B));
263: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", NULL));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode MatUpdate_LMVMDQN(Mat B, Vec X, Vec F)
268: {
269: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
270: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
271: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
272: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
274: PetscBool is_ddfp, is_dbfgs, is_dqn;
275: PetscDeviceContext dctx;
277: PetscFunctionBegin;
278: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
279: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
280: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
281: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
282: PetscCall(PetscDeviceContextGetCurrentContext(&dctx));
283: if (lmvm->prev_set) {
284: Vec FX[2];
285: PetscScalar dotFX[2];
286: PetscScalar stFprev;
287: PetscScalar curvature, yTy;
288: PetscReal curvtol;
290: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
291: PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
292: /* Test if the updates can be accepted */
293: FX[0] = lmvm->Fprev; /* dotFX[0] = s^T Fprev */
294: FX[1] = F; /* dotFX[1] = s^T F */
295: PetscCall(VecMDot(lmvm->Xprev, 2, FX, dotFX));
296: PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
297: PetscCall(VecDot(lmvm->Fprev, lmvm->Fprev, &yTy));
298: stFprev = PetscConj(dotFX[0]);
299: curvature = PetscConj(dotFX[1] - dotFX[0]); /* s^T y */
300: if (PetscRealPart(yTy) < lmvm->eps) {
301: curvtol = 0.0;
302: } else {
303: curvtol = lmvm->eps * PetscRealPart(yTy);
304: }
305: if (PetscRealPart(curvature) > curvtol) {
306: PetscInt m = lmvm->m;
307: PetscInt k = lmvm->k;
308: PetscInt h_old = k - oldest_update(m, k);
309: PetscInt h_new = k + 1 - oldest_update(m, k + 1);
310: PetscInt idx = recycle_index(m, k);
312: /* Update is good, accept it */
313: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
314: lqn->num_updates++;
315: lqn->watchdog = 0;
316: lqn->needPQ = PETSC_TRUE;
318: if (h_old == m && lqn->strategy == MAT_LMVM_DENSE_REORDER) {
319: if (is_dqn) {
320: PetscCall(MatMove_LR3(B, lqn->StY_triu, m - 1));
321: PetscCall(MatMove_LR3(B, lqn->YtS_triu, m - 1));
322: } else if (is_dbfgs) {
323: PetscCall(MatMove_LR3(B, lqn->StY_triu, m - 1));
324: } else if (is_ddfp) {
325: PetscCall(MatMove_LR3(B, lqn->YtS_triu, m - 1));
326: } else {
327: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "MatUpdate_LMVMDQN is only available for dense derived types. (DBFGS, DDFP, DQN");
328: }
329: }
331: if (lqn->use_recursive && (is_dbfgs || is_ddfp)) lqn->yts[idx] = PetscRealPart(curvature);
333: if (is_dqn || is_dbfgs) { /* implement the scheme of Byrd, Nocedal, and Schnabel to save a MatMultTranspose call in the common case the *
334: * H_k is immediately applied to F after begin updated. The S^T y computation can be split up as S^T (F - F_prev) */
335: PetscInt local_n;
336: PetscScalar *StFprev;
337: PetscMemType memtype;
338: PetscInt StYidx;
340: StYidx = (lqn->strategy == MAT_LMVM_DENSE_REORDER) ? history_index(m, lqn->num_updates, k) : idx;
341: if (!lqn->StFprev) PetscCall(VecDuplicate(lqn->rwork1, &lqn->StFprev));
342: PetscCall(VecGetLocalSize(lqn->StFprev, &local_n));
343: PetscCall(VecGetArrayAndMemType(lqn->StFprev, &StFprev, &memtype));
344: if (local_n) {
345: if (PetscMemTypeHost(memtype)) {
346: StFprev[idx] = stFprev;
347: } else {
348: PetscCall(PetscDeviceRegisterMemory(&stFprev, PETSC_MEMTYPE_HOST, 1 * sizeof(stFprev)));
349: PetscCall(PetscDeviceRegisterMemory(StFprev, memtype, local_n * sizeof(*StFprev)));
350: PetscCall(PetscDeviceArrayCopy(dctx, &StFprev[idx], &stFprev, 1));
351: PetscCall(PetscDeviceContextSynchronize(dctx));
352: }
353: }
354: PetscCall(VecRestoreArrayAndMemType(lqn->StFprev, &StFprev));
356: {
357: Vec this_sy_col;
358: /* Now StFprev is updated for the new S vector. Write -StFprev into the appropriate row */
359: PetscCall(MatDenseGetColumnVecWrite(lqn->StY_triu, StYidx, &this_sy_col));
360: PetscCall(VecAXPBY(this_sy_col, -1.0, 0.0, lqn->StFprev));
362: /* Now compute the new StFprev */
363: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, lqn->StFprev, 0, h_new));
364: lqn->St_count++;
366: /* Now add StFprev: this_sy_col == S^T (F - Fprev) == S^T y */
367: PetscCall(VecAXPY(this_sy_col, 1.0, lqn->StFprev));
369: if (lqn->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, this_sy_col, lqn->num_updates, lqn->cyclic_work_vec));
370: PetscCall(MatDenseRestoreColumnVecWrite(lqn->StY_triu, StYidx, &this_sy_col));
371: }
372: }
374: if (is_ddfp || is_dqn) {
375: PetscInt YtSidx;
377: YtSidx = (lqn->strategy == MAT_LMVM_DENSE_REORDER) ? history_index(m, lqn->num_updates, k) : idx;
379: {
380: Vec this_ys_col;
382: PetscCall(MatDenseGetColumnVecWrite(lqn->YtS_triu, YtSidx, &this_ys_col));
383: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, lmvm->Xprev, this_ys_col, 0, h_new));
384: lqn->Yt_count++;
386: if (lqn->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, this_ys_col, lqn->num_updates, lqn->cyclic_work_vec));
387: PetscCall(MatDenseRestoreColumnVecWrite(lqn->YtS_triu, YtSidx, &this_ys_col));
388: }
389: }
391: if (is_dbfgs || is_dqn) {
392: PetscCall(MatGetDiagonal(lqn->StY_triu, lqn->diag_vec));
393: } else if (is_ddfp) {
394: PetscCall(MatGetDiagonal(lqn->YtS_triu, lqn->diag_vec));
395: } else {
396: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "MatUpdate_LMVMDQN is only available for dense derived types. (DBFGS, DDFP, DQN");
397: }
399: if (lqn->strategy == MAT_LMVM_DENSE_REORDER) {
400: if (!lqn->diag_vec_recycle_order) PetscCall(VecDuplicate(lqn->diag_vec, &lqn->diag_vec_recycle_order));
401: PetscCall(VecCopy(lqn->diag_vec, lqn->diag_vec_recycle_order));
402: PetscCall(VecHistoryOrderToRecycleOrder(B, lqn->diag_vec_recycle_order, lqn->num_updates, lqn->cyclic_work_vec));
403: } else {
404: if (!lqn->diag_vec_recycle_order) {
405: PetscCall(PetscObjectReference((PetscObject)lqn->diag_vec));
406: lqn->diag_vec_recycle_order = lqn->diag_vec;
407: }
408: }
410: PetscCall(SymBroydenRescaleUpdate(B, lqn->rescale));
411: } else {
412: /* Update is bad, skip it */
413: ++lmvm->nrejects;
414: ++lqn->watchdog;
415: PetscInt m = lmvm->m;
416: PetscInt k = lmvm->k;
417: PetscInt h = k - oldest_update(m, k);
419: /* we still have to maintain StFprev */
420: if (!lqn->StFprev) PetscCall(VecDuplicate(lqn->rwork1, &lqn->StFprev));
421: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, lqn->StFprev, 0, h));
422: lqn->St_count++;
423: }
424: }
426: if (lqn->watchdog > lqn->max_seq_rejects) PetscCall(MatLMVMReset(B, PETSC_FALSE));
428: /* Save the solution and function to be used in the next update */
429: PetscCall(VecCopy(X, lmvm->Xprev));
430: PetscCall(VecCopy(F, lmvm->Fprev));
431: PetscCall(PetscObjectReference((PetscObject)F));
432: PetscCall(VecDestroy(&lqn->Fprev_ref));
433: lqn->Fprev_ref = F;
434: PetscCall(PetscObjectStateGet((PetscObject)F, &lqn->Fprev_state));
435: lmvm->prev_set = PETSC_TRUE;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: static PetscErrorCode MatDestroyThenCopy(Mat src, Mat *dst)
440: {
441: PetscFunctionBegin;
442: PetscCall(MatDestroy(dst));
443: if (src) PetscCall(MatDuplicate(src, MAT_COPY_VALUES, dst));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode VecDestroyThenCopy(Vec src, Vec *dst)
448: {
449: PetscFunctionBegin;
450: PetscCall(VecDestroy(dst));
451: if (src) {
452: PetscCall(VecDuplicate(src, dst));
453: PetscCall(VecCopy(src, *dst));
454: }
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: static PetscErrorCode MatCopy_LMVMDQN(Mat B, Mat M, MatStructure str)
459: {
460: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
461: Mat_DQN *blqn = (Mat_DQN *)bdata->ctx;
462: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
463: Mat_DQN *mlqn = (Mat_DQN *)mdata->ctx;
464: PetscBool is_dbfgs, is_ddfp, is_dqn;
466: PetscFunctionBegin;
467: PetscCall(SymBroydenRescaleCopy(blqn->rescale, mlqn->rescale));
468: mlqn->num_updates = blqn->num_updates;
469: mlqn->num_mult_updates = blqn->num_mult_updates;
470: mlqn->dense_type = blqn->dense_type;
471: mlqn->strategy = blqn->strategy;
472: mlqn->S_count = 0;
473: mlqn->St_count = 0;
474: mlqn->Y_count = 0;
475: mlqn->Yt_count = 0;
476: mlqn->watchdog = blqn->watchdog;
477: mlqn->max_seq_rejects = blqn->max_seq_rejects;
478: mlqn->use_recursive = blqn->use_recursive;
479: mlqn->needPQ = blqn->needPQ;
480: if (blqn->allocated) {
481: PetscCall(MatAllocate_LMVMDQN_Internal(M));
482: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
483: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
484: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
485: PetscCall(MatDestroyThenCopy(blqn->HY, &mlqn->BS));
486: PetscCall(VecDestroyThenCopy(blqn->StFprev, &mlqn->StFprev));
487: PetscCall(MatDestroyThenCopy(blqn->StY_triu, &mlqn->StY_triu));
488: PetscCall(MatDestroyThenCopy(blqn->StY_triu_strict, &mlqn->StY_triu_strict));
489: PetscCall(MatDestroyThenCopy(blqn->YtS_triu, &mlqn->YtS_triu));
490: PetscCall(MatDestroyThenCopy(blqn->YtS_triu_strict, &mlqn->YtS_triu_strict));
491: PetscCall(MatDestroyThenCopy(blqn->YtHY, &mlqn->YtHY));
492: PetscCall(MatDestroyThenCopy(blqn->StBS, &mlqn->StBS));
493: PetscCall(MatDestroyThenCopy(blqn->J, &mlqn->J));
494: PetscCall(VecDestroyThenCopy(blqn->diag_vec, &mlqn->diag_vec));
495: PetscCall(VecDestroyThenCopy(blqn->diag_vec_recycle_order, &mlqn->diag_vec_recycle_order));
496: PetscCall(VecDestroyThenCopy(blqn->inv_diag_vec, &mlqn->inv_diag_vec));
497: if (blqn->use_recursive && (is_dbfgs || is_ddfp)) {
498: for (PetscInt i = 0; i < bdata->m; i++) {
499: PetscCall(VecDestroyThenCopy(blqn->PQ[i], &mlqn->PQ[i]));
500: mlqn->yts[i] = blqn->yts[i];
501: if (is_dbfgs) {
502: mlqn->stp[i] = blqn->stp[i];
503: } else if (is_ddfp) {
504: mlqn->ytq[i] = blqn->ytq[i];
505: }
506: }
507: }
508: }
509: PetscCall(PetscObjectReference((PetscObject)blqn->Fprev_ref));
510: PetscCall(VecDestroy(&mlqn->Fprev_ref));
511: mlqn->Fprev_ref = blqn->Fprev_ref;
512: mlqn->Fprev_state = blqn->Fprev_state;
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: static PetscErrorCode MatMult_LMVMDQN(Mat B, Vec X, Vec Z)
517: {
518: PetscFunctionBegin;
519: PetscCall(MatMult_LMVMDDFP(B, X, Z));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode MatSolve_LMVMDQN(Mat H, Vec F, Vec dX)
524: {
525: PetscFunctionBegin;
526: PetscCall(MatSolve_LMVMDBFGS(H, F, dX));
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: static PetscErrorCode MatLMVMSymBroydenSetDelta_LMVMDQN(Mat B, PetscScalar delta)
531: {
532: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
533: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
535: PetscFunctionBegin;
536: PetscCall(SymBroydenRescaleSetDelta(B, lqn->rescale, PetscAbsReal(PetscRealPart(delta))));
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: /*
541: This dense representation uses Davidon-Fletcher-Powell (DFP) for MatMult,
542: and Broyden-Fletcher-Goldfarb-Shanno (BFGS) for MatSolve. This implementation
543: results in avoiding costly Cholesky factorization, at the cost of duality cap.
544: Please refer to MatLMVMDDFP and MatLMVMDBFGS for more information.
545: */
546: PetscErrorCode MatCreate_LMVMDQN(Mat B)
547: {
548: Mat_LMVM *lmvm;
549: Mat_DQN *lqn;
551: PetscFunctionBegin;
552: PetscCall(MatCreate_LMVM(B));
553: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDQN));
554: PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
555: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
556: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
557: B->ops->view = MatView_LMVMDQN;
558: B->ops->setup = MatSetUp_LMVMDQN;
559: B->ops->setfromoptions = MatSetFromOptions_LMVMDQN;
560: B->ops->destroy = MatDestroy_LMVMDQN;
562: lmvm = (Mat_LMVM *)B->data;
563: lmvm->ops->reset = MatReset_LMVMDQN;
564: lmvm->ops->update = MatUpdate_LMVMDQN;
565: lmvm->ops->mult = MatMult_LMVMDQN;
566: lmvm->ops->solve = MatSolve_LMVMDQN;
567: lmvm->ops->copy = MatCopy_LMVMDQN;
569: lmvm->ops->multht = lmvm->ops->mult;
570: lmvm->ops->solveht = lmvm->ops->solve;
572: PetscCall(PetscNew(&lqn));
573: lmvm->ctx = (void *)lqn;
574: lqn->allocated = PETSC_FALSE;
575: lqn->use_recursive = PETSC_FALSE;
576: lqn->needPQ = PETSC_FALSE;
577: lqn->watchdog = 0;
578: lqn->max_seq_rejects = lmvm->m / 2;
579: lqn->strategy = MAT_LMVM_DENSE_INPLACE;
581: PetscCall(SymBroydenRescaleCreate(&lqn->rescale));
582: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_LMVMDQN));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: MatCreateLMVMDQN - Creates a dense representation of the limited-memory
588: Quasi-Newton approximation to a Hessian.
590: Collective
592: Input Parameters:
593: + comm - MPI communicator
594: . n - number of local rows for storage vectors
595: - N - global size of the storage vectors
597: Output Parameter:
598: . B - the matrix
600: Level: advanced
602: Note:
603: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
604: paradigm instead of this routine directly.
606: .seealso: `MatCreate()`, `MATLMVM`, `MATLMVMDBFGS`, `MATLMVMDDFP`, `MatCreateLMVMDDFP()`, `MatCreateLMVMDBFGS()`
607: @*/
608: PetscErrorCode MatCreateLMVMDQN(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
609: {
610: PetscFunctionBegin;
611: PetscCall(KSPInitializePackage());
612: PetscCall(MatCreate(comm, B));
613: PetscCall(MatSetSizes(*B, n, n, N, N));
614: PetscCall(MatSetType(*B, MATLMVMDQN));
615: PetscCall(MatSetUp(*B));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: static PetscErrorCode MatDQNApplyJ0Fwd(Mat B, Vec X, Vec Z)
620: {
621: PetscFunctionBegin;
622: PetscCall(MatLMVMApplyJ0Fwd(B, X, Z));
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: static PetscErrorCode MatDQNApplyJ0Inv(Mat B, Vec F, Vec dX)
627: {
628: PetscFunctionBegin;
629: PetscCall(MatLMVMApplyJ0Inv(B, F, dX));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: /* This is not Bunch-Kaufman LDLT: here L is strictly lower triangular part of STY */
634: static PetscErrorCode MatGetLDLT(Mat B, Mat result)
635: {
636: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
637: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
638: PetscInt m_local;
640: PetscFunctionBegin;
641: if (!lbfgs->temp_mat) PetscCall(MatDuplicate(lbfgs->YtS_triu_strict, MAT_SHARE_NONZERO_PATTERN, &lbfgs->temp_mat));
642: PetscCall(MatCopy(lbfgs->YtS_triu_strict, lbfgs->temp_mat, SAME_NONZERO_PATTERN));
643: PetscCall(MatDiagonalScale(lbfgs->temp_mat, lbfgs->inv_diag_vec, NULL));
644: PetscCall(MatGetLocalSize(result, &m_local, NULL));
645: // need to conjugate and conjugate again because we have MatTransposeMatMult but not MatHermitianTransposeMatMult()
646: PetscCall(MatConjugate(lbfgs->temp_mat));
647: if (m_local) {
648: Mat temp_local, YtS_local, result_local;
649: PetscCall(MatDenseGetLocalMatrix(lbfgs->YtS_triu_strict, &YtS_local));
650: PetscCall(MatDenseGetLocalMatrix(lbfgs->temp_mat, &temp_local));
651: PetscCall(MatDenseGetLocalMatrix(result, &result_local));
652: PetscCall(MatTransposeMatMult(YtS_local, temp_local, MAT_REUSE_MATRIX, PETSC_DETERMINE, &result_local));
653: }
654: PetscCall(MatConjugate(result));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: static PetscErrorCode MatLMVMDBFGSUpdateMultData(Mat B)
659: {
660: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
661: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
662: PetscInt m = lmvm->m, m_local;
663: PetscInt k = lmvm->k;
664: PetscInt h = k - oldest_update(m, k);
665: PetscInt j_0;
666: PetscInt prev_oldest;
667: Mat J_local;
668: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
669: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
671: PetscFunctionBegin;
672: if (!lbfgs->YtS_triu_strict) {
673: PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->YtS_triu_strict));
674: PetscCall(MatDestroy(&lbfgs->StBS));
675: PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->StBS));
676: PetscCall(MatDestroy(&lbfgs->J));
677: PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->J));
678: PetscCall(MatDestroy(&lbfgs->BS));
679: PetscCall(MatDuplicate(Yfull, MAT_SHARE_NONZERO_PATTERN, &lbfgs->BS));
680: PetscCall(MatShift(lbfgs->StBS, 1.0));
681: lbfgs->num_mult_updates = oldest_update(m, k);
682: }
683: if (lbfgs->num_mult_updates == k) PetscFunctionReturn(PETSC_SUCCESS);
685: /* B_0 may have been updated, we must recompute B_0 S and S^T B_0 S */
686: for (PetscInt j = oldest_update(m, k); j < k; j++) {
687: Vec s_j;
688: Vec Bs_j;
689: Vec StBs_j;
690: PetscInt S_idx = recycle_index(m, j);
691: PetscInt StBS_idx = lbfgs->strategy == MAT_LMVM_DENSE_INPLACE ? S_idx : history_index(m, k, j);
693: PetscCall(MatDenseGetColumnVecWrite(lbfgs->BS, S_idx, &Bs_j));
694: PetscCall(MatDenseGetColumnVecRead(Sfull, S_idx, &s_j));
695: PetscCall(MatDQNApplyJ0Fwd(B, s_j, Bs_j));
696: PetscCall(MatDenseRestoreColumnVecRead(Sfull, S_idx, &s_j));
697: PetscCall(MatDenseGetColumnVecWrite(lbfgs->StBS, StBS_idx, &StBs_j));
698: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, Bs_j, StBs_j, 0, h));
699: lbfgs->St_count++;
700: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, StBs_j, lbfgs->num_updates, lbfgs->cyclic_work_vec));
701: PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->StBS, StBS_idx, &StBs_j));
702: PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->BS, S_idx, &Bs_j));
703: }
704: prev_oldest = oldest_update(m, lbfgs->num_mult_updates);
705: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER && prev_oldest < oldest_update(m, k)) {
706: /* move the YtS entries that have been computed and need to be kept back up */
707: PetscInt m_keep = m - (oldest_update(m, k) - prev_oldest);
709: PetscCall(MatMove_LR3(B, lbfgs->YtS_triu_strict, m_keep));
710: }
711: PetscCall(MatGetLocalSize(lbfgs->YtS_triu_strict, &m_local, NULL));
712: j_0 = PetscMax(lbfgs->num_mult_updates, oldest_update(m, k));
713: for (PetscInt j = j_0; j < k; j++) {
714: PetscInt S_idx = recycle_index(m, j);
715: PetscInt YtS_idx = lbfgs->strategy == MAT_LMVM_DENSE_INPLACE ? S_idx : history_index(m, k, j);
716: Vec s_j, Yts_j;
718: PetscCall(MatDenseGetColumnVecRead(Sfull, S_idx, &s_j));
719: PetscCall(MatDenseGetColumnVecWrite(lbfgs->YtS_triu_strict, YtS_idx, &Yts_j));
720: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, s_j, Yts_j, 0, h));
721: lbfgs->Yt_count++;
722: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, Yts_j, lbfgs->num_updates, lbfgs->cyclic_work_vec));
723: PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->YtS_triu_strict, YtS_idx, &Yts_j));
724: PetscCall(MatDenseRestoreColumnVecRead(Sfull, S_idx, &s_j));
725: /* zero the corresponding row */
726: if (m_local > 0) {
727: Mat YtS_local, YtS_row;
729: PetscCall(MatDenseGetLocalMatrix(lbfgs->YtS_triu_strict, &YtS_local));
730: PetscCall(MatDenseGetSubMatrix(YtS_local, YtS_idx, YtS_idx + 1, PETSC_DECIDE, PETSC_DECIDE, &YtS_row));
731: PetscCall(MatZeroEntries(YtS_row));
732: PetscCall(MatDenseRestoreSubMatrix(YtS_local, &YtS_row));
733: }
734: }
735: if (!lbfgs->inv_diag_vec) PetscCall(VecDuplicate(lbfgs->diag_vec, &lbfgs->inv_diag_vec));
736: PetscCall(VecCopy(lbfgs->diag_vec, lbfgs->inv_diag_vec));
737: PetscCall(VecReciprocal(lbfgs->inv_diag_vec));
738: PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
739: PetscCall(MatSetFactorType(J_local, MAT_FACTOR_NONE));
740: PetscCall(MatGetLDLT(B, lbfgs->J));
741: PetscCall(MatAXPY(lbfgs->J, 1.0, lbfgs->StBS, SAME_NONZERO_PATTERN));
742: if (m_local) {
743: PetscCall(MatSetOption(J_local, MAT_SPD, PETSC_TRUE));
744: PetscCall(MatCholeskyFactor(J_local, NULL, NULL));
745: }
746: lbfgs->num_mult_updates = lbfgs->num_updates;
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
750: /* Solves for
751: * [ I | -S R^{-T} ] [ I | 0 ] [ H_0 | 0 ] [ I | Y ] [ I ]
752: * [-----+---] [-----+---] [---+---] [-------------]
753: * [ Y^T | I ] [ 0 | D ] [ 0 | I ] [ -R^{-1} S^T ] */
755: static PetscErrorCode MatSolve_LMVMDBFGS(Mat H, Vec F, Vec dX)
756: {
757: Mat_LMVM *lmvm = (Mat_LMVM *)H->data;
758: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
759: Vec rwork1 = lbfgs->rwork1;
760: PetscInt m = lmvm->m;
761: PetscInt k = lmvm->k;
762: PetscInt h = k - oldest_update(m, k);
763: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
764: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
765: PetscObjectState Fstate;
767: PetscFunctionBegin;
768: VecCheckSameSize(F, 2, dX, 3);
769: VecCheckMatCompatible(H, dX, 3, F, 2);
771: /* Block Version */
772: if (!lbfgs->num_updates) {
773: PetscCall(MatDQNApplyJ0Inv(H, F, dX));
774: PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
775: }
777: PetscCall(PetscObjectStateGet((PetscObject)F, &Fstate));
778: if (F == lbfgs->Fprev_ref && Fstate == lbfgs->Fprev_state) {
779: PetscCall(VecCopy(lbfgs->StFprev, rwork1));
780: } else {
781: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, rwork1, 0, h));
782: lbfgs->St_count++;
783: }
785: /* Reordering rwork1, as STY is in history order, while S is in recycled order */
786: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
787: PetscCall(MatUpperTriangularSolveInPlace(H, lbfgs->StY_triu, rwork1, PETSC_FALSE, lbfgs->num_updates, lbfgs->strategy));
788: PetscCall(VecScale(rwork1, -1.0));
789: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
791: PetscCall(VecCopy(F, lbfgs->column_work));
792: PetscCall(MatMultAddColumnRange(Yfull, rwork1, lbfgs->column_work, lbfgs->column_work, 0, h));
793: lbfgs->Y_count++;
795: PetscCall(VecPointwiseMult(rwork1, lbfgs->diag_vec_recycle_order, rwork1));
796: PetscCall(MatDQNApplyJ0Inv(H, lbfgs->column_work, dX));
798: PetscCall(MatMultHermitianTransposeAddColumnRange(Yfull, dX, rwork1, rwork1, 0, h));
799: lbfgs->Yt_count++;
801: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
802: PetscCall(MatUpperTriangularSolveInPlace(H, lbfgs->StY_triu, rwork1, PETSC_TRUE, lbfgs->num_updates, lbfgs->strategy));
803: PetscCall(VecScale(rwork1, -1.0));
804: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
806: PetscCall(MatMultAddColumnRange(Sfull, rwork1, dX, dX, 0, h));
807: lbfgs->S_count++;
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: /* Solves for
812: B_0 - [ Y | B_0 S] [ -D | L^T ]^-1 [ Y^T ]
813: [-----+-----------] [---------]
814: [ L | S^T B_0 S ] [ S^T B_0 ]
816: Above is equivalent to
818: B_0 - [ Y | B_0 S] [[ I | 0 ][ -D | 0 ][ I | -D^{-1} L^T ]]^-1 [ Y^T ]
819: [[-----------+---][-----+---][---+-------------]] [---------]
820: [[ -L D^{-1} | I ][ 0 | J ][ 0 | I ]] [ S^T B_0 ]
822: where J = S^T B_0 S + L D^{-1} L^T
824: becomes
826: B_0 - [ Y | B_0 S] [ I | D^{-1} L^T ][ -D^{-1} | 0 ][ I | 0 ] [ Y^T ]
827: [---+------------][----------+--------][----------+---] [---------]
828: [ 0 | I ][ 0 | J^{-1} ][ L D^{-1} | I ] [ S^T B_0 ]
830: =
832: B_0 + [ Y | B_0 S] [ D^{-1} | 0 ][ I | L^T ][ I | 0 ][ I | 0 ] [ Y^T ]
833: [--------+---][---+-----][---+---------][----------+---] [---------]
834: [ 0 | I ][ 0 | I ][ 0 | -J^{-1} ][ L D^{-1} | I ] [ S^T B_0 ]
836: (Note that YtS_triu_strict is L^T)
837: Byrd, Nocedal, Schnabel 1994
839: Alternative approach: considering the fact that DFP is dual to BFGS, use MatMult of DPF:
840: (See ddfp.c's MatMult_LMVMDDFP)
842: */
843: static PetscErrorCode MatMult_LMVMDBFGS(Mat B, Vec X, Vec Z)
844: {
845: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
846: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
847: Mat J_local;
848: PetscInt idx, i, j, m_local, local_n;
849: PetscInt m = lmvm->m;
850: PetscInt k = lmvm->k;
851: PetscInt h = k - oldest_update(m, k);
852: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
853: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
855: PetscFunctionBegin;
856: VecCheckSameSize(X, 2, Z, 3);
857: VecCheckMatCompatible(B, X, 2, Z, 3);
859: /* Cholesky Version */
860: /* Start with the B0 term */
861: PetscCall(MatDQNApplyJ0Fwd(B, X, Z));
862: if (!lbfgs->num_updates) PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
864: if (lbfgs->use_recursive) {
865: PetscDeviceContext dctx;
866: PetscMemType memtype;
867: PetscScalar stz, ytx, stp, sjtpi, yjtsi, *workscalar;
868: PetscInt oldest = oldest_update(m, k);
870: PetscCall(PetscDeviceContextGetCurrentContext(&dctx));
871: /* Recursive formulation to avoid Cholesky. Not a dense formulation */
872: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, lbfgs->rwork1, 0, h));
873: lbfgs->Yt_count++;
875: PetscCall(VecGetLocalSize(lbfgs->rwork1, &local_n));
877: if (lbfgs->needPQ) {
878: PetscInt oldest = oldest_update(m, k);
879: for (i = oldest; i < k; ++i) {
880: idx = recycle_index(m, i);
881: /* column_work = S[idx] */
882: PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work, idx));
883: PetscCall(MatDQNApplyJ0Fwd(B, lbfgs->column_work, lbfgs->PQ[idx]));
884: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, lbfgs->column_work, lbfgs->rwork3, 0, h));
885: PetscCall(VecGetArrayAndMemType(lbfgs->rwork3, &workscalar, &memtype));
886: for (j = oldest; j < i; ++j) {
887: PetscInt idx_j = recycle_index(m, j);
888: /* Copy yjtsi in device-aware manner */
889: if (local_n) {
890: if (PetscMemTypeHost(memtype)) {
891: yjtsi = workscalar[idx_j];
892: } else {
893: PetscCall(PetscDeviceRegisterMemory(&yjtsi, PETSC_MEMTYPE_HOST, sizeof(yjtsi)));
894: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
895: PetscCall(PetscDeviceArrayCopy(dctx, &yjtsi, &workscalar[idx_j], 1));
896: PetscCall(PetscDeviceContextSynchronize(dctx));
897: }
898: }
899: PetscCallMPI(MPI_Bcast(&yjtsi, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)B)));
900: /* column_work2 = S[j] */
901: PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work2, idx_j));
902: PetscCall(VecDot(lbfgs->PQ[idx], lbfgs->column_work2, &sjtpi));
903: /* column_work2 = Y[j] */
904: PetscCall(MatGetColumnVector(Yfull, lbfgs->column_work2, idx_j));
905: /* Compute the pure BFGS component of the forward product */
906: PetscCall(VecAXPBYPCZ(lbfgs->PQ[idx], -sjtpi / lbfgs->stp[idx_j], yjtsi / lbfgs->yts[idx_j], 1.0, lbfgs->PQ[idx_j], lbfgs->column_work2));
907: }
908: PetscCall(VecRestoreArrayAndMemType(lbfgs->rwork3, &workscalar));
909: PetscCall(VecDot(lbfgs->PQ[idx], lbfgs->column_work, &stp));
910: lbfgs->stp[idx] = PetscRealPart(stp);
911: }
912: lbfgs->needPQ = PETSC_FALSE;
913: }
915: PetscCall(VecGetArrayAndMemType(lbfgs->rwork1, &workscalar, &memtype));
916: for (i = oldest; i < k; ++i) {
917: idx = recycle_index(m, i);
918: /* Copy stz[i], ytx[i] in device-aware manner */
919: if (local_n) {
920: if (PetscMemTypeHost(memtype)) {
921: ytx = workscalar[idx];
922: } else {
923: PetscCall(PetscDeviceRegisterMemory(&ytx, PETSC_MEMTYPE_HOST, 1 * sizeof(ytx)));
924: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
925: PetscCall(PetscDeviceArrayCopy(dctx, &ytx, &workscalar[idx], 1));
926: PetscCall(PetscDeviceContextSynchronize(dctx));
927: }
928: }
929: PetscCallMPI(MPI_Bcast(&ytx, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)B)));
930: /* column_work : S[i], column_work2 : Y[i] */
931: PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work, idx));
932: PetscCall(MatGetColumnVector(Yfull, lbfgs->column_work2, idx));
933: PetscCall(VecDot(Z, lbfgs->column_work, &stz));
934: PetscCall(VecAXPBYPCZ(Z, -stz / lbfgs->stp[idx], ytx / lbfgs->yts[idx], 1.0, lbfgs->PQ[idx], lbfgs->column_work2));
935: }
936: PetscCall(VecRestoreArrayAndMemType(lbfgs->rwork1, &workscalar));
937: } else {
938: PetscCall(MatLMVMDBFGSUpdateMultData(B));
939: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, lbfgs->rwork1, 0, h));
940: lbfgs->Yt_count++;
941: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, Z, lbfgs->rwork2, 0, h));
942: lbfgs->St_count++;
943: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) {
944: PetscCall(VecRecycleOrderToHistoryOrder(B, lbfgs->rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
945: PetscCall(VecRecycleOrderToHistoryOrder(B, lbfgs->rwork2, lbfgs->num_updates, lbfgs->cyclic_work_vec));
946: }
948: PetscCall(VecPointwiseMult(lbfgs->rwork1, lbfgs->rwork1, lbfgs->inv_diag_vec));
949: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(lbfgs->YtS_triu_strict));
950: PetscCall(MatMultTransposeAdd(lbfgs->YtS_triu_strict, lbfgs->rwork1, lbfgs->rwork2, lbfgs->rwork2));
951: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(lbfgs->YtS_triu_strict));
953: if (!lbfgs->rwork2_local) PetscCall(VecCreateLocalVector(lbfgs->rwork2, &lbfgs->rwork2_local));
954: if (!lbfgs->rwork3_local) PetscCall(VecCreateLocalVector(lbfgs->rwork3, &lbfgs->rwork3_local));
955: PetscCall(VecGetLocalVectorRead(lbfgs->rwork2, lbfgs->rwork2_local));
956: PetscCall(VecGetLocalVector(lbfgs->rwork3, lbfgs->rwork3_local));
957: PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
958: PetscCall(VecGetSize(lbfgs->rwork2_local, &m_local));
959: if (m_local) {
960: PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
961: PetscCall(MatSolve(J_local, lbfgs->rwork2_local, lbfgs->rwork3_local));
962: }
963: PetscCall(VecRestoreLocalVector(lbfgs->rwork3, lbfgs->rwork3_local));
964: PetscCall(VecRestoreLocalVectorRead(lbfgs->rwork2, lbfgs->rwork2_local));
965: PetscCall(VecScale(lbfgs->rwork3, -1.0));
967: PetscCall(MatMult(lbfgs->YtS_triu_strict, lbfgs->rwork3, lbfgs->rwork2));
968: PetscCall(VecPointwiseMult(lbfgs->rwork2, lbfgs->rwork2, lbfgs->inv_diag_vec));
969: PetscCall(VecAXPY(lbfgs->rwork1, 1.0, lbfgs->rwork2));
971: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) {
972: PetscCall(VecHistoryOrderToRecycleOrder(B, lbfgs->rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
973: PetscCall(VecHistoryOrderToRecycleOrder(B, lbfgs->rwork3, lbfgs->num_updates, lbfgs->cyclic_work_vec));
974: }
976: PetscCall(MatMultAddColumnRange(Yfull, lbfgs->rwork1, Z, Z, 0, h));
977: lbfgs->Y_count++;
978: PetscCall(MatMultAddColumnRange(lbfgs->BS, lbfgs->rwork3, Z, Z, 0, h));
979: lbfgs->S_count++;
980: }
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: /*
985: This dense representation reduces the L-BFGS update to a series of
986: matrix-vector products with dense matrices in lieu of the conventional matrix-free
987: two-loop algorithm.
988: */
989: PetscErrorCode MatCreate_LMVMDBFGS(Mat B)
990: {
991: Mat_LMVM *lmvm;
992: Mat_DQN *lbfgs;
994: PetscFunctionBegin;
995: PetscCall(MatCreate_LMVM(B));
996: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDBFGS));
997: PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
998: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
999: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
1000: B->ops->view = MatView_LMVMDQN;
1001: B->ops->setup = MatSetUp_LMVMDQN;
1002: B->ops->setfromoptions = MatSetFromOptions_LMVMDQN;
1003: B->ops->destroy = MatDestroy_LMVMDQN;
1005: lmvm = (Mat_LMVM *)B->data;
1006: lmvm->ops->reset = MatReset_LMVMDQN;
1007: lmvm->ops->update = MatUpdate_LMVMDQN;
1008: lmvm->ops->mult = MatMult_LMVMDBFGS;
1009: lmvm->ops->solve = MatSolve_LMVMDBFGS;
1010: lmvm->ops->copy = MatCopy_LMVMDQN;
1012: lmvm->ops->multht = lmvm->ops->mult;
1013: lmvm->ops->solveht = lmvm->ops->solve;
1015: PetscCall(PetscNew(&lbfgs));
1016: lmvm->ctx = (void *)lbfgs;
1017: lbfgs->allocated = PETSC_FALSE;
1018: lbfgs->use_recursive = PETSC_TRUE;
1019: lbfgs->needPQ = PETSC_TRUE;
1020: lbfgs->watchdog = 0;
1021: lbfgs->max_seq_rejects = lmvm->m / 2;
1022: lbfgs->strategy = MAT_LMVM_DENSE_INPLACE;
1024: PetscCall(SymBroydenRescaleCreate(&lbfgs->rescale));
1025: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_LMVMDQN));
1026: PetscFunctionReturn(PETSC_SUCCESS);
1027: }
1029: /*@
1030: MatCreateLMVMDBFGS - Creates a dense representation of the limited-memory
1031: Broyden-Fletcher-Goldfarb-Shanno (BFGS) approximation to a Hessian.
1033: Collective
1035: Input Parameters:
1036: + comm - MPI communicator
1037: . n - number of local rows for storage vectors
1038: - N - global size of the storage vectors
1040: Output Parameter:
1041: . B - the matrix
1043: Level: advanced
1045: Note:
1046: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
1047: paradigm instead of this routine directly.
1049: .seealso: `MatCreate()`, `MATLMVM`, `MATLMVMDBFGS`, `MatCreateLMVMBFGS()`
1050: @*/
1051: PetscErrorCode MatCreateLMVMDBFGS(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
1052: {
1053: PetscFunctionBegin;
1054: PetscCall(KSPInitializePackage());
1055: PetscCall(MatCreate(comm, B));
1056: PetscCall(MatSetSizes(*B, n, n, N, N));
1057: PetscCall(MatSetType(*B, MATLMVMDBFGS));
1058: PetscCall(MatSetUp(*B));
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: /* here R is strictly upper triangular part of STY */
1063: static PetscErrorCode MatGetRTDR(Mat B, Mat result)
1064: {
1065: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1066: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1067: PetscInt m_local;
1069: PetscFunctionBegin;
1070: if (!ldfp->temp_mat) PetscCall(MatDuplicate(ldfp->StY_triu_strict, MAT_SHARE_NONZERO_PATTERN, &ldfp->temp_mat));
1071: PetscCall(MatCopy(ldfp->StY_triu_strict, ldfp->temp_mat, SAME_NONZERO_PATTERN));
1072: PetscCall(MatDiagonalScale(ldfp->temp_mat, ldfp->inv_diag_vec, NULL));
1073: PetscCall(MatGetLocalSize(result, &m_local, NULL));
1074: // need to conjugate and conjugate again because we have MatTransposeMatMult but not MatHermitianTransposeMatMult()
1075: PetscCall(MatConjugate(ldfp->temp_mat));
1076: if (m_local) {
1077: Mat temp_local, StY_local, result_local;
1078: PetscCall(MatDenseGetLocalMatrix(ldfp->StY_triu_strict, &StY_local));
1079: PetscCall(MatDenseGetLocalMatrix(ldfp->temp_mat, &temp_local));
1080: PetscCall(MatDenseGetLocalMatrix(result, &result_local));
1081: PetscCall(MatTransposeMatMult(StY_local, temp_local, MAT_REUSE_MATRIX, PETSC_DETERMINE, &result_local));
1082: }
1083: PetscCall(MatConjugate(result));
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: static PetscErrorCode MatLMVMDDFPUpdateSolveData(Mat B)
1088: {
1089: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1090: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1091: PetscInt m = lmvm->m, m_local;
1092: PetscInt k = lmvm->k;
1093: PetscInt h = k - oldest_update(m, k);
1094: PetscInt j_0;
1095: PetscInt prev_oldest;
1096: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
1097: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
1098: Mat J_local;
1100: PetscFunctionBegin;
1101: if (!ldfp->StY_triu_strict) {
1102: PetscCall(MatDuplicate(ldfp->YtS_triu, MAT_SHARE_NONZERO_PATTERN, &ldfp->StY_triu_strict));
1103: PetscCall(MatDestroy(&ldfp->YtHY));
1104: PetscCall(MatDuplicate(ldfp->YtS_triu, MAT_SHARE_NONZERO_PATTERN, &ldfp->YtHY));
1105: PetscCall(MatDestroy(&ldfp->J));
1106: PetscCall(MatDuplicate(ldfp->YtS_triu, MAT_SHARE_NONZERO_PATTERN, &ldfp->J));
1107: PetscCall(MatDestroy(&ldfp->HY));
1108: PetscCall(MatDuplicate(Yfull, MAT_SHARE_NONZERO_PATTERN, &ldfp->HY));
1109: PetscCall(MatShift(ldfp->YtHY, 1.0));
1110: ldfp->num_mult_updates = oldest_update(m, k);
1111: }
1112: if (ldfp->num_mult_updates == k) PetscFunctionReturn(PETSC_SUCCESS);
1114: /* H_0 may have been updated, we must recompute H_0 Y and Y^T H_0 Y */
1115: for (PetscInt j = oldest_update(m, k); j < k; j++) {
1116: Vec y_j;
1117: Vec Hy_j;
1118: Vec YtHy_j;
1119: PetscInt Y_idx = recycle_index(m, j);
1120: PetscInt YtHY_idx = ldfp->strategy == MAT_LMVM_DENSE_INPLACE ? Y_idx : history_index(m, k, j);
1122: PetscCall(MatDenseGetColumnVecWrite(ldfp->HY, Y_idx, &Hy_j));
1123: PetscCall(MatDenseGetColumnVecRead(Yfull, Y_idx, &y_j));
1124: PetscCall(MatDQNApplyJ0Inv(B, y_j, Hy_j));
1125: PetscCall(MatDenseRestoreColumnVecRead(Yfull, Y_idx, &y_j));
1126: PetscCall(MatDenseGetColumnVecWrite(ldfp->YtHY, YtHY_idx, &YtHy_j));
1127: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, Hy_j, YtHy_j, 0, h));
1128: ldfp->Yt_count++;
1129: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, YtHy_j, ldfp->num_updates, ldfp->cyclic_work_vec));
1130: PetscCall(MatDenseRestoreColumnVecWrite(ldfp->YtHY, YtHY_idx, &YtHy_j));
1131: PetscCall(MatDenseRestoreColumnVecWrite(ldfp->HY, Y_idx, &Hy_j));
1132: }
1133: prev_oldest = oldest_update(m, ldfp->num_mult_updates);
1134: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER && prev_oldest < oldest_update(m, k)) {
1135: /* move the YtS entries that have been computed and need to be kept back up */
1136: PetscInt m_keep = m - (oldest_update(m, k) - prev_oldest);
1138: PetscCall(MatMove_LR3(B, ldfp->StY_triu_strict, m_keep));
1139: }
1140: PetscCall(MatGetLocalSize(ldfp->StY_triu_strict, &m_local, NULL));
1141: j_0 = PetscMax(ldfp->num_mult_updates, oldest_update(m, k));
1142: for (PetscInt j = j_0; j < k; j++) {
1143: PetscInt Y_idx = recycle_index(m, j);
1144: PetscInt StY_idx = ldfp->strategy == MAT_LMVM_DENSE_INPLACE ? Y_idx : history_index(m, k, j);
1145: Vec y_j, Sty_j;
1147: PetscCall(MatDenseGetColumnVecRead(Yfull, Y_idx, &y_j));
1148: PetscCall(MatDenseGetColumnVecWrite(ldfp->StY_triu_strict, StY_idx, &Sty_j));
1149: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, y_j, Sty_j, 0, h));
1150: ldfp->St_count++;
1151: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, Sty_j, ldfp->num_updates, ldfp->cyclic_work_vec));
1152: PetscCall(MatDenseRestoreColumnVecWrite(ldfp->StY_triu_strict, StY_idx, &Sty_j));
1153: PetscCall(MatDenseRestoreColumnVecRead(Yfull, Y_idx, &y_j));
1154: /* zero the corresponding row */
1155: if (m_local > 0) {
1156: Mat StY_local, StY_row;
1158: PetscCall(MatDenseGetLocalMatrix(ldfp->StY_triu_strict, &StY_local));
1159: PetscCall(MatDenseGetSubMatrix(StY_local, StY_idx, StY_idx + 1, PETSC_DECIDE, PETSC_DECIDE, &StY_row));
1160: PetscCall(MatZeroEntries(StY_row));
1161: PetscCall(MatDenseRestoreSubMatrix(StY_local, &StY_row));
1162: }
1163: }
1164: if (!ldfp->inv_diag_vec) PetscCall(VecDuplicate(ldfp->diag_vec, &ldfp->inv_diag_vec));
1165: PetscCall(VecCopy(ldfp->diag_vec, ldfp->inv_diag_vec));
1166: PetscCall(VecReciprocal(ldfp->inv_diag_vec));
1167: PetscCall(MatDenseGetLocalMatrix(ldfp->J, &J_local));
1168: PetscCall(MatSetFactorType(J_local, MAT_FACTOR_NONE));
1169: PetscCall(MatGetRTDR(B, ldfp->J));
1170: PetscCall(MatAXPY(ldfp->J, 1.0, ldfp->YtHY, SAME_NONZERO_PATTERN));
1171: if (m_local) {
1172: PetscCall(MatSetOption(J_local, MAT_SPD, PETSC_TRUE));
1173: PetscCall(MatCholeskyFactor(J_local, NULL, NULL));
1174: }
1175: ldfp->num_mult_updates = ldfp->num_updates;
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: /* Solves for
1181: H_0 - [ S | H_0 Y] [ -D | R.T ]^-1 [ S^T ]
1182: [-----+-----------] [---------]
1183: [ R | Y^T H_0 Y ] [ Y^T H_0 ]
1185: Above is equivalent to
1187: H_0 - [ S | H_0 Y] [[ I | 0 ][ -D | 0 ][ I | -D^{-1} R^T ]]^-1 [ S^T ]
1188: [[-----------+---][----+---][---+-------------]] [---------]
1189: [[ -R D^{-1} | I ][ 0 | J ][ 0 | I ]] [ Y^T H_0 ]
1191: where J = Y^T H_0 Y + R D^{-1} R.T
1193: becomes
1195: H_0 - [ S | H_0 Y] [ I | D^{-1} R^T ][ -D^{-1} | 0 ][ I | 0 ] [ S^T ]
1196: [---+------------][----------+--------][----------+---] [---------]
1197: [ 0 | I ][ 0 | J^{-1} ][ R D^{-1} | I ] [ Y^T H_0 ]
1199: =
1201: H_0 + [ S | H_0 Y] [ D^{-1} | 0 ][ I | R^T ][ I | 0 ][ I | 0 ] [ S^T ]
1202: [--------+---][---+-----][---+---------][----------+---] [---------]
1203: [ 0 | I ][ 0 | I ][ 0 | -J^{-1} ][ R D^{-1} | I ] [ Y^T H_0 ]
1205: (Note that StY_triu_strict is R)
1206: Byrd, Nocedal, Schnabel 1994
1208: */
1209: static PetscErrorCode MatSolve_LMVMDDFP(Mat H, Vec F, Vec dX)
1210: {
1211: Mat_LMVM *lmvm = (Mat_LMVM *)H->data;
1212: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1213: PetscInt m = lmvm->m;
1214: PetscInt k = lmvm->k;
1215: PetscInt h = k - oldest_update(m, k);
1216: PetscInt idx, i, j, local_n;
1217: PetscInt m_local;
1218: Mat J_local;
1219: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
1220: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
1222: PetscFunctionBegin;
1223: VecCheckSameSize(F, 2, dX, 3);
1224: VecCheckMatCompatible(H, dX, 3, F, 2);
1226: /* Cholesky Version */
1227: /* Start with the B0 term */
1228: PetscCall(MatDQNApplyJ0Inv(H, F, dX));
1229: if (!ldfp->num_updates) PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
1231: if (ldfp->use_recursive) {
1232: PetscDeviceContext dctx;
1233: PetscMemType memtype;
1234: PetscScalar stf, ytx, ytq, yjtqi, sjtyi, *workscalar;
1236: PetscCall(PetscDeviceContextGetCurrentContext(&dctx));
1237: /* Recursive formulation to avoid Cholesky. Not a dense formulation */
1238: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, ldfp->rwork1, 0, h));
1239: ldfp->Yt_count++;
1241: PetscCall(VecGetLocalSize(ldfp->rwork1, &local_n));
1243: PetscInt oldest = oldest_update(m, k);
1245: if (ldfp->needPQ) {
1246: PetscInt oldest = oldest_update(m, k);
1247: for (i = oldest; i < k; ++i) {
1248: idx = recycle_index(m, i);
1249: /* column_work = S[idx] */
1250: PetscCall(MatGetColumnVector(Yfull, ldfp->column_work, idx));
1251: PetscCall(MatDQNApplyJ0Inv(H, ldfp->column_work, ldfp->PQ[idx]));
1252: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, ldfp->column_work, ldfp->rwork3, 0, h));
1253: PetscCall(VecGetArrayAndMemType(ldfp->rwork3, &workscalar, &memtype));
1254: for (j = oldest; j < i; ++j) {
1255: PetscInt idx_j = recycle_index(m, j);
1256: /* Copy sjtyi in device-aware manner */
1257: if (local_n) {
1258: if (PetscMemTypeHost(memtype)) {
1259: sjtyi = workscalar[idx_j];
1260: } else {
1261: PetscCall(PetscDeviceRegisterMemory(&sjtyi, PETSC_MEMTYPE_HOST, 1 * sizeof(sjtyi)));
1262: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
1263: PetscCall(PetscDeviceArrayCopy(dctx, &sjtyi, &workscalar[idx_j], 1));
1264: PetscCall(PetscDeviceContextSynchronize(dctx));
1265: }
1266: }
1267: PetscCallMPI(MPI_Bcast(&sjtyi, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)H)));
1268: /* column_work2 = Y[j] */
1269: PetscCall(MatGetColumnVector(Yfull, ldfp->column_work2, idx_j));
1270: PetscCall(VecDot(ldfp->PQ[idx], ldfp->column_work2, &yjtqi));
1271: /* column_work2 = Y[j] */
1272: PetscCall(MatGetColumnVector(Sfull, ldfp->column_work2, idx_j));
1273: /* Compute the pure BFGS component of the forward product */
1274: PetscCall(VecAXPBYPCZ(ldfp->PQ[idx], -yjtqi / ldfp->ytq[idx_j], sjtyi / ldfp->yts[idx_j], 1.0, ldfp->PQ[idx_j], ldfp->column_work2));
1275: }
1276: PetscCall(VecRestoreArrayAndMemType(ldfp->rwork3, &workscalar));
1277: PetscCall(VecDot(ldfp->PQ[idx], ldfp->column_work, &ytq));
1278: ldfp->ytq[idx] = PetscRealPart(ytq);
1279: }
1280: ldfp->needPQ = PETSC_FALSE;
1281: }
1283: PetscCall(VecGetArrayAndMemType(ldfp->rwork1, &workscalar, &memtype));
1284: for (i = oldest; i < k; ++i) {
1285: idx = recycle_index(m, i);
1286: /* Copy stz[i], ytx[i] in device-aware manner */
1287: if (local_n) {
1288: if (PetscMemTypeHost(memtype)) {
1289: stf = workscalar[idx];
1290: } else {
1291: PetscCall(PetscDeviceRegisterMemory(&stf, PETSC_MEMTYPE_HOST, sizeof(stf)));
1292: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
1293: PetscCall(PetscDeviceArrayCopy(dctx, &stf, &workscalar[idx], 1));
1294: PetscCall(PetscDeviceContextSynchronize(dctx));
1295: }
1296: }
1297: PetscCallMPI(MPI_Bcast(&stf, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)H)));
1298: /* column_work : S[i], column_work2 : Y[i] */
1299: PetscCall(MatGetColumnVector(Sfull, ldfp->column_work, idx));
1300: PetscCall(MatGetColumnVector(Yfull, ldfp->column_work2, idx));
1301: PetscCall(VecDot(dX, ldfp->column_work2, &ytx));
1302: PetscCall(VecAXPBYPCZ(dX, -ytx / ldfp->ytq[idx], stf / ldfp->yts[idx], 1.0, ldfp->PQ[idx], ldfp->column_work));
1303: }
1304: PetscCall(VecRestoreArrayAndMemType(ldfp->rwork1, &workscalar));
1305: } else {
1306: PetscCall(MatLMVMDDFPUpdateSolveData(H));
1307: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, ldfp->rwork1, 0, h));
1308: ldfp->St_count++;
1309: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, dX, ldfp->rwork2, 0, h));
1310: ldfp->Yt_count++;
1311: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) {
1312: PetscCall(VecRecycleOrderToHistoryOrder(H, ldfp->rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1313: PetscCall(VecRecycleOrderToHistoryOrder(H, ldfp->rwork2, ldfp->num_updates, ldfp->cyclic_work_vec));
1314: }
1316: PetscCall(VecPointwiseMult(ldfp->rwork3, ldfp->rwork1, ldfp->inv_diag_vec));
1317: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(ldfp->StY_triu_strict));
1318: PetscCall(MatMultTransposeAdd(ldfp->StY_triu_strict, ldfp->rwork3, ldfp->rwork2, ldfp->rwork2));
1319: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(ldfp->StY_triu_strict));
1321: if (!ldfp->rwork2_local) PetscCall(VecCreateLocalVector(ldfp->rwork2, &ldfp->rwork2_local));
1322: if (!ldfp->rwork3_local) PetscCall(VecCreateLocalVector(ldfp->rwork3, &ldfp->rwork3_local));
1323: PetscCall(VecGetLocalVectorRead(ldfp->rwork2, ldfp->rwork2_local));
1324: PetscCall(VecGetLocalVector(ldfp->rwork3, ldfp->rwork3_local));
1325: PetscCall(MatDenseGetLocalMatrix(ldfp->J, &J_local));
1326: PetscCall(VecGetSize(ldfp->rwork2_local, &m_local));
1327: if (m_local) {
1328: Mat J_local;
1330: PetscCall(MatDenseGetLocalMatrix(ldfp->J, &J_local));
1331: PetscCall(MatSolve(J_local, ldfp->rwork2_local, ldfp->rwork3_local));
1332: }
1333: PetscCall(VecRestoreLocalVector(ldfp->rwork3, ldfp->rwork3_local));
1334: PetscCall(VecRestoreLocalVectorRead(ldfp->rwork2, ldfp->rwork2_local));
1335: PetscCall(VecScale(ldfp->rwork3, -1.0));
1337: PetscCall(MatMultAdd(ldfp->StY_triu_strict, ldfp->rwork3, ldfp->rwork1, ldfp->rwork1));
1338: PetscCall(VecPointwiseMult(ldfp->rwork1, ldfp->rwork1, ldfp->inv_diag_vec));
1340: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) {
1341: PetscCall(VecHistoryOrderToRecycleOrder(H, ldfp->rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1342: PetscCall(VecHistoryOrderToRecycleOrder(H, ldfp->rwork3, ldfp->num_updates, ldfp->cyclic_work_vec));
1343: }
1345: PetscCall(MatMultAddColumnRange(Sfull, ldfp->rwork1, dX, dX, 0, h));
1346: ldfp->S_count++;
1347: PetscCall(MatMultAddColumnRange(ldfp->HY, ldfp->rwork3, dX, dX, 0, h));
1348: ldfp->Y_count++;
1349: }
1350: PetscFunctionReturn(PETSC_SUCCESS);
1351: }
1353: /* Solves for
1354: (Theorem 1, Erway, Jain, and Marcia, 2013)
1356: B_0 - [ Y | B_0 S] [ -R^{-T} (D + S^T B_0 S) R^{-1} | R^{-T} ] [ Y^T ]
1357: ---------------------------------+--------] [---------]
1358: [ R^{-1} | 0 ] [ S^T B_0 ]
1360: (Note: R above is right triangular part of YTS)
1361: which becomes,
1363: [ I | -Y L^{-T} ] [ I | 0 ] [ B_0 | 0 ] [ I | S ] [ I ]
1364: [-----+---] [-----+---] [---+---] [-------------]
1365: [ S^T | I ] [ 0 | D ] [ 0 | I ] [ -L^{-1} Y^T ]
1367: (Note: L above is right triangular part of STY)
1369: */
1370: static PetscErrorCode MatMult_LMVMDDFP(Mat B, Vec X, Vec Z)
1371: {
1372: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1373: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1374: Vec rwork1 = ldfp->rwork1;
1375: PetscInt m = lmvm->m;
1376: PetscInt k = lmvm->k;
1377: PetscInt h = k - oldest_update(m, k);
1378: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
1379: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
1380: PetscObjectState Xstate;
1382: PetscFunctionBegin;
1383: VecCheckSameSize(X, 2, Z, 3);
1384: VecCheckMatCompatible(B, X, 2, Z, 3);
1386: /* DFP Version. Erway, Jain, Marcia, 2013, Theorem 1 */
1387: /* Block Version */
1388: if (!ldfp->num_updates) {
1389: PetscCall(MatDQNApplyJ0Fwd(B, X, Z));
1390: PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
1391: }
1393: PetscCall(PetscObjectStateGet((PetscObject)X, &Xstate));
1394: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, rwork1, 0, h));
1396: /* Reordering rwork1, as STY is in history order, while Y is in recycled order */
1397: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1398: PetscCall(MatUpperTriangularSolveInPlace(B, ldfp->YtS_triu, rwork1, PETSC_FALSE, ldfp->num_updates, ldfp->strategy));
1399: PetscCall(VecScale(rwork1, -1.0));
1400: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1402: PetscCall(VecCopy(X, ldfp->column_work));
1403: PetscCall(MatMultAddColumnRange(Sfull, rwork1, ldfp->column_work, ldfp->column_work, 0, h));
1404: ldfp->S_count++;
1406: PetscCall(VecPointwiseMult(rwork1, ldfp->diag_vec_recycle_order, rwork1));
1407: PetscCall(MatDQNApplyJ0Fwd(B, ldfp->column_work, Z));
1409: PetscCall(MatMultHermitianTransposeAddColumnRange(Sfull, Z, rwork1, rwork1, 0, h));
1410: ldfp->St_count++;
1412: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1413: PetscCall(MatUpperTriangularSolveInPlace(B, ldfp->YtS_triu, rwork1, PETSC_TRUE, ldfp->num_updates, ldfp->strategy));
1414: PetscCall(VecScale(rwork1, -1.0));
1415: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1417: PetscCall(MatMultAddColumnRange(Yfull, rwork1, Z, Z, 0, h));
1418: ldfp->Y_count++;
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: /*
1423: This dense representation reduces the L-DFP update to a series of
1424: matrix-vector products with dense matrices in lieu of the conventional
1425: matrix-free two-loop algorithm.
1426: */
1427: PetscErrorCode MatCreate_LMVMDDFP(Mat B)
1428: {
1429: Mat_LMVM *lmvm;
1430: Mat_DQN *ldfp;
1432: PetscFunctionBegin;
1433: PetscCall(MatCreate_LMVM(B));
1434: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDDFP));
1435: PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
1436: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
1437: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
1438: B->ops->view = MatView_LMVMDQN;
1439: B->ops->setup = MatSetUp_LMVMDQN;
1440: B->ops->setfromoptions = MatSetFromOptions_LMVMDQN;
1441: B->ops->destroy = MatDestroy_LMVMDQN;
1443: lmvm = (Mat_LMVM *)B->data;
1444: lmvm->ops->reset = MatReset_LMVMDQN;
1445: lmvm->ops->update = MatUpdate_LMVMDQN;
1446: lmvm->ops->mult = MatMult_LMVMDDFP;
1447: lmvm->ops->solve = MatSolve_LMVMDDFP;
1448: lmvm->ops->copy = MatCopy_LMVMDQN;
1450: lmvm->ops->multht = lmvm->ops->mult;
1451: lmvm->ops->solveht = lmvm->ops->solve;
1453: PetscCall(PetscNew(&ldfp));
1454: lmvm->ctx = (void *)ldfp;
1455: ldfp->allocated = PETSC_FALSE;
1456: ldfp->watchdog = 0;
1457: ldfp->max_seq_rejects = lmvm->m / 2;
1458: ldfp->strategy = MAT_LMVM_DENSE_INPLACE;
1459: ldfp->use_recursive = PETSC_TRUE;
1460: ldfp->needPQ = PETSC_TRUE;
1462: PetscCall(SymBroydenRescaleCreate(&ldfp->rescale));
1463: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_LMVMDQN));
1464: PetscFunctionReturn(PETSC_SUCCESS);
1465: }
1467: /*@
1468: MatCreateLMVMDDFP - Creates a dense representation of the limited-memory
1469: Davidon-Fletcher-Powell (DFP) approximation to a Hessian.
1471: Collective
1473: Input Parameters:
1474: + comm - MPI communicator
1475: . n - number of local rows for storage vectors
1476: - N - global size of the storage vectors
1478: Output Parameter:
1479: . B - the matrix
1481: Level: advanced
1483: Note:
1484: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
1485: paradigm instead of this routine directly.
1487: .seealso: `MatCreate()`, `MATLMVM`, `MATLMVMDDFP`, `MatCreateLMVMDFP()`
1488: @*/
1489: PetscErrorCode MatCreateLMVMDDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
1490: {
1491: PetscFunctionBegin;
1492: PetscCall(KSPInitializePackage());
1493: PetscCall(MatCreate(comm, B));
1494: PetscCall(MatSetSizes(*B, n, n, N, N));
1495: PetscCall(MatSetType(*B, MATLMVMDDFP));
1496: PetscCall(MatSetUp(*B));
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1500: /*@
1501: MatLMVMDenseSetType - Sets the memory storage type for dense `MATLMVM`
1503: Input Parameters:
1504: + B - the `MATLMVM` matrix
1505: - type - scale type, see `MatLMVMDenseSetType`
1507: Options Database Keys:
1508: + -mat_lqn_type (reorder|inplace) - set the strategy
1509: . -mat_lbfgs_type (reorder|inplace) - set the strategy
1510: - -mat_ldfp_type (reorder|inplace) - set the strategy
1512: Level: intermediate
1514: MatLMVMDenseTypes\:
1515: + `MAT_LMVM_DENSE_REORDER` - reorders memory to minimize kernel launch
1516: - `MAT_LMVM_DENSE_INPLACE` - launches kernel inplace to minimize memory movement
1518: .seealso: [](ch_ksp), `MATLMVMDQN`, `MATLMVMDBFGS`, `MATLMVMDDFP`, `MatLMVMDenseType`
1519: @*/
1520: PetscErrorCode MatLMVMDenseSetType(Mat B, MatLMVMDenseType type)
1521: {
1522: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1523: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
1525: PetscFunctionBegin;
1527: lqn->strategy = type;
1528: PetscFunctionReturn(PETSC_SUCCESS);
1529: }