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_LEFTMOST_H 26 #define TRLIB_LEFTMOST_H 27 28 #define TRLIB_LMR_CONV (0) 29 #define TRLIB_LMR_ITMAX (-1) 30 #define TRLIB_LMR_NEWTON_BREAK (3) 31 32 /** Computes smallest eigenvalue of symmetric tridiagonal matrix 33 * :math:`T \in \mathbb R^{n\times n}`, 34 * using a iteration based on last-pivot function of Parlett and Reid. 35 * 36 * Let :math:`T = \begin{pmatrix} T_1 & & \\ & \ddots & \\ & & T_\ell \end{pmatrix}` 37 * be composed into irreducible blocks :math:`T_i`. 38 * 39 * Calls :c:func:`trlib_leftmost_irreducible` on every irreducible block in case of coldstart, 40 * in case of warmstart just updates information on :math:`T_\ell`. 41 * 42 * :param nirblk: number of irreducible blocks :math:`\ell`, ensure :math:`\ell > 0` 43 * :type nirblk: trlib_int_t, input 44 * :param irblk: pointer to indices of irreducible blocks, length :c:data:`nirblk+1`: 45 * 46 * - :c:data:`irblk[i]` is the start index of block :math:`i` in :c:data:`diag` and :c:data:`offdiag` 47 * - :c:data:`irblk[i+1] - 1` is the stop index of block :math:`i` 48 * - :c:data:`irblk[i+1] - irred[i]` the dimension :math:`n_\ell` of block :math:`i` 49 * - :c:data:`irblk[nirred]` the dimension of :math:`T` 50 * 51 * :type irblk: trlib_int_t, input 52 * :param diag: pointer to array holding diagonal of :math:`T`, length :c:data:`irblk[nirblk]` 53 * :type diag: trlib_flt_t, input 54 * :param offdiag: pointer to array holding offdiagonal of :math:`T`, length :c:data:`irblk[nirblk]` 55 * :type offdiag: trlib_flt_t, input 56 * :param warm: set :math:`\ge 1` if you want to update information about block :math:`\ell`, provide values in :c:data:`leftmost_minor`, :c:data:`ileftmost`, :c:data:`leftmost`; else ``0`` 57 * :type warm: trlib_int_t, input 58 * :param leftmost_minor: smallest eigenvalue of principal :math:`(n_\ell-1)\times (n_\ell-1)` submatrix of :math:`T_\ell` 59 * :type leftmost_minor: trlib_flt_t, input 60 * :param itmax: maximum number of iterations 61 * :type itmax: trlib_int_t, input 62 * :param tol_abs: absolute stopping tolerance in Reid-Parlett zero finding, good default may be :math:`\sqrt{\texttt{macheps}}^{3/4}` (:c:macro:`TRLIB_EPS_POW_75`) 63 * :type tol_abs: trlib_flt_t, input 64 * :param verbose: determines the verbosity level of output that is written to :c:data:`fout` 65 * :type verbose: trlib_int_t, input 66 * :param unicode: set to ``1`` if :c:data:`fout` can handle unicode, otherwise to ``0`` 67 * :type unicode: trlib_int_t, input 68 * :param prefix: string that is printed before iteration output 69 * :type prefix: char, input 70 * :param fout: output stream 71 * :type fout: FILE, input 72 * :param timing: gives timing details, provide allocated zero initialized memory of length :c:func:`trlib_leftmost_timing_size` 73 * 74 * ====== ================================ 75 * block description 76 * ====== ================================ 77 * 0 total duration 78 * ====== ================================ 79 * 80 * :type timing: trlib_int_t, input/output 81 * :param ileftmost: index of block that corresponds to absolute smallest eigenvalue 82 * :type ileftmost: trlib_int_t, input/output 83 * :param leftmost: smallest eigenvalue of :math:`T`, length :math:`\ell` 84 * 85 * - on entry: allocated memory 86 * - on exit: :c:data:`leftmost[i]` smallest eigenvalue of :math:`T_i` 87 * 88 * :type leftmost: trlib_flt_t, input/output 89 * :returns: status 90 * 91 * - :c:macro:`TRLIB_LMR_CONV` success 92 * - :c:macro:`TRLIB_LMR_ITMAX` iteration limit exceeded 93 * 94 * :rtype: trlib_int_t 95 */ 96 97 trlib_int_t trlib_leftmost( 98 trlib_int_t nirblk, trlib_int_t *irblk, trlib_flt_t *diag, trlib_flt_t *offdiag, 99 trlib_int_t warm, trlib_flt_t leftmost_minor, trlib_int_t itmax, trlib_flt_t tol_abs, 100 trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout, 101 trlib_int_t *timing, trlib_int_t *ileftmost, trlib_flt_t *leftmost); 102 103 /** Computes smallest eigenvalue of irreducible symmetric tridiagonal matrix 104 * :math:`T \in \mathbb R^{n\times n}`, 105 * using a iteration based on last-pivot function of Parlett and Reid. 106 * 107 * Method is sketched on p. 516 in [Gould1999]_. 108 * 109 * Note that this function most likely will fail in the case of a reducible matrix 110 * (:c:data:`offdiag` contains 0). 111 * 112 * **Convergence** 113 * 114 * Convergence is reported if :math:`\texttt{up}-\texttt{low} \le \texttt{tol}\_\texttt{abs} * \max\{1, \vert \texttt{low} \vert, \vert \texttt{up} \vert \}` or :math:`\texttt{prlp} \le \texttt{tol}\_\texttt{abs}`, :math:`\texttt{low}` and :math:`\texttt{up}` denote bracket values enclosing the leftmost eigenvalue and :math:`\texttt{prlp}` denotes the last-pivot function value used in root finding. 115 * 116 * :param n: dimension, ensure :math:`n > 0` 117 * :type n: trlib_int_t, input 118 * :param diag: pointer to array holding diagonal of :math:`T`, length :c:data:`n` 119 * :type diag: trlib_flt_t, input 120 * :param offdiag: pointer to array holding offdiagonal of :math:`T`, length :c:data:`n-1` 121 * :type offdiag: trlib_flt_t, input 122 * :param warm: set :math:`\ge 1` if you provide a valid value in :c:data:`leftmost_minor`, else ``0``. Exact value determines which model will be used for zero-finding 123 * 124 * ===== ====== 125 * warm model 126 * ===== ====== 127 * 0 linear model for rational function 128 * 1 sensible heuristic model choice for lifted rational function 129 * 2 asymptotic quadratic model :math:`\theta^2 + b \theta + c` for lifted rational function 130 * 3 taylor quadratic model :math:`a \theta^2 + b \theta + c` for lifted rational function 131 * 4 linear model :math:`b \theta + c` for lifted rational function 132 * ===== ====== 133 * 134 * :type warm: trlib_int_t, input 135 * :param leftmost_minor: smallest eigenvalue of principal :math:`(n-1)\times (n-1)` submatrix of :math:`T` 136 * :type leftmost_minor: trlib_flt_t, input 137 * :param itmax: maximum number of iterations 138 * :type itmax: trlib_int_t, input 139 * :param tol_abs: absolute stopping tolerance in Reid-Parlett zero finding, good default may be :math:`\sqrt{\texttt{macheps}}^{3/4}` (:c:macro:`TRLIB_EPS_POW_75`) 140 * :type tol_abs: trlib_flt_t, input 141 * :param verbose: determines the verbosity level of output that is written to :c:data:`fout` 142 * :type verbose: trlib_int_t, input 143 * :param unicode: set to ``1`` if :c:data:`fout` can handle unicode, otherwise to ``0`` 144 * :type unicode: trlib_int_t, input 145 * :param prefix: string that is printed before iteration output 146 * :type prefix: char, input 147 * :param fout: output stream 148 * :type fout: FILE, input 149 * :param timing: gives timing details, provide allocated zero initialized memory of length :c:func:`trlib_leftmost_timing_size` 150 * 151 * ====== ================================ 152 * block description 153 * ====== ================================ 154 * 0 total duration 155 * ====== ================================ 156 * 157 * :type timing: trlib_int_t, input/output 158 * :param leftmost: smallest eigenvalue of :math:`T` 159 * :type leftmost: trlib_flt_t, output 160 * :param iter_pr: number of Parlett-Reid iterations 161 * :type iter_pr: trlib_int_t, output 162 * 163 * :returns: status 164 * 165 * - :c:macro:`TRLIB_LMR_CONV` success 166 * - :c:macro:`TRLIB_LMR_ITMAX` iteration limit exceeded 167 * 168 * :rtype: trlib_int_t 169 */ 170 171 trlib_int_t trlib_leftmost_irreducible( 172 trlib_int_t n, trlib_flt_t *diag, trlib_flt_t *offdiag, 173 trlib_int_t warm, trlib_flt_t leftmost_minor, trlib_int_t itmax, trlib_flt_t tol_abs, 174 trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout, 175 trlib_int_t *timing, trlib_flt_t *leftmost, trlib_int_t *iter_pr); 176 177 /** size that has to be allocated for :c:data:`timing` in :c:func:`trlib_leftmost_irreducible` and :c:func:`trlib_leftmost` 178 */ 179 trlib_int_t trlib_leftmost_timing_size(void); 180 181 #endif 182