LOOS  v2.3.2
MatrixOps.cpp
1 /*
2  MatrixOps.cpp
3 */
4 
5 
6 
7 /*
8  This file is part of LOOS.
9 
10  LOOS (Lightweight Object-Oriented Structure library)
11  Copyright (c) 2009, Tod D. Romo, Alan Grossfield
12  Department of Biochemistry and Biophysics
13  School of Medicine & Dentistry, University of Rochester
14 
15  This package (LOOS) is free software: you can redistribute it and/or modify
16  it under the terms of the GNU General Public License as published by
17  the Free Software Foundation under version 3 of the License.
18 
19  This package is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  GNU General Public License for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with this program. If not, see <http://www.gnu.org/licenses/>.
26 */
27 
28 
29 
30 #include <MatrixOps.hpp>
31 
32 
33 namespace loos {
34  namespace Math {
35 
37  {
38  // Compute [U,D] = eig(M)
39 
40  f77int n = M.rows();
41  char jobz = 'V';
42  char uplo = 'L';
43  f77int lda = n;
44  double dummy;
45  DoubleMatrix W(n, 1);
46  f77int lwork = -1;
47  f77int info;
48 
49  dsyev_(&jobz, &uplo, &n, M.get(), &lda, W.get(), &dummy, &lwork, &info);
50  if (info != 0)
51  throw(NumericalError("DSYEV failed to give an estimate of space required"), info);
52 
53  lwork = static_cast<f77int>(dummy);
54  double* work = new double[lwork+1];
55 
56  dsyev_(&jobz, &uplo, &n, M.get(), &lda, W.get(), work, &lwork, &info);
57  if (info != 0)
58  throw(NumericalError("DSYEV reported an error"), info);
59 
60  return(W);
61  }
62 
63 
65  {
66  // Compute [U,D] = eig(M)
67 
68  f77int n = M.rows();
69  char jobz = 'V';
70  char uplo = 'L';
71  f77int lda = n;
72  float dummy;
73  RealMatrix W(n, 1);
74  f77int lwork = -1;
75  f77int info;
76 
77  ssyev_(&jobz, &uplo, &n, M.get(), &lda, W.get(), &dummy, &lwork, &info);
78  if (info != 0)
79  throw(NumericalError("SSYEV failed to give an estimate of space required"), info);
80 
81  lwork = static_cast<f77int>(dummy);
82  float* work = new float[lwork+1];
83 
84  ssyev_(&jobz, &uplo, &n, M.get(), &lda, W.get(), work, &lwork, &info);
85  if (info != 0)
86  throw(NumericalError("SSYEV reported an error"), info);
87 
88  return(W);
89  }
90 
91 
92 
93 
94  boost::tuple<RealMatrix, RealMatrix, RealMatrix> svd(RealMatrix& M) {
95  f77int m = M.rows();
96  f77int n = M.cols();
97  f77int sn = m<n ? m : n;
98  long estimate = (m*m + n*n + m*n + sn) * sizeof(float);
99  // This is somewhat vestigial... We continue to estimate the
100  // total storage required for the SVD in case we want to do
101  // something with it in the future...
102 
103 
104 
105  char jobu = 'A', jobvt = 'A';
106  f77int lda = m, ldu = m, ldvt = n, lwork= -1, info;
107  float prework[10], *work;
108 
109  RealMatrix U(m, m);
110  RealMatrix S(sn, 1);
111  RealMatrix Vt(n, n);
112 
113  sgesvd_(&jobu, &jobvt, &m, &n, M.get(), &lda, S.get(), U.get(), &ldu, Vt.get(), &ldvt, prework, &lwork, &info);
114  if (info != 0)
115  throw(NumericalError("SGESVD estimate reported an error"), info);
116 
117  lwork = (f77int)prework[0];
118  estimate += lwork * sizeof(double);
119  work = new float[lwork];
120 
121  sgesvd_(&jobu, &jobvt, &m, &n, M.get(), &lda, S.get(), U.get(), &ldu, Vt.get(), &ldvt, work, &lwork, &info);
122  if (info != 0)
123  throw(NumericalError("SGESVD reported an error"), info);
124 
125 
126  delete[] work;
127 
128  boost::tuple<RealMatrix, RealMatrix, RealMatrix> result(U, S, Vt);
129  return(result);
130  }
131 
132 
133 
134  boost::tuple<DoubleMatrix, DoubleMatrix, DoubleMatrix> svd(DoubleMatrix& M) {
135  f77int m = M.rows();
136  f77int n = M.cols();
137  f77int sn = m<n ? m : n;
138  long estimate = (m*m + n*n + m*n + sn) * sizeof(double);
139 
140 
141  char jobu = 'A', jobvt = 'A';
142  f77int lda = m, ldu = m, ldvt = n, lwork= -1, info;
143  double prework[10], *work;
144 
145  DoubleMatrix U(m, m);
146  DoubleMatrix S(sn, 1);
147  DoubleMatrix Vt(n, n);
148 
149  dgesvd_(&jobu, &jobvt, &m, &n, M.get(), &lda, S.get(), U.get(), &ldu, Vt.get(), &ldvt, prework, &lwork, &info);
150  if (info != 0)
151  throw(NumericalError("DGESVD estimate reported an error", info));
152 
153 
154  lwork = (f77int)prework[0];
155  estimate += lwork * sizeof(double);
156  work = new double[lwork];
157 
158  dgesvd_(&jobu, &jobvt, &m, &n, M.get(), &lda, S.get(), U.get(), &ldu, Vt.get(), &ldvt, work, &lwork, &info);
159  if (info != 0)
160  throw(NumericalError("DGESVD reported an error", info));
161 
162  delete[] work;
163 
164  boost::tuple<DoubleMatrix, DoubleMatrix, DoubleMatrix> result(U, S, Vt);
165  return(result);
166  }
167 
168  // Multiply two matrices using BLAS
169 
170  RealMatrix MMMultiply(const RealMatrix& A, const RealMatrix& B, const bool transa, const bool transb) {
171 
172  f77int m = transa ? A.cols() : A.rows();
173  f77int n = transb ? B.rows() : B.cols();
174  f77int k = transa ? A.rows() : A.cols();
175  float alpha = 1.0;
176  float beta = 0.0;
177 
178  f77int lda = transa ? k : m;
179  f77int ldb = transb ? n : k;
180  f77int ldc = m;
181 
182 
183  RealMatrix C(m, n);
184 
185 #if defined(__linux__) | defined(__CYGWIN__) || defined(__FreeBSD__)
186  char ta = (transa ? 'T' : 'N');
187  char tb = (transb ? 'T' : 'N');
188 
189  sgemm_(&ta, &tb, &m, &n, &k, &alpha, A.get(), &lda, B.get(), &ldb, &beta, C.get(), &ldc);
190 #else
191  cblas_sgemm(CblasColMajor, transa ? CblasTrans : CblasNoTrans, transb ? CblasTrans : CblasNoTrans,
192  m, n, k, alpha, A.get(), lda, B.get(), ldb, beta, C.get(), ldc);
193 #endif
194 
195  return(C);
196  }
197 
198 
199  DoubleMatrix MMMultiply(const DoubleMatrix& A, const DoubleMatrix& B, const bool transa, const bool transb) {
200 
201  f77int m = transa ? A.cols() : A.rows();
202  f77int n = transb ? B.rows() : B.cols();
203  f77int k = transa ? A.rows() : A.cols();
204  double alpha = 1.0;
205  double beta = 0.0;
206 
207  f77int lda = transa ? k : m;
208  f77int ldb = transb ? n : k;
209  f77int ldc = m;
210 
211 
212  DoubleMatrix C(m, n);
213 
214 #if defined(__linux__) || defined(__CYGWIN__) || defined(__FreeBSD__)
215  char ta = (transa ? 'T' : 'N');
216  char tb = (transb ? 'T' : 'N');
217 
218  dgemm_(&ta, &tb, &m, &n, &k, &alpha, A.get(), &lda, B.get(), &ldb, &beta, C.get(), &ldc);
219 #else
220  cblas_dgemm(CblasColMajor, transa ? CblasTrans : CblasNoTrans, transb ? CblasTrans : CblasNoTrans,
221  m, n, k, alpha, A.get(), lda, B.get(), ldb, beta, C.get(), ldc);
222 #endif
223 
224  return(C);
225  }
226 
227 
228  // Pseudo-inverse of a matrix using the SVD
229 
230  RealMatrix invert(RealMatrix& A, const float eps) {
231 
232  // The SVD (at least dgesvd) will destroy the source matrix, so we
233  // need to make an explicit copy...
234 
235  RealMatrix B(A.copy());
236  boost::tuple<RealMatrix, RealMatrix, RealMatrix> res = svd(B);
237  RealMatrix U(boost::get<0>(res));
238  RealMatrix S(boost::get<1>(res));
239  RealMatrix Vt(boost::get<2>(res));
240 
241 
242  for (uint i=0; i<Vt.rows(); ++i) {
243  double sv = S[i];
244  if (sv < eps)
245  sv = 0.0;
246  else
247  sv = 1.0 / sv;
248 
249  for (uint j=0; j<Vt.cols(); ++j)
250  Vt(i,j) *= sv;
251  }
252 
253  RealMatrix Ai = MMMultiply(Vt, U, true, true);
254  return(Ai);
255  }
256 
257 
258  DoubleMatrix invert(DoubleMatrix& A, const double eps) {
259 
260  // The SVD (at least dgesvd) will destroy the source matrix, so we
261  // need to make an explicit copy...
262 
263  DoubleMatrix B(A.copy());
264  boost::tuple<DoubleMatrix, DoubleMatrix, DoubleMatrix> res = svd(B);
265  DoubleMatrix U(boost::get<0>(res));
266  DoubleMatrix S(boost::get<1>(res));
267  DoubleMatrix Vt(boost::get<2>(res));
268 
269 
270  for (uint i=0; i<Vt.rows(); ++i) {
271  double sv = S[i];
272  if (sv < eps)
273  sv = 0.0;
274  else
275  sv = 1.0 / sv;
276 
277  for (uint j=0; j<Vt.cols(); ++j)
278  Vt(i,j) *= sv;
279  }
280 
281  DoubleMatrix Ai = MMMultiply(Vt, U, true, true);
282  return(Ai);
283  }
284 
285 
286  void operator+=(RealMatrix& A, const RealMatrix& B) {
287  if (A.rows() != B.rows() || A.cols() != B.cols())
288  throw(std::logic_error("Matrices are not the same size"));
289 
290  for (uint i=0; i<A.rows() * A.cols(); ++i)
291  A[i] += B[i];
292  }
293 
294  RealMatrix operator+(const RealMatrix& A, const RealMatrix& B) {
295  RealMatrix C(A.copy());
296  C += B;
297  return(C);
298  }
299 
300 
301  void operator+=(DoubleMatrix& A, const DoubleMatrix& B) {
302  if (A.rows() != B.rows() || A.cols() != B.cols())
303  throw(std::logic_error("Matrices are not the same size"));
304 
305  for (uint i=0; i<A.rows() * A.cols(); ++i)
306  A[i] += B[i];
307  }
308 
309  DoubleMatrix operator+(const DoubleMatrix& A, const DoubleMatrix& B) {
310  DoubleMatrix C(A.copy());
311  C += B;
312  return(C);
313  }
314 
315  // ----------------------------------------------------------
316 
317 
318  void operator-=(RealMatrix& A, const RealMatrix& B) {
319  if (A.rows() != B.rows() || A.cols() != B.cols())
320  throw(std::logic_error("Matrices are not the same size"));
321 
322  for (uint i=0; i<A.rows() * A.cols(); ++i)
323  A[i] -= B[i];
324  }
325 
326  RealMatrix operator-(const RealMatrix& A, const RealMatrix& B) {
327  RealMatrix C(A.copy());
328  C -= B;
329  return(C);
330  }
331 
332 
333 
334  void operator-=(DoubleMatrix& A, const DoubleMatrix& B) {
335  if (A.rows() != B.rows() || A.cols() != B.cols())
336  throw(std::logic_error("Matrices are not the same size"));
337 
338  for (uint i=0; i<A.rows() * A.cols(); ++i)
339  A[i] -= B[i];
340  }
341 
342  DoubleMatrix operator-(const DoubleMatrix& A, const DoubleMatrix& B) {
343  DoubleMatrix C(A.copy());
344  C -= B;
345  return(C);
346  }
347 
348 
349  // ----------------------------------------------------------
350 
351  void operator*=(RealMatrix& A, const float d) {
352  for (uint i=0; i<A.rows() * A.cols(); ++i)
353  A[i] *= d;
354  }
355 
356  RealMatrix operator*(const RealMatrix& A, const float d) {
357  RealMatrix B(A.copy());
358  B *= d;
359  return(B);
360  }
361 
362  void operator*=(RealMatrix& A, const RealMatrix& B) {
363  RealMatrix C = MMMultiply(A, B);
364  A=C;
365  }
366 
367  RealMatrix operator*(const RealMatrix& A, const RealMatrix& B) {
368  RealMatrix C = MMMultiply(A, B);
369  return(C);
370  }
371 
372 
373 
374  void operator*=(DoubleMatrix& A, const double d) {
375  for (uint i=0; i<A.rows() * A.cols(); ++i)
376  A[i] *= d;
377  }
378 
379  DoubleMatrix operator*(const DoubleMatrix& A, const double d) {
380  DoubleMatrix B(A.copy());
381  B *= d;
382  return(B);
383  }
384 
385  void operator*=(DoubleMatrix& A, const DoubleMatrix& B) {
386  DoubleMatrix C = MMMultiply(A, B);
387  A=C;
388  }
389 
390  DoubleMatrix operator*(const DoubleMatrix& A, const DoubleMatrix& B) {
391  DoubleMatrix C = MMMultiply(A, B);
392  return(C);
393  }
394 
395 
396  // ----------------------------------------------------------
397 
398  RealMatrix operator-(RealMatrix& A) {
399  RealMatrix B(A.copy());
400  for (uint i=0; i<B.rows() * B.cols(); ++i)
401  B[i] = -B[i];
402  return(B);
403  }
404 
405 
406  DoubleMatrix operator-(DoubleMatrix& A) {
407  DoubleMatrix B(A.copy());
408  for (uint i=0; i<B.rows() * B.cols(); ++i)
409  B[i] = -B[i];
410  return(B);
411  }
412 
413 
414 
415 
416 
417  // Matrix holds column vectors. Each vector is normalized...
419  for (uint i=0; i<A.cols(); ++i) {
420  double sum = 0.0;
421  for (uint j=0; j<A.rows(); ++j)
422  sum += A(j, i) * A(j, i);
423 
424  if (sum <= 0) {
425  for (uint j=0; j<A.rows(); ++j)
426  A(j, i) = 0.0;
427  } else {
428  sum = sqrt(sum);
429  for (uint j=0; j<A.rows(); ++j)
430  A(j, i) /= sum;
431  }
432  }
433  }
434 
435 
436  // ***DEPRECATED***
437  DoubleMatrix deye(const uint n) {
438  return(eye<DoubleMatrix>(n));
439  }
440 
441 
442  }
443 
444 }
445 
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
boost::tuple< RealMatrix, RealMatrix, RealMatrix > svd(RealMatrix &M)
Compute the SVD of a single precision matrix.
Definition: MatrixOps.cpp:94
DoubleMatrix eigenDecomp(DoubleMatrix &M)
Compute eigendecomposition of M.
Definition: MatrixOps.cpp:36
RealMatrix invert(RealMatrix &A, const float eps)
Pseudo-inverse of a matrix using the SVD.
Definition: MatrixOps.cpp:230
void operator+=(RealMatrix &A, const RealMatrix &B)
Overloaded operators for RealMatrix matrices (see important note below)
Definition: MatrixOps.cpp:286
Coord< T > operator*(const Matrix44< T > &, const Coord< T > &)
Definition: Matrix44.hpp:244
Exception caused by a blas/atlas error.
Definition: exceptions.hpp:91
Namespace for most things not already encapsulated within a class.
RealMatrix MMMultiply(const RealMatrix &A, const RealMatrix &B, const bool transa, const bool transb)
Matrix-matrix multiply (using BLAS)
Definition: MatrixOps.cpp:170
Matrix< T, OrderPolicy, StoragePolicy > copy(void) const
Deep copy...
Definition: MatrixImpl.hpp:170
void normalizeColumns(DoubleMatrix &A)
Normalizes each column as a column-vector.
Definition: MatrixOps.cpp:418