Actual source code: bqnk.c
1: #include <../src/tao/bound/impls/bqnk/bqnk.h>
2: #include <petscksp.h>
4: static PetscErrorCode TaoBQNKComputeHessian(Tao tao)
5: {
6: TAO_BNK *bnk = (TAO_BNK *)tao->data;
7: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
8: PetscReal gnorm2, delta;
10: PetscFunctionBegin;
11: /* Alias the LMVM matrix into the TAO hessian */
12: if (tao->hessian) PetscCall(MatDestroy(&tao->hessian));
13: if (tao->hessian_pre) PetscCall(MatDestroy(&tao->hessian_pre));
14: PetscCall(PetscObjectReference((PetscObject)bqnk->B));
15: tao->hessian = bqnk->B;
16: PetscCall(PetscObjectReference((PetscObject)bqnk->B));
17: tao->hessian_pre = bqnk->B;
18: /* Update the Hessian with the latest solution */
19: if (bqnk->is_spd) {
20: gnorm2 = bnk->gnorm * bnk->gnorm;
21: if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
22: if (bnk->f == 0.0) {
23: delta = 2.0 / gnorm2;
24: } else {
25: delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2;
26: }
27: PetscCall(MatLMVMSymBroydenSetDelta(bqnk->B, delta));
28: }
29: PetscCall(MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient));
30: PetscCall(MatLMVMResetShift(tao->hessian));
31: /* Prepare the reduced sub-matrices for the inactive set */
32: PetscCall(MatDestroy(&bnk->H_inactive));
33: if (bnk->active_idx) {
34: PetscCall(MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive));
35: PetscCall(PCLMVMSetIS(bqnk->pc, bnk->inactive_idx));
36: } else {
37: PetscCall(PetscObjectReference((PetscObject)tao->hessian));
38: bnk->H_inactive = tao->hessian;
39: PetscCall(PCLMVMClearIS(bqnk->pc));
40: }
41: PetscCall(MatDestroy(&bnk->Hpre_inactive));
42: PetscCall(PetscObjectReference((PetscObject)bnk->H_inactive));
43: bnk->Hpre_inactive = bnk->H_inactive;
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type)
48: {
49: TAO_BNK *bnk = (TAO_BNK *)tao->data;
50: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
52: PetscFunctionBegin;
53: PetscCall(TaoBNKComputeStep(tao, shift, ksp_reason, step_type));
54: if (*ksp_reason < 0) {
55: /* Krylov solver failed to converge so reset the LMVM matrix */
56: PetscCall(MatLMVMReset(bqnk->B, PETSC_FALSE));
57: PetscCall(MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient));
58: }
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: PetscErrorCode TaoSolve_BQNK(Tao tao)
63: {
64: TAO_BNK *bnk = (TAO_BNK *)tao->data;
65: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
66: Mat_LMVM *lmvm = (Mat_LMVM *)bqnk->B->data;
67: Mat_LMVM *J0;
68: PetscBool flg = PETSC_FALSE;
70: PetscFunctionBegin;
71: if (!tao->recycle) {
72: PetscCall(MatLMVMReset(bqnk->B, PETSC_FALSE));
73: lmvm->nresets = 0;
74: if (lmvm->J0) {
75: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &flg));
76: if (flg) {
77: J0 = (Mat_LMVM *)lmvm->J0->data;
78: J0->nresets = 0;
79: }
80: }
81: }
82: PetscCall((*bqnk->solve)(tao));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PetscErrorCode TaoSetUp_BQNK(Tao tao)
87: {
88: TAO_BNK *bnk = (TAO_BNK *)tao->data;
89: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
90: PetscInt n, N;
91: PetscBool is_lmvm, is_set, is_sym;
93: PetscFunctionBegin;
94: PetscCall(TaoSetUp_BNK(tao));
95: PetscCall(VecGetLocalSize(tao->solution, &n));
96: PetscCall(VecGetSize(tao->solution, &N));
97: PetscCall(MatSetSizes(bqnk->B, n, n, N, N));
98: PetscCall(MatLMVMAllocate(bqnk->B, tao->solution, bnk->unprojected_gradient));
99: PetscCall(PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm));
100: PetscCheck(is_lmvm, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type");
101: PetscCall(MatIsSymmetricKnown(bqnk->B, &is_set, &is_sym));
102: PetscCheck(is_set && is_sym, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric");
103: PetscCall(KSPGetPC(tao->ksp, &bqnk->pc));
104: PetscCall(PCSetType(bqnk->pc, PCLMVM));
105: PetscCall(PCLMVMSetMatLMVM(bqnk->pc, bqnk->B));
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: static PetscErrorCode TaoSetFromOptions_BQNK(Tao tao, PetscOptionItems PetscOptionsObject)
110: {
111: TAO_BNK *bnk = (TAO_BNK *)tao->data;
112: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
113: PetscBool is_set;
115: PetscFunctionBegin;
116: PetscCall(TaoSetFromOptions_BNK(tao, PetscOptionsObject));
117: if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION;
118: PetscCall(MatSetOptionsPrefix(bqnk->B, ((PetscObject)tao)->prefix));
119: PetscCall(MatAppendOptionsPrefix(bqnk->B, "tao_bqnk_"));
120: PetscCall(MatSetFromOptions(bqnk->B));
121: PetscCall(MatIsSPDKnown(bqnk->B, &is_set, &bqnk->is_spd));
122: if (!is_set) bqnk->is_spd = PETSC_FALSE;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer)
127: {
128: TAO_BNK *bnk = (TAO_BNK *)tao->data;
129: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
130: PetscBool isascii;
132: PetscFunctionBegin;
133: PetscCall(TaoView_BNK(tao, viewer));
134: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
135: if (isascii) {
136: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
137: PetscCall(MatView(bqnk->B, viewer));
138: PetscCall(PetscViewerPopFormat(viewer));
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static PetscErrorCode TaoDestroy_BQNK(Tao tao)
144: {
145: TAO_BNK *bnk = (TAO_BNK *)tao->data;
146: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
148: PetscFunctionBegin;
149: PetscCall(MatDestroy(&bnk->Hpre_inactive));
150: PetscCall(MatDestroy(&bnk->H_inactive));
151: PetscCall(MatDestroy(&bqnk->B));
152: PetscCall(PetscFree(bnk->ctx));
153: PetscCall(TaoDestroy_BNK(tao));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao)
158: {
159: TAO_BNK *bnk;
160: TAO_BQNK *bqnk;
162: PetscFunctionBegin;
163: PetscCall(TaoCreate_BNK(tao));
164: tao->ops->solve = TaoSolve_BQNK;
165: tao->ops->setfromoptions = TaoSetFromOptions_BQNK;
166: tao->ops->destroy = TaoDestroy_BQNK;
167: tao->ops->view = TaoView_BQNK;
168: tao->ops->setup = TaoSetUp_BQNK;
169: tao->uses_hessian_matrices = PETSC_FALSE;
171: bnk = (TAO_BNK *)tao->data;
172: bnk->computehessian = TaoBQNKComputeHessian;
173: bnk->computestep = TaoBQNKComputeStep;
174: bnk->init_type = BNK_INIT_DIRECTION;
176: PetscCall(PetscNew(&bqnk));
177: bnk->ctx = (void *)bqnk;
178: bqnk->is_spd = PETSC_TRUE;
180: PetscCall(MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B));
181: PetscCall(PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1));
182: PetscCall(MatSetType(bqnk->B, MATLMVMSR1));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: TaoGetLMVMMatrix - Returns a pointer to the internal LMVM matrix. Valid
188: only for quasi-Newton family of methods.
190: Input Parameter:
191: . tao - `Tao` solver context
193: Output Parameter:
194: . B - LMVM matrix
196: Level: advanced
198: .seealso: `TAOBQNLS`, `TAOBQNKLS`, `TAOBQNKTL`, `TAOBQNKTR`, `MATLMVM`, `TaoSetLMVMMatrix()`
199: @*/
200: PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B)
201: {
202: TAO_BNK *bnk = (TAO_BNK *)tao->data;
203: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
204: PetscBool flg = PETSC_FALSE;
206: PetscFunctionBegin;
207: PetscCall(PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, ""));
208: PetscCheck(flg, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
209: *B = bqnk->B;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@
214: TaoSetLMVMMatrix - Sets an external LMVM matrix into the Tao solver. Valid
215: only for quasi-Newton family of methods.
217: QN family of methods create their own LMVM matrices and users who wish to
218: manipulate this matrix should use TaoGetLMVMMatrix() instead.
220: Input Parameters:
221: + tao - Tao solver context
222: - B - LMVM matrix
224: Level: advanced
226: .seealso: `TAOBQNLS`, `TAOBQNKLS`, `TAOBQNKTL`, `TAOBQNKTR`, `MATLMVM`, `TaoGetLMVMMatrix()`
227: @*/
228: PetscErrorCode TaoSetLMVMMatrix(Tao tao, Mat B)
229: {
230: TAO_BNK *bnk = (TAO_BNK *)tao->data;
231: TAO_BQNK *bqnk = (TAO_BQNK *)bnk->ctx;
232: PetscBool flg = PETSC_FALSE;
234: PetscFunctionBegin;
235: PetscCall(PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, ""));
236: PetscCheck(flg, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
237: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &flg));
238: PetscCheck(flg, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given matrix is not an LMVM matrix");
239: if (bqnk->B) PetscCall(MatDestroy(&bqnk->B));
240: PetscCall(PetscObjectReference((PetscObject)B));
241: bqnk->B = B;
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }