IT++ Logo
svd.cpp
Go to the documentation of this file.
1 
29 #ifndef _MSC_VER
30 # include <itpp/config.h>
31 #else
32 # include <itpp/config_msvc.h>
33 #endif
34 
35 #if defined(HAVE_LAPACK)
36 # include <itpp/base/algebra/lapack.h>
37 #endif
38 
39 #include <itpp/base/algebra/svd.h>
40 
41 
42 namespace itpp
43 {
44 
45 #if defined(HAVE_LAPACK)
46 
47 bool svd(const mat &A, vec &S)
48 {
49  char jobu = 'N', jobvt = 'N';
50  int m, n, lda, ldu, ldvt, lwork, info;
51  m = lda = ldu = A.rows();
52  n = ldvt = A.cols();
53  lwork = std::max(3 * std::min(m, n) + std::max(m, n), 5 * std::min(m, n));
54 
55  mat U, V;
56  S.set_size(std::min(m, n), false);
57  vec work(lwork);
58 
59  mat B(A);
60 
61  // The theoretical calculation of lwork above results in the minimum size
62  // needed for dgesvd_ to run to completion without having memory errors.
63  // For speed improvement it is best to set lwork=-1 and have dgesvd_
64  // calculate the best workspace requirement.
65  int lwork_tmp = -1;
66  dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
67  V._data(), &ldvt, work._data(), &lwork_tmp, &info);
68  if (info == 0) {
69  lwork = static_cast<int>(work(0));
70  work.set_size(lwork, false);
71  }
72 
73  dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
74  V._data(), &ldvt, work._data(), &lwork, &info);
75 
76  return (info == 0);
77 }
78 
79 bool svd(const cmat &A, vec &S)
80 {
81  char jobu = 'N', jobvt = 'N';
82  int m, n, lda, ldu, ldvt, lwork, info;
83  m = lda = ldu = A.rows();
84  n = ldvt = A.cols();
85  lwork = 2 * std::min(m, n) + std::max(m, n);
86 
87  cvec U, V;
88  S.set_size(std::min(m, n), false);
89  cvec work(lwork);
90  vec rwork(5*std::min(m, n));
91 
92  cmat B(A);
93 
94  // The theoretical calculation of lwork above results in the minimum size
95  // needed for zgesvd_ to run to completion without having memory errors.
96  // For speed improvement it is best to set lwork=-1 and have zgesvd_
97  // calculate the best workspace requirement.
98  int lwork_tmp = -1;
99  zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
100  V._data(), &ldvt, work._data(), &lwork_tmp, rwork._data(), &info);
101  if (info == 0) {
102  lwork = static_cast<int>(real(work(0)));
103  work.set_size(lwork, false);
104  }
105 
106  zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
107  V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info);
108 
109  return (info == 0);
110 }
111 
112 bool svd(const mat &A, mat &U, vec &S, mat &V)
113 {
114  char jobu = 'A', jobvt = 'A';
115  int m, n, lda, ldu, ldvt, lwork, info;
116  m = lda = ldu = A.rows();
117  n = ldvt = A.cols();
118  lwork = std::max(3 * std::min(m, n) + std::max(m, n), 5 * std::min(m, n));
119 
120  U.set_size(m, m, false);
121  V.set_size(n, n, false);
122  S.set_size(std::min(m, n), false);
123  vec work(lwork);
124 
125  mat B(A);
126 
127  // The theoretical calculation of lwork above results in the minimum size
128  // needed for dgesvd_ to run to completion without having memory errors.
129  // For speed improvement it is best to set lwork=-1 and have dgesvd_
130  // calculate the best workspace requirement.
131  int lwork_tmp = -1;
132  dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
133  V._data(), &ldvt, work._data(), &lwork_tmp, &info);
134  if (info == 0) {
135  lwork = static_cast<int>(work(0));
136  work.set_size(lwork, false);
137  }
138 
139  dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
140  V._data(), &ldvt, work._data(), &lwork, &info);
141 
142  V = V.T(); // This is probably slow!!!
143 
144  return (info == 0);
145 }
146 
147 bool svd(const cmat &A, cmat &U, vec &S, cmat &V)
148 {
149  char jobu = 'A', jobvt = 'A';
150  int m, n, lda, ldu, ldvt, lwork, info;
151  m = lda = ldu = A.rows();
152  n = ldvt = A.cols();
153  lwork = 2 * std::min(m, n) + std::max(m, n);
154 
155  U.set_size(m, m, false);
156  V.set_size(n, n, false);
157  S.set_size(std::min(m, n), false);
158  cvec work(lwork);
159  vec rwork(5 * std::min(m, n));
160 
161  cmat B(A);
162 
163  // The theoretical calculation of lwork above results in the minimum size
164  // needed for zgesvd_ to run to completion without having memory errors.
165  // For speed improvement it is best to set lwork=-1 and have zgesvd_
166  // calculate the best workspace requirement.
167  int lwork_tmp = -1;
168  zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
169  V._data(), &ldvt, work._data(), &lwork_tmp, rwork._data(), &info);
170  if (info == 0) {
171  lwork = static_cast<int>(real(work(0)));
172  work.set_size(lwork, false);
173  }
174 
175  zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
176  V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info);
177 
178  V = V.H(); // This is slow!!!
179 
180  return (info == 0);
181 }
182 
183 #else
184 
185 bool svd(const mat &A, vec &S)
186 {
187  it_error("LAPACK library is needed to use svd() function");
188  return false;
189 }
190 
191 bool svd(const cmat &A, vec &S)
192 {
193  it_error("LAPACK library is needed to use svd() function");
194  return false;
195 }
196 
197 bool svd(const mat &A, mat &U, vec &S, mat &V)
198 {
199  it_error("LAPACK library is needed to use svd() function");
200  return false;
201 }
202 
203 bool svd(const cmat &A, cmat &U, vec &S, cmat &V)
204 {
205  it_error("LAPACK library is needed to use svd() function");
206  return false;
207 }
208 
209 #endif // HAVE_LAPACK
210 
211 vec svd(const mat &A)
212 {
213  vec S;
214  svd(A, S);
215  return S;
216 }
217 
218 vec svd(const cmat &A)
219 {
220  vec S;
221  svd(A, S);
222  return S;
223 }
224 
225 } //namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:19 for IT++ by Doxygen 1.8.2