Revision 795 branches/replacement/src/linad99/dmat42.cpp
dmat42.cpp (revision 795)  

2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
5 
* Copyright (c) 2009 ADMB Foundation 

5 
* Copyright (c) 2009, 2010 ADMB Foundation


6  6 
*/ 
7 
/** 

8 
* \file 

9 
* This file deals with the Singular Value Decomposition 

10 
* of a matrix. The format of the decomposition follows 

11 
* the format given in "Numerical Recipes in C", 2nd edition, 

12 
* Press, Teukolsky, Vetterling, Flannery, section 2.6. 

13 
*/ 

14  
7  15 
#include <fvar.hpp> 
8  16 
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : fabs(a)) 
9  17  
...  ...  
49  57 
a(_a), w(_w), v(_v) 
50  58 
{} 
51  59  
52 
/** Singular value decomposition. 

53 
*/ 

60 
/** 

61 
* Singular value decomposition. 

62 
*/ 

54  63 
sing_val_decomp singval_decomp(const dmatrix &_a) 
55  64 
{ 
56  65 
if (_a.indexmin() !=1 ) 
...  ...  
58  67 
cerr << "index error in singval_decomp" << endl; 
59  68 
ad_exit(1); 
60  69 
} 
61 
int m=_a.indexmax();


62 
int n=_a(1).indexmax();


70 
int m = _a.indexmax();


71 
int n = _a(1).indexmax();


63  72 
dmatrix a(1,m,1,n); 
64  73 
a=_a; 
65  74 
dvector w(1,n); 
66  75 
dmatrix u(1,m,1,n); 
67  76 
dmatrix v(1,n,1,n); 
68  77  
69 
int k; 

70 
double eps, tol; 

71 
eps = tol = 1.e12; 

72 
k = svd(m,n,1,1,eps,tol,a,w,u,v); 

78 
double eps = 1.e12; 

79 
double tol = eps; 

80 
int k = svd(m,n,1,1,eps,tol,a,w,u,v); 

73  81 
if(k!=0) 
74  82 
{ 
75  83 
cerr << "Error in singval_decomp in iteration " << k << endl; 
...  ...  
78  86 
return sing_val_decomp(u,w,v); 
79  87 
} 
80  88  
81  
82  
83  
84  
85  
86  
87 
/** Singular value decomposition. 

89 
/** 

90 
* Singular value decomposition. 

88  91 
* \ingroup svd 
89  92 
* 
90  93 
* \param m \f$m\f$ the number of rows of \f$A\f$ 
91  94 
* \param n \f$n\f$ the number of columns of \f$A\f$ 
92 
* \param withu true if the \f$U\f$part is wanted 

93 
* \param withv true if the \f$V\f$part is wanted 

94 
* \param eps 

95 
* \param tol 

95 
* \param withu true if the \f$U\f$part is wanted (true=1, false=0)


96 
* \param withv true if the \f$V\f$part is wanted (true=1, false=0)


97 
* \param eps i.e Epsilon


98 
* \param tol the tolerance used


96  99 
* \param aa \f$A\f$ 
97  100 
* \param _q \f$q\f$ 
98  101 
* \param _u \f$U\f$ 
...  ...  
145  148 
return k; 
146  149 
} 
147  150  
148 
/** Singular value decomposition. 

151 
/** 

152 
* Singular value decomposition. 

149  153 
* \ingroup svd 
150  154 
* 
151  155 
* Used to find the svd of a matrix when \f$m<n\f$. 
...  ...  
422  426 
return retval; 
423  427 
} 
424  428  
425 
/** \ingroup svd 
Also available in: Unified diff