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
40
41
42namespace itpp
43{
44
45#if defined(HAVE_LAPACK)
46
47bool 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
79bool 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
112bool 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
147bool 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
185bool svd(const mat &A, vec &S)
186{
187 it_error("LAPACK library is needed to use svd() function");
188 return false;
189}
190
191bool svd(const cmat &A, vec &S)
192{
193 it_error("LAPACK library is needed to use svd() function");
194 return false;
195}
196
197bool 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
203bool 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
211vec svd(const mat &A)
212{
213 vec S;
214 svd(A, S);
215 return S;
216}
217
218vec svd(const cmat &A)
219{
220 vec S;
221 svd(A, S);
222 return S;
223}
224
225} //namespace itpp
#define it_error(s)
Abort unconditionally.
Definition: itassert.h:126
bool svd(const mat &A, vec &S)
Get singular values s of a real matrix A using SVD.
Definition: svd.cpp:185
vec real(const cvec &data)
Real part of complex values.
Definition: elem_math.cpp:157
itpp namespace
Definition: itmex.h:37
Definitions of Singular Value Decompositions.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.5