~ubuntu-branches/ubuntu/saucy/rheolef/saucy-proposed : revision 9

1

#ifndef _RHEOLEF_TENSOR_SVD_ALGO_H

2

#define _RHEOLEF_TENSOR_SVD_ALGO_H

3

///

4

/// This file is part of Rheolef.

5

///

6

7

///

8

/// Rheolef is free software; you can redistribute it and/or modify

9

/// it under the terms of the GNU General Public License as published by

10

/// the Free Software Foundation; either version 2 of the License, or

11

/// (at your option) any later version.

12

///

13

/// Rheolef is distributed in the hope that it will be useful,

14

/// but WITHOUT ANY WARRANTY; without even the implied warranty of

15

/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

16

/// GNU General Public License for more details.

17

///

18

/// You should have received a copy of the GNU General Public License

19

/// along with Rheolef; if not, write to the Free Software

20

/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

21

///

22

/// =========================================================================

23

#include <cmath>

24

namespace rheolef {

25

26

template<class T>

27

inline const T MIN(const T &a, const T &b)

28

{return b < a ? (b) : (a);}

29

template<class T>

30

inline const T MAX(const T &a, const T &b)

31

{return b > a ? (b) : (a);}

32

template<class T>

33

inline const T SIGN(const T &a, const T &b)

34

{return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}

35

template<class T>

36

inline const T SQR(const T a) {return a*a;}

37

38

template<class T>

39

T pythag(const T& a, const T& b)

40

{

41

T absa=std::fabs(a);

42

T absb=std::fabs(b);

43

if (absa > absb) return absa*::sqrt(1.0+SQR(absb/absa));

44

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

45

}

46

47

template <class Matrix1, class Matrix2, class Vector, class T, class Size>

48

void

49

svdcmp (Matrix1& a, Vector& w, Matrix2& v, Size anrow, Size ancol, T dummy)

50

{

51

bool flag;

52

int i,its,j,jj,k,l,nm;

53

T anorm,c,f,g,h,s,scale,x,y,z;

54

55

int m=anrow;

56

int n=ancol;

57

T rv1 [n]; // automatically allocated and destroyed

58

g=scale=anorm=0.0;

59

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

60

l=i+2;

61

rv1[i]=scale*g;

62

g=s=scale=0.0;

63

if (i < m) {

64

for (k=i;k<m;k++) scale += std::fabs(a[k][i]);

65

if (scale != 0.0) {

66

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

67

a[k][i] /= scale;

68

s += a[k][i]*a[k][i];

69

}

70

f=a[i][i];

71

g = -SIGN(::sqrt(s),f);

72

h=f*g-s;

73

a[i][i]=f-g;

74

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

75

for (s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];

76

f=s/h;

77

for (k=i;k<m;k++) a[k][j] += f*a[k][i];

78

}

79

for (k=i;k<m;k++) a[k][i] *= scale;

80

}

81

}

82

w[i]=scale *g;

83

g=s=scale=0.0;

84

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

85

for (k=l-1;k<n;k++) scale += std::fabs(a[i][k]);

86

if (scale != 0.0) {

87

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

88

a[i][k] /= scale;

89

s += a[i][k]*a[i][k];

90

}

91

f=a[i][l-1];

92

g = -SIGN(::sqrt(s),f);

93

h=f*g-s;

94

a[i][l-1]=f-g;

95

for (k=l-1;k<n;k++) rv1[k]=a[i][k]/h;

96

for (j=l-1;j<m;j++) {

97

for (s=0.0,k=l-1;k<n;k++) s += a[j][k]*a[i][k];

98

for (k=l-1;k<n;k++) a[j][k] += s*rv1[k];

99

}

100

for (k=l-1;k<n;k++) a[i][k] *= scale;

101

}

102

}

103

anorm=MAX(anorm,(std::fabs(w[i])+std::fabs(rv1[i])));

104

}

105

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

106

if (i < n-1) {

107

if (g != 0.0) {

108

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

109

v[j][i]=(a[i][j]/a[i][l])/g;

110

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

111

for (s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];

112

for (k=l;k<n;k++) v[k][j] += s*v[k][i];

113

}

114

}

115

for (j=l;j<n;j++) v[i][j]=v[j][i]=0.0;

116

}

117

v[i][i]=1.0;

118

g=rv1[i];

119

l=i;

120

}

121

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

122

l=i+1;

123

g=w[i];

124

for (j=l;j<n;j++) a[i][j]=0.0;

125

if (g != 0.0) {

126

g=1.0/g;

127

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

128

for (s=0.0,k=l;k<m;k++) s += a[k][i]*a[k][j];

129

f=(s/a[i][i])*g;

130

for (k=i;k<m;k++) a[k][j] += f*a[k][i];

131

}

132

for (j=i;j<m;j++) a[j][i] *= g;

133

} else for (j=i;j<m;j++) a[j][i]=0.0;

134

++a[i][i];

135

}

136

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

137

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

138

flag=true;

139

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

140

nm=l-1;

141

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

142

flag=false;

143

break;

144

}

145

if (std::fabs(w[nm])+anorm == anorm) break;

146

}

147

if (flag) {

148

c=0.0;

149

s=1.0;

150

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

151

f=s*rv1[i];

152

rv1[i]=c*rv1[i];

153

if (std::fabs(f)+anorm == anorm) break;

154

g=w[i];

155

h=pythag(f,g);

156

w[i]=h;

157

h=1.0/h;

158

c=g*h;

159

s = -f*h;

160

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

161

y=a[j][nm];

162

z=a[j][i];

163

a[j][nm]=y*c+z*s;

164

a[j][i]=z*c-y*s;

165

}

166

}

167

}

168

z=w[k];

169

if (l == k) {

170

if (z < 0.0) {

171

w[k] = -z;

172

for (j=0;j<n;j++) v[j][k] = -v[j][k];

173

}

174

break;

175

}

176

if (its == 29) error_macro("no convergence in 30 svdcmp iterations");

177

x=w[l];

178

nm=k-1;

179

y=w[nm];

180

g=rv1[nm];

181

h=rv1[k];

182

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

183

g=pythag(f,1.0);

184

f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;

185

c=s=1.0;

186

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

187

i=j+1;

188

g=rv1[i];

189

y=w[i];

190

h=s*g;

191

g=c*g;

192

z=pythag(f,h);

193

rv1[j]=z;

194

c=f/z;

195

s=h/z;

196

f=x*c+g*s;

197

g=g*c-x*s;

198

h=y*s;

199

y *= c;

200

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

201

x=v[jj][j];

202

z=v[jj][i];

203

v[jj][j]=x*c+z*s;

204

v[jj][i]=z*c-x*s;

205

}

206

z=pythag(f,h);

207

w[j]=z;

208

if (z) {

209

z=1.0/z;

210

c=f*z;

211

s=h*z;

212

}

213

f=c*g+s*y;

214

x=c*y-s*g;

215

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

216

y=a[jj][j];

217

z=a[jj][i];

218

a[jj][j]=y*c+z*s;

219

a[jj][i]=z*c-y*s;

220

}

221

}

222

rv1[l]=0.0;

223

rv1[k]=f;

224

w[k]=x;

225

}

226

}

227

}

228

}// namespace rheolef

229

#endif // _RHEOLEF_TENSOR_SVD_ALGO_H