TBCI Numerical high perf. C++ Library 2.8.0
svd_solver.h
Go to the documentation of this file.
1
9
10
11#ifndef TBCI_SOLVER_SVD_SOLVER_H
12#define TBCI_SOLVER_SVD_SOLVER_H
13
14#include "tbci/matrix.h"
15#ifdef TYP
16# define SVDFLOAT TYP
17#else
18# define SVDFLOAT typename MatrixType::value_type
19#endif
20
21NAMESPACE_TBCI
22
23#ifdef NEED_BCPLX_WORKAROUND
24# define SVDREAL(x) CPLX__ real(x)
25#else
26# define SVDREAL(x) (x)
27#endif
28
30template <typename ARG>
31inline double sv_decomp_pythag(ARG a, ARG b)
32{
33 const double absa = MATH__ fabs(a);
34 const double absb = MATH__ fabs(b);
35 if (absa > absb)
36 return absa * MATH__ sqrt(1.0 + (absb/absa)*(absb/absa));
37 else
38 return (absb==0.0? 0.0: absb * MATH__ sqrt(1.0 + (absa/absb)*(absa/absb)));
39}
40
41
42//TODO: Check which of the variables used should be double rather
43// then SVDFLOAT (which can be complex, e.g.)
44
45INST2(template <typename MAT, typename VEC > class MAT friend void sv_decomp (MAT &, MAT &, VEC &);)
46template <typename MatrixType /*= TBCI__ Matrix<double>*/,
47 typename VectorType /*= TBCI__ Vector<double>*/ >
53void sv_decomp (MatrixType & A, MatrixType & V, VectorType & W)
54{
55 const int m = A.rows();
56 const int n = A.columns();
57 int l = 0, i, k; // M$ annoyance (for scoping bug)
58 VectorType rv1(n);
59 SVDFLOAT g = 0.0, scale = 0.0;
60 double anorm = 0.0;
61 for (i = 0; i < n; ++i) {
62 SVDFLOAT s = 0.0;
63 l = i+1;
64 rv1[i] = scale*g;
65 g = scale = 0.0;
66 if (i < m) {
67 for (k = i; k < m; ++k)
68 scale += MATH__ fabs(A(k,i));
69 if (scale != (SVDFLOAT)0) {
70 for (k = i; k < m; ++k) {
71 A(k,i) /= scale;
72 s += A(k,i)*A(k,i);
73 }
74 SVDFLOAT f = A(i,i);
75 g = SVDREAL(f)>SVDREAL((SVDFLOAT)0) ? - MATH__ sqrt(s) : MATH__ sqrt(s);
76 SVDFLOAT h = f*g - s;
77 A(i,i) = f - g;
78 for (int j = l; j < n; ++j) {
79 s = (SVDFLOAT)0;
80 for (k = i; k < m; ++k)
81 s += A(k,i)*A(k,j);
82 SVDFLOAT f = s/h;
83 for (k = i; k < m; ++k)
84 A(k,j) += f*A(k,i);
85 }
86 for (k = i; k < m; ++k)
87 A(k,i) *= scale;
88 }
89 }
90 W[i] = scale * g;
91 g = s = scale = 0.0;
92 if (i < m && i != (n-1)) {
93 for (k = l; k < n; ++k)
94 scale += MATH__ fabs(A(i,k));
95 if (scale != (SVDFLOAT)0) {
96 for (k = l; k < n; ++k) {
97 A(i,k) /= scale;
98 s += A(i,k)*A(i,k);
99 }
100 SVDFLOAT f = A(i,l);
101 g = SVDREAL(f)>SVDREAL((SVDFLOAT)0) ? - MATH__ sqrt(s) : MATH__ sqrt(s);
102 SVDFLOAT h = f*g - s;
103 A(i,l) = f - g;
104 for (k=l; k < n; ++k)
105 rv1[k] = A(i,k)/h;
106 for (int j = l; j < m; ++j) {
107 s = 0.0;
108 for (k = l; k < n; ++k)
109 s += A(j,k)*A(i,k);
110 for (k = l; k < n; ++k)
111 A(j,k) += s*rv1[k];
112 }
113 for(k = l; k < n; ++k)
114 A(i,k) *= scale;
115 }
116 }
117 anorm = MAX(anorm, MATH__ fabs(W[i]) + MATH__ fabs(rv1[i]));
118 }
119 for (i = n-1; i >= 0; --i) {
120 if (i < (n-1)) {
121 if (g != (SVDFLOAT)0) {
122 int j;
123 for (j = l; j < n; ++j )
124 V(j,i) = (A(i,j)/A(i,l))/g;
125 for (j = l; j < n; ++j) {
126 SVDFLOAT s = 0.0;
127 for (k = l; k < n; ++k )
128 s += A(i,k)*V(k,j);
129 for (k = l; k < n; ++k )
130 V(k,j) += s*V(k,i);
131 }
132 }
133 for (int j = l; j < n; ++j)
134 V(i,j) = V(j,i) = 0.0;
135 }
136 V(i,i) = 1.0;
137 g = rv1[i];
138 l = i;
139 }
140 for (i = MIN(m,n)-1; i >= 0; --i) {
141 int j;
142 int l = i+1;
143 g = W[i];
144 for (j = l; j < n; ++j)
145 A(i,j) = 0.0;
146 if (g != (SVDFLOAT)0) {
147 g = 1.0/g;
148 for (j = l; j < n; ++j) {
149 SVDFLOAT s = 0.0;
150 for (k = l; k < m; ++k)
151 s += A(k,i)*A(k,j);
152 SVDFLOAT f = (s/A(i,i))*g;
153 for (k = i; k < m; ++k)
154 A(k,j) += f*A(k,i);
155 }
156 for (j = i; j < m; ++j)
157 A(j,i) *= g;
158 } else
159 for (j = i; j < m; ++j)
160 A(j,i) = 0.0;
161 A(i,i) += 1.0;
162 }
163 for (k = n-1; k >= 0; --k) {
164 for (int its = 0; its < 30; ++its) {
165 bool flag = true;
166 int nm=0;
167 for (l = k; l >= 0; --l) {
168 nm = l-1; // note that rv1[0] is always zero
169 if ((MATH__ fabs(rv1[l]) + anorm) == anorm) {
170 flag = false;
171 break;
172 }
173 if ((MATH__ fabs(W[nm]) + anorm) == anorm)
174 break;
175 }
176 if (flag) {
177 SVDFLOAT c = 0.0;
178 SVDFLOAT s = 1.0;
179 for (i = l; i < k; ++i) {
180 SVDFLOAT f = s*rv1[i];
181 rv1[i] *= c;
182 if ((MATH__ fabs(f) + anorm) == anorm)
183 break;
184 SVDFLOAT g = W[i];
185 SVDFLOAT h = sv_decomp_pythag(f, g);
186 W[i] = h;
187 h = 1.0/h;
188 c = g*h;
189 s = -f*h;
190 for (int j = 0; j < m; ++j) {
191 SVDFLOAT y = A(j,nm);
192 SVDFLOAT z = A(j,i);
193 A(j,nm) = y*c + z*s;
194 A(j,i) = z*c - y*s;
195 }
196 }
197 }
198 SVDFLOAT z = W[k];
199 if (l == k) { // convergence
200 if (SVDREAL(z) < SVDREAL((SVDFLOAT)0)) {
201 W[k] = -z;
202 for (int j = 0; j < n; ++j)
203 V(j,k) = -V(j,k);
204 }
205 break;
206 }
207 if (its == 29)
208 STD__ cerr << "SVD: No convergence in 30 iterations.\n";
209 SVDFLOAT x = W[l];
210 nm = k-1;
211 SVDFLOAT y = W[nm];
212 SVDFLOAT g = rv1[nm];
213 SVDFLOAT h = rv1[k];
214 SVDFLOAT f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
215 g = sv_decomp_pythag(f, (SVDFLOAT)1.0);
216 f = ((x-z)*(x+z) +
217 h*((y / (f + (SVDREAL(f)>SVDREAL((SVDFLOAT)0) ? MATH__ fabs(g) : - MATH__ fabs(g))))-h))/x;
218 SVDFLOAT c = 1.0;
219 SVDFLOAT s = 1.0;
220 for (int j = l; j <= nm; ++j) {
221 int jj; // You know who annoyance
222 i = j+1;
223 SVDFLOAT g = rv1[i];
224 SVDFLOAT y = W[i];
225 SVDFLOAT h = s*g;
226 g*=c;
227 SVDFLOAT z = sv_decomp_pythag(f, h);
228 rv1[j] = z;
229 c = f/z;
230 s = h/z;
231 f = x*c + g*s;
232 g = g*c - x*s;
233 h = y*s;
234 y *= c;
235 for (jj = 0; jj < n; ++jj) {
236 SVDFLOAT x = V(jj,j);
237 SVDFLOAT z = V(jj,i);
238 V(jj,j) = x*c + z*s;
239 V(jj,i) = z*c - x*s;
240 }
241 z = sv_decomp_pythag(f, h);
242 W[j] = z;
243 if (z != (SVDFLOAT)0) {
244 z = 1.0/z;
245 c = f*z;
246 s = h*z;
247 }
248 f = c*g + s*y;
249 x = c*y - s*g;
250 for (jj = 0; jj < m; ++jj) {
251 SVDFLOAT y = A(jj,j);
252 SVDFLOAT z = A(jj,i);
253 A(jj,j) = y*c + z*s;
254 A(jj,i) = z*c - y*s;
255 }
256 }
257 rv1[l] = 0.0;
258 rv1[k] = f;
259 W[k] = x;
260 }
261 }
262}
263
264
265INST2(template <typename MAT, typename VEC > class MAT friend void sv_decomp_backsub (const MAT &, const MAT &, const VEC &, const VEC &, VEC &);)
266template <typename MatrixType /*= TBCI__ Matrix<double>*/,
267 typename VectorType /*= TBCI__ Vector<double>*/ >
268void sv_decomp_backsub (const MatrixType & U, const MatrixType & V,
269 const VectorType & W, const VectorType & b,
270 VectorType & x)
271{
272 const int n = W.size();
273 const int m = b.size();
274 VectorType tmp(n);
275 int j;
276 for (j = 0; j < n; ++j) {
277 SVDFLOAT s = 0.0;
278 if (W[j] != (SVDFLOAT)0) {
279 for (int i = 0; i < m; ++i)
280 s += U(i,j)*b[i];
281 s /= W[j];
282 }
283 tmp[j] = s;
284 }
285 for (j = 0; j < n; ++j) {
286 SVDFLOAT s = 0.0;
287 for (int i = 0; i < n; ++i)
288 s += V(j,i)*tmp[i];
289 x[j] = s;
290 }
291}
292
293INST(template <typename VEC > class VEC friend /*static*/ int fix_condition (VEC& vec, const double cndno);)
294/* Fix those singular values that are close to zero */
295template <typename VectorType >
296static int fix_condition (VectorType& vec, const double cndno = 1e-12)
297{
298 int fixed = 0; unsigned i;
299 double maxsv = 0.0;
300 const unsigned size = vec.size();
301 //typename VectorType::value_type * RESTRICT const VEC = &vec.setval(0);
302 for (i = 0; i < size; ++i)
303 maxsv = MAX(maxsv, MATH__ fabs(vec[i]));
304 double minsv = maxsv * cndno;
305 for (i = 0; i < size; ++i)
306 if (MATH__ fabs(vec[i]) < minsv) {
307 fixed++;
308 vec[i] = (typename VectorType::value_type)0.0;
309 }
310 return fixed;
311}
312
313
314
315INST2(template <typename MAT, typename VEC > class MAT friend VEC svd_solve (const MAT& mat, const VEC& b, const double cndno);)
323template <typename MatrixType, typename VectorType>
324VectorType svd_solve (const MatrixType& mat, const VectorType& b, const double cndno = 1e-12)
325{
326 MatrixType A(mat);
327 MatrixType V(mat.columns(), mat.columns());
328 VectorType W(mat.columns());
329 sv_decomp <MatrixType, VectorType > (A, V, W);
330 fix_condition(W, cndno);
331 // we could calc mat^-1 = V * diag (1/W) * A * A^T now
332 VectorType x (mat.columns());
333 sv_decomp_backsub(A, V, W, b, x);
334 return x;
335}
336
337
338// Optimized version for our own matrix & vector classes, using TBCI
339INST2(template <typename MAT, typename T > class MAT friend TVector<T> sv_decomp_backsub (const MAT &, const MAT &, const Vector<T> &, const Vector<T> &);)
340template <typename MatrixType, typename T >
341TVector<T> sv_decomp_backsub (const MatrixType & U, const MatrixType & V,
342 const Vector<T> & W, const Vector<T> & b)
343{
344 const int n = W.size();
345 const int m = b.size();
346 TVector<T> x (n);
347 Vector<T> tmp(n);
348 int j;
349 for (j = 0; j < n; ++j) {
350 SVDFLOAT s = 0.0;
351 if (W[j] != (SVDFLOAT)0) {
352 for (int i = 0; i < m; ++i)
353 s += U(i,j)*b[i];
354 s /= W[j];
355 }
356 tmp[j] = s;
357 }
358 for (j = 0; j < n; ++j) {
359 SVDFLOAT s = 0.0;
360 for (int i = 0; i < n; ++i)
361 s += V(j,i)*tmp[i];
362 x.set (s, j);
363 }
364 return x;
365}
366
367INST2(template <typename MAT, typename T > class MAT friend TVector<T> svd_solve (const MAT& mat, const Vector<T>& b, const double cndno);)
378template <typename MatrixType, typename T >
379TVector<T> svd_solve (const MatrixType& mat, const Vector<T>& b, const double cndno = 1e-12)
380{
381 MatrixType A(mat);
382 MatrixType V(mat.columns(), mat.columns());
383 Vector<T> W(mat.columns());
384 sv_decomp <MatrixType, Vector<T> > (A, V, W);
385 fix_condition (W, cndno);
386 return sv_decomp_backsub (A, V, W, b);
387}
388
389
393
394#undef SVDFLOAT
395
396NAMESPACE_END
397#endif /* TBCI_SOLVER_SVD_SOLVER_H */
x
const Vector< T > const Vector< T > & x
Definition LM_fit.h:97
h
const Vector< T > const Vector< T > const Vector< T > int T h
Definition LM_fit.h:99
i
int i
Definition LM_fit.h:71
y
doublereal y
Definition TOMS_707.C:27
STD__
#define STD__
Definition basics.h:338
MIN
#define MIN(a, b)
Definition basics.h:655
NAMESPACE_END
#define NAMESPACE_END
Definition basics.h:323
INST
#define INST(x)
Definition basics.h:238
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition basics.h:317
MATH__
#define MATH__
Definition basics.h:339
INST2
#define INST2(x, y)
Definition basics.h:239
MAX
#define MAX(a, b)
Definition basics.h:656
U
#define U
Definition bdmatlib.cc:21
MAT
#define MAT
Definition bdmatsolverlib.cc:27
VEC
#define VEC
Definition bdmatsolverlib.cc:28
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition vector.h:73
TVector::size
unsigned long size() const
Definition vector.h:104
fabs
double fabs(const TBCI__ cplx< T > &c)
Definition cplx.h:746
sqrt
NAMESPACE_END NAMESPACE_CSTD TBCI__ cplx< T > sqrt(const TBCI__ cplx< T > &z)
Definition cplx.h:751
flag
long int flag
Definition f2c.h:66
c
return c
Definition f_matrix.h:760
b
F_TMatrix< T > b
Definition f_matrix.h:736
MatrixType
#define MatrixType
Definition itersolverlib.cc:25
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition matrix_kernels.h:26
fix_condition
static int fix_condition(VectorType &vec, const double cndno=1e-12)
Definition svd_solver.h:296
svd_solve
VectorType svd_solve(const MatrixType &mat, const VectorType &b, const double cndno=1e-12)
The interface to use the SVD solver: Solve the equation mat * x = b.
Definition svd_solver.h:324
sv_decomp_backsub
void sv_decomp_backsub(const MatrixType &U, const MatrixType &V, const VectorType &W, const VectorType &b, VectorType &x)
Definition svd_solver.h:268
SVDREAL
#define SVDREAL(x)
Definition svd_solver.h:26
SVDFLOAT
#define SVDFLOAT
svd_solver.h
Definition svd_solver.h:18
sv_decomp
void sv_decomp(MatrixType &A, MatrixType &V, VectorType &W)
Singular Value Decomposition Decomposes matrix A into U * diag(W) * V^T A is overwritten with U.
Definition svd_solver.h:53
sv_decomp_pythag
double sv_decomp_pythag(ARG a, ARG b)
Calculate sqrt (a^2 + b^2) avoiding under/overflows.
Definition svd_solver.h:31
Generated by doxygen 1.16.1