1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
/* MIT License
*
* Copyright (c) 2016--2017 Felix Lenders
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*
*/
#ifndef TRLIB_LEFTMOST_H
#define TRLIB_LEFTMOST_H
#define TRLIB_LMR_CONV (0)
#define TRLIB_LMR_ITMAX (-1)
#define TRLIB_LMR_NEWTON_BREAK (3)
/** Computes smallest eigenvalue of symmetric tridiagonal matrix
* :math:`T \in \mathbb R^{n\times n}`,
* using a iteration based on last-pivot function of Parlett and Reid.
*
* Let :math:`T = \begin{pmatrix} T_1 & & \\ & \ddots & \\ & & T_\ell \end{pmatrix}`
* be composed into irreducible blocks :math:`T_i`.
*
* Calls :c:func:`trlib_leftmost_irreducible` on every irreducible block in case of coldstart,
* in case of warmstart just updates information on :math:`T_\ell`.
*
* :param nirblk: number of irreducible blocks :math:`\ell`, ensure :math:`\ell > 0`
* :type nirblk: trlib_int_t, input
* :param irblk: pointer to indices of irreducible blocks, length :c:data:`nirblk+1`:
*
* - :c:data:`irblk[i]` is the start index of block :math:`i` in :c:data:`diag` and :c:data:`offdiag`
* - :c:data:`irblk[i+1] - 1` is the stop index of block :math:`i`
* - :c:data:`irblk[i+1] - irred[i]` the dimension :math:`n_\ell` of block :math:`i`
* - :c:data:`irblk[nirred]` the dimension of :math:`T`
*
* :type irblk: trlib_int_t, input
* :param diag: pointer to array holding diagonal of :math:`T`, length :c:data:`irblk[nirblk]`
* :type diag: trlib_flt_t, input
* :param offdiag: pointer to array holding offdiagonal of :math:`T`, length :c:data:`irblk[nirblk]`
* :type offdiag: trlib_flt_t, input
* :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``
* :type warm: trlib_int_t, input
* :param leftmost_minor: smallest eigenvalue of principal :math:`(n_\ell-1)\times (n_\ell-1)` submatrix of :math:`T_\ell`
* :type leftmost_minor: trlib_flt_t, input
* :param itmax: maximum number of iterations
* :type itmax: trlib_int_t, input
* :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`)
* :type tol_abs: trlib_flt_t, input
* :param verbose: determines the verbosity level of output that is written to :c:data:`fout`
* :type verbose: trlib_int_t, input
* :param unicode: set to ``1`` if :c:data:`fout` can handle unicode, otherwise to ``0``
* :type unicode: trlib_int_t, input
* :param prefix: string that is printed before iteration output
* :type prefix: char, input
* :param fout: output stream
* :type fout: FILE, input
* :param timing: gives timing details, provide allocated zero initialized memory of length :c:func:`trlib_leftmost_timing_size`
*
* ====== ================================
* block description
* ====== ================================
* 0 total duration
* ====== ================================
*
* :type timing: trlib_int_t, input/output
* :param ileftmost: index of block that corresponds to absolute smallest eigenvalue
* :type ileftmost: trlib_int_t, input/output
* :param leftmost: smallest eigenvalue of :math:`T`, length :math:`\ell`
*
* - on entry: allocated memory
* - on exit: :c:data:`leftmost[i]` smallest eigenvalue of :math:`T_i`
*
* :type leftmost: trlib_flt_t, input/output
* :returns: status
*
* - :c:macro:`TRLIB_LMR_CONV` success
* - :c:macro:`TRLIB_LMR_ITMAX` iteration limit exceeded
*
* :rtype: trlib_int_t
*/
trlib_int_t trlib_leftmost(
trlib_int_t nirblk, trlib_int_t *irblk, trlib_flt_t *diag, trlib_flt_t *offdiag,
trlib_int_t warm, trlib_flt_t leftmost_minor, trlib_int_t itmax, trlib_flt_t tol_abs,
trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout,
trlib_int_t *timing, trlib_int_t *ileftmost, trlib_flt_t *leftmost);
/** Computes smallest eigenvalue of irreducible symmetric tridiagonal matrix
* :math:`T \in \mathbb R^{n\times n}`,
* using a iteration based on last-pivot function of Parlett and Reid.
*
* Method is sketched on p. 516 in [Gould1999]_.
*
* Note that this function most likely will fail in the case of a reducible matrix
* (:c:data:`offdiag` contains 0).
*
* **Convergence**
*
* 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.
*
* :param n: dimension, ensure :math:`n > 0`
* :type n: trlib_int_t, input
* :param diag: pointer to array holding diagonal of :math:`T`, length :c:data:`n`
* :type diag: trlib_flt_t, input
* :param offdiag: pointer to array holding offdiagonal of :math:`T`, length :c:data:`n-1`
* :type offdiag: trlib_flt_t, input
* :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
*
* ===== ======
* warm model
* ===== ======
* 0 linear model for rational function
* 1 sensible heuristic model choice for lifted rational function
* 2 asymptotic quadratic model :math:`\theta^2 + b \theta + c` for lifted rational function
* 3 taylor quadratic model :math:`a \theta^2 + b \theta + c` for lifted rational function
* 4 linear model :math:`b \theta + c` for lifted rational function
* ===== ======
*
* :type warm: trlib_int_t, input
* :param leftmost_minor: smallest eigenvalue of principal :math:`(n-1)\times (n-1)` submatrix of :math:`T`
* :type leftmost_minor: trlib_flt_t, input
* :param itmax: maximum number of iterations
* :type itmax: trlib_int_t, input
* :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`)
* :type tol_abs: trlib_flt_t, input
* :param verbose: determines the verbosity level of output that is written to :c:data:`fout`
* :type verbose: trlib_int_t, input
* :param unicode: set to ``1`` if :c:data:`fout` can handle unicode, otherwise to ``0``
* :type unicode: trlib_int_t, input
* :param prefix: string that is printed before iteration output
* :type prefix: char, input
* :param fout: output stream
* :type fout: FILE, input
* :param timing: gives timing details, provide allocated zero initialized memory of length :c:func:`trlib_leftmost_timing_size`
*
* ====== ================================
* block description
* ====== ================================
* 0 total duration
* ====== ================================
*
* :type timing: trlib_int_t, input/output
* :param leftmost: smallest eigenvalue of :math:`T`
* :type leftmost: trlib_flt_t, output
* :param iter_pr: number of Parlett-Reid iterations
* :type iter_pr: trlib_int_t, output
*
* :returns: status
*
* - :c:macro:`TRLIB_LMR_CONV` success
* - :c:macro:`TRLIB_LMR_ITMAX` iteration limit exceeded
*
* :rtype: trlib_int_t
*/
trlib_int_t trlib_leftmost_irreducible(
trlib_int_t n, trlib_flt_t *diag, trlib_flt_t *offdiag,
trlib_int_t warm, trlib_flt_t leftmost_minor, trlib_int_t itmax, trlib_flt_t tol_abs,
trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout,
trlib_int_t *timing, trlib_flt_t *leftmost, trlib_int_t *iter_pr);
/** size that has to be allocated for :c:data:`timing` in :c:func:`trlib_leftmost_irreducible` and :c:func:`trlib_leftmost`
*/
trlib_int_t trlib_leftmost_timing_size(void);
#endif
|
