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: }
352: }
353: PetscCall(VecRestoreArrayAndMemType(lqn->StFprev, &StFprev));
355: {
356: Vec this_sy_col;
357: /* Now StFprev is updated for the new S vector. Write -StFprev into the appropriate row */
358: PetscCall(MatDenseGetColumnVecWrite(lqn->StY_triu, StYidx, &this_sy_col));
359: PetscCall(VecAXPBY(this_sy_col, -1.0, 0.0, lqn->StFprev));
361: /* Now compute the new StFprev */
362: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, lqn->StFprev, 0, h_new));
363: lqn->St_count++;
365: /* Now add StFprev: this_sy_col == S^T (F - Fprev) == S^T y */
366: PetscCall(VecAXPY(this_sy_col, 1.0, lqn->StFprev));
368: if (lqn->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, this_sy_col, lqn->num_updates, lqn->cyclic_work_vec));
369: PetscCall(MatDenseRestoreColumnVecWrite(lqn->StY_triu, StYidx, &this_sy_col));
370: }
371: }
373: if (is_ddfp || is_dqn) {
374: PetscInt YtSidx;
376: YtSidx = (lqn->strategy == MAT_LMVM_DENSE_REORDER) ? history_index(m, lqn->num_updates, k) : idx;
378: {
379: Vec this_ys_col;
381: PetscCall(MatDenseGetColumnVecWrite(lqn->YtS_triu, YtSidx, &this_ys_col));
382: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, lmvm->Xprev, this_ys_col, 0, h_new));
383: lqn->Yt_count++;
385: if (lqn->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, this_ys_col, lqn->num_updates, lqn->cyclic_work_vec));
386: PetscCall(MatDenseRestoreColumnVecWrite(lqn->YtS_triu, YtSidx, &this_ys_col));
387: }
388: }
390: if (is_dbfgs || is_dqn) {
391: PetscCall(MatGetDiagonal(lqn->StY_triu, lqn->diag_vec));
392: } else if (is_ddfp) {
393: PetscCall(MatGetDiagonal(lqn->YtS_triu, lqn->diag_vec));
394: } else {
395: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "MatUpdate_LMVMDQN is only available for dense derived types. (DBFGS, DDFP, DQN");
396: }
398: if (lqn->strategy == MAT_LMVM_DENSE_REORDER) {
399: if (!lqn->diag_vec_recycle_order) PetscCall(VecDuplicate(lqn->diag_vec, &lqn->diag_vec_recycle_order));
400: PetscCall(VecCopy(lqn->diag_vec, lqn->diag_vec_recycle_order));
401: PetscCall(VecHistoryOrderToRecycleOrder(B, lqn->diag_vec_recycle_order, lqn->num_updates, lqn->cyclic_work_vec));
402: } else {
403: if (!lqn->diag_vec_recycle_order) {
404: PetscCall(PetscObjectReference((PetscObject)lqn->diag_vec));
405: lqn->diag_vec_recycle_order = lqn->diag_vec;
406: }
407: }
409: PetscCall(SymBroydenRescaleUpdate(B, lqn->rescale));
410: } else {
411: /* Update is bad, skip it */
412: ++lmvm->nrejects;
413: ++lqn->watchdog;
414: PetscInt m = lmvm->m;
415: PetscInt k = lmvm->k;
416: PetscInt h = k - oldest_update(m, k);
418: /* we still have to maintain StFprev */
419: if (!lqn->StFprev) PetscCall(VecDuplicate(lqn->rwork1, &lqn->StFprev));
420: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, lqn->StFprev, 0, h));
421: lqn->St_count++;
422: }
423: }
425: if (lqn->watchdog > lqn->max_seq_rejects) PetscCall(MatLMVMReset(B, PETSC_FALSE));
427: /* Save the solution and function to be used in the next update */
428: PetscCall(VecCopy(X, lmvm->Xprev));
429: PetscCall(VecCopy(F, lmvm->Fprev));
430: PetscCall(PetscObjectReference((PetscObject)F));
431: PetscCall(VecDestroy(&lqn->Fprev_ref));
432: lqn->Fprev_ref = F;
433: PetscCall(PetscObjectStateGet((PetscObject)F, &lqn->Fprev_state));
434: lmvm->prev_set = PETSC_TRUE;
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode MatDestroyThenCopy(Mat src, Mat *dst)
439: {
440: PetscFunctionBegin;
441: PetscCall(MatDestroy(dst));
442: if (src) PetscCall(MatDuplicate(src, MAT_COPY_VALUES, dst));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: static PetscErrorCode VecDestroyThenCopy(Vec src, Vec *dst)
447: {
448: PetscFunctionBegin;
449: PetscCall(VecDestroy(dst));
450: if (src) {
451: PetscCall(VecDuplicate(src, dst));
452: PetscCall(VecCopy(src, *dst));
453: }
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode MatCopy_LMVMDQN(Mat B, Mat M, MatStructure str)
458: {
459: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
460: Mat_DQN *blqn = (Mat_DQN *)bdata->ctx;
461: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
462: Mat_DQN *mlqn = (Mat_DQN *)mdata->ctx;
463: PetscBool is_dbfgs, is_ddfp, is_dqn;
465: PetscFunctionBegin;
466: PetscCall(SymBroydenRescaleCopy(blqn->rescale, mlqn->rescale));
467: mlqn->num_updates = blqn->num_updates;
468: mlqn->num_mult_updates = blqn->num_mult_updates;
469: mlqn->dense_type = blqn->dense_type;
470: mlqn->strategy = blqn->strategy;
471: mlqn->S_count = 0;
472: mlqn->St_count = 0;
473: mlqn->Y_count = 0;
474: mlqn->Yt_count = 0;
475: mlqn->watchdog = blqn->watchdog;
476: mlqn->max_seq_rejects = blqn->max_seq_rejects;
477: mlqn->use_recursive = blqn->use_recursive;
478: mlqn->needPQ = blqn->needPQ;
479: if (blqn->allocated) {
480: PetscCall(MatAllocate_LMVMDQN_Internal(M));
481: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
482: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
483: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
484: PetscCall(MatDestroyThenCopy(blqn->HY, &mlqn->BS));
485: PetscCall(VecDestroyThenCopy(blqn->StFprev, &mlqn->StFprev));
486: PetscCall(MatDestroyThenCopy(blqn->StY_triu, &mlqn->StY_triu));
487: PetscCall(MatDestroyThenCopy(blqn->StY_triu_strict, &mlqn->StY_triu_strict));
488: PetscCall(MatDestroyThenCopy(blqn->YtS_triu, &mlqn->YtS_triu));
489: PetscCall(MatDestroyThenCopy(blqn->YtS_triu_strict, &mlqn->YtS_triu_strict));
490: PetscCall(MatDestroyThenCopy(blqn->YtHY, &mlqn->YtHY));
491: PetscCall(MatDestroyThenCopy(blqn->StBS, &mlqn->StBS));
492: PetscCall(MatDestroyThenCopy(blqn->J, &mlqn->J));
493: PetscCall(VecDestroyThenCopy(blqn->diag_vec, &mlqn->diag_vec));
494: PetscCall(VecDestroyThenCopy(blqn->diag_vec_recycle_order, &mlqn->diag_vec_recycle_order));
495: PetscCall(VecDestroyThenCopy(blqn->inv_diag_vec, &mlqn->inv_diag_vec));
496: if (blqn->use_recursive && (is_dbfgs || is_ddfp)) {
497: for (PetscInt i = 0; i < bdata->m; i++) {
498: PetscCall(VecDestroyThenCopy(blqn->PQ[i], &mlqn->PQ[i]));
499: mlqn->yts[i] = blqn->yts[i];
500: if (is_dbfgs) {
501: mlqn->stp[i] = blqn->stp[i];
502: } else if (is_ddfp) {
503: mlqn->ytq[i] = blqn->ytq[i];
504: }
505: }
506: }
507: }
508: PetscCall(PetscObjectReference((PetscObject)blqn->Fprev_ref));
509: PetscCall(VecDestroy(&mlqn->Fprev_ref));
510: mlqn->Fprev_ref = blqn->Fprev_ref;
511: mlqn->Fprev_state = blqn->Fprev_state;
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: static PetscErrorCode MatMult_LMVMDQN(Mat B, Vec X, Vec Z)
516: {
517: PetscFunctionBegin;
518: PetscCall(MatMult_LMVMDDFP(B, X, Z));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: static PetscErrorCode MatSolve_LMVMDQN(Mat H, Vec F, Vec dX)
523: {
524: PetscFunctionBegin;
525: PetscCall(MatSolve_LMVMDBFGS(H, F, dX));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: static PetscErrorCode MatLMVMSymBroydenSetDelta_LMVMDQN(Mat B, PetscScalar delta)
530: {
531: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
532: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
534: PetscFunctionBegin;
535: PetscCall(SymBroydenRescaleSetDelta(B, lqn->rescale, PetscAbsReal(PetscRealPart(delta))));
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*
540: This dense representation uses Davidon-Fletcher-Powell (DFP) for MatMult,
541: and Broyden-Fletcher-Goldfarb-Shanno (BFGS) for MatSolve. This implementation
542: results in avoiding costly Cholesky factorization, at the cost of duality cap.
543: Please refer to MatLMVMDDFP and MatLMVMDBFGS for more information.
544: */
545: PetscErrorCode MatCreate_LMVMDQN(Mat B)
546: {
547: Mat_LMVM *lmvm;
548: Mat_DQN *lqn;
550: PetscFunctionBegin;
551: PetscCall(MatCreate_LMVM(B));
552: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDQN));
553: PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
554: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
555: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
556: B->ops->view = MatView_LMVMDQN;
557: B->ops->setup = MatSetUp_LMVMDQN;
558: B->ops->setfromoptions = MatSetFromOptions_LMVMDQN;
559: B->ops->destroy = MatDestroy_LMVMDQN;
561: lmvm = (Mat_LMVM *)B->data;
562: lmvm->ops->reset = MatReset_LMVMDQN;
563: lmvm->ops->update = MatUpdate_LMVMDQN;
564: lmvm->ops->mult = MatMult_LMVMDQN;
565: lmvm->ops->solve = MatSolve_LMVMDQN;
566: lmvm->ops->copy = MatCopy_LMVMDQN;
568: lmvm->ops->multht = lmvm->ops->mult;
569: lmvm->ops->solveht = lmvm->ops->solve;
571: PetscCall(PetscNew(&lqn));
572: lmvm->ctx = (void *)lqn;
573: lqn->allocated = PETSC_FALSE;
574: lqn->use_recursive = PETSC_FALSE;
575: lqn->needPQ = PETSC_FALSE;
576: lqn->watchdog = 0;
577: lqn->max_seq_rejects = lmvm->m / 2;
578: lqn->strategy = MAT_LMVM_DENSE_INPLACE;
580: PetscCall(SymBroydenRescaleCreate(&lqn->rescale));
581: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_LMVMDQN));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: /*@
586: MatCreateLMVMDQN - Creates a dense representation of the limited-memory
587: Quasi-Newton approximation to a Hessian.
589: Collective
591: Input Parameters:
592: + comm - MPI communicator
593: . n - number of local rows for storage vectors
594: - N - global size of the storage vectors
596: Output Parameter:
597: . B - the matrix
599: Level: advanced
601: Note:
602: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
603: paradigm instead of this routine directly.
605: .seealso: `MatCreate()`, `MATLMVM`, `MATLMVMDBFGS`, `MATLMVMDDFP`, `MatCreateLMVMDDFP()`, `MatCreateLMVMDBFGS()`
606: @*/
607: PetscErrorCode MatCreateLMVMDQN(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
608: {
609: PetscFunctionBegin;
610: PetscCall(KSPInitializePackage());
611: PetscCall(MatCreate(comm, B));
612: PetscCall(MatSetSizes(*B, n, n, N, N));
613: PetscCall(MatSetType(*B, MATLMVMDQN));
614: PetscCall(MatSetUp(*B));
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: static PetscErrorCode MatDQNApplyJ0Fwd(Mat B, Vec X, Vec Z)
619: {
620: PetscFunctionBegin;
621: PetscCall(MatLMVMApplyJ0Fwd(B, X, Z));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static PetscErrorCode MatDQNApplyJ0Inv(Mat B, Vec F, Vec dX)
626: {
627: PetscFunctionBegin;
628: PetscCall(MatLMVMApplyJ0Inv(B, F, dX));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /* This is not Bunch-Kaufman LDLT: here L is strictly lower triangular part of STY */
633: static PetscErrorCode MatGetLDLT(Mat B, Mat result)
634: {
635: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
636: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
637: PetscInt m_local;
639: PetscFunctionBegin;
640: if (!lbfgs->temp_mat) PetscCall(MatDuplicate(lbfgs->YtS_triu_strict, MAT_SHARE_NONZERO_PATTERN, &lbfgs->temp_mat));
641: PetscCall(MatCopy(lbfgs->YtS_triu_strict, lbfgs->temp_mat, SAME_NONZERO_PATTERN));
642: PetscCall(MatDiagonalScale(lbfgs->temp_mat, lbfgs->inv_diag_vec, NULL));
643: PetscCall(MatGetLocalSize(result, &m_local, NULL));
644: // need to conjugate and conjugate again because we have MatTransposeMatMult but not MatHermitianTransposeMatMult()
645: PetscCall(MatConjugate(lbfgs->temp_mat));
646: if (m_local) {
647: Mat temp_local, YtS_local, result_local;
648: PetscCall(MatDenseGetLocalMatrix(lbfgs->YtS_triu_strict, &YtS_local));
649: PetscCall(MatDenseGetLocalMatrix(lbfgs->temp_mat, &temp_local));
650: PetscCall(MatDenseGetLocalMatrix(result, &result_local));
651: PetscCall(MatTransposeMatMult(YtS_local, temp_local, MAT_REUSE_MATRIX, PETSC_DETERMINE, &result_local));
652: }
653: PetscCall(MatConjugate(result));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: static PetscErrorCode MatLMVMDBFGSUpdateMultData(Mat B)
658: {
659: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
660: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
661: PetscInt m = lmvm->m, m_local;
662: PetscInt k = lmvm->k;
663: PetscInt h = k - oldest_update(m, k);
664: PetscInt j_0;
665: PetscInt prev_oldest;
666: Mat J_local;
667: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
668: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
670: PetscFunctionBegin;
671: if (!lbfgs->YtS_triu_strict) {
672: PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->YtS_triu_strict));
673: PetscCall(MatDestroy(&lbfgs->StBS));
674: PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->StBS));
675: PetscCall(MatDestroy(&lbfgs->J));
676: PetscCall(MatDuplicate(lbfgs->StY_triu, MAT_SHARE_NONZERO_PATTERN, &lbfgs->J));
677: PetscCall(MatDestroy(&lbfgs->BS));
678: PetscCall(MatDuplicate(Yfull, MAT_SHARE_NONZERO_PATTERN, &lbfgs->BS));
679: PetscCall(MatShift(lbfgs->StBS, 1.0));
680: lbfgs->num_mult_updates = oldest_update(m, k);
681: }
682: if (lbfgs->num_mult_updates == k) PetscFunctionReturn(PETSC_SUCCESS);
684: /* B_0 may have been updated, we must recompute B_0 S and S^T B_0 S */
685: for (PetscInt j = oldest_update(m, k); j < k; j++) {
686: Vec s_j;
687: Vec Bs_j;
688: Vec StBs_j;
689: PetscInt S_idx = recycle_index(m, j);
690: PetscInt StBS_idx = lbfgs->strategy == MAT_LMVM_DENSE_INPLACE ? S_idx : history_index(m, k, j);
692: PetscCall(MatDenseGetColumnVecWrite(lbfgs->BS, S_idx, &Bs_j));
693: PetscCall(MatDenseGetColumnVecRead(Sfull, S_idx, &s_j));
694: PetscCall(MatDQNApplyJ0Fwd(B, s_j, Bs_j));
695: PetscCall(MatDenseRestoreColumnVecRead(Sfull, S_idx, &s_j));
696: PetscCall(MatDenseGetColumnVecWrite(lbfgs->StBS, StBS_idx, &StBs_j));
697: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, Bs_j, StBs_j, 0, h));
698: lbfgs->St_count++;
699: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, StBs_j, lbfgs->num_updates, lbfgs->cyclic_work_vec));
700: PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->StBS, StBS_idx, &StBs_j));
701: PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->BS, S_idx, &Bs_j));
702: }
703: prev_oldest = oldest_update(m, lbfgs->num_mult_updates);
704: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER && prev_oldest < oldest_update(m, k)) {
705: /* move the YtS entries that have been computed and need to be kept back up */
706: PetscInt m_keep = m - (oldest_update(m, k) - prev_oldest);
708: PetscCall(MatMove_LR3(B, lbfgs->YtS_triu_strict, m_keep));
709: }
710: PetscCall(MatGetLocalSize(lbfgs->YtS_triu_strict, &m_local, NULL));
711: j_0 = PetscMax(lbfgs->num_mult_updates, oldest_update(m, k));
712: for (PetscInt j = j_0; j < k; j++) {
713: PetscInt S_idx = recycle_index(m, j);
714: PetscInt YtS_idx = lbfgs->strategy == MAT_LMVM_DENSE_INPLACE ? S_idx : history_index(m, k, j);
715: Vec s_j, Yts_j;
717: PetscCall(MatDenseGetColumnVecRead(Sfull, S_idx, &s_j));
718: PetscCall(MatDenseGetColumnVecWrite(lbfgs->YtS_triu_strict, YtS_idx, &Yts_j));
719: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, s_j, Yts_j, 0, h));
720: lbfgs->Yt_count++;
721: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, Yts_j, lbfgs->num_updates, lbfgs->cyclic_work_vec));
722: PetscCall(MatDenseRestoreColumnVecWrite(lbfgs->YtS_triu_strict, YtS_idx, &Yts_j));
723: PetscCall(MatDenseRestoreColumnVecRead(Sfull, S_idx, &s_j));
724: /* zero the corresponding row */
725: if (m_local > 0) {
726: Mat YtS_local, YtS_row;
728: PetscCall(MatDenseGetLocalMatrix(lbfgs->YtS_triu_strict, &YtS_local));
729: PetscCall(MatDenseGetSubMatrix(YtS_local, YtS_idx, YtS_idx + 1, PETSC_DECIDE, PETSC_DECIDE, &YtS_row));
730: PetscCall(MatZeroEntries(YtS_row));
731: PetscCall(MatDenseRestoreSubMatrix(YtS_local, &YtS_row));
732: }
733: }
734: if (!lbfgs->inv_diag_vec) PetscCall(VecDuplicate(lbfgs->diag_vec, &lbfgs->inv_diag_vec));
735: PetscCall(VecCopy(lbfgs->diag_vec, lbfgs->inv_diag_vec));
736: PetscCall(VecReciprocal(lbfgs->inv_diag_vec));
737: PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
738: PetscCall(MatSetFactorType(J_local, MAT_FACTOR_NONE));
739: PetscCall(MatGetLDLT(B, lbfgs->J));
740: PetscCall(MatAXPY(lbfgs->J, 1.0, lbfgs->StBS, SAME_NONZERO_PATTERN));
741: if (m_local) {
742: PetscCall(MatSetOption(J_local, MAT_SPD, PETSC_TRUE));
743: PetscCall(MatCholeskyFactor(J_local, NULL, NULL));
744: }
745: lbfgs->num_mult_updates = lbfgs->num_updates;
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /* Solves for
750: * [ I | -S R^{-T} ] [ I | 0 ] [ H_0 | 0 ] [ I | Y ] [ I ]
751: * [-----+---] [-----+---] [---+---] [-------------]
752: * [ Y^T | I ] [ 0 | D ] [ 0 | I ] [ -R^{-1} S^T ] */
754: static PetscErrorCode MatSolve_LMVMDBFGS(Mat H, Vec F, Vec dX)
755: {
756: Mat_LMVM *lmvm = (Mat_LMVM *)H->data;
757: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
758: Vec rwork1 = lbfgs->rwork1;
759: PetscInt m = lmvm->m;
760: PetscInt k = lmvm->k;
761: PetscInt h = k - oldest_update(m, k);
762: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
763: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
764: PetscObjectState Fstate;
766: PetscFunctionBegin;
767: VecCheckSameSize(F, 2, dX, 3);
768: VecCheckMatCompatible(H, dX, 3, F, 2);
770: /* Block Version */
771: if (!lbfgs->num_updates) {
772: PetscCall(MatDQNApplyJ0Inv(H, F, dX));
773: PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
774: }
776: PetscCall(PetscObjectStateGet((PetscObject)F, &Fstate));
777: if (F == lbfgs->Fprev_ref && Fstate == lbfgs->Fprev_state) {
778: PetscCall(VecCopy(lbfgs->StFprev, rwork1));
779: } else {
780: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, rwork1, 0, h));
781: lbfgs->St_count++;
782: }
784: /* Reordering rwork1, as STY is in history order, while S is in recycled order */
785: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
786: PetscCall(MatUpperTriangularSolveInPlace(H, lbfgs->StY_triu, rwork1, PETSC_FALSE, lbfgs->num_updates, lbfgs->strategy));
787: PetscCall(VecScale(rwork1, -1.0));
788: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
790: PetscCall(VecCopy(F, lbfgs->column_work));
791: PetscCall(MatMultAddColumnRange(Yfull, rwork1, lbfgs->column_work, lbfgs->column_work, 0, h));
792: lbfgs->Y_count++;
794: PetscCall(VecPointwiseMult(rwork1, lbfgs->diag_vec_recycle_order, rwork1));
795: PetscCall(MatDQNApplyJ0Inv(H, lbfgs->column_work, dX));
797: PetscCall(MatMultHermitianTransposeAddColumnRange(Yfull, dX, rwork1, rwork1, 0, h));
798: lbfgs->Yt_count++;
800: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
801: PetscCall(MatUpperTriangularSolveInPlace(H, lbfgs->StY_triu, rwork1, PETSC_TRUE, lbfgs->num_updates, lbfgs->strategy));
802: PetscCall(VecScale(rwork1, -1.0));
803: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(H, rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
805: PetscCall(MatMultAddColumnRange(Sfull, rwork1, dX, dX, 0, h));
806: lbfgs->S_count++;
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /* Solves for
811: B_0 - [ Y | B_0 S] [ -D | L^T ]^-1 [ Y^T ]
812: [-----+-----------] [---------]
813: [ L | S^T B_0 S ] [ S^T B_0 ]
815: Above is equivalent to
817: B_0 - [ Y | B_0 S] [[ I | 0 ][ -D | 0 ][ I | -D^{-1} L^T ]]^-1 [ Y^T ]
818: [[-----------+---][-----+---][---+-------------]] [---------]
819: [[ -L D^{-1} | I ][ 0 | J ][ 0 | I ]] [ S^T B_0 ]
821: where J = S^T B_0 S + L D^{-1} L^T
823: becomes
825: B_0 - [ Y | B_0 S] [ I | D^{-1} L^T ][ -D^{-1} | 0 ][ I | 0 ] [ Y^T ]
826: [---+------------][----------+--------][----------+---] [---------]
827: [ 0 | I ][ 0 | J^{-1} ][ L D^{-1} | I ] [ S^T B_0 ]
829: =
831: B_0 + [ Y | B_0 S] [ D^{-1} | 0 ][ I | L^T ][ I | 0 ][ I | 0 ] [ Y^T ]
832: [--------+---][---+-----][---+---------][----------+---] [---------]
833: [ 0 | I ][ 0 | I ][ 0 | -J^{-1} ][ L D^{-1} | I ] [ S^T B_0 ]
835: (Note that YtS_triu_strict is L^T)
836: Byrd, Nocedal, Schnabel 1994
838: Alternative approach: considering the fact that DFP is dual to BFGS, use MatMult of DPF:
839: (See ddfp.c's MatMult_LMVMDDFP)
841: */
842: static PetscErrorCode MatMult_LMVMDBFGS(Mat B, Vec X, Vec Z)
843: {
844: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
845: Mat_DQN *lbfgs = (Mat_DQN *)lmvm->ctx;
846: Mat J_local;
847: PetscInt idx, i, j, m_local, local_n;
848: PetscInt m = lmvm->m;
849: PetscInt k = lmvm->k;
850: PetscInt h = k - oldest_update(m, k);
851: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
852: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
854: PetscFunctionBegin;
855: VecCheckSameSize(X, 2, Z, 3);
856: VecCheckMatCompatible(B, X, 2, Z, 3);
858: /* Cholesky Version */
859: /* Start with the B0 term */
860: PetscCall(MatDQNApplyJ0Fwd(B, X, Z));
861: if (!lbfgs->num_updates) PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
863: if (lbfgs->use_recursive) {
864: PetscDeviceContext dctx;
865: PetscMemType memtype;
866: PetscScalar stz, ytx, stp, sjtpi, yjtsi, *workscalar;
867: PetscInt oldest = oldest_update(m, k);
869: PetscCall(PetscDeviceContextGetCurrentContext(&dctx));
870: /* Recursive formulation to avoid Cholesky. Not a dense formulation */
871: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, lbfgs->rwork1, 0, h));
872: lbfgs->Yt_count++;
874: PetscCall(VecGetLocalSize(lbfgs->rwork1, &local_n));
876: if (lbfgs->needPQ) {
877: PetscInt oldest = oldest_update(m, k);
878: for (i = oldest; i < k; ++i) {
879: idx = recycle_index(m, i);
880: /* column_work = S[idx] */
881: PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work, idx));
882: PetscCall(MatDQNApplyJ0Fwd(B, lbfgs->column_work, lbfgs->PQ[idx]));
883: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, lbfgs->column_work, lbfgs->rwork3, 0, h));
884: PetscCall(VecGetArrayAndMemType(lbfgs->rwork3, &workscalar, &memtype));
885: for (j = oldest; j < i; ++j) {
886: PetscInt idx_j = recycle_index(m, j);
887: /* Copy yjtsi in device-aware manner */
888: if (local_n) {
889: if (PetscMemTypeHost(memtype)) {
890: yjtsi = workscalar[idx_j];
891: } else {
892: PetscCall(PetscDeviceRegisterMemory(&yjtsi, PETSC_MEMTYPE_HOST, sizeof(yjtsi)));
893: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
894: PetscCall(PetscDeviceArrayCopy(dctx, &yjtsi, &workscalar[idx_j], 1));
895: }
896: }
897: PetscCallMPI(MPI_Bcast(&yjtsi, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)B)));
898: /* column_work2 = S[j] */
899: PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work2, idx_j));
900: PetscCall(VecDot(lbfgs->PQ[idx], lbfgs->column_work2, &sjtpi));
901: /* column_work2 = Y[j] */
902: PetscCall(MatGetColumnVector(Yfull, lbfgs->column_work2, idx_j));
903: /* Compute the pure BFGS component of the forward product */
904: PetscCall(VecAXPBYPCZ(lbfgs->PQ[idx], -sjtpi / lbfgs->stp[idx_j], yjtsi / lbfgs->yts[idx_j], 1.0, lbfgs->PQ[idx_j], lbfgs->column_work2));
905: }
906: PetscCall(VecDot(lbfgs->PQ[idx], lbfgs->column_work, &stp));
907: lbfgs->stp[idx] = PetscRealPart(stp);
908: }
909: lbfgs->needPQ = PETSC_FALSE;
910: }
912: PetscCall(VecGetArrayAndMemType(lbfgs->rwork1, &workscalar, &memtype));
913: for (i = oldest; i < k; ++i) {
914: idx = recycle_index(m, i);
915: /* Copy stz[i], ytx[i] in device-aware manner */
916: if (local_n) {
917: if (PetscMemTypeHost(memtype)) {
918: ytx = workscalar[idx];
919: } else {
920: PetscCall(PetscDeviceRegisterMemory(&ytx, PETSC_MEMTYPE_HOST, 1 * sizeof(ytx)));
921: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
922: PetscCall(PetscDeviceArrayCopy(dctx, &ytx, &workscalar[idx], 1));
923: }
924: }
925: PetscCallMPI(MPI_Bcast(&ytx, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)B)));
926: /* column_work : S[i], column_work2 : Y[i] */
927: PetscCall(MatGetColumnVector(Sfull, lbfgs->column_work, idx));
928: PetscCall(MatGetColumnVector(Yfull, lbfgs->column_work2, idx));
929: PetscCall(VecDot(Z, lbfgs->column_work, &stz));
930: PetscCall(VecAXPBYPCZ(Z, -stz / lbfgs->stp[idx], ytx / lbfgs->yts[idx], 1.0, lbfgs->PQ[idx], lbfgs->column_work2));
931: }
932: PetscCall(VecRestoreArrayAndMemType(lbfgs->rwork1, &workscalar));
933: } else {
934: PetscCall(MatLMVMDBFGSUpdateMultData(B));
935: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, lbfgs->rwork1, 0, h));
936: lbfgs->Yt_count++;
937: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, Z, lbfgs->rwork2, 0, h));
938: lbfgs->St_count++;
939: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) {
940: PetscCall(VecRecycleOrderToHistoryOrder(B, lbfgs->rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
941: PetscCall(VecRecycleOrderToHistoryOrder(B, lbfgs->rwork2, lbfgs->num_updates, lbfgs->cyclic_work_vec));
942: }
944: PetscCall(VecPointwiseMult(lbfgs->rwork1, lbfgs->rwork1, lbfgs->inv_diag_vec));
945: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(lbfgs->YtS_triu_strict));
946: PetscCall(MatMultTransposeAdd(lbfgs->YtS_triu_strict, lbfgs->rwork1, lbfgs->rwork2, lbfgs->rwork2));
947: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(lbfgs->YtS_triu_strict));
949: if (!lbfgs->rwork2_local) PetscCall(VecCreateLocalVector(lbfgs->rwork2, &lbfgs->rwork2_local));
950: if (!lbfgs->rwork3_local) PetscCall(VecCreateLocalVector(lbfgs->rwork3, &lbfgs->rwork3_local));
951: PetscCall(VecGetLocalVectorRead(lbfgs->rwork2, lbfgs->rwork2_local));
952: PetscCall(VecGetLocalVector(lbfgs->rwork3, lbfgs->rwork3_local));
953: PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
954: PetscCall(VecGetSize(lbfgs->rwork2_local, &m_local));
955: if (m_local) {
956: PetscCall(MatDenseGetLocalMatrix(lbfgs->J, &J_local));
957: PetscCall(MatSolve(J_local, lbfgs->rwork2_local, lbfgs->rwork3_local));
958: }
959: PetscCall(VecRestoreLocalVector(lbfgs->rwork3, lbfgs->rwork3_local));
960: PetscCall(VecRestoreLocalVectorRead(lbfgs->rwork2, lbfgs->rwork2_local));
961: PetscCall(VecScale(lbfgs->rwork3, -1.0));
963: PetscCall(MatMult(lbfgs->YtS_triu_strict, lbfgs->rwork3, lbfgs->rwork2));
964: PetscCall(VecPointwiseMult(lbfgs->rwork2, lbfgs->rwork2, lbfgs->inv_diag_vec));
965: PetscCall(VecAXPY(lbfgs->rwork1, 1.0, lbfgs->rwork2));
967: if (lbfgs->strategy == MAT_LMVM_DENSE_REORDER) {
968: PetscCall(VecHistoryOrderToRecycleOrder(B, lbfgs->rwork1, lbfgs->num_updates, lbfgs->cyclic_work_vec));
969: PetscCall(VecHistoryOrderToRecycleOrder(B, lbfgs->rwork3, lbfgs->num_updates, lbfgs->cyclic_work_vec));
970: }
972: PetscCall(MatMultAddColumnRange(Yfull, lbfgs->rwork1, Z, Z, 0, h));
973: lbfgs->Y_count++;
974: PetscCall(MatMultAddColumnRange(lbfgs->BS, lbfgs->rwork3, Z, Z, 0, h));
975: lbfgs->S_count++;
976: }
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*
981: This dense representation reduces the L-BFGS update to a series of
982: matrix-vector products with dense matrices in lieu of the conventional matrix-free
983: two-loop algorithm.
984: */
985: PetscErrorCode MatCreate_LMVMDBFGS(Mat B)
986: {
987: Mat_LMVM *lmvm;
988: Mat_DQN *lbfgs;
990: PetscFunctionBegin;
991: PetscCall(MatCreate_LMVM(B));
992: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDBFGS));
993: PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
994: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
995: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
996: B->ops->view = MatView_LMVMDQN;
997: B->ops->setup = MatSetUp_LMVMDQN;
998: B->ops->setfromoptions = MatSetFromOptions_LMVMDQN;
999: B->ops->destroy = MatDestroy_LMVMDQN;
1001: lmvm = (Mat_LMVM *)B->data;
1002: lmvm->ops->reset = MatReset_LMVMDQN;
1003: lmvm->ops->update = MatUpdate_LMVMDQN;
1004: lmvm->ops->mult = MatMult_LMVMDBFGS;
1005: lmvm->ops->solve = MatSolve_LMVMDBFGS;
1006: lmvm->ops->copy = MatCopy_LMVMDQN;
1008: lmvm->ops->multht = lmvm->ops->mult;
1009: lmvm->ops->solveht = lmvm->ops->solve;
1011: PetscCall(PetscNew(&lbfgs));
1012: lmvm->ctx = (void *)lbfgs;
1013: lbfgs->allocated = PETSC_FALSE;
1014: lbfgs->use_recursive = PETSC_TRUE;
1015: lbfgs->needPQ = PETSC_TRUE;
1016: lbfgs->watchdog = 0;
1017: lbfgs->max_seq_rejects = lmvm->m / 2;
1018: lbfgs->strategy = MAT_LMVM_DENSE_INPLACE;
1020: PetscCall(SymBroydenRescaleCreate(&lbfgs->rescale));
1021: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_LMVMDQN));
1022: PetscFunctionReturn(PETSC_SUCCESS);
1023: }
1025: /*@
1026: MatCreateLMVMDBFGS - Creates a dense representation of the limited-memory
1027: Broyden-Fletcher-Goldfarb-Shanno (BFGS) approximation to a Hessian.
1029: Collective
1031: Input Parameters:
1032: + comm - MPI communicator
1033: . n - number of local rows for storage vectors
1034: - N - global size of the storage vectors
1036: Output Parameter:
1037: . B - the matrix
1039: Level: advanced
1041: Note:
1042: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
1043: paradigm instead of this routine directly.
1045: .seealso: `MatCreate()`, `MATLMVM`, `MATLMVMDBFGS`, `MatCreateLMVMBFGS()`
1046: @*/
1047: PetscErrorCode MatCreateLMVMDBFGS(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
1048: {
1049: PetscFunctionBegin;
1050: PetscCall(KSPInitializePackage());
1051: PetscCall(MatCreate(comm, B));
1052: PetscCall(MatSetSizes(*B, n, n, N, N));
1053: PetscCall(MatSetType(*B, MATLMVMDBFGS));
1054: PetscCall(MatSetUp(*B));
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: /* here R is strictly upper triangular part of STY */
1059: static PetscErrorCode MatGetRTDR(Mat B, Mat result)
1060: {
1061: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1062: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1063: PetscInt m_local;
1065: PetscFunctionBegin;
1066: if (!ldfp->temp_mat) PetscCall(MatDuplicate(ldfp->StY_triu_strict, MAT_SHARE_NONZERO_PATTERN, &ldfp->temp_mat));
1067: PetscCall(MatCopy(ldfp->StY_triu_strict, ldfp->temp_mat, SAME_NONZERO_PATTERN));
1068: PetscCall(MatDiagonalScale(ldfp->temp_mat, ldfp->inv_diag_vec, NULL));
1069: PetscCall(MatGetLocalSize(result, &m_local, NULL));
1070: // need to conjugate and conjugate again because we have MatTransposeMatMult but not MatHermitianTransposeMatMult()
1071: PetscCall(MatConjugate(ldfp->temp_mat));
1072: if (m_local) {
1073: Mat temp_local, StY_local, result_local;
1074: PetscCall(MatDenseGetLocalMatrix(ldfp->StY_triu_strict, &StY_local));
1075: PetscCall(MatDenseGetLocalMatrix(ldfp->temp_mat, &temp_local));
1076: PetscCall(MatDenseGetLocalMatrix(result, &result_local));
1077: PetscCall(MatTransposeMatMult(StY_local, temp_local, MAT_REUSE_MATRIX, PETSC_DETERMINE, &result_local));
1078: }
1079: PetscCall(MatConjugate(result));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: static PetscErrorCode MatLMVMDDFPUpdateSolveData(Mat B)
1084: {
1085: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1086: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1087: PetscInt m = lmvm->m, m_local;
1088: PetscInt k = lmvm->k;
1089: PetscInt h = k - oldest_update(m, k);
1090: PetscInt j_0;
1091: PetscInt prev_oldest;
1092: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
1093: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
1094: Mat J_local;
1096: PetscFunctionBegin;
1097: if (!ldfp->StY_triu_strict) {
1098: PetscCall(MatDuplicate(ldfp->YtS_triu, MAT_SHARE_NONZERO_PATTERN, &ldfp->StY_triu_strict));
1099: PetscCall(MatDestroy(&ldfp->YtHY));
1100: PetscCall(MatDuplicate(ldfp->YtS_triu, MAT_SHARE_NONZERO_PATTERN, &ldfp->YtHY));
1101: PetscCall(MatDestroy(&ldfp->J));
1102: PetscCall(MatDuplicate(ldfp->YtS_triu, MAT_SHARE_NONZERO_PATTERN, &ldfp->J));
1103: PetscCall(MatDestroy(&ldfp->HY));
1104: PetscCall(MatDuplicate(Yfull, MAT_SHARE_NONZERO_PATTERN, &ldfp->HY));
1105: PetscCall(MatShift(ldfp->YtHY, 1.0));
1106: ldfp->num_mult_updates = oldest_update(m, k);
1107: }
1108: if (ldfp->num_mult_updates == k) PetscFunctionReturn(PETSC_SUCCESS);
1110: /* H_0 may have been updated, we must recompute H_0 Y and Y^T H_0 Y */
1111: for (PetscInt j = oldest_update(m, k); j < k; j++) {
1112: Vec y_j;
1113: Vec Hy_j;
1114: Vec YtHy_j;
1115: PetscInt Y_idx = recycle_index(m, j);
1116: PetscInt YtHY_idx = ldfp->strategy == MAT_LMVM_DENSE_INPLACE ? Y_idx : history_index(m, k, j);
1118: PetscCall(MatDenseGetColumnVecWrite(ldfp->HY, Y_idx, &Hy_j));
1119: PetscCall(MatDenseGetColumnVecRead(Yfull, Y_idx, &y_j));
1120: PetscCall(MatDQNApplyJ0Inv(B, y_j, Hy_j));
1121: PetscCall(MatDenseRestoreColumnVecRead(Yfull, Y_idx, &y_j));
1122: PetscCall(MatDenseGetColumnVecWrite(ldfp->YtHY, YtHY_idx, &YtHy_j));
1123: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, Hy_j, YtHy_j, 0, h));
1124: ldfp->Yt_count++;
1125: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, YtHy_j, ldfp->num_updates, ldfp->cyclic_work_vec));
1126: PetscCall(MatDenseRestoreColumnVecWrite(ldfp->YtHY, YtHY_idx, &YtHy_j));
1127: PetscCall(MatDenseRestoreColumnVecWrite(ldfp->HY, Y_idx, &Hy_j));
1128: }
1129: prev_oldest = oldest_update(m, ldfp->num_mult_updates);
1130: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER && prev_oldest < oldest_update(m, k)) {
1131: /* move the YtS entries that have been computed and need to be kept back up */
1132: PetscInt m_keep = m - (oldest_update(m, k) - prev_oldest);
1134: PetscCall(MatMove_LR3(B, ldfp->StY_triu_strict, m_keep));
1135: }
1136: PetscCall(MatGetLocalSize(ldfp->StY_triu_strict, &m_local, NULL));
1137: j_0 = PetscMax(ldfp->num_mult_updates, oldest_update(m, k));
1138: for (PetscInt j = j_0; j < k; j++) {
1139: PetscInt Y_idx = recycle_index(m, j);
1140: PetscInt StY_idx = ldfp->strategy == MAT_LMVM_DENSE_INPLACE ? Y_idx : history_index(m, k, j);
1141: Vec y_j, Sty_j;
1143: PetscCall(MatDenseGetColumnVecRead(Yfull, Y_idx, &y_j));
1144: PetscCall(MatDenseGetColumnVecWrite(ldfp->StY_triu_strict, StY_idx, &Sty_j));
1145: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, y_j, Sty_j, 0, h));
1146: ldfp->St_count++;
1147: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, Sty_j, ldfp->num_updates, ldfp->cyclic_work_vec));
1148: PetscCall(MatDenseRestoreColumnVecWrite(ldfp->StY_triu_strict, StY_idx, &Sty_j));
1149: PetscCall(MatDenseRestoreColumnVecRead(Yfull, Y_idx, &y_j));
1150: /* zero the corresponding row */
1151: if (m_local > 0) {
1152: Mat StY_local, StY_row;
1154: PetscCall(MatDenseGetLocalMatrix(ldfp->StY_triu_strict, &StY_local));
1155: PetscCall(MatDenseGetSubMatrix(StY_local, StY_idx, StY_idx + 1, PETSC_DECIDE, PETSC_DECIDE, &StY_row));
1156: PetscCall(MatZeroEntries(StY_row));
1157: PetscCall(MatDenseRestoreSubMatrix(StY_local, &StY_row));
1158: }
1159: }
1160: if (!ldfp->inv_diag_vec) PetscCall(VecDuplicate(ldfp->diag_vec, &ldfp->inv_diag_vec));
1161: PetscCall(VecCopy(ldfp->diag_vec, ldfp->inv_diag_vec));
1162: PetscCall(VecReciprocal(ldfp->inv_diag_vec));
1163: PetscCall(MatDenseGetLocalMatrix(ldfp->J, &J_local));
1164: PetscCall(MatSetFactorType(J_local, MAT_FACTOR_NONE));
1165: PetscCall(MatGetRTDR(B, ldfp->J));
1166: PetscCall(MatAXPY(ldfp->J, 1.0, ldfp->YtHY, SAME_NONZERO_PATTERN));
1167: if (m_local) {
1168: PetscCall(MatSetOption(J_local, MAT_SPD, PETSC_TRUE));
1169: PetscCall(MatCholeskyFactor(J_local, NULL, NULL));
1170: }
1171: ldfp->num_mult_updates = ldfp->num_updates;
1172: PetscFunctionReturn(PETSC_SUCCESS);
1173: }
1175: /* Solves for
1177: H_0 - [ S | H_0 Y] [ -D | R.T ]^-1 [ S^T ]
1178: [-----+-----------] [---------]
1179: [ R | Y^T H_0 Y ] [ Y^T H_0 ]
1181: Above is equivalent to
1183: H_0 - [ S | H_0 Y] [[ I | 0 ][ -D | 0 ][ I | -D^{-1} R^T ]]^-1 [ S^T ]
1184: [[-----------+---][----+---][---+-------------]] [---------]
1185: [[ -R D^{-1} | I ][ 0 | J ][ 0 | I ]] [ Y^T H_0 ]
1187: where J = Y^T H_0 Y + R D^{-1} R.T
1189: becomes
1191: H_0 - [ S | H_0 Y] [ I | D^{-1} R^T ][ -D^{-1} | 0 ][ I | 0 ] [ S^T ]
1192: [---+------------][----------+--------][----------+---] [---------]
1193: [ 0 | I ][ 0 | J^{-1} ][ R D^{-1} | I ] [ Y^T H_0 ]
1195: =
1197: H_0 + [ S | H_0 Y] [ D^{-1} | 0 ][ I | R^T ][ I | 0 ][ I | 0 ] [ S^T ]
1198: [--------+---][---+-----][---+---------][----------+---] [---------]
1199: [ 0 | I ][ 0 | I ][ 0 | -J^{-1} ][ R D^{-1} | I ] [ Y^T H_0 ]
1201: (Note that StY_triu_strict is R)
1202: Byrd, Nocedal, Schnabel 1994
1204: */
1205: static PetscErrorCode MatSolve_LMVMDDFP(Mat H, Vec F, Vec dX)
1206: {
1207: Mat_LMVM *lmvm = (Mat_LMVM *)H->data;
1208: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1209: PetscInt m = lmvm->m;
1210: PetscInt k = lmvm->k;
1211: PetscInt h = k - oldest_update(m, k);
1212: PetscInt idx, i, j, local_n;
1213: PetscInt m_local;
1214: Mat J_local;
1215: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
1216: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
1218: PetscFunctionBegin;
1219: VecCheckSameSize(F, 2, dX, 3);
1220: VecCheckMatCompatible(H, dX, 3, F, 2);
1222: /* Cholesky Version */
1223: /* Start with the B0 term */
1224: PetscCall(MatDQNApplyJ0Inv(H, F, dX));
1225: if (!ldfp->num_updates) PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
1227: if (ldfp->use_recursive) {
1228: PetscDeviceContext dctx;
1229: PetscMemType memtype;
1230: PetscScalar stf, ytx, ytq, yjtqi, sjtyi, *workscalar;
1232: PetscCall(PetscDeviceContextGetCurrentContext(&dctx));
1233: /* Recursive formulation to avoid Cholesky. Not a dense formulation */
1234: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, ldfp->rwork1, 0, h));
1235: ldfp->Yt_count++;
1237: PetscCall(VecGetLocalSize(ldfp->rwork1, &local_n));
1239: PetscInt oldest = oldest_update(m, k);
1241: if (ldfp->needPQ) {
1242: PetscInt oldest = oldest_update(m, k);
1243: for (i = oldest; i < k; ++i) {
1244: idx = recycle_index(m, i);
1245: /* column_work = S[idx] */
1246: PetscCall(MatGetColumnVector(Yfull, ldfp->column_work, idx));
1247: PetscCall(MatDQNApplyJ0Inv(H, ldfp->column_work, ldfp->PQ[idx]));
1248: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, ldfp->column_work, ldfp->rwork3, 0, h));
1249: PetscCall(VecGetArrayAndMemType(ldfp->rwork3, &workscalar, &memtype));
1250: for (j = oldest; j < i; ++j) {
1251: PetscInt idx_j = recycle_index(m, j);
1252: /* Copy sjtyi in device-aware manner */
1253: if (local_n) {
1254: if (PetscMemTypeHost(memtype)) {
1255: sjtyi = workscalar[idx_j];
1256: } else {
1257: PetscCall(PetscDeviceRegisterMemory(&sjtyi, PETSC_MEMTYPE_HOST, 1 * sizeof(sjtyi)));
1258: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
1259: PetscCall(PetscDeviceArrayCopy(dctx, &sjtyi, &workscalar[idx_j], 1));
1260: }
1261: }
1262: PetscCallMPI(MPI_Bcast(&sjtyi, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)H)));
1263: /* column_work2 = Y[j] */
1264: PetscCall(MatGetColumnVector(Yfull, ldfp->column_work2, idx_j));
1265: PetscCall(VecDot(ldfp->PQ[idx], ldfp->column_work2, &yjtqi));
1266: /* column_work2 = Y[j] */
1267: PetscCall(MatGetColumnVector(Sfull, ldfp->column_work2, idx_j));
1268: /* Compute the pure BFGS component of the forward product */
1269: PetscCall(VecAXPBYPCZ(ldfp->PQ[idx], -yjtqi / ldfp->ytq[idx_j], sjtyi / ldfp->yts[idx_j], 1.0, ldfp->PQ[idx_j], ldfp->column_work2));
1270: }
1271: PetscCall(VecDot(ldfp->PQ[idx], ldfp->column_work, &ytq));
1272: ldfp->ytq[idx] = PetscRealPart(ytq);
1273: }
1274: ldfp->needPQ = PETSC_FALSE;
1275: }
1277: PetscCall(VecGetArrayAndMemType(ldfp->rwork1, &workscalar, &memtype));
1278: for (i = oldest; i < k; ++i) {
1279: idx = recycle_index(m, i);
1280: /* Copy stz[i], ytx[i] in device-aware manner */
1281: if (local_n) {
1282: if (PetscMemTypeHost(memtype)) {
1283: stf = workscalar[idx];
1284: } else {
1285: PetscCall(PetscDeviceRegisterMemory(&stf, PETSC_MEMTYPE_HOST, sizeof(stf)));
1286: PetscCall(PetscDeviceRegisterMemory(workscalar, memtype, local_n * sizeof(*workscalar)));
1287: PetscCall(PetscDeviceArrayCopy(dctx, &stf, &workscalar[idx], 1));
1288: }
1289: }
1290: PetscCallMPI(MPI_Bcast(&stf, 1, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)H)));
1291: /* column_work : S[i], column_work2 : Y[i] */
1292: PetscCall(MatGetColumnVector(Sfull, ldfp->column_work, idx));
1293: PetscCall(MatGetColumnVector(Yfull, ldfp->column_work2, idx));
1294: PetscCall(VecDot(dX, ldfp->column_work2, &ytx));
1295: PetscCall(VecAXPBYPCZ(dX, -ytx / ldfp->ytq[idx], stf / ldfp->yts[idx], 1.0, ldfp->PQ[idx], ldfp->column_work));
1296: }
1297: PetscCall(VecRestoreArrayAndMemType(ldfp->rwork1, &workscalar));
1298: } else {
1299: PetscCall(MatLMVMDDFPUpdateSolveData(H));
1300: PetscCall(MatMultHermitianTransposeColumnRange(Sfull, F, ldfp->rwork1, 0, h));
1301: ldfp->St_count++;
1302: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, dX, ldfp->rwork2, 0, h));
1303: ldfp->Yt_count++;
1304: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) {
1305: PetscCall(VecRecycleOrderToHistoryOrder(H, ldfp->rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1306: PetscCall(VecRecycleOrderToHistoryOrder(H, ldfp->rwork2, ldfp->num_updates, ldfp->cyclic_work_vec));
1307: }
1309: PetscCall(VecPointwiseMult(ldfp->rwork3, ldfp->rwork1, ldfp->inv_diag_vec));
1310: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(ldfp->StY_triu_strict));
1311: PetscCall(MatMultTransposeAdd(ldfp->StY_triu_strict, ldfp->rwork3, ldfp->rwork2, ldfp->rwork2));
1312: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(ldfp->StY_triu_strict));
1314: if (!ldfp->rwork2_local) PetscCall(VecCreateLocalVector(ldfp->rwork2, &ldfp->rwork2_local));
1315: if (!ldfp->rwork3_local) PetscCall(VecCreateLocalVector(ldfp->rwork3, &ldfp->rwork3_local));
1316: PetscCall(VecGetLocalVectorRead(ldfp->rwork2, ldfp->rwork2_local));
1317: PetscCall(VecGetLocalVector(ldfp->rwork3, ldfp->rwork3_local));
1318: PetscCall(MatDenseGetLocalMatrix(ldfp->J, &J_local));
1319: PetscCall(VecGetSize(ldfp->rwork2_local, &m_local));
1320: if (m_local) {
1321: Mat J_local;
1323: PetscCall(MatDenseGetLocalMatrix(ldfp->J, &J_local));
1324: PetscCall(MatSolve(J_local, ldfp->rwork2_local, ldfp->rwork3_local));
1325: }
1326: PetscCall(VecRestoreLocalVector(ldfp->rwork3, ldfp->rwork3_local));
1327: PetscCall(VecRestoreLocalVectorRead(ldfp->rwork2, ldfp->rwork2_local));
1328: PetscCall(VecScale(ldfp->rwork3, -1.0));
1330: PetscCall(MatMultAdd(ldfp->StY_triu_strict, ldfp->rwork3, ldfp->rwork1, ldfp->rwork1));
1331: PetscCall(VecPointwiseMult(ldfp->rwork1, ldfp->rwork1, ldfp->inv_diag_vec));
1333: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) {
1334: PetscCall(VecHistoryOrderToRecycleOrder(H, ldfp->rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1335: PetscCall(VecHistoryOrderToRecycleOrder(H, ldfp->rwork3, ldfp->num_updates, ldfp->cyclic_work_vec));
1336: }
1338: PetscCall(MatMultAddColumnRange(Sfull, ldfp->rwork1, dX, dX, 0, h));
1339: ldfp->S_count++;
1340: PetscCall(MatMultAddColumnRange(ldfp->HY, ldfp->rwork3, dX, dX, 0, h));
1341: ldfp->Y_count++;
1342: }
1343: PetscFunctionReturn(PETSC_SUCCESS);
1344: }
1346: /* Solves for
1347: (Theorem 1, Erway, Jain, and Marcia, 2013)
1349: B_0 - [ Y | B_0 S] [ -R^{-T} (D + S^T B_0 S) R^{-1} | R^{-T} ] [ Y^T ]
1350: ---------------------------------+--------] [---------]
1351: [ R^{-1} | 0 ] [ S^T B_0 ]
1353: (Note: R above is right triangular part of YTS)
1354: which becomes,
1356: [ I | -Y L^{-T} ] [ I | 0 ] [ B_0 | 0 ] [ I | S ] [ I ]
1357: [-----+---] [-----+---] [---+---] [-------------]
1358: [ S^T | I ] [ 0 | D ] [ 0 | I ] [ -L^{-1} Y^T ]
1360: (Note: L above is right triangular part of STY)
1362: */
1363: static PetscErrorCode MatMult_LMVMDDFP(Mat B, Vec X, Vec Z)
1364: {
1365: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1366: Mat_DQN *ldfp = (Mat_DQN *)lmvm->ctx;
1367: Vec rwork1 = ldfp->rwork1;
1368: PetscInt m = lmvm->m;
1369: PetscInt k = lmvm->k;
1370: PetscInt h = k - oldest_update(m, k);
1371: Mat Sfull = lmvm->basis[LMBASIS_S]->vecs;
1372: Mat Yfull = lmvm->basis[LMBASIS_Y]->vecs;
1373: PetscObjectState Xstate;
1375: PetscFunctionBegin;
1376: VecCheckSameSize(X, 2, Z, 3);
1377: VecCheckMatCompatible(B, X, 2, Z, 3);
1379: /* DFP Version. Erway, Jain, Marcia, 2013, Theorem 1 */
1380: /* Block Version */
1381: if (!ldfp->num_updates) {
1382: PetscCall(MatDQNApplyJ0Fwd(B, X, Z));
1383: PetscFunctionReturn(PETSC_SUCCESS); /* No updates stored yet */
1384: }
1386: PetscCall(PetscObjectStateGet((PetscObject)X, &Xstate));
1387: PetscCall(MatMultHermitianTransposeColumnRange(Yfull, X, rwork1, 0, h));
1389: /* Reordering rwork1, as STY is in history order, while Y is in recycled order */
1390: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1391: PetscCall(MatUpperTriangularSolveInPlace(B, ldfp->YtS_triu, rwork1, PETSC_FALSE, ldfp->num_updates, ldfp->strategy));
1392: PetscCall(VecScale(rwork1, -1.0));
1393: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1395: PetscCall(VecCopy(X, ldfp->column_work));
1396: PetscCall(MatMultAddColumnRange(Sfull, rwork1, ldfp->column_work, ldfp->column_work, 0, h));
1397: ldfp->S_count++;
1399: PetscCall(VecPointwiseMult(rwork1, ldfp->diag_vec_recycle_order, rwork1));
1400: PetscCall(MatDQNApplyJ0Fwd(B, ldfp->column_work, Z));
1402: PetscCall(MatMultHermitianTransposeAddColumnRange(Sfull, Z, rwork1, rwork1, 0, h));
1403: ldfp->St_count++;
1405: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecRecycleOrderToHistoryOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1406: PetscCall(MatUpperTriangularSolveInPlace(B, ldfp->YtS_triu, rwork1, PETSC_TRUE, ldfp->num_updates, ldfp->strategy));
1407: PetscCall(VecScale(rwork1, -1.0));
1408: if (ldfp->strategy == MAT_LMVM_DENSE_REORDER) PetscCall(VecHistoryOrderToRecycleOrder(B, rwork1, ldfp->num_updates, ldfp->cyclic_work_vec));
1410: PetscCall(MatMultAddColumnRange(Yfull, rwork1, Z, Z, 0, h));
1411: ldfp->Y_count++;
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: /*
1416: This dense representation reduces the L-DFP update to a series of
1417: matrix-vector products with dense matrices in lieu of the conventional
1418: matrix-free two-loop algorithm.
1419: */
1420: PetscErrorCode MatCreate_LMVMDDFP(Mat B)
1421: {
1422: Mat_LMVM *lmvm;
1423: Mat_DQN *ldfp;
1425: PetscFunctionBegin;
1426: PetscCall(MatCreate_LMVM(B));
1427: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDDFP));
1428: PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
1429: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
1430: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
1431: B->ops->view = MatView_LMVMDQN;
1432: B->ops->setup = MatSetUp_LMVMDQN;
1433: B->ops->setfromoptions = MatSetFromOptions_LMVMDQN;
1434: B->ops->destroy = MatDestroy_LMVMDQN;
1436: lmvm = (Mat_LMVM *)B->data;
1437: lmvm->ops->reset = MatReset_LMVMDQN;
1438: lmvm->ops->update = MatUpdate_LMVMDQN;
1439: lmvm->ops->mult = MatMult_LMVMDDFP;
1440: lmvm->ops->solve = MatSolve_LMVMDDFP;
1441: lmvm->ops->copy = MatCopy_LMVMDQN;
1443: lmvm->ops->multht = lmvm->ops->mult;
1444: lmvm->ops->solveht = lmvm->ops->solve;
1446: PetscCall(PetscNew(&ldfp));
1447: lmvm->ctx = (void *)ldfp;
1448: ldfp->allocated = PETSC_FALSE;
1449: ldfp->watchdog = 0;
1450: ldfp->max_seq_rejects = lmvm->m / 2;
1451: ldfp->strategy = MAT_LMVM_DENSE_INPLACE;
1452: ldfp->use_recursive = PETSC_TRUE;
1453: ldfp->needPQ = PETSC_TRUE;
1455: PetscCall(SymBroydenRescaleCreate(&ldfp->rescale));
1456: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_LMVMDQN));
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@
1461: MatCreateLMVMDDFP - Creates a dense representation of the limited-memory
1462: Davidon-Fletcher-Powell (DFP) approximation to a Hessian.
1464: Collective
1466: Input Parameters:
1467: + comm - MPI communicator
1468: . n - number of local rows for storage vectors
1469: - N - global size of the storage vectors
1471: Output Parameter:
1472: . B - the matrix
1474: Level: advanced
1476: Note:
1477: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
1478: paradigm instead of this routine directly.
1480: .seealso: `MatCreate()`, `MATLMVM`, `MATLMVMDDFP`, `MatCreateLMVMDFP()`
1481: @*/
1482: PetscErrorCode MatCreateLMVMDDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
1483: {
1484: PetscFunctionBegin;
1485: PetscCall(KSPInitializePackage());
1486: PetscCall(MatCreate(comm, B));
1487: PetscCall(MatSetSizes(*B, n, n, N, N));
1488: PetscCall(MatSetType(*B, MATLMVMDDFP));
1489: PetscCall(MatSetUp(*B));
1490: PetscFunctionReturn(PETSC_SUCCESS);
1491: }
1493: /*@
1494: MatLMVMDenseSetType - Sets the memory storage type for dense `MATLMVM`
1496: Input Parameters:
1497: + B - the `MATLMVM` matrix
1498: - type - scale type, see `MatLMVMDenseSetType`
1500: Options Database Keys:
1501: + -mat_lqn_type (reorder|inplace) - set the strategy
1502: . -mat_lbfgs_type (reorder|inplace) - set the strategy
1503: - -mat_ldfp_type (reorder|inplace) - set the strategy
1505: Level: intermediate
1507: MatLMVMDenseTypes\:
1508: + `MAT_LMVM_DENSE_REORDER` - reorders memory to minimize kernel launch
1509: - `MAT_LMVM_DENSE_INPLACE` - launches kernel inplace to minimize memory movement
1511: .seealso: [](ch_ksp), `MATLMVMDQN`, `MATLMVMDBFGS`, `MATLMVMDDFP`, `MatLMVMDenseType`
1512: @*/
1513: PetscErrorCode MatLMVMDenseSetType(Mat B, MatLMVMDenseType type)
1514: {
1515: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1516: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
1518: PetscFunctionBegin;
1520: lqn->strategy = type;
1521: PetscFunctionReturn(PETSC_SUCCESS);
1522: }