libtbcinumlib/html/svd__solver_8h_source.html

 #ifndef TBCI_SOLVER_SVD_SOLVER_H

 #define TBCI_SOLVER_SVD_SOLVER_H


 #include "tbci/matrix.h"

 #ifdef TYP

 # define SVDFLOAT TYP

 #else

 # define SVDFLOAT typename MatrixType::value_type

 #endif


 NAMESPACE_TBCI


 #ifdef NEED_BCPLX_WORKAROUND

 # define SVDREAL(x) CPLX__ real(x)

 #else

 # define SVDREAL(x) (x)

 #endif


 template <typename ARG>

 inline double sv_decomp_pythag(ARG a, ARG b)

 {

         const double absa = MATH__ fabs(a);

         const double absb = MATH__ fabs(b);

         if (absa > absb)

                 return absa * MATH__ sqrt(1.0 + (absb/absa)*(absb/absa));

         else

                 return (absb==0.0? 0.0: absb * MATH__ sqrt(1.0 + (absa/absb)*(absa/absb)));

 }


 //TODO: Check which of the variables used should be double rather

 // then SVDFLOAT (which can be complex, e.g.)


 INST2(template <typename MAT, typename VEC > class MAT friend void sv_decomp (MAT &, MAT &, VEC &);)

 template <typename MatrixType /*= TBCI__ Matrix<double>*/,

           typename VectorType /*= TBCI__ Vector<double>*/ >

 void sv_decomp (MatrixType & A, MatrixType & V, VectorType & W)

 {

         const int m = A.rows();

         const int n = A.columns();

         int l = 0, i, k;        // M$ annoyance (for scoping bug)

         VectorType rv1(n);

         SVDFLOAT g = 0.0, scale = 0.0;

         double anorm = 0.0;

         for (i = 0; i < n; ++i) {

                 SVDFLOAT s = 0.0;

                 l = i+1;

                 rv1[i] = scale*g;

                 g = scale = 0.0;

                 if (i < m) {

                         for (k = i; k < m; ++k)

                                 scale += MATH__ fabs(A(k,i));

                         if (scale != (SVDFLOAT)0) {

                                 for (k = i; k < m; ++k) {

                                         A(k,i) /= scale;

                                         s += A(k,i)*A(k,i);

                                 }

                                 SVDFLOAT f = A(i,i);

                                 g = SVDREAL(f)>SVDREAL((SVDFLOAT)0) ? - MATH__ sqrt(s) : MATH__ sqrt(s);

                                 SVDFLOAT h = f*g - s;

                                 A(i,i) = f - g;

                                 for (int j = l; j < n; ++j) {

                                         s = (SVDFLOAT)0;

                                         for (k = i; k < m; ++k)

                                                 s += A(k,i)*A(k,j);

                                         SVDFLOAT f = s/h;

                                         for (k = i; k < m; ++k)

                                                 A(k,j) += f*A(k,i);

                                 }

                                 for (k = i; k < m; ++k)

                                         A(k,i) *= scale;

                         }

                 }

                 W[i] = scale * g;

                 g = s = scale = 0.0;

                 if (i < m && i != (n-1)) {

                         for (k = l; k < n; ++k)

                                 scale += MATH__ fabs(A(i,k));

                         if (scale != (SVDFLOAT)0) {

                                 for (k = l; k < n; ++k) {

                                         A(i,k) /= scale;

                                         s += A(i,k)*A(i,k);

                                 }

                                 SVDFLOAT f = A(i,l);

                                 g = SVDREAL(f)>SVDREAL((SVDFLOAT)0) ? - MATH__ sqrt(s) : MATH__ sqrt(s);

                                 SVDFLOAT h = f*g - s;

                                 A(i,l) = f - g;

                                 for (k=l; k < n; ++k)

                                         rv1[k] = A(i,k)/h;

                                 for (int j = l; j < m; ++j) {

                                         s = 0.0;

                                         for (k = l; k < n; ++k)

                                                 s += A(j,k)*A(i,k);

                                         for (k = l; k < n; ++k)

                                                 A(j,k) += s*rv1[k];

                                 }

                                 for(k = l; k < n; ++k)

                                         A(i,k) *= scale;

                         }

                 }

                 anorm = MAX(anorm, MATH__ fabs(W[i]) + MATH__ fabs(rv1[i]));

         }

         for (i = n-1; i >= 0; --i) {

                 if (i < (n-1)) {

                         if (g != (SVDFLOAT)0) {

                                 int j;

                                 for (j = l; j < n; ++j )

                                         V(j,i) = (A(i,j)/A(i,l))/g;

                                 for (j = l; j < n; ++j) {

                                         SVDFLOAT s = 0.0;

                                         for (k = l; k < n; ++k )

                                                 s += A(i,k)*V(k,j);

                                         for (k = l; k < n; ++k )

                                                 V(k,j) += s*V(k,i);

                                 }

                         }

                         for (int j = l; j < n; ++j)

                                 V(i,j) = V(j,i) = 0.0;

                 }

                 V(i,i) = 1.0;

                 g = rv1[i];

                 l = i;

         }

         for (i = MIN(m,n)-1; i >= 0; --i) {

                 int j;

                 int l = i+1;

                 g = W[i];

                 for (j = l; j < n; ++j)

                         A(i,j) = 0.0;

                 if (g != (SVDFLOAT)0) {

                         g = 1.0/g;

                         for (j = l; j < n; ++j) {

                                 SVDFLOAT s = 0.0;

                                 for (k = l; k < m; ++k)

                                         s += A(k,i)*A(k,j);

                                 SVDFLOAT f = (s/A(i,i))*g;

                                 for (k = i; k < m; ++k)

                                         A(k,j) += f*A(k,i);

                         }

                         for (j = i; j < m; ++j)

                                 A(j,i) *= g;

                 } else

                         for (j = i; j < m; ++j)

                                 A(j,i) = 0.0;

                 A(i,i) += 1.0;

         }

         for (k = n-1; k >= 0; --k) {

                 for (int its = 0; its < 30; ++its) {

                         bool flag = true;

                         int nm=0;

                         for (l = k; l >= 0; --l) {

                                 nm = l-1; // note that rv1[0] is always zero

                                 if ((MATH__ fabs(rv1[l]) + anorm) == anorm) {

                                         flag = false;

                                         break;

                                 }

                                 if ((MATH__ fabs(W[nm]) + anorm) == anorm)

                                         break;

                         }

                         if (flag) {

                                 SVDFLOAT c = 0.0;

                                 SVDFLOAT s = 1.0;

                                 for (i = l; i < k; ++i) {

                                         SVDFLOAT f = s*rv1[i];

                                         rv1[i] *= c;

                                         if ((MATH__ fabs(f) + anorm) == anorm)

                                                 break;

                                         SVDFLOAT g = W[i];

                                         SVDFLOAT h = sv_decomp_pythag(f, g);

                                         W[i] = h;

                                         h = 1.0/h;

                                         c =  g*h;

                                         s = -f*h;

                                         for (int j = 0; j < m; ++j) {

                                                 SVDFLOAT y = A(j,nm);

                                                 SVDFLOAT z = A(j,i);

                                                 A(j,nm) = y*c + z*s;

                                                 A(j,i)  = z*c - y*s;

                                         }

                                 }

                         }

                         SVDFLOAT z = W[k];

                         if (l == k) {   // convergence

                                 if (SVDREAL(z) < SVDREAL((SVDFLOAT)0)) {

                                         W[k] = -z;

                                         for (int j = 0; j < n; ++j)

                                                 V(j,k) = -V(j,k);

                                 }

                                 break;

                         }

                         if (its == 29)

                                 STD__ cerr << "SVD: No convergence in 30 iterations.\n";

                         SVDFLOAT x = W[l];

                         nm = k-1;

                         SVDFLOAT y = W[nm];

                         SVDFLOAT g = rv1[nm];

                         SVDFLOAT h = rv1[k];

                         SVDFLOAT f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);

                         g = sv_decomp_pythag(f, (SVDFLOAT)1.0);

                         f = ((x-z)*(x+z) +

                              h*((y / (f + (SVDREAL(f)>SVDREAL((SVDFLOAT)0) ? MATH__ fabs(g) : - MATH__ fabs(g))))-h))/x;

                         SVDFLOAT c = 1.0;

                         SVDFLOAT s = 1.0;

                         for (int j = l; j <= nm; ++j) {

                                 int jj; // You know who annoyance

                                 i = j+1;

                                 SVDFLOAT g = rv1[i];

                                 SVDFLOAT y = W[i];

                                 SVDFLOAT h = s*g;

                                 g*=c;

                                 SVDFLOAT z = sv_decomp_pythag(f, h);

                                 rv1[j] = z;

                                 c = f/z;

                                 s = h/z;

                                 f = x*c + g*s;

                                 g = g*c - x*s;

                                 h = y*s;

                                 y *= c;

                                 for (jj = 0; jj < n; ++jj) {

                                         SVDFLOAT x = V(jj,j);

                                         SVDFLOAT z = V(jj,i);

                                         V(jj,j) = x*c + z*s;

                                         V(jj,i) = z*c - x*s;

                                 }

                                 z = sv_decomp_pythag(f, h);

                                 W[j] = z;

                                 if (z != (SVDFLOAT)0) {

                                         z = 1.0/z;

                                         c = f*z;

                                         s = h*z;

                                 }

                                 f = c*g + s*y;

                                 x = c*y - s*g;

                                 for (jj = 0; jj < m; ++jj) {

                                         SVDFLOAT y = A(jj,j);

                                         SVDFLOAT z = A(jj,i);

                                         A(jj,j) = y*c + z*s;

                                         A(jj,i) = z*c - y*s;

                                 }

                         }

                         rv1[l] = 0.0;

                         rv1[k] = f;

                         W[k] = x;

                 }

         }

 }


 INST2(template <typename MAT, typename VEC > class MAT friend void sv_decomp_backsub (const MAT &, const MAT &, const VEC &, const VEC &, VEC &);)

 template <typename MatrixType /*= TBCI__ Matrix<double>*/,

           typename VectorType /*= TBCI__ Vector<double>*/ >

 void sv_decomp_backsub (const MatrixType & U, const MatrixType & V,

                         const VectorType & W, const VectorType & b,

                         VectorType & x)

 {

         const int n = W.size();

         const int m = b.size();

         VectorType tmp(n);

         int j;

         for (j = 0; j < n; ++j) {

                 SVDFLOAT s = 0.0;

                 if (W[j] != (SVDFLOAT)0) {

                         for (int i = 0; i < m; ++i)

                                 s += U(i,j)*b[i];

                         s /= W[j];

                 }

                 tmp[j] = s;

         }

         for (j = 0; j < n; ++j) {

                 SVDFLOAT s = 0.0;

                 for (int i = 0; i < n; ++i)

                         s += V(j,i)*tmp[i];

                 x[j] = s;

         }

 }


 INST(template <typename VEC > class VEC friend /*static*/ int fix_condition (VEC& vec, const double cndno);)

 /* Fix those singular values that are close to zero */

 template <typename VectorType >

 static int fix_condition (VectorType& vec, const double cndno = 1e-12)

 {

         int fixed = 0; unsigned i;

         double maxsv = 0.0;

         const unsigned size = vec.size();

         //typename VectorType::value_type * RESTRICT const VEC = &vec.setval(0);

         for (i = 0; i < size; ++i)

                 maxsv = MAX(maxsv, MATH__ fabs(vec[i]));

         double minsv = maxsv * cndno;

         for (i = 0; i < size; ++i)

                 if (MATH__ fabs(vec[i]) < minsv) {

                         fixed++;

                         vec[i] = (typename VectorType::value_type)0.0;

                 }

         return fixed;

 }


 INST2(template <typename MAT, typename VEC > class MAT friend VEC svd_solve (const MAT& mat, const VEC& b, const double cndno);)

 template <typename MatrixType, typename VectorType>

 VectorType svd_solve (const MatrixType& mat, const VectorType& b, const double cndno = 1e-12)

 {

         MatrixType A(mat);

         MatrixType V(mat.columns(), mat.columns());

         VectorType W(mat.columns());

         sv_decomp <MatrixType, VectorType > (A, V, W);

         fix_condition(W, cndno);

         // we could calc mat^-1 = V * diag (1/W) * A * A^T now

         VectorType x (mat.columns());

         sv_decomp_backsub(A, V, W, b, x);

         return x;

 }


 // Optimized version for our own matrix & vector classes, using TBCI

 INST2(template <typename MAT, typename T > class MAT friend TVector<T> sv_decomp_backsub (const MAT &, const MAT &, const Vector<T> &, const Vector<T> &);)

 template <typename MatrixType, typename T >

 TVector<T> sv_decomp_backsub (const MatrixType & U, const MatrixType & V,

                               const Vector<T> & W, const Vector<T> & b)

 {

         const int n = W.size();

         const int m = b.size();

         TVector<T> x (n);

         Vector<T> tmp(n);

         int j;

         for (j = 0; j < n; ++j) {

                 SVDFLOAT s = 0.0;

                 if (W[j] != (SVDFLOAT)0) {

                         for (int i = 0; i < m; ++i)

                                 s += U(i,j)*b[i];

                         s /= W[j];

                 }

                 tmp[j] = s;

         }

         for (j = 0; j < n; ++j) {

                 SVDFLOAT s = 0.0;

                 for (int i = 0; i < n; ++i)

                         s += V(j,i)*tmp[i];

                 x.set (s, j);

         }

         return x;

 }


 INST2(template <typename MAT, typename T > class MAT friend TVector<T> svd_solve (const MAT& mat, const Vector<T>& b, const double cndno);)

 template <typename MatrixType, typename T >

 TVector<T> svd_solve (const MatrixType& mat, const Vector<T>& b, const double cndno = 1e-12)

 {

         MatrixType A(mat);

         MatrixType V(mat.columns(), mat.columns());

         Vector<T>  W(mat.columns());

         sv_decomp <MatrixType, Vector<T> > (A, V, W);

         fix_condition (W, cndno);

         return sv_decomp_backsub (A, V, W, b);

 }


 #undef SVDFLOAT


 NAMESPACE_END

 #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:97

c
return c
Definition: f_matrix.h:760

std::fabs
double fabs(const int a)
Definition: basics.h:1215

MIN
#define MIN(a, b)
Definition: basics.h:655

NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition: basics.h:317

VEC
#define VEC
Definition: bdmatsolverlib.cc:28

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_backsub
void sv_decomp_backsub(const MatrixType &U, const MatrixType &V, const VectorType &W, const VectorType &b, VectorType &x)
Definition: svd_solver.h:268

TVector::set
T & set(const T &val, const unsigned long i) const
Definition: vector.h:194

MAT
#define MAT
Definition: bdmatsolverlib.cc:27

z
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > & z
Definition: LM_fit.h:172

std::sqrt
double sqrt(const int a)
Definition: basics.h:1216

U
#define U
Definition: bdmatlib.cc:21

MatrixType
#define MatrixType
Definition: itersolverlib.cc:39

fix_condition
static int fix_condition(VectorType &vec, const double cndno=1e-12)
Definition: svd_solver.h:296

SVDFLOAT
#define SVDFLOAT
svd_solver.h
Definition: svd_solver.h:18

b
F_TMatrix< T > b
Definition: f_matrix.h:736

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

INST2
#define INST2(x, y)
Definition: basics.h:239

TVector::size
unsigned long size() const
Definition: vector.h:104

flag
long int flag
Definition: f2c.h:66

INST
#define INST(x)
Definition: basics.h:238

i
int i
Definition: LM_fit.h:71

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

SVDREAL
#define SVDREAL(x)
Definition: svd_solver.h:26

STD__
#define STD__
Definition: basics.h:338

TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition: bvector.h:52

MAX
#define MAX(a, b)
Definition: basics.h:656

NAMESPACE_END
#define NAMESPACE_END
Definition: basics.h:323

Vector
Definition: bvector.h:54

y
const Vector< T > Vector< T > Vector< T > Vector< T > & y
Definition: LM_fit.h:172

a
const unsigned TMatrix< T > const Matrix< T > * a
Definition: matrix_kernels.h:25

MATH__
#define MATH__
Definition: basics.h:339