1 /* MIT License 2 * 3 * Copyright (c) 2016--2017 Felix Lenders 4 * 5 * Permission is hereby granted, free of charge, to any person obtaining a copy 6 * of this software and associated documentation files (the "Software"), to deal 7 * in the Software without restriction, including without limitation the rights 8 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 9 * copies of the Software, and to permit persons to whom the Software is 10 * furnished to do so, subject to the following conditions: 11 * 12 * The above copyright notice and this permission notice shall be included in all 13 * copies or substantial portions of the Software. 14 * 15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 16 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 17 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 18 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 19 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 20 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 21 * SOFTWARE. 22 * 23 */ 24 25 #ifndef TRLIB_KRYLOV_H 26 #define TRLIB_KRYLOV_H 27 28 #define TRLIB_CLR_CONV_BOUND (0) 29 #define TRLIB_CLR_CONV_INTERIOR (1) 30 #define TRLIB_CLR_APPROX_HARD (2) 31 #define TRLIB_CLR_NEWTON_BREAK (3) 32 #define TRLIB_CLR_HARD_INIT_LAM (4) 33 #define TRLIB_CLR_FAIL_HARD (5) 34 #define TRLIB_CLR_UNBDBEL (6) 35 #define TRLIB_CLR_UNLIKE_CONV (7) 36 #define TRLIB_CLR_CONTINUE (10) 37 #define TRLIB_CLR_ITMAX (-1) 38 #define TRLIB_CLR_FAIL_FACTOR (-3) 39 #define TRLIB_CLR_FAIL_LINSOLVE (-4) 40 #define TRLIB_CLR_FAIL_NUMERIC (-5) 41 #define TRLIB_CLR_FAIL_TTR (-7) 42 #define TRLIB_CLR_PCINDEF (-8) 43 #define TRLIB_CLR_UNEXPECT_INT (-9) 44 45 #define TRLIB_CLT_CG (1) 46 #define TRLIB_CLT_L (2) 47 48 #define TRLIB_CLA_TRIVIAL (0) 49 #define TRLIB_CLA_INIT (1) 50 #define TRLIB_CLA_RETRANSF (2) 51 #define TRLIB_CLA_UPDATE_STATIO (3) 52 #define TRLIB_CLA_UPDATE_GRAD (4) 53 #define TRLIB_CLA_UPDATE_DIR (5) 54 #define TRLIB_CLA_NEW_KRYLOV (6) 55 #define TRLIB_CLA_CONV_HARD (7) 56 #define TRLIB_CLA_OBJVAL (8) 57 58 #define TRLIB_CLS_INIT (1) 59 #define TRLIB_CLS_HOTSTART (2) 60 #define TRLIB_CLS_HOTSTART_G (3) 61 #define TRLIB_CLS_HOTSTART_P (4) 62 #define TRLIB_CLS_HOTSTART_R (5) 63 #define TRLIB_CLS_HOTSTART_T (6) 64 #define TRLIB_CLS_VEC_INIT (7) 65 #define TRLIB_CLS_CG_NEW_ITER (8) 66 #define TRLIB_CLS_CG_UPDATE_S (9) 67 #define TRLIB_CLS_CG_UPDATE_GV (10) 68 #define TRLIB_CLS_CG_UPDATE_P (11) 69 #define TRLIB_CLS_LANCZOS_SWT (12) 70 #define TRLIB_CLS_L_UPDATE_P (13) 71 #define TRLIB_CLS_L_CMP_CONV (14) 72 #define TRLIB_CLS_L_CMP_CONV_RT (15) 73 #define TRLIB_CLS_L_CHK_CONV (16) 74 #define TRLIB_CLS_L_NEW_ITER (17) 75 #define TRLIB_CLS_CG_IF_IRBLK_P (18) 76 #define TRLIB_CLS_CG_IF_IRBLK_C (19) 77 #define TRLIB_CLS_CG_IF_IRBLK_N (20) 78 79 #define TRLIB_CLC_NO_EXP_INV (0) 80 #define TRLIB_CLC_EXP_INV_LOC (1) 81 #define TRLIB_CLC_EXP_INV_GLO (2) 82 83 #define TRLIB_CLT_CG_INT (0) 84 #define TRLIB_CLT_CG_BOUND (1) 85 #define TRLIB_CLT_LANCZOS (2) 86 #define TRLIB_CLT_HOTSTART (3) 87 88 /** 89 * Solves trust region subproblem: Computes minimizer to 90 * 91 * :math:`\min \frac 12 \langle s, H s \rangle + \langle g_0, s \rangle` 92 * subject to the trust region constraint :math:`\Vert s \Vert_M \le r`, 93 * where 94 * 95 * - :math:`H` is available via matrix-vector products :math:`v \mapsto Hv`, 96 * - :math:`\Vert s \Vert_M = \sqrt{\langle s, Ms \rangle}` with :math:`M` positive definite, 97 * - :math:`M` is available via matrix-vector products of its inverse, :math:`v \mapsto M^{-1}v`. 98 * 99 * The minimizer is a global minimizer (modulo floating point), 100 * as long as the hard case does not occur. 101 * The hard case is characterized by the fact that the eigenspace 102 * corresponding to the smallest eigenvalue is degenerate. 103 * 104 * In that case a global minimizer in the Krylov subspaces sampled so far is returned, 105 * sampling the complete search space can be forced by setting :c:data:`ctl_invariant` to :c:macro:`TRLIB_CLC_EXP_INV_GLO`, 106 * but if this is desired factorization based method will most likely be much more efficient. 107 * 108 * A preconditioned Conjugate Gradient / Lanczos method is used, 109 * that possibly employs a reduction to a tridiagnoal subproblem 110 * that is solved using Moré-Sorensens method by explicitly 111 * factorizing the tridiagonal matrix, calling :c:func:`trlib_tri_factor_min`. 112 * 113 * The method builds upon algorithm 5.1 in [Gould1999]_, details of the implementation can be found [Lenders2016]_: 114 * 115 * **Convergence** 116 * 117 * The stopping criterion is based on the gradient of the Lagrangian. Convergence in iteration :math:`i` is reported as soon as 118 * 119 * - interior case: :math:`\Vert g_{i+1} \Vert_{M^{-1}} \le \max\{ \texttt{tol}\_\texttt{abs}\_\texttt{i}, \eta_i \, \Vert g_0 \Vert_{M^{-1}} \}`. 120 * - boundary case: :math:`\Vert g_{i+1} \Vert_{M^{-1}} \le \max\{ \texttt{tol}\_\texttt{abs}\_\texttt{b}, \eta_b \, \Vert g_0 \Vert_{M^{-1}} \}`. 121 * - hard case: :c:data:`ctl_invariant` determines if this is the sole stopping criterion or if further invariant subspaces are going to be explored. See the documentation of this option. 122 * 123 * Here 124 * :math:`\eta_i = \begin{cases} \texttt{tol}\_\texttt{rel}\_\texttt{i} & \texttt{tol}\_\texttt{rel}\_\texttt{i} > 0 \\ \min\{ 0.5, \sqrt{\Vert g_{i+1} \Vert_{M^{-1}}} \} & \texttt{tol}\_\texttt{rel}\_\texttt{i} = -1 \\ \min\{ 0.5, \Vert g_{i+1} \Vert_{M^{-1}} \} & \texttt{tol}\_\texttt{rel}\_\texttt{i} = -2 \end{cases},` 125 * :math:`\eta_b = \begin{cases} \texttt{tol}\_\texttt{rel}\_\texttt{b} & \texttt{tol}\_\texttt{rel}\_\texttt{b} > 0 \\ \min\{ 0.5, \sqrt{\Vert g_{i+1} \Vert_{M^{-1}}} \} & \texttt{tol}\_\texttt{rel}\_\texttt{b} = -1 \\ \min\{ 0.5, \Vert g_{i+1} \Vert_{M^{-1}} \} & \texttt{tol}\_\texttt{rel}\_\texttt{b} = -2, \\ \max\{10^{-6}, \min\{ 0.5, \sqrt{\Vert g_{i+1} \Vert_{M^{-1}}} \}\} & \texttt{tol}\_\texttt{rel}\_\texttt{b} = -3 \\ \max\{10^{-6}, \min\{ 0.5, \Vert g_{i+1} \Vert_{M^{-1}}\} \} & \texttt{tol}\_\texttt{rel}\_\texttt{b} = -4 \end{cases}.` 126 * 127 * Choice of :math:`\eta_i` and :math:`\eta_b` depend on the application, the author has good overall experience in unconstrained optimization with :math:`\texttt{tol}\_\texttt{rel}\_\texttt{i} = -2` and :math:`\texttt{tol}\_\texttt{rel}\_\texttt{b} = -3`. 128 * Some remarks to keep in mind when choosing the tolerances: 129 * 130 * - Choosing fixed small values for :math:`\eta_i` and :math:`\eta_b`, for example :math:`\texttt{tol}\_\texttt{rel}\_\texttt{i} = 10^{-8}` and :math:`\texttt{tol}\_\texttt{rel}\_\texttt{b} = 10^{-6}` lead to quick convergence in a NLP algorithm with lowest number of function evaluations but at a possible excessive amount of matrix-vector products as this also solves problems accurately far away from the solution. 131 * - Choosing :math:`\eta \sim O(\sqrt{\Vert g \Vert_{M^{-1}}})` leads to superlinear convergence in a nonlinear programming algorithm. 132 * - Choosing :math:`\eta \sim O(\Vert g \Vert_{M^{-1}})` leads to quadratic convergence in a nonlinear programming algorithm. 133 * - It is questionable whether it makes sense to solve the boundary with very high accuracy, as the final solution of a nonlinear program should be interior. 134 * 135 * 136 * **Calling Scheme** 137 * 138 * This function employs a reverse communication paradigm. 139 * The functions exits whenever there is an action to be performed 140 * by the user as indicated by the :c:data:`action`. 141 * The user should perform this action and continue with the other 142 * values unchanged as long as the return value is positive. 143 * 144 * **User provided storage** 145 * 146 * The user has to manage 5/6 vectors referred by the names 147 * :math:`g`, :math:`g_-`, :math:`v`, :math:`s`, :math:`p`, :math:`\textit{Hp}` 148 * and a matrix :math:`Q` with :c:data:`itmax` columns to store the 149 * Lanczos directions :math:`p_i`. 150 * The user must perform the actions on those as indicated by the return value. 151 * 152 * In the case of trivial preconditioner, :math:`M = I`, it always 153 * holds :math:`v = g` and the vector :math:`v` is not necessary. 154 * 155 * :math:`s` holds the solution candidate. 156 * 157 * Note that the action :math:`s \leftarrow Q h` will only sometimes be requested 158 * in the very final iteration before convergence. If memory and not computational 159 * load is an issue, you may very well save :c:data:`iter`, :c:data:`ityp`, :c:data:`flt1`, :c:data:`flt2` 160 * :c:data:`flt2` instead of :math:`q_j` and when :math:`s \leftarrow Q s` is requested 161 * simultaniously recompute the directions :math:`q_j` and update the direction 162 * :math:`s` once :math:`q_j` has been computed as :math:`s \leftarrow h_j q_j` 163 * with initialization :math:`s \leftarrow 0`. 164 * 165 * Furthermore the user has to provide a integer and floating point workspace, details 166 * are described in :c:data:`iwork` and :c:data:`fwork`. 167 * 168 * **Resolves** 169 * 170 * *Reentry with new radius* 171 * 172 * You can efficiently hotstart from old results if you have a new problem 173 * with *decreased* trust region radius. Just set `status` to :c:macro:`TRLIB_CLS_HOTSTART`. 174 * Furthermore hotstarting with increased trust region radius should be 175 * trivial as you should be able to just increase the radius and take off 176 * the computation from the previous stage. However as this is an untypical 177 * scenario, it has not been tested at all. 178 * 179 * *Reentry to get minimizer of regularized problem* 180 * 181 * You can reenter to compute the solution of a convexified trust region problem. 182 * You may be interested to this once you see the effect of bad conditioning or also as a cheap retry in a trust region algorithm 183 * before continuing to the next point after a discarded step. 184 * 185 * In this case, call the function with `status` set to :c:macro:`TRLIB_CLS_HOTSTART_P`. 186 * 187 * *Reentry to get unconstrained minimizer of constant regularized problem* 188 * 189 * After a successful solve you can compute the norm of the unconstrained minimizer to the problem with regularized 190 * hessian :math:`H + \beta M` (:math:`\beta` large enough such that :math:`H + \beta M` is positive definite) 191 * in the previously expanded Krylov subspace. You will certainly be interested in this if you want to implement 192 * the TRACE algorithm described in [Curtis2016]_. 193 * 194 * In this case, call the function with `status` set to :c:macro:`TRLIB_CLS_HOTSTART_R` and ``radius`` set to :math:`\beta`. 195 * 196 * On exit, the ``obj`` field of :c:data:`fwork` yields the norm of the solution vector. 197 * 198 * If you not only want the norm but also the unconstrained minimizer, you can get it by one step of :c:func:`TRLIB_CLA_RETRANSF`. 199 * However since this is usually not the case, you are not asked to do this in the reverse communication scheme. 200 * 201 * *Reentry to find suitable TRACE regularization parameter* 202 * 203 * For the aforementioned TRACE algorithm, there is also the need to compute :math:`\beta` such that :math:`\sigma_{\text l} \le \frac{\beta}{\Vert s(\beta) \Vert} \le \sigma_{\text u}`, where :math:`\Vert s(\beta) \Vert` denotes the norm of the regularized unconstrained minimizer. 204 * 205 * To get this, set :c:data:`status` to :c:macro:`TRLIB_CLS_HOTSTART_T` and :c:data:`radius` to an initial guess for :math:`\beta`, :c:data:`tol_rel_i` to :math:`\sigma_{\text l}` and :c:data:`tol_rel_b` to :math:`\sigma_{\text u}`. The field ``obj`` of :c:data:`fwork` contains the desired solution norm on exit and :c:data:`flt1` is the desired regularization parameter :math:`\beta`. 206 * 207 * If you not only want the norm but also the unconstrained minimizer, you can get it by one step of :c:macro:`TRLIB_CLA_RETRANSF`. 208 * However since this is usually not the case, you are not asked to do this in the reverse communication scheme. 209 * 210 * *Reentry with new gradient* 211 * 212 * You can efficiently hotstart from old results if you have a new problem 213 * with changed :math:`g_0` where the previous sampled Krylov subspace is 214 * large enough to contain the solution to the new problem, convergence will 215 * *not* be checked: 216 * 217 * - get pointer ``gt`` to negative gradient of tridiagonal problem in `fwork` with :c:func:`trlib_krylov_gt` 218 * - store ``gt[0]`` and overwrite :math:`\texttt{gt} \leftarrow - Q_i^T \texttt{grad}` 219 * - set :c:data:`status` to :c:macro:`TRLIB_CLS_HOTSTART_G` and start reverse communication process 220 * - to reset data for new problem make sure that you restore ``gt[0]`` and set ``gt[1:] = 0`` for those elements previously overwritten 221 * 222 * **Hard case** 223 * 224 * If you want to investigate the problem also if the hard case, you can either sample through the invariant subspaces (:c:data:`ctl_invariant` set to :c:macro:`TRLIB_CLC_EXP_INV_GLO`) or solve the problem with a gradient for which the hard does not occur and then hotstart with the actual gradient. 225 226 * 227 * :param init: set to :c:macro:`TRLIB_CLS_INIT` on first call, to :c:macro:`TRLIB_CLS_HOTSTART` on hotstart with smaller radius and otherwise to ``0``. 228 * :type init: trlib_int_t, input 229 * :param radius: trust region radius 230 * :type radius: trlib_flt_t, input 231 * :param equality: set to ``1`` if trust region constraint should be enforced as equality 232 * :type equality: trlib_int_t, input 233 * :param itmax: maximum number of iterations 234 * :type itmax: trlib_int_t, input 235 * :param itmax_lanczos: maximum number of Lanczos type iterations. 236 * You are strongly advised to set this to a small number (say 25) unless you know better. 237 * Keep in mind that Lanczos type iteration are only performed when curvature 238 * information is flat and Lanczos may amplify rounding errors without 239 * reorthogonalization. If you allow a lot of Lanczos type iterations 240 * consider reorthogonalizing the new direction against the vector storage. 241 * 242 * :type itmax_lanczos: trlib_int_t, input 243 * :param tol_rel_i: relative stopping tolerance for interior solution 244 * :type tol_rel_i: trlib_flt_t, input 245 * :param tol_abs_i: absolute stopping tolerance for interior solution 246 * :type tol_abs_i: trlib_flt_t, input 247 * :param tol_rel_b: relative stopping tolerance for boundary solution 248 * :type tol_rel_b: trlib_flt_t, input 249 * :param tol_abs_b: absolute stopping tolerance for boundary solution 250 * :type tol_abs_b: trlib_flt_t, input 251 * :param zero: threshold that determines when :math:`\langle p, Hp \rangle` is considered to be zero and thus control eventual Lanczos switch 252 * :type zero: trlib_flt_t, input 253 * :param obj_lo: lower bound on objective, returns if a point is found with function value :math:`\le \texttt{obj}\_\texttt{lo}`. Set to a positive value to ignore this bound. 254 * :type obj_lo: trlib_flt_t, input 255 * :param ctl_invariant: 256 * 257 * - set to :c:macro:`TRLIB_CLC_NO_EXP_INV` if you want to search only in the first invariant Krylov subspace 258 * - set to :c:macro:`TRLIB_CLC_EXP_INV_LOC` if you want to continue to expand the Krylov subspaces but terminate if there is convergence indication in the subspaces sampled so far. 259 * - set to :c:macro:`TRLIB_CLC_EXP_INV_GLO` if you want to continue to expand the Krylov subspaces even in the case of convergence to a local minimizer within in the subspaces sampled so far. Reverse communication continues as long as ctl_invariant is set to :c:macro:`TRLIB_CLC_EXP_INV_GLO`, so you should reset :c:data:`ctl_invariant` to either :c:macro:`TRLIB_CLC_EXP_INV_LOC` or :c:macro:`TRLIB_CLC_EXP_INV_LOC` if there is no reason to continue, for example because you cannot find a new nonzero random vector orthogonal to the sampled directions if :c:data:`action` is :c:macro:`TRLIB_CLA_NEW_KRYLOV`. 260 * 261 * :type ctl_invariant: trlib_int_t, input 262 * :param convexify: set to `1` if you like to monitor if the tridiagonal solution and the backtransformed solution match and if not resolve with a convexified problem, else set to `0` 263 * :type convexify: trlib_int_t, input 264 * :param earlyterm: set to `1` if you like to terminate in the boundary case if it unlikely that much progress will be made fast but no convergence is reached, else set to `0` 265 * :type earlyterm: trlib_int_t, input 266 * :param g_dot_g: dot product :math:`\langle g, g \rangle` 267 * :type g_dot_g: trlib_flt_t, input 268 * :param v_dot_g: dot product :math:`\langle v, g \rangle` 269 * :type v_dot_g: trlib_flt_t, input 270 * :param p_dot_Hp: dot product :math:`\langle p, Hp \rangle` 271 * :type p_dot_Hp: trlib_flt_t, input 272 * :param iwork: integer workspace for problem, internal memory layout described in table below 273 * 274 * - on first call, provide allocated memory for :c:data:`iwork_size` entries provided by :c:func:`trlib_krylov_memory_size` 275 * - leave untouched between iterations 276 * 277 * ====== ================= ========================== 278 * start end(exclusive) description 279 * ====== ================= ========================== 280 * 0 1 internal status flag 281 * 1 2 iteration counter 282 * 2 3 flag that indicates if :math:`H` is positive definite in sampled Krylov subspace 283 * 3 4 flag that indicates if interior solution is expected 284 * 4 5 flag that indicates if warmstart information to :c:data:`leftmost` is available 285 * 5 6 index to block with smallest leftmost 286 * 6 7 flag that indicates if warmstart information to :math:`\lambda_0` is available 287 * 7 8 flag that indicates if warmstart information to :math:`\lambda` is available 288 * 8 9 iteration in which switched to Lanczos iteration, ``-1``: no switch occurred 289 * 9 10 return code from :c:func:`trlib_tri_factor_min` 290 * 10 11 :c:data:`sub_fail` exit code from subrotines called in :c:func:`trlib_tri_factor_min` 291 * 11 12 number of newton iterations needed in :c:func:`trlib_tri_factor_min` 292 * 12 13 last iteration in which a headline has been printed 293 * 13 14 kind of iteration headline that has been printed last 294 * 14 15 status variable if convexified version is going to be resolved 295 * 15 16 number of irreducible blocks 296 * 16 17 + ``itmax`` decomposition into irredecubile block, :c:data:`irblk` for :c:func:`trlib_tri_factor_min` 297 * ====== ================= ========================== 298 * 299 * :type iwork: trlib_int_t, input/output 300 * :param fwork: floating point workspace for problem, internal memory layout described in table below 301 * 302 * - on very first call, provide allocated memory for at least :c:data:`fwork_size` entries that has been initialized using :c:func:`trlib_prepare_memory`, :c:data:`fwork_size` can be obtained by calling :c:func:`trlib_krylov_memory_size` 303 * - can be used for different problem instances with matching dimensions and itmax without reinitialization 304 * - leave untouched between iterations 305 * 306 * ======================== ================================= ============================================= 307 * start end (exclusive) description 308 * ======================== ================================= ============================================= 309 * 0 1 stopping tolerance in case of interior solution 310 * 1 2 stopping tolerance in case of boundary solution 311 * 2 3 dot product :math:`\langle v, g \rangle` in current iteration 312 * 3 4 dot product :math:`\langle p, Hp \rangle` in current iteration 313 * 4 5 ratio between projected CG gradient and Lanczos direction in current iteration 314 * 5 6 ratio between projected CG gradient and Lanczos direction in previous iteration 315 * 6 7 Lagrange multiplier :math:`\lambda_0` for trust region constraint 316 * 7 8 Lagrange multiplier :math:`\lambda` for trust region constraint 317 * 8 9 objective function value in current iterate 318 * 9 10 :math:`\langle s_i, Mp_i \rangle` 319 * 10 11 :math:`\langle p_i, Mp_i \rangle` 320 * 11 12 :math:`\langle s_i, Ms_i \rangle` 321 * 12 13 :math:`\sigma_i` 322 * 13 14 max Rayleigh quotient, :math:`\max_i \frac{\langle p, Hp \rangle}{\langle p, M p \rangle}` 323 * 14 15 min Rayleigh quotient, :math:`\min_i \frac{\langle p, Hp \rangle}{\langle p, M p \rangle}` 324 * 15 16 + :c:data:`itmax` :math:`\alpha_i, i \ge 0`, step length CG 325 * 16 + :c:data:`itmax` 17 + 2 :c:data:`itmax` :math:`\beta_i, i \ge 0`, step update factor CG 326 * 17 + 2 :c:data:`itmax` 18 + 3 :c:data:`itmax` :c:data:`neglin` for :c:func:`trlib_tri_factor_min`, just given by :math:`- \gamma_0 e_1` 327 * 18 + 3 :c:data:`itmax` 19 + 4 :c:data:`itmax` solution :math:`h_0` of tridiagonal subproblem provided as :c:data:`sol` by :c:func:`trlib_tri_factor_min` 328 * 19 + 4 :c:data:`itmax` 20 + 5 :c:data:`itmax` solution :math:`h` of tridiagonal subproblem provided as :c:data:`sol` by :c:func:`trlib_tri_factor_min` 329 * 20 + 5 :c:data:`itmax` 21 + 6 :c:data:`itmax` :math:`\delta_i, i \ge 0`, curvature in Lanczos, diagonal of :math:`T` in Lanczos tridiagonalization process 330 * 21 + 6 :c:data:`itmax` 22 + 7 :c:data:`itmax` diagonal of Cholesky of :math:`T_0 + \lambda_0 I` 331 * 22 + 7 :c:data:`itmax` 23 + 8 :c:data:`itmax` diagonal of Cholesky of :math:`T + \lambda I` 332 * 23 + 8 :c:data:`itmax` 23 + 9 :c:data:`itmax` :math:`\gamma_i, i \ge 1`, norm of gradients in Lanczos; provides offdiagonal of :math:`T` in Lanczos tridiagonalization process 333 * 23 + 9 :c:data:`itmax` 23 + 10 :c:data:`itmax` offdiagonal of Cholesky factorization of :math:`T_0 + \lambda_0 I` 334 * 23 + 10 :c:data:`itmax` 23 + 11 :c:data:`itmax` offdiagonal of Cholesky factorization of :math:`T + \lambda I` 335 * 23 + 11 :c:data:`itmax` 24 + 12 :c:data:`itmax` :c:data:`ones` for :c:func:`trlib_tri_factor_min` and :c:func:`trlib_eigen_inverse` 336 * 24 + 12 :c:data:`itmax` 25 + 13 :c:data:`itmax` :c:data:`leftmost` for :c:func:`trlib_tri_factor_min` 337 * 25 + 13 :c:data:`itmax` 26 + 14 :c:data:`itmax` :c:data:`regdiag` for :c:func:`trlib_tri_factor_regularize_posdef` 338 * 26 + 14 :c:data:`itmax` 27 + 15 :c:data:`itmax` history of convergence criteria values 339 * 27 + 15 :c:data:`itmax` 27 + 15 :c:data:`itmax` + ``fmz`` :c:data:`fwork` for :c:func:`trlib_tri_factor_min`, ``fmz`` is given by :c:func:`trlib_tri_factor_memory_size` with argument ``itmax+1`` 340 * ======================== ================================= ============================================= 341 * 342 * 343 * :type fwork: trlib_flt_t, input/output 344 * :param refine: set to ``1`` if iterative refinement should be used on solving linear systems, otherwise to ``0`` 345 * :type refine: trlib_int_t, input 346 * :param verbose: determines the verbosity level of output that is written to :c:data:`fout` 347 * :type verbose: trlib_int_t, input 348 * :param unicode: set to ``1`` if :c:data:`fout` can handle unicode, otherwise to ``0`` 349 * :type unicode: trlib_int_t, input 350 * :param prefix: string that is printed before iteration output 351 * :type prefix: char, input 352 * :param fout: output stream 353 * :type fout: FILE, input 354 * :param timing: gives timing details, all values are multiples of nanoseconds, provide zero allocated memory of length :c:func:`trlib_krylov_timing_size` 355 * 356 * ====== ============ 357 * block description 358 * ====== ============ 359 * 0 total duration 360 * 1 :c:data:`timing` of :c:func:`trlib_tri_factor_min` 361 * ====== ============ 362 * 363 * :type timing: trlib_int_t, input/output 364 * :param action: The user should perform the following action depending on :c:data:`action` and :c:data:`ityp` on the vectors he manages, see the table below. 365 * The table makes use of the notation explained in the *User provided storage* above and the following: 366 * 367 * - :math:`i`: :c:data:`iter` 368 * - :math:`q_j`: :math:`j`-th column of :math:`Q` 369 * - :math:`Q_i`: matrix consisting of the first :math:`i+1` columns of :math:`Q`, :math:`Q_i = (q_0, \ldots, q_i)` 370 * - :math:`h_i`: vector of length :math:`i+1` stored in :c:data:`fwork` at start position :c:data:`h_pointer` provided by :c:func:`trlib_krylov_memory_size` 371 * - :math:`p \leftarrow \perp_M Q_j`: optionally :math:`M`-reorthonormalize :math:`p` against :math:`Q_j` 372 * - :math:`g \leftarrow \texttt{rand}\perp Q_j` find nonzero random vector that is orthogonal to :math:`Q_j` 373 * - Note that :c:macro:`TRLIB_CLA_NEW_KRYLOV` is unlikely and only occurs on problems that employ the hard case and only if :c:data:`ctl_invariant` :math:`\neq` :c:macro:`TRLIB_CLC_NO_EXP_INV`. 374 * If you want to safe yourself from the trouble implementing this and are confident that you don't need to expand several invariant Krylov subspaces, just ensure that :c:data:`ctl_invariant` is set to :c:macro:`TRLIB_CLC_NO_EXP_INV`. 375 * 376 * =================================== ================================================= ============================= 377 * action ityp command 378 * =================================== ================================================= ============================= 379 * :c:macro:`TRLIB_CLA_TRIVIAL` :c:macro:`TRLIB_CLT_CG`, :c:macro:`TRLIB_CLT_L` do nothing 380 * :c:macro:`TRLIB_CLA_RETRANSF` :c:macro:`TRLIB_CLT_CG`, :c:macro:`TRLIB_CLT_L` :math:`s \leftarrow Q_i h_i` 381 * :c:macro:`TRLIB_CLA_INIT` :c:macro:`TRLIB_CLT_CG`, :c:macro:`TRLIB_CLT_L` :math:`s \leftarrow 0`, :math:`g \leftarrow g_0`, :math:`g_- \leftarrow 0`, :math:`v \leftarrow M^{-1} g`, :math:`p \leftarrow -v`, :math:`\textit{Hp} \leftarrow Hp`, :math:`\texttt{g}\_\texttt{dot}\_\texttt{g} \leftarrow \langle g, g \rangle`, :math:`\texttt{v}\_\texttt{dot}\_\texttt{g} \leftarrow \langle v, g \rangle`, :math:`\texttt{p}\_\texttt{dot}\_\texttt{Hp} \leftarrow \langle p, \textit{Hp} \rangle`, :math:`q_0 \leftarrow \frac{1}{\sqrt{\texttt{v}\_\texttt{dot}\_\texttt{g}}} v` 382 * :c:macro:`TRLIB_CLA_UPDATE_STATIO` :c:macro:`TRLIB_CLT_CG` :math:`s \leftarrow s + \texttt{flt1} \, p` 383 * :c:macro:`TRLIB_CLA_UPDATE_STATIO` :c:macro:`TRLIB_CLT_L` do nothing 384 * :c:macro:`TRLIB_CLA_UPDATE_GRAD` :c:macro:`TRLIB_CLT_CG` :math:`q_i \leftarrow \texttt{flt2} \, v`, :math:`g_- \leftarrow g`, :math:`g \leftarrow g + \texttt{flt1} \, \textit{Hp}`, :math:`v \leftarrow M^{-1} g`, :math:`\texttt{g}\_\texttt{dot}\_\texttt{g} \leftarrow \langle g, g \rangle`, :math:`\texttt{v}\_\texttt{dot}\_\texttt{g} \leftarrow \langle v, g \rangle` 385 * :c:macro:`TRLIB_CLA_UPDATE_GRAD` :c:macro:`TRLIB_CLT_L` :math:`s \leftarrow \textit{Hp} + \texttt{flt1}\, g + \texttt{flt2}\, g_-`, :math:`g_- \leftarrow \texttt{flt3}\, g`, :math:`g \leftarrow s`, :math:`v \leftarrow M^{-1} g`, :math:`\texttt{v}\_\texttt{dot}\_\texttt{g} \leftarrow \langle v, g \rangle` 386 * :c:macro:`TRLIB_CLA_UPDATE_DIR` :c:macro:`TRLIB_CLT_CG` :math:`p \leftarrow \texttt{flt1} \, v + \texttt{flt2} \, p` with :math:`\texttt{flt1} = -1`, :math:`\textit{Hp} \leftarrow Hp`, :math:`\texttt{p}\_\texttt{dot}\_\texttt{Hp} \leftarrow \langle p, \textit{Hp} \rangle` 387 * :c:macro:`TRLIB_CLA_UPDATE_DIR` :c:macro:`TRLIB_CLT_L` :math:`p \leftarrow \texttt{flt1} \, v + \texttt{flt2} \, p` with :math:`\texttt{flt2} = 0`, :math:`p \leftarrow \perp_M Q_{i-1}`, :math:`\textit{Hp} \leftarrow Hp`, :math:`\texttt{p}\_\texttt{dot}\_\texttt{Hp} \leftarrow \langle p, \textit{Hp} \rangle`, :math:`q_i \leftarrow p` 388 * :c:macro:`TRLIB_CLA_NEW_KRYLOV` :c:macro:`TRLIB_CLT_CG`, :c:macro:`TRLIB_CLT_L` :math:`g \leftarrow \texttt{rand}\perp Q_{i-1}`, :math:`g_- \leftarrow 0`, :math:`v \leftarrow M^{-1} g`, :math:`\texttt{v}\_\texttt{dot}\_\texttt{g} \leftarrow \langle v, g \rangle`, :math:`p \leftarrow \frac{1}{\sqrt{\texttt{v}\_\texttt{dot}\_\texttt{g}}} v`, :math:`\texttt{p}\_\texttt{dot}\_\texttt{Hp} \leftarrow \langle p, \textit{Hp} \rangle`, :math:`q_{i+1} \leftarrow p` 389 * :c:macro:`TRLIB_CLA_CONV_HARD` :c:macro:`TRLIB_CLT_CG`, :c:macro:`TRLIB_CLT_L` :math:`\texttt{v}\_\texttt{dot}\_\texttt{g} \leftarrow \langle Hs+g_0+\texttt{flt1}\, Ms, M^{-1}(Hs+g_0) + \texttt{flt1} s \rangle` 390 * :c:macro:`TRLIB_CLA_OBJVAL` :c:macro:`TRLIB_CLT_CG`, :c:macro:`TRLIB_CLT_L` :math:`\texttt{g}\_\texttt{dot}\_\texttt{g} \leftarrow \tfrac 12 \langle s, Hs \rangle + \langle g, s \rangle` 391 * =================================== ================================================= ============================= 392 * 393 * :type action: trlib_int_t, output 394 * :param iter: iteration counter to tell user position in vector storage 395 * :type iter: trlib_int_t, output 396 * :param ityp: iteration type, see :c:data:`action` 397 * :type ityp: trlib_int_t, output 398 * :param flt1: floating point value that user needs for actions 399 * :type flt1: trlib_flt_t, output 400 * :param flt2: floating point value that user needs for actions 401 * :type flt2: trlib_flt_t, output 402 * :param flt3: floating point value that user needs for actions 403 * :type flt3: trlib_flt_t, output 404 * 405 * :returns: status flag with following meaning: 406 * 407 * - :c:macro:`TRLIB_CLR_CONTINUE` no convergence yet, continue in reverse communication 408 * - :c:macro:`TRLIB_CLR_CONV_BOUND` successful exit with converged solution on boundary, end reverse communication process 409 * - :c:macro:`TRLIB_CLR_CONV_INTERIOR` successful exit with converged interior solution, end reverse communication process 410 * - :c:macro:`TRLIB_CLR_APPROX_HARD` successful exit with approximate solution, hard case occurred, end reverse communication process 411 * - :c:macro:`TRLIB_CLR_NEWTON_BREAK` exit with breakdown in Newton iteration in :c:func:`trlib_tri_factor_min`, most likely converged to boundary solution 412 * - :c:macro:`TRLIB_TTR_HARD_INIT_LAM` hard case encountered without being able to find suitable initial :math:`\lambda` for Newton iteration, returned approximate stationary point that maybe suboptimal 413 * - :c:macro:`TRLIB_CLR_ITMAX` iteration limit exceeded, end reverse communication process 414 * - :c:macro:`TRLIB_CLR_UNBDBEL` problem seems to be unbounded from below, end reverse communication process 415 * - :c:macro:`TRLIB_CLR_FAIL_FACTOR` failure on factorization, end reverse communication process 416 * - :c:macro:`TRLIB_CLR_FAIL_LINSOLVE` failure on backsolve, end reverse communication process 417 * - :c:macro:`TRLIB_CLR_FAIL_NUMERIC` failure as one of the values :c:data:`v_dot_g` or :c:data:`p_dot_Hp` are not a floating point number 418 * - :c:macro:`TRLIB_CLR_UNLIKE_CONV` exit early as convergence seems to be unlikely 419 * - :c:macro:`TRLIB_CLR_PCINDEF` preconditioner apperas to be indefinite, end reverse communication process 420 * - :c:macro:`TRLIB_CLR_UNEXPECT_INT` unexpected interior solution found, expected boundary solution, end reverse communication process 421 * - :c:macro:`TRLIB_CLR_FAIL_TTR` failure occurred in :c:func:`trlib_tri_factor_min`, check :c:data:`iwork[7]` and :c:data:`iwork[8]` for details 422 * - :c:macro:`TRLIB_CLR_FAIL_HARD` failure due to occurrence of hard case: invariant subspace encountered without local convergence and request for early termination without exploring further invariant subspaces 423 * 424 * :rtype: trlib_int_t 425 * 426 * 427 */ 428 429 trlib_int_t trlib_krylov_min( 430 trlib_int_t init, trlib_flt_t radius, trlib_int_t equality, trlib_int_t itmax, trlib_int_t itmax_lanczos, 431 trlib_flt_t tol_rel_i, trlib_flt_t tol_abs_i, 432 trlib_flt_t tol_rel_b, trlib_flt_t tol_abs_b, trlib_flt_t zero, trlib_flt_t obj_lo, 433 trlib_int_t ctl_invariant, trlib_int_t convexify, trlib_int_t earlyterm, 434 trlib_flt_t g_dot_g, trlib_flt_t v_dot_g, trlib_flt_t p_dot_Hp, 435 trlib_int_t *iwork, trlib_flt_t *fwork, trlib_int_t refine, 436 trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout, trlib_int_t *timing, 437 trlib_int_t *action, trlib_int_t *iter, trlib_int_t *ityp, 438 trlib_flt_t *flt1, trlib_flt_t *flt2, trlib_flt_t *flt3 439 ); 440 441 /** Prepares floating point workspace for :c:func::`trlib_krylov_min` 442 * 443 * Initializes floating point workspace :c:data:`fwork` for :c:func:`trlib_krylov_min` 444 * 445 * :param itmax: maximum number of iterations 446 * :type itmax: trlib_int_t, input 447 * :param fwork: floating point workspace to be used by :c:func:`trlib_krylov_min` 448 * :type fwork: trlib_flt_t, input/output 449 * 450 * :returns: ``0`` 451 * :rtype: trlib_int_t 452 */ 453 454 trlib_int_t trlib_krylov_prepare_memory(trlib_int_t itmax, trlib_flt_t *fwork); 455 456 /** Gives information on memory that has to be allocated for :c:func::`trlib_krylov_min` 457 * 458 * :param itmax: maximum number of iterations 459 * :type itmax: trlib_int_t, input 460 * :param iwork_size: size of integer workspace iwork that has to be allocated for :c:func:`trlib_krylov_min` 461 * :type iwork_size: trlib_int_t, output 462 * :param fwork_size: size of floating point workspace fwork that has to be allocated for :c:func:`trlib_krylov_min` 463 * :type fwork_size: trlib_int_t, output 464 * :param h_pointer: start index of vector :math:`h` that has to be used in reverse communication in action :c:macro:`TRLIB_CLA_RETRANSF` 465 * :type h_pointer: trlib_int_t, output 466 * 467 * :returns: ``0`` 468 * :rtype: trlib_int_t 469 */ 470 471 trlib_int_t trlib_krylov_memory_size(trlib_int_t itmax, trlib_int_t *iwork_size, trlib_int_t *fwork_size, trlib_int_t *h_pointer); 472 473 /** size that has to be allocated for :c:data:`timing` in :c:func:`trlib_krylov_min` 474 * 475 * :returns: ``0`` 476 * :rtype: trlib_int_t 477 */ 478 trlib_int_t trlib_krylov_timing_size(void); 479 480 /** Gives pointer to negative gradient of tridiagonal problem 481 * 482 * :param itmax: itmax maximum number of iterations 483 * :type itmax: trlib_int_t, input 484 * :param gt_pointer: pointer to negative gradient of tridiagonal subproblem 485 * :type gt_pointer: trlib_int_t, output 486 * 487 * :returns: ``0`` 488 * :rtype: trlib_int_t 489 */ 490 491 trlib_int_t trlib_krylov_gt(trlib_int_t itmax, trlib_int_t *gt_pointer); 492 493 #endif 494