Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SVD.cpp
1 /*
2 This file is part of the BIAS library (Basic ImageAlgorithmS).
3 
4 Copyright (C) 2003-2009 (see file CONTACT for details)
5  Multimediale Systeme der Informationsverarbeitung
6  Institut fuer Informatik
7  Christian-Albrechts-Universitaet Kiel
8 
9 
10 BIAS is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 2.1 of the License, or
13 (at your option) any later version.
14 
15 BIAS is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Lesser General Public License for more details.
19 
20 You should have received a copy of the GNU Lesser General Public License
21 along with BIAS; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
24 
25 
26 #include "SVD.hh"
27 
28 #include <Base/Debug/Exception.hh>
29 #include <Base/Common/CompareFloatingPoint.hh>
30 
31 #include <algorithm>
32 #include <cmath>
33 #include <vector>
34 
35 #include "Lapack.hh"
36 #ifdef BIAS_HAVE_OPENCV
37 #include <cv.h>
38 #include <cxcore.h>
39 #endif
40 using namespace BIAS;
41 using namespace std;
42 
43 SVD::SVD(const Matrix<double> &M, double ZeroThreshold,
44  bool AbsoluteZeroThreshold) {
45  Decomposed_=false;
46  AbsoluteZeroThreshold_ = AbsoluteZeroThreshold;
47  Compute( M, ZeroThreshold);
48 }
49 
50 #ifdef BIAS_HAVE_OPENCV
51 int SVD::ComputeOpenCV(const Matrix<double>& M, double ZeroThreshold){
52  ZeroThreshold_ = ZeroThreshold;
53  A_ = M;
54  unsigned cols = M.GetCols();
55  unsigned rows = M.GetRows();
56  CvMat* mat = cvCreateMat(rows,cols,CV_32FC1);
57  cvZero(mat);
58  //copy BIAS to OpenCV matrix
59  for(unsigned i=0;i<rows;i++){
60  for(unsigned j=0;j<cols;j++){
61  CV_MAT_ELEM( *mat, float, i, j ) = (float)M[i][j];
62  }
63  }
64 
65  CvMat* U = cvCreateMat(rows,rows,CV_32FC1); cvZero(U);
66  CvMat* D = cvCreateMat(rows,cols,CV_32FC1); cvZero(D);
67  CvMat* V = cvCreateMat(cols,cols,CV_32FC1); cvZero(V);
68  cvSVD(mat, D, U, V, CV_SVD_V_T); // A = U D V^T
69 
70  U_.newsize(rows,rows);
71  VT_.newsize(cols,cols);
72  S_.newsize(rows);
73  for(unsigned i=0;i<rows;i++){
74  for(unsigned j=0;j<rows;j++){
75  U_[i][j] = (double)CV_MAT_ELEM(*U,float,i,j);
76  }
77  if(i< rows && i < cols){
78  //singular values are on diagonal of matrix D
79  S_[i] = (double)CV_MAT_ELEM(*D,float,i,i);
80  }
81  }
82  for(unsigned i=0;i<cols;i++){
83  for(unsigned j=0;j<cols;j++){
84  VT_[i][j] = (double)CV_MAT_ELEM(*V,float,i,j);
85  }
86  }
87 
88  cvReleaseMat(&U);
89  cvReleaseMat(&D);
90  cvReleaseMat(&V);
91  cvReleaseMat(&mat);
92 
93  //need extra if -> ensure S_.size != 0
94  if (S_[0] < numeric_limits<double>::epsilon())
95  return -1;
96  else
97  Decomposed_=true;
98 
99  return 0;
100 }
101 #endif
102 int SVD::Compute(const Matrix<double> & M, double ZeroThreshold)
103 {
104  Decomposed_=false;
105  ZeroThreshold_ = ZeroThreshold;
106  return General_Eigenproblem_GeneralMatrix_Lapack( M );
107 }
108 
109 
111 {
112  // S_ will contain the singular values of A sorted with s(i) <= s(i+1)
113  // U_ will contain m x m left orthogonal/unitary matrix with
114  // eigenvectors in columns
115  // VT_ will contain n x n right _transposed_ orthogonal/unitary matrix
116  // with eigenvectors in rows (because of transpose)
117  A_=M;
118  int res=0;
119  res = General_singular_value_decomposition(A_, S_, U_, VT_);
120 
121  if ( S_.size() != 0 && U_.num_rows() != 0 && U_.num_cols() != 0 &&
122  VT_.num_rows() != 0 && VT_.num_cols() != 0 && res == 0) {
123  //need extra if -> ensure S_.size != 0
124  if (S_[0] < numeric_limits<double>::epsilon()) {
125  res = -1;
126  } else Decomposed_=true;
127  } else {
128  res = -1;
129  }
130 
131  return res;
132 }
133 
134 // -- Solve the matrix-vector system M x = y, returning x.
136 {
137  /// @bug bug in this function (woelk)
138  BEXCEPTION("Error, there is a bug in this function, use Solve(Matrix, "
139  "vector, Vector instead");
140  Vector<double> x(VT_.num_cols()); // solution matrix vector
141  Vector<double> UTy;
142 
143  //cout << "U_=" << U_ << " VT_= " << VT_ << " S_=" << S_;
144 
145  if(Decomposed_) {
146  if (y.size() <= U_.num_rows()){
147  Matrix<double> UT=U_.Transpose();
148  // at first multiply with U^T
149  if (U_.num_rows() < U_.num_cols()) {// augment y with extra rows of
150  Vector<double> yy(U_.num_rows(), 0.0); // zeros, so that it match
151  for (int k=0; k< y.size();k++)
152  yy[k] = y[k]; // cols of u.transpose.
153  UTy = UT * yy;
154  } else {
155  UTy = UT * y;
156  }
157 
158  // multiply with diagonal 1/S_
159  int i=0;
160  for (; i < UTy.size(); i++)
161  x[i]=(fabs(S_[i]) >= ZeroThreshold_) ? UTy[i]/S_[i] : 0;//UTy[i]*S_[i];
162 
163  for (;i<x.size();i++) x[i]=0;
164 
165  return VT_.Transpose() * x; // premultiply with v.
166 
167  } else {
168  BIASERR("size of right hand side of Ax = y is too big. Length(y) = "
169  << y.size());
170  x.newsize(0);
171  return x;
172  }
173  } else {
174  BIASERR("SVD holds no decomposed matrix.");
175  x.newsize(0);
176  return x;
177  }
178 }
179 
180 #ifdef BIAS_HAVE_OPENCV
181  /** \brief Calculates new inverse with OpenCV cvInvert*/
183 {
184  unsigned cols = A.GetCols();
185  unsigned rows = A.GetRows();
186  Matrix<double> B(rows,cols,0);
187  CvMat* a = cvCreateMat(rows,cols,CV_32FC1);
188  cvZero(a);
189  //copy BIAS to OpenCV matrix
190  for(unsigned i=0;i<rows;i++){
191  for(unsigned j=0;j<cols;j++){
192  CV_MAT_ELEM( *a, float, i, j ) = (float)A[i][j];
193  }
194  }
195  CvMat* b = cvCreateMat(rows,cols,CV_32FC1);
196  cvZero(b);
197  //CV_LU - Gaussian elimination with optimal pivot element
198  //CV_SVD - Singular decomposition method
199  double ret = cvInvert(a,b,CV_SVD);
200  if(ret==0.0) BIASWARN("Could not invert Matrix");
201  //copy OpenCV to BIAS matrix
202  for(unsigned i=0;i<rows;i++){
203  for(unsigned j=0;j<cols;j++){
204  B[i][j] = (double)CV_MAT_ELEM( *a, float, i, j );
205  }
206  }
207  cvReleaseMat(&a);
208  cvReleaseMat(&b);
209 
210  return B;
211 }
212 #endif
213 
215 {
216  if (!Decomposed_) {//no decomposition, return 0x0-matrix
217  Matrix<double> res(0,0);
218  return res;
219  }
220 
221  Matrix<double> V(VT_.num_cols(), VT_.num_rows());
222  Matrix<double> UT(U_.num_cols(), U_.num_rows());
223  Matrix<double> SdT(VT_.num_rows(), U_.num_cols());
224  SdT.SetZero();
225 
226  if (AbsoluteZeroThreshold_) {
227  for (register int i = 0; i < S_.size(); i++){
228  SdT[i][i] = (S_[i] < ZeroThreshold_ && S_[i] > -ZeroThreshold_)
229  ? (0.0) : (1.0/S_[i]);
230  }
231  } else {
232  for (register int i = 0; i < S_.size(); i++){
233  SdT[i][i]=(S_[i]/S_[0] < ZeroThreshold_ && S_[i]/S_[0] > -ZeroThreshold_)
234  ? (0.0) : (1.0/S_[i]);
235  }
236  }
237  V = VT_.Transpose();
238  UT = U_.Transpose();
239 
240  return V * SdT * UT;
241 }
242 
244 {
245  if (!Decomposed_) {//no decomposition, return 0x0-matrix
246  Matrix<double> res(0,0);
247  return res;
248  }
249 
250  Matrix<double> V(VT_.num_cols(), VT_.num_rows());
251  Matrix<double> UT(U_.num_cols(), U_.num_rows());
252  Matrix<double> SdT(VT_.num_rows(), U_.num_cols());
253  SdT.SetZero();
254 
255  for (register int i = 0; i < rank; i++){
256  SdT[i][i] = (1.0/S_[i]);
257  }
258  for (register int i = rank; i < S_.size(); i++){
259  SdT[i][i] = 0;
260  }
261  V = VT_.Transpose();
262  UT = U_.Transpose();
263 
264  return V * SdT * UT;
265 }
266 
267 
269 {
270  if (Compute(A)==0)
271  return Invert();
272  else return Matrix<double>(0,0);
273 }
274 
275 
277 {
278 #ifdef BIAS_DEBUG
279  // check symmetry and symmetric positive definit'ness
280  // 1. symmetry
281  bool valid = (A_.GetCols() == A_.GetRows());
282  if (valid){
283  const int num = A_.GetCols();
284  for (int r=0; r<num; r++){
285  for (int c=r; c<num; c++){
286  // dont use double epsilon but 1e-10 for numeric comparison
287  valid = valid && Equal(A_[r][c], A_[c][r], 1e-10);
288  if (!valid) {
289  cout << setprecision(20) << A_[r][c] << "\t"<< setprecision(20) << A_[c][r]<<endl;
290  break;
291  break;
292  }
293  }
294  }
295  }
296  if (!valid){
297  //BIASERR("matrix is not symmetric "<<A_);
298  BEXCEPTION("matrix is not symmetric "/*<<A_*/);
299  }
300  // 2. positive definit = all eigenvalues are positive
302  for (register int i=0; i<eig.size(); i++){
303  valid = valid && (GreaterEqual(eig[i], 0.0, 1e-10));
304  if (!valid){
305  cout << "eigenvalue "<<setprecision(20)<<eig[i]<<endl;
306  break;
307  }
308  }
309  if (!valid){
310  //BIASERR("matrix is not positive definit "<<A_);
311  BEXCEPTION("matrix is not positive definit "/*<<A_*/);
312  }
313 #endif
314  if (!Decomposed_){
315  int res = 0;
316  if ((res = Compute(A_))!=0) {
317  BIASERR("SVD failed with result "<<res);
318  return Matrix<double>(VT_.num_rows(), U_.num_cols(), MatrixZero);
319  }
320  }
321  // [U S V] = svd(A); sqrtm(A) = U * sqrt(S) * V'
322  Matrix<double> SqrtS(VT_.num_rows(), U_.num_cols());
323  SqrtS.SetZero();
324  if (AbsoluteZeroThreshold_) {
325  for (register int i=0; i<S_.size(); i++){
326  SqrtS[i][i] = (fabs(S_[i])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
327  }
328  } else {
329  for (register int i=0; i<S_.size(); i++){
330  SqrtS[i][i] = (fabs(S_[i]/S_[0])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
331  }
332  }
333 
334  return U_ * SqrtS * VT_;
335 }
336 
338 {
339  if (Compute(A)==0)
340  return Sqrt();
341  else return Matrix<double>(0,0);
342 }
343 
345 {
346 #ifdef BIAS_DEBUG
347  // check symmetry and symmetric positive definit'ness
348  // 1. symmetry
349  bool valid = (A_.GetCols() == A_.GetRows());
350  if (valid){
351  const int num = A_.GetCols();
352  for (int r=0; r<num; r++){
353  for (int c=r; c<num; c++){
354  valid = valid && Equal(A_[r][c], A_[c][r], 1e-12);
355  if (!valid) {
356  cout<<" Invalid matrix ("<<num<<"x"<<num<<") element at "
357  <<r<<" "<<c<<endl;
358  cout << setprecision(20) << A_[r][c] << "\t"<< setprecision(20)
359  << A_[c][r]<<endl;
360  break;
361  }
362  }
363  if (!valid) { break; }
364  }
365  }
366  if (!valid){
367  BEXCEPTION("matrix is not symmetric "<<A_<<endl
368  <<"Fix this by calling A_.MakeSymmetric() in CALLING class");
369  }
370 
371  // 2. positive definit = all eigenvalues are positive
373  for (register int i=0; i<eig.size(); i++){
374  valid = valid && (GreaterEqual(eig[i], 0.0, 1e-12));
375  if (!valid){
376  //cout << "eigenvalue "<<setprecision(20)<<eig[i]<<endl;
377  break;
378  }
379  }
380  if (!valid){
381  //BIASERR("matrix is not positive definit "<<A_);
382  //BEXCEPTION("matrix is not positive definit "/*<<A_*/);
383  }
384 #endif
385 
386  if (!Decomposed_){
387  if (Compute(A_)!=0) {
388  // error in decomposition
389  return Matrix<double>(0,0); // empty matrix
390  }
391  }
392  // [U S V] = svd(A); sqrtm(A) = U * sqrt(S)
393  Matrix<double> SqrtS(VT_.num_rows(), U_.num_cols(), MatrixZero);
394  if (AbsoluteZeroThreshold_) {
395  for (register int i=0; i<S_.size(); i++){
396  SqrtS[i][i] = (fabs(S_[i])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
397  }
398  } else {
399  for (register int i=0; i<S_.size(); i++){
400  SqrtS[i][i] = (fabs(S_[i]/S_[0])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
401  }
402  }
403 
404  return U_ * SqrtS;
405 }
406 
408 {
409  if (Compute(A)==0)
410  return SqrtT();
411  else return Matrix<double>(0,0);
412 }
413 
414 
415 const int SVD::AbsNullspaceDim() const
416 {
417  int dim = 0;
418  int index = S_.size()-1;
419 
420  //cerr << "ZeroThreshold_: "<<ZeroThreshold_<<endl;
421  while (index >= 0 && fabs(S_[index]) < ZeroThreshold_ ) {
422  // singular values are sorted in descending order is still (around zero
423  index--;
424  dim++;
425  }
426 
427  if (A_.num_cols()>A_.num_rows()){
428  dim+=A_.num_cols()-A_.num_rows();
429  }
430  return dim;
431 }
432 
433 const int SVD::RelNullspaceDim() const
434 {
435  int dim = 0;
436  int index = S_.size()-1;
437 
438  // use greatest singular value as reference
439  const double RefSingValue(S_[0]);
440  if (RefSingValue<=DBL_EPSILON) {
441  // the greatest singular value is zero => all are zero
442  return 0;
443  }
444 
445  while ((S_[index]/RefSingValue)<ZeroThreshold_){
446  dim++;
447  if (--index<1) break; // S[0]/S[0]==1 always => dont need to check that
448  }
449 
450  // check if we have an "underdetermined equation system" by dimensionality
451  if (A_.num_cols()>A_.num_rows()){
452  dim += A_.num_cols()-A_.num_rows();
453  }
454 
455  return dim;
456 }
457 
458 
459 const int SVD::AbsRightNullspaceDim() const {
460  return AbsNullspaceDim();
461 }
462 
463 const int SVD::RelRightNullspaceDim() const {
464  return RelNullspaceDim();
465 }
466 
467 const int SVD::AbsLeftNullspaceDim() const
468 {
469  int dim = 0;
470  int index = S_.size()-1;
471 
472  //cerr << "ZeroThreshold_: "<<ZeroThreshold_<<endl;
473  while (index >= 0 && fabs(S_[index]) < ZeroThreshold_ ) {
474  // singular values are sorted in descending order
475  index--;
476  dim++;
477  }
478  // compensate null space for rectangular matrices
479  if (A_.num_rows()>A_.num_cols()){
480  dim += A_.num_rows()-A_.num_cols();
481  }
482  return dim;
483 }
484 
485 const int SVD::RelLeftNullspaceDim() const
486 {
487  int dim = 0;
488  int index = S_.size()-1;
489 
490  if (GetSingularValue(0)==0){
491  dim=NullspaceDim();
492  } else {
493  double maxs=S_[0];
494  while (index>=0 && fabs(S_[index]/maxs)<ZeroThreshold_){
495  index--;
496  dim++;
497  }
498  // compensate null space for rectangular matrices
499  if (A_.num_rows()>A_.num_cols()){
500  dim+=A_.num_rows()-A_.num_cols();
501  }
502  }
503  return dim;
504 }
505 
506 unsigned int SVD::Rank()
507 {
508  BIASDOUT(D_SVD_RANK, "S_.Size(): "<<S_.Size()<<" NullspaceDim(): "
509  << NullspaceDim()<<" A_.num_cols(): "<<A_.num_cols()
510  <<" A_.num_rows(): "<<A_.num_rows()<<" S_: "<<S_);
511  unsigned int lowdim=A_.num_rows()<A_.num_cols()?A_.num_rows():A_.num_cols();
512  return (lowdim-NullspaceDim());
513 }
514 
515 
517 {
518  int res = 0;
519  // make squared symmetric matrix A^T * A
520  Matrix<double> ATA(A.num_cols(), A.num_cols());
521  //ATA = A.Transpose()*A;
522  A.GetSystemMatrix(ATA);
523 
524  // compute ATA = U * S * V^T
525  Compute(ATA);
526  Vector<double> d(A.num_rows()), z(A.num_cols());
527  Vector<double> S = GetS();
528  BIASCDOUT(D_SVD_SOLVE, "S: "<<S<<endl);
529 
530  // A * x = b => x = (A^T * A )^-1 * A^T * b
531  // = (U * S * V^T)^-1 * A^T * b
532  // by symmetry we know: A^T * A = (A^T * A)^T
533  // = (V * S * U^T)^-1 * A^T * b
534  // since U and V^T are square orthogonal matrices here
535  // = U^-T * S^-1 * V^T * A^T * b
536  // = U * S^-1 * d
537  d = GetVT() * A.Transpose() * b;
538 
539  // compute z = S^-1 * d
540  BIASCDOUT(D_SVD_SOLVE, "d: "<<d<<endl);
541  if (AbsoluteZeroThreshold_) {
542  for (int i=0; i<A.num_cols(); i++) {
543  if (fabs(S[i])>ZeroThreshold_) {
544  z[i] = d[i]/S[i];
545  } else {
546  z[i] = 0.0;
547  res--;
548  }
549  }
550  } else {
551  for (int i=0; i<A.num_cols(); i++) {
552  if (fabs(S[i]/S_[0])>ZeroThreshold_) {
553  z[i] = d[i]/S[i];
554  } else {
555  z[i] = 0.0;
556  res--;
557  }
558  }
559  }
560 
561  if (res<0) {
562  BIASDOUT(D_SVD_SOLVE,"Error in SVD, degenerate case in Solve ("
563  << -res <<" dim. solution space), S="
564  << S);
565  }
566 
567  // get final solution x = U * z
568  x = GetU() * z;
569  return res;
570 }
571 
572 
574 {
575  std::vector<double> eigValues(VT_.GetRows());
576  // now we need the eigenvalues (with sign !)
577  // compute this by projecting the (normalized) eigenvectors onto their images
578  Vector<double> eigenvector, eigenimage;
579  for (unsigned int i=0; i<VT_.GetRows(); i++) {
580  // get ith eigenvector and its image
581  eigenvector=VT_.GetRow(i);
582  eigenimage = A_ * eigenvector;
583  // compute scalar product
584  eigValues[i] = eigenimage * eigenvector;
585  }
586  std::sort(eigValues.begin(),eigValues.end());
587  BIAS::Vector<double> result(VT_.GetRows(),&eigValues[0]);
588 
589  return result;
590 }
591 
int General_Eigenproblem_GeneralMatrix_Lapack(const Matrix< double > &M)
solve the general (non-special) eigenvalue/eigenvector problem of a general (non-symmetric) matrix M ...
Definition: SVD.cpp:110
const int AbsNullspaceDim() const
returns dim of nullspace using the absolute threshold criterion
Definition: SVD.cpp:415
const int AbsRightNullspaceDim() const
returns the dim of the matrix&#39;s right nullspace, using absolute threshold criterion ...
Definition: SVD.cpp:459
Subscript num_cols() const
Definition: cmat.h:320
Vector< double > Solve(const Vector< double > &y) const
Definition: SVD.cpp:135
Matrix< T > Transpose() const
transpose function, storing data destination matrix
Definition: Matrix.hh:823
unsigned int Rank()
returns the rank of A_
Definition: SVD.cpp:506
Matrix< double > SqrtT()
returns the square root of a symmetric positive definite matrix M A = sqrt(M) = U * sqrt(S)...
Definition: SVD.cpp:344
BIAS::Vector< double > GetEigenValues()
Call this after Compute(), returns the eigenvalues of the matrix in ascending order (smallest first)...
Definition: SVD.cpp:573
int Compute(const Matrix< double > &M, double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD)
set a new matrix and compute its decomposition.
Definition: SVD.cpp:102
const int RelNullspaceDim() const
compare singular values against greatest, not absolute
Definition: SVD.cpp:433
void SetZero()
Sets all values to zero.
Definition: Matrix.hh:856
unsigned int GetRows() const
Definition: Matrix.hh:202
bool GreaterEqual(const T left, const T right, const T eps=std::numeric_limits< T >::epsilon())
comparison function for floating point values
const int RelLeftNullspaceDim() const
Definition: SVD.cpp:485
Matrix< double > InvertOpenCV(const Matrix< double > &A)
Calculates new inverse with OpenCV cvInvert.
Definition: SVD.cpp:182
SVD(const Matrix< double > &M, double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD, bool AbsoluteZeroThreshold=true)
solve the general eigenproblem of the rectangular matrix M or with other words, make the U*S*V^T deco...
Definition: SVD.cpp:43
int General_singular_value_decomposition(char jobu, char jobvt, const BIAS::Matrix< double > &A, BIAS::Vector< double > &ret_S, BIAS::Matrix< double > &ret_U, BIAS::Matrix< double > &ret_VT)
solve general eigenproblem.
Definition: Lapack.cpp:688
bool Equal(const T left, const T right, const T eps)
comparison function for floating point values See http://www.boost.org/libs/test/doc/components/test_...
const int AbsLeftNullspaceDim() const
returns the dim of the matrix&#39;s left nullspace, using absolute threshold criterion ...
Definition: SVD.cpp:467
unsigned int GetCols() const
Definition: Matrix.hh:204
int ComputeOpenCV(const Matrix< double > &M, double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD)
Use OpenCV for decomposition as Lapack frequently leads to crashes under windows. ...
Definition: SVD.cpp:51
Matrix< double > Invert()
returns pseudoinverse of A = U * S * V^T A^+ = V * S^+ * U^T
Definition: SVD.cpp:214
Matrix< double > Sqrt()
returns the square root of a symmetric positive definite matrix M A = sqrt(M) = U * sqrt(S) * V_T...
Definition: SVD.cpp:276
Subscript num_rows() const
Definition: cmat.h:319
BIAS::Vector< double > Upper_symmetric_eigenvalue_solve(const BIAS::Matrix< double > &A)
Solve symmetric eigenvalue problem (eigenvalues only)
Definition: Lapack.cpp:438
void GetSystemMatrix(Matrix< T > &dest) const
compute square system matrix dest = A^T * A
Definition: Matrix.hh:1075
Subscript size() const
Definition: vec.h:262
const int RelRightNullspaceDim() const
Definition: SVD.cpp:463