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 #include "trlib.h"
26 #include "trlib_private.h"
27
28 #include "_c99compat.h"
29
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)30 trlib_int_t trlib_leftmost(
31 trlib_int_t nirblk, trlib_int_t *irblk, trlib_flt_t *diag, trlib_flt_t *offdiag,
32 trlib_int_t warm, trlib_flt_t leftmost_minor, trlib_int_t itmax, trlib_flt_t tol_abs,
33 trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout,
34 trlib_int_t *timing, trlib_int_t *ileftmost, trlib_flt_t *leftmost) {
35 trlib_int_t ret = 0, curit = 0;
36 if(! warm) {
37 trlib_int_t curret = 0;
38 trlib_int_t ii = 0;
39 ret = 0;
40 for(ii = 0; ii < nirblk; ++ii) {
41 curret = trlib_leftmost_irreducible(irblk[ii+1]-irblk[ii], diag+irblk[ii], offdiag+irblk[ii], 0, 0.0, itmax,
42 tol_abs, verbose, unicode, prefix, fout, timing, leftmost+ii, &curit);
43 if (curret == 0) { ret = curret; }
44 }
45 *ileftmost = 0;
46 for(ii = 1; ii < nirblk; ++ii) {
47 if (leftmost[ii] < leftmost[*ileftmost]) { *ileftmost = ii; }
48 }
49 }
50 else {
51 ret = trlib_leftmost_irreducible(irblk[nirblk] - irblk[nirblk-1], diag+irblk[nirblk-1], offdiag+irblk[nirblk-1],
52 1, leftmost_minor, itmax, tol_abs, verbose, unicode, prefix, fout, timing, leftmost+nirblk-1, &curit);
53 if (leftmost[nirblk-1] < leftmost[*ileftmost]) { *ileftmost = nirblk-1; }
54 }
55 return ret;
56 }
57
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)58 trlib_int_t trlib_leftmost_irreducible(
59 trlib_int_t n, trlib_flt_t *diag, trlib_flt_t *offdiag,
60 trlib_int_t warm, trlib_flt_t leftmost_minor, trlib_int_t itmax, trlib_flt_t tol_abs,
61 trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout,
62 trlib_int_t *timing, trlib_flt_t *leftmost, trlib_int_t *iter_pr) {
63 // Local variables
64 #if TRLIB_MEASURE_TIME
65 struct timespec verystart, start, end;
66 TRLIB_TIC(verystart)
67 #endif
68 trlib_int_t jj = 0; // local counter variable
69 trlib_flt_t low = 0.0; // lower bracket variable: low <= leftmost for desired value
70 trlib_flt_t up = 0.0; // upper bracket variable: leftmost <= up for desired value
71 trlib_flt_t leftmost_attempt = 0.0; // trial step for leftmost eigenvalue
72 trlib_flt_t dleftmost = 0.0; // increment
73 trlib_flt_t prlp = 0.0; // value of Parlett-Reid-Last-Pivot function
74 trlib_flt_t obyprlp = 0.0; // quotient used in Cholesky computation
75 trlib_flt_t dprlp = 0.0; // derivative of Parlett-Reid-Last-Pivot function wrt to leftmost
76 trlib_flt_t ddprlp = 0.0; // second derivative of Parlett-Reid-Last-Pivot function wrt to leftmost
77 trlib_int_t n_neg_piv = 0; // number of negative pivots in factorization
78 trlib_flt_t quad_abs = 0.0; // absolute coefficient in quadratic model
79 trlib_flt_t quad_lin = 0.0; // linear coefficient in quadratic model
80 trlib_flt_t quad_qua = 0.0; // quadratic coefficient in quadratic model
81 trlib_flt_t zerodum = 0.0; // dummy return variables from quadratic equation
82 trlib_flt_t oabs0 = 0.0, oabs1 = 0.0; // temporaries in Gershgorin limit computation
83
84 trlib_int_t continue_outer_loop = 0; // local spaghetti code control variable
85 trlib_int_t model_type = 0;
86 trlib_int_t ii = 0;
87 *leftmost = 0.0; // estimation of desired leftmost eigenvalue
88 *iter_pr = 0; // iteration counter
89
90 // trivial case: one-dimensional. return diagonal value
91 if (n == 1) { *leftmost = diag[0]; TRLIB_RETURN(TRLIB_LMR_CONV) }
92
93 /* set bracket interval derived from Gershgorin circles
94 Gershgorin:
95 eigenvalues are contained in the union of balls centered at
96 diag_i with radius sum of absolute values in column i, except diagonal element
97 this estimation is rough and could be improved by circle component analysis
98 determine if worth doing */
99
100 oabs0 = fabs(offdiag[0]); oabs1 = fabs(offdiag[n-2]);
101 low = fmin( diag[0] - oabs0, diag[n-1] - oabs1 );
102 up = fmax( diag[0] + oabs0, diag[n-1] - oabs1 );
103 for(ii = 1; ii < n-1; ++ii ) {
104 oabs1 = fabs(offdiag[ii]);
105 low = fmin( low, diag[ii] - oabs0 - oabs1 );
106 up = fmax( up, diag[ii] + oabs0 + oabs1 );
107 oabs0 = oabs1;
108 }
109
110 /* set leftmost to sensible initialization
111 on warmstart, provided leftmost is eigenvalue of principal (n-1) * (n-1) submatrix
112 by eigenvalue interlacing theorem desired value <= provided leftmost
113 on coldstart, start close lower bound as hopefully this is a good estimation */
114 if ( warm ) {
115 // provided leftmost is an upper bound and a pole of Parlett-Reid Value, thus pertub a bit
116 up = fmin(up, leftmost_minor); *leftmost = leftmost_minor - .1*(up-low); //*leftmost = leftmost_minor - TRLIB_EPS_POW_4;
117 }
118 else { leftmost_minor = 0.0; *leftmost = low + .1*(up-low); }; // ensure sanity on leftmost_minor and start with lower bound
119 // Parlett-Reid Iteration, note we can assume n > 1
120 itmax = itmax*n;
121
122 while (1) {
123 /* iterate to obtain Parlett-Reid last pivot value of -leftmost == 0.0
124 this iteration uses a safeguard bracket [low, up] such that always low <= leftmost <= up
125 note that T - t*I is positive definite for t <= desired leftmost
126 steps of iteration:
127
128 (1) compute Parlett-Reid last pivot value which is D_n in a LDL^T factorization of T
129 obtain derivative d D_n / d leftmost as byproduct in the very same recursion
130 track if breakdown would occur in factorization, happens if either
131 (a) a pivot become zero
132 (b) more than one negative pivot present
133 if breakdown would occurs this means that Parlett-Reid value is infinite
134 end iteration at premature point and restart with adapted bounds and estimation:
135 (a) a pivot became zero:
136 if last pivot zero --> goal reached, exit
137 if previous zero --> T - leftmost I not positive definite, thus desired value <= leftmost
138 (b) multiple neg privots --> T - leftmost I indefinite, thus desired value <= leftmost
139 (2) compute a trial update for leftmost. Possibilities
140 (a) Newton update
141 (b) zero of suitable model of analytic expression,
142 analytic expression is given by prlp(t) = det(T-t*I)/det(U-t*I) with U principal (n-1)*(n-1) submatrix
143 prlp(t) has a pole at leftmost_minor, so better use lifted
144 lprlp(t) = (leftmost_minor-t)*prlp(t)
145 Gould proposes model m(t) = a+bt+t^2 (correct asymptotic, very good for t far away)
146 Other choices would be m(t) = a+bt (Newton on lifted model)
147 m(t) = a+bt+ct^2 (Taylor, very good close to zero, possibly really bad far away)
148 do (b) if warmstart where user provided leftmost(U), otherwise go route (a)
149
150 (3) take trial step if inside bracket, otherwise bisect
151
152 stop iteration if either bracket is sufficiently small or Parlett-Reid value is close enough to zero */
153
154 *iter_pr += 1;
155
156 // test if iteration limit exceeded
157 if ( *iter_pr > itmax ) { TRLIB_RETURN(TRLIB_LMR_ITMAX) }
158
159 // initialize: no negative pivots so far
160 n_neg_piv = 0;
161
162 // print iteration headline every 10 iterations
163 if (*iter_pr % 10 == 1) {
164 TRLIB_PRINTLN_1("%6s%8s%14s%14s%14s%14s%14s%6s%6s", " it ", " action ", " low ", " leftmost ", " up ", " dleftmost ", " prlp ", " nneg ", " br ")
165 }
166 TRLIB_PRINTLN_1("%6ld%8s%14e%14e%14e", *iter_pr, " entry ", low, *leftmost, up)
167
168 // compute pivot and derivative of LDL^T factorization of T - leftmost I
169 continue_outer_loop = 0;
170 for( jj = 0; jj < n; ++jj ) {
171 /* compute jjth pivot
172 special case for jj == 0 since offdiagonal is missing */
173 if (jj == 0) { prlp = diag[0] - *leftmost; dprlp = -1.0; ddprlp = 0.0; }
174 else{
175 // update pivot as pivot = d_j - leftmost - o_{j-1}^2/pivot
176 // thus dpivot/dleftmost = - 1.0 + o_{j-1}^2/pivot^2 * dpivot
177 // d^2 pivot/dleftmost^2 = (o_{j-1}^2/pivot^2)(ddpivot - 2 dpivot^2/pivot)
178 obyprlp = offdiag[jj-1]/prlp;
179 //dprlp = -1.0 + obyprlp*dprlp*obyprlp;// * dprlp;
180 //prlp = diag[jj] - offdiag[jj-1]*offdiag[jj-1]/prlp - *leftmost;
181 ddprlp = obyprlp*obyprlp*(ddprlp - 2.0*dprlp*dprlp/prlp);
182 dprlp = -1.0 + offdiag[jj-1]*offdiag[jj-1]*dprlp / (prlp*prlp);
183 prlp = diag[jj] - *leftmost - offdiag[jj-1]*obyprlp;
184 }
185
186 // check for breakdown
187 if (prlp == 0.0) {
188 // if last pivot and no negative pivots encountered --> finished
189 if (n_neg_piv == 0 && jj+1 == n) { TRLIB_RETURN(TRLIB_LMR_CONV) }
190 else{
191 /* if not last pivot or negative pivots encountered:
192 estimation provides a new upper bound; reset estimation */
193 up = *leftmost;
194 *leftmost = 0.5 * (low+up);
195 continue_outer_loop = 1;
196 break; // continue outer loop
197 }
198 }
199 else if ( prlp < 0.0 ) {
200 n_neg_piv += 1;
201 up = *leftmost;
202 if (n_neg_piv > 1) {
203 // more than one negative pivot: factorization would fail, to the right of singularity!
204 *leftmost = 0.5 * (low+up);
205 continue_outer_loop = 1;
206 break; // continue outer loop
207 }
208 }
209 }
210
211 if (continue_outer_loop) {
212 TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14s%14e%6ld%6ld", "", " bisecp ", low, *leftmost, up, "", prlp, n_neg_piv, jj)
213 continue;
214 }
215
216 // we have survived computing the Last-Pivot value without finding a zero pivot and at most one negative pivot
217
218 // adapt bracket, no negative pivots encountered: leftmost provides new lower bound, otherwise upper bound
219 if (n_neg_piv == 0) { low = *leftmost; }
220 else { up = *leftmost; }
221
222 // test if bracket interval is small or last pivot has converged to zero
223 if (up-low <= tol_abs * fmax(1.0, fmax(fabs(low), fabs(up))) || fabs(prlp) <= tol_abs) {
224 TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14s%14e%6ld%6ld", "", " conv ", low, *leftmost, up, "", prlp, n_neg_piv, jj)
225 TRLIB_RETURN(TRLIB_LMR_CONV)
226 }
227
228 /* compute trial step for new leftmost
229 on coldstart do Newton iteration, on warmstart find zero of model of analytic expression */
230 // select suitable model for analytic expression in dependence of warm
231 // warm: 1 --> heuristic depending on estimation if already close to zero
232 // warm: 2 --> use asymptotic quadratic model of lifted prlp
233 // warm: 3 --> use Taylor quadratic model of lifted prlp
234 // warm: 4 --> use linear (newton) model of lifted prlp
235 if (warm) {
236 if ( warm == 2 || (warm == 1 && up-low >= .1*fmax(1.0, fabs(*leftmost))) ) {
237 /* use analytic model m(t) = (t-a)(t-b)/(t-leftmost_minor) for prlp(t)
238 fit a, b such that m matches function value and derivative
239 at current estimation and compute left zero of numerator */
240 quad_lin = -(2.0*(*leftmost)+prlp+((*leftmost)-leftmost_minor)*dprlp);
241 quad_abs = -(((*leftmost)-leftmost_minor)*prlp+(*leftmost)*(quad_lin+(*leftmost)));
242 trlib_quadratic_zero(quad_abs, quad_lin, TRLIB_EPS_POW_75, 0, 0, "", NULL, &leftmost_attempt, &zerodum);
243 model_type = 2; dleftmost = leftmost_attempt - *leftmost;
244 }
245 if( warm > 2 || (warm == 1 && up-low <.1*fmax(1.0, fabs(*leftmost))) ) {
246 /* use quadratic Taylor model for pole lifted function (leftmost_minor-t)*prlp(t) */
247 quad_qua = -(dprlp + .5*((*leftmost)-leftmost_minor)*ddprlp);
248 quad_lin = -(prlp + ((*leftmost)-leftmost_minor)*dprlp);
249 quad_abs = (leftmost_minor-(*leftmost))*prlp;
250 if ( warm == 4 || fabs(quad_qua) <= TRLIB_EPS * TRLIB_EPS || quad_lin*quad_lin - 4.0*quad_qua*quad_abs <= 0.0 ) { // tiny curvature, resort to Newton step of lifted function
251 model_type = 3; dleftmost = -quad_abs/quad_lin; leftmost_attempt = *leftmost + dleftmost;
252 }
253 else { // scale coefficients to normalize equation
254 quad_lin = quad_lin/quad_qua;
255 quad_abs = quad_abs/quad_qua;
256 trlib_quadratic_zero(quad_abs, quad_lin, TRLIB_EPS_POW_75, 0, 0, "", NULL, &dleftmost, &zerodum);
257 model_type = 4; leftmost_attempt = *leftmost + dleftmost;
258 }
259 }
260 }
261 else { model_type = 1; dleftmost = -prlp/dprlp; leftmost_attempt = *leftmost + dleftmost; } // Newton step
262
263 // assess if we can use trial step
264 if (low <= leftmost_attempt && leftmost_attempt <= up) {
265 if( fabs(dleftmost) <= tol_abs * fmax(1.0, fmax(fabs(low), fabs(up))) ) { TRLIB_RETURN(TRLIB_LMR_NEWTON_BREAK) }
266 *leftmost = leftmost_attempt;
267 }
268 else {
269 // if warmstart information available, lifted newton step may still be feasible
270 if(warm) {
271 quad_lin = -(prlp + ((*leftmost)-leftmost_minor)*dprlp);
272 quad_abs = (leftmost_minor-(*leftmost))*prlp;
273 model_type = 3; dleftmost = -quad_abs/quad_lin; leftmost_attempt = *leftmost + dleftmost;
274 }
275 // if that fails, newton step may still be feasible
276 if(low > leftmost_attempt || up < leftmost_attempt) {
277 model_type = 1; dleftmost = -prlp/dprlp; leftmost_attempt = *leftmost + dleftmost;
278 }
279 // now out of options, do bisection
280 if(low > leftmost_attempt || up < leftmost_attempt) {
281 model_type = 0; *leftmost = .5*(low+up);
282 }
283 }
284 if ( verbose > 0 ) {
285 if ( model_type == 0 ) { TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14e%14e%6ld%6ld", "", " bisecs ", low, *leftmost, up, .5*(up-low), prlp, n_neg_piv, jj) }
286 if ( model_type == 1 ) { TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14e%14e%6ld%6ld", "", " piv 1 ", low, *leftmost, up, dleftmost, prlp, n_neg_piv, jj) }
287 if ( model_type == 2 ) { TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14e%14e%6ld%6ld", "", " lpiv q ", low, *leftmost, up, dleftmost, prlp, n_neg_piv, jj) }
288 if ( model_type == 3 ) { TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14e%14e%6ld%6ld", "", " lpiv 1 ", low, *leftmost, up, dleftmost, prlp, n_neg_piv, jj) }
289 if ( model_type == 4 ) { TRLIB_PRINTLN_1("%6s%8s%14e%14e%14e%14e%14e%6ld%6ld", "", " lpiv 2 ", low, *leftmost, up, dleftmost, prlp, n_neg_piv, jj) }
290 }
291
292 }
293 }
294
trlib_leftmost_timing_size()295 trlib_int_t trlib_leftmost_timing_size() {
296 #if TRLIB_MEASURE_TIME
297 return 1;
298 #endif
299 return 0;
300 }
301