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_eigen_inverse(trlib_int_t n,trlib_flt_t * diag,trlib_flt_t * offdiag,trlib_flt_t lam_init,trlib_int_t itmax,trlib_flt_t tol_abs,trlib_flt_t * ones,trlib_flt_t * diag_fac,trlib_flt_t * offdiag_fac,trlib_flt_t * eig,trlib_int_t verbose,trlib_int_t unicode,char * prefix,FILE * fout,trlib_int_t * timing,trlib_flt_t * lam_pert,trlib_flt_t * pert,trlib_int_t * iter_inv)30 trlib_int_t trlib_eigen_inverse(
31 trlib_int_t n, trlib_flt_t *diag, trlib_flt_t *offdiag,
32 trlib_flt_t lam_init, trlib_int_t itmax, trlib_flt_t tol_abs,
33 trlib_flt_t *ones, trlib_flt_t *diag_fac, trlib_flt_t *offdiag_fac,
34 trlib_flt_t *eig, trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout,
35 trlib_int_t *timing, trlib_flt_t *lam_pert, trlib_flt_t *pert, trlib_int_t *iter_inv) {
36 // Local variables
37 #if TRLIB_MEASURE_TIME
38 struct timespec verystart, start, end;
39 TRLIB_TIC(verystart)
40 #endif
41 trlib_int_t info_fac = 0; // status variable for factorization
42 trlib_flt_t invnorm = 0.0; // 1/norm of eig before normalization
43 trlib_flt_t minuslam = - lam_init; // negative of current estimation of eigenvalue
44 trlib_int_t inc = 1; trlib_int_t nm = n-1;
45
46 trlib_int_t seeds[TRLIB_EIR_N_STARTVEC];
47 trlib_flt_t residuals[TRLIB_EIR_N_STARTVEC];
48
49 trlib_int_t jj = 0;
50 trlib_int_t kk = 0;
51 trlib_int_t seedpivot = 0;
52
53 *iter_inv = 0; // iteration counter
54 *pert = 0.0; // perturbation factor to update lam until factorization is possible
55
56 *iter_inv = 0;
57 *pert = 0.0;
58 info_fac = 0;
59 invnorm = 0.0;
60 minuslam = - lam_init;
61
62 // obtain factorization of T - lam*I, perturb until possible
63 // iter_inv is misused in this loop as flag if we can find a suitable lambda to start with
64 *iter_inv = TRLIB_EIR_FAIL_FACTOR;
65 while (*pert <= 1.0/TRLIB_EPS) {
66 // set diag_fac to diag - lam
67 TRLIB_DCOPY(&n, diag, &inc, diag_fac, &inc) // diag_fac <-- diag
68 TRLIB_DAXPY(&n, &minuslam, ones, &inc, diag_fac, &inc) // diag_fac <-- diag_fac - lam
69 TRLIB_DCOPY(&nm, offdiag, &inc, offdiag_fac, &inc) // offdiag_fac <-- offdiag
70 TRLIB_DPTTRF(&n, diag_fac, offdiag_fac, &info_fac); // compute factorization
71 if (info_fac == 0) { *iter_inv = 0; break; }
72 if (*pert == 0.0) {
73 *pert = TRLIB_EPS_POW_4 * fmax(1.0, -lam_init);
74 }
75 else {
76 *pert = 10.0*(*pert);
77 }
78 minuslam = *pert - lam_init;
79 }
80 *lam_pert = -minuslam;
81
82 if ( *iter_inv == TRLIB_EIR_FAIL_FACTOR ) { TRLIB_PRINTLN_2("Failure on factorizing in inverse correction!") TRLIB_RETURN(TRLIB_EIR_FAIL_FACTOR) }
83
84 // try with TRLIB_EIR_N_STARTVEC different start vectors and hope that it converges for one
85 seeds[0] = time(NULL);
86 for(jj = 1; jj < TRLIB_EIR_N_STARTVEC; ++jj ) { seeds[jj] = rand(); }
87 for(jj = 0; jj < TRLIB_EIR_N_STARTVEC; ++jj ) {
88 *iter_inv = 0;
89 srand((unsigned) seeds[jj]);
90 for(kk = 0; kk < n; ++kk ) { eig[kk] = ((trlib_flt_t)rand()/(trlib_flt_t)RAND_MAX); }
91
92 TRLIB_DNRM2(invnorm, &n, eig, &inc) invnorm = 1.0/invnorm;
93 TRLIB_DSCAL(&n, &invnorm, eig, &inc) // normalize eig
94 // perform inverse iteration
95 while (1) {
96 *iter_inv += 1;
97
98 if ( *iter_inv > itmax ) { break; }
99
100 // solve (T - lam*I)*eig_new = eig_old
101 TRLIB_DPTTRS(&n, &inc, diag_fac, offdiag_fac, eig, &n, &info_fac)
102 if( info_fac != 0 ) { TRLIB_PRINTLN_2("Failure on solving inverse correction!") TRLIB_RETURN(TRLIB_EIR_FAIL_LINSOLVE) }
103
104 // normalize eig
105 TRLIB_DNRM2(invnorm, &n, eig, &inc) invnorm = 1.0/invnorm;
106 TRLIB_DSCAL(&n, &invnorm, eig, &inc)
107
108 residuals[jj] = fabs(invnorm - *pert);
109
110 // check for convergence
111 if (residuals[jj] <= tol_abs ) { TRLIB_RETURN(TRLIB_EIR_CONV) }
112 }
113 }
114
115 // no convergence with any of the starting values.
116 // take the seed with least residual and redo computation
117 for(jj = 0; jj < TRLIB_EIR_N_STARTVEC; ++jj) { if (residuals[jj] < residuals[seedpivot]) { seedpivot = jj; } }
118
119 *iter_inv = 0;
120 srand((unsigned) seeds[seedpivot]);
121 for(kk = 0; kk < n; ++kk ) { eig[kk] = ((trlib_flt_t)rand()/(trlib_flt_t)RAND_MAX); }
122
123 TRLIB_DNRM2(invnorm, &n, eig, &inc) invnorm = 1.0/invnorm;
124 TRLIB_DSCAL(&n, &invnorm, eig, &inc) // normalize eig
125 // perform inverse iteration
126 while (1) {
127 *iter_inv += 1;
128
129 if ( *iter_inv > itmax ) { break; }
130
131 // solve (T - lam*I)*eig_new = eig_old
132 TRLIB_DPTTRS(&n, &inc, diag_fac, offdiag_fac, eig, &n, &info_fac)
133 if( info_fac != 0 ) { TRLIB_PRINTLN_2("Failure on solving inverse correction!") TRLIB_RETURN(TRLIB_EIR_FAIL_LINSOLVE) }
134
135 // normalize eig
136 TRLIB_DNRM2(invnorm, &n, eig, &inc) invnorm = 1.0/invnorm;
137 TRLIB_DSCAL(&n, &invnorm, eig, &inc)
138
139 residuals[seedpivot] = fabs(invnorm - *pert);
140
141 // check for convergence
142 if (residuals[seedpivot] <= tol_abs ) { TRLIB_RETURN(TRLIB_EIR_CONV) }
143 }
144
145 TRLIB_RETURN(TRLIB_EIR_ITMAX)
146 }
147
trlib_eigen_timing_size()148 trlib_int_t trlib_eigen_timing_size() {
149 #if TRLIB_MEASURE_TIME
150 return 1 + TRLIB_SIZE_TIMING_LINALG;
151 #endif
152 return 0;
153 }
154
155