Imagine++
Matrix.h
1  // ===========================================================================
2 // Imagine++ Libraries
3 // Copyright (C) Imagine
4 // For detailed information: http://imagine.enpc.fr/software
5 // ===========================================================================
6 
7 
8 namespace Imagine {
11 
12  template<typename T> class SymMatrix;
13 
22  template <typename T>
23  class Matrix : public MultiArray<T,2> {
24  typedef MultiArray<T,2> Base;
25  public:
27  typedef typename Base::iterator iterator;
30 
36  Matrix() : Base() {}
44  Matrix(int M,int N) : Base(M,N) {}
57  Matrix(T* t,int M,int N,bool handleDelete=false) : Base(t,M,N,handleDelete) {}
64  Matrix(const Base& A) : Base(A) {}
71  Matrix(const SymMatrix<T>& A);
74  ~Matrix() {}
82  Matrix& operator=(const Matrix& A) { Base::operator=(A); return *this; }
89  Matrix clone() const
90  {
91  Matrix A(nrow(),ncol());
92  equalize((int)this->totalSize(),this->data(),1,A.data(),1);
93  return A;
94  }
101  int nrow() const { return this->width(); }
108  int ncol() const { return this->height(); }
116  Matrix& fill(T x) { Base::fill(x); return *this; }
124  static Matrix Zero(int M,int N){
125  Matrix Z(M,N);
126  memset(Z.data(),0,Z.totalSize()*sizeof(T));
127  return Z;
128  }
136  static Matrix Identity(int N){
137  Matrix I=Matrix::Zero(N,N);
138  for (int i=0;i<N;i++)
139  I(i,i)=T(1);
140  return I;
141  }
149  Vector<T> getCol(int j) const {
150  assert(j>=0 && j<ncol());
151  Vector<T> v(nrow());
152  equalize(nrow(),this->data()+nrow()*j,1,v.data(),1);
153  return v;
154  }
162  Vector<T> getRow(int i) const {
163  assert(i>=0 && i<nrow());
164  Vector<T> v(ncol());
165  equalize(ncol(),this->data()+i,nrow(),v.data(),1);
166  return v;
167  }
175  int n=std::min(nrow(),ncol());
176  Vector<T> v(n);
177  for (int i=0;i<n;i++)
178  v[i]=(*this)(i,i);
179  return v;
180  }
189  Matrix& setCol(int j, const Vector<T>& v) {
190  assert(v.size()==(size_t)nrow() && j>=0 && j<ncol());
191  equalize(nrow(),v.data(),1,this->data()+nrow()*j,1);
192  return *this;
193  }
202  Matrix& setRow(int i, const Vector<T>& v) {
203  assert(v.size()==(size_t)ncol() && i>=0 && i<nrow());
204  equalize(ncol(),v.data(),1,this->data()+i,nrow());
205  return *this;
206  }
215  int n=int(v.size());
216  assert(n==std::min(nrow(),ncol()));
217  for (int i=0;i<n;i++)
218  (*this)(i,i)=v[i];
219  return *this;
220  }
230  assert (j>=0 && j<ncol());
231  return Vector<T>(this->data()+(size_t)nrow()*j,nrow());
232  }
243  Vector<T> getSubColRef(int j, int i0, int m) {
244  assert (j>=0 && j<ncol() && i0>=0 && m>=0 && i0+m<=nrow());
245  return Vector<T>(this->data()+i0+(size_t)+nrow()*j,m);
246  }
256  Matrix getSubMat(int i0, int m, int j0, int n) const {
257  assert (i0>=0 && i0+m<=nrow() && j0>=0 && j0+n<=ncol());
258  Matrix retMat(m,n);
259  for(int j=0;j<n;j++)
260  equalize(m,this->data()+(j+j0)*nrow()+i0, 1, retMat.data()+j*m, 1);
261  return retMat;
262  }
272  Matrix getColsRef(int j0,int n) {
273  assert (j0>=0 && n>=0 && j0+n<=ncol());
274  return Matrix(this->data()+(size_t)nrow()*j0, nrow(), n);
275  }
283  Vector<T> operator *(const Vector<T>& v) const {
284  assert((size_t)ncol()==v.size());
285  Vector<T> y(nrow());
286  y.fill(T(0));
287  affinity(nrow(), ncol(), T(1), this->data(), v.data(), T(0), y.data(), '0');
288  return y;
289  }
298  friend Vector<T> tmult(const Matrix& A,const Vector<T>& v) {
299  assert((size_t)A.nrow()==v.size());
300  Vector<T> y(A.ncol());
301  y.fill(T(0));
302  affinity(A.nrow(), A.ncol(), T(1), A.data(), v.data(), T(0), y.data(), 'T');
303  return y;
304  }
312  friend Matrix transpose(const Matrix& A) {
313  Matrix r(A.ncol(),A.nrow());
314  for(int i=0;i<A.nrow();i++) for(int j=0;j<A.ncol();j++) r(j,i)=A(i,j);
315  return r;
316  }
324  friend Matrix inverse(const Matrix& A) {
325  assert(A.nrow()==A.ncol());
326  Matrix invA=A.clone();
327  int *pivots=new int[A.nrow()];
328  int info,n=A.ncol();
329  // LU Inverse
330  matrixInverse(&n,invA.data(),&n,pivots,&info);
331  if (info) {
332  std::cerr<<"Cannot invert matrix"<<std::endl;
333  invA.fill(T(0));
334  delete [] pivots;
335  return invA;
336  }
337  delete[] pivots;
338  return invA;
339  }
347  Matrix operator *(const Matrix& B) const {
348  assert(ncol()==B.nrow());
349  int p=ncol();
350  Matrix C(nrow(),B.ncol());
351  C.fill(T(0));
352  matrixProduct(C.nrow(),C.ncol(),p,T(1),this->data(),nrow(),B.data(),B.nrow(),T(0),C.data(), 'N', 'N');
353  return C;
354  }
363  friend Matrix tmult(const Matrix& A,const Matrix& B) {
364  assert(A.nrow()==B.nrow());
365  int p=A.nrow();
366  Matrix C(A.ncol(),B.ncol());
367  C.fill(T(0));
368  matrixProduct(C.nrow(),C.ncol(),p,T(1),A.data(),A.nrow(),B.data(),B.nrow(),T(0),C.data(),'T','N');
369  return C;
370  }
379  friend Matrix multt(const Matrix& A,const Matrix& B) {
380  assert(A.ncol()==B.ncol());
381  int p=A.ncol();
382  Matrix C(A.nrow(),B.nrow());
383  C.fill(T(0));
384  matrixProduct(C.nrow(),C.ncol(),p,T(1),A.data(),A.nrow(),B.data(),B.nrow(),T(0),C.data(),'N','T');
385  return C;
386  }
395  friend Matrix tmultt(const Matrix& A,const Matrix& B) {
396  assert(A.nrow()==B.ncol());
397  int p=A.nrow();
398  Matrix C(A.ncol(),B.nrow());
399  C.fill(T(0));
400  matrixProduct(C.nrow(),C.ncol(),p,T(1),A.data(),A.nrow(),B.data(),B.nrow(),T(0),C.data(),'T','T');
401  return C;
402  }
410  Matrix operator *(const SymMatrix<T>& B) const;
418  Matrix operator +(const Matrix& B) const {
419  Matrix C=clone();
420  return (C+=B);
421  }
429  Matrix operator -(const Matrix& B) const {
430  Matrix C=clone();
431  return (C-=B);
432  }
440  Matrix& operator +=(const Matrix& B) {
441  assert(ncol()==B.ncol());
442  assert(nrow()==B.nrow());
443  combine((int)this->totalSize(), 1, B.data(), this->data());
444  return *this;
445  }
453  Matrix& operator -=(const Matrix& B) {
454  assert(ncol()==B.ncol());
455  assert(nrow()==B.nrow());
456  combine((int)this->totalSize(), -1, B.data(), this->data());
457  return *this;
458  }
465  Matrix operator - () const {
466  return (*this)*T(-1);
467  }
475  Matrix operator *(T x) const {
476  Matrix C=clone();
477  return (C*=x);
478  }
486  Matrix operator /(T x) const {
487  Matrix C=clone();
488  return (C/=x);
489  }
498  multiply((int)this->totalSize(), x, this->data());
499  return *this;
500  }
509  multiply((int)this->totalSize(), 1 / x, this->data());
510  return *this;
511  }
520  friend inline Matrix operator *(T x, const Matrix& B) {return B*x;}
529  friend T dot(const Matrix& a,const Matrix& b) {
530  assert(a.nrow()==b.nrow() && a.ncol()==b.ncol());
531  return scalarProduct((int)a.totalSize(), a.data(), b.data());
532  }
546  friend T norm(const Matrix& A,char type='F') {
547  int nr=A.nrow(),nc=A.ncol();
548  T n = matrixNorm(&type,&nr,&nc,(T*)A.data(),&nr);
549  return n;
550  }
558  friend T rcond(const Matrix& A)
559  {
560  if (A.nrow()<A.ncol())
561  return rcond(transpose(A));
562  // QR decomposition
563  Matrix tmpA=A.clone();
564  int info,m=A.nrow(),n=A.ncol();
565  // Estimation of condition number
566  T rcond;
567  conditionNumber(&n,tmpA.data(),&m,&rcond,&info);
568  if (info!=0) {
569  std::cerr << "Bad cond number " << std::endl;
570  rcond = T(-1);
571  }
572  return rcond;
573  }
589  friend void svd(const Matrix& A,Matrix &U,Vector<T> &S, Matrix &Vt,bool all=false) {
590  Matrix cpy=A.clone();
591  U=Matrix(A.nrow(),A.nrow());
592  Vt=Matrix(A.ncol(),A.ncol());
593  S=Vector<T>(std::min(A.nrow(),A.ncol()));
594  int nr=A.nrow(),nc=A.ncol();
595  int info;
596  singularValuesDecomposition(all?"A":"S",&nr,&nc,cpy.data(),S.data(),U.data(),Vt.data(),&info);
597  if (info) {
598  std::cerr<<"Cannot compute SVD"<<std::endl;
599  return;
600  }
601  }
615  friend Vector<T> linSolve(const Matrix& A,const Vector<T>& b)
616  {
617  assert(b.size() == (size_t)A.nrow());
618  if(A.nrow()==A.ncol())
619  {
620  // Linear system is not over- or under- determined
621  Matrix tmpA=A.clone();
622  Vector<T> tmpb=b.clone();
623  // LU factorization
624  int *pivots=new int[A.nrow()];
625  int info,n=A.nrow();
626  LUSystemSolve(&n,tmpA.data(),pivots,tmpb.data(),&info);
627  if (info) {
628  std::cerr<<"Cannot solve linear system"<<std::endl;
629  tmpb.fill(T(0));
630  delete [] pivots;
631  return tmpb;
632  tmpb.fill(T(0));
633  std::cerr<<"Cannot solve linear system"<<std::endl;
634  }
635  delete[] pivots;
636  return tmpb;
637  } else {
638  // Least-squares solution in the over-determined case, minimum norm solution in the under-determined one
639  int info,m=A.nrow(),n=A.ncol(),o=std::max(m,n);
640  Matrix tmpA=A.clone();
641  Vector<T> tmpb(o);
642  for (int i=0;i<m;i++) tmpb[i] = b[i];
643  QRSystemSolve(&m,&n,tmpA.data(),tmpb.data(),&info);
644  if (info) {
645  tmpb.fill(T(0));
646  std::cerr<<"Cannot solve linear system"<<std::endl;
647  }
648  return tmpb.getSubVect(0,A.ncol());
649  }
650  }
660  friend Matrix pseudoInverse(const Matrix& A,T tolrel=0) {
661  int m=A.nrow(),n=A.ncol();
662  if (n > m)
663  return transpose(pseudoInverse(transpose(A)));
664  Matrix U,Vt;
665  Vector<T> S;
666  svd(A,U,S,Vt);
667  if (tolrel==0)
668  tolrel=std::numeric_limits<T>::epsilon();
669  T tol = S[0] * tolrel;
670  for (int i=0;i<n;i++) {
671  T s=(S[i]>tol)?T(1)/S[i]:T(0);
672  for (int j=0;j<n;j++)
673  Vt(i,j)*=s;
674  }
675  return transpose(U.getColsRef(0,n)*Vt);
676  }
688  friend bool QR(const Matrix& A,Matrix& Q, Matrix& R, bool all=false){
689  assert(A.nrow()>=A.ncol());
690  int info,m=A.nrow(),n=A.ncol(),qn=(all?m:n);
691  if(all) {
692  // Copy A in beginning of Q
693  Q = Matrix(A.nrow(),A.nrow());
694  equalize((int)A.totalSize(),A.data(),1,Q.data(),1);
695  R=Matrix(A.nrow(),A.ncol());
696  QRFactorization(&m,&qn,&n,Q.data(),R.data(),&info);
697  } else
698  {
699  Q = A.clone();
700  R=Matrix(A.ncol(),A.ncol());
701  QRFactorization(&m,&qn,&n,Q.data(),R.data(),&info);
702  }
703  if(info!=1) {
704  std::cerr << "Failed QR" << std::endl;
705  return false;
706  }
707  return true;
708  }
720  friend Matrix cholesky(const Matrix& A,bool low=true) {
721  assert(A.nrow()==A.ncol());
722  Matrix b=A.clone();
723  int n=A.nrow();
724  int info;
725  choleskyDecomposition(low?"L":"U",&n,b.data(),&info);
726  if (info) {
727  std::cerr << "Bad cholesky" << std::endl;
728  b.fill(T(0));
729  return b;
730  }
731  if (low) {
732  for (int j=1; j<n; j++)
733  for (int i=0; i<j; i++)
734  b(i,j)=0;
735  } else {
736  for (int j=1; j<n; j++)
737  for (int i=0; i<j; i++)
738  b(j,i)=0;
739  }
740  return b;
741  }
742 
743 
752  friend T det(const Matrix& A) {
753  assert(A.nrow()==A.ncol());
754  Matrix b=A.clone();
755  // LU decomposition
756  int info,m=A.nrow(),n=A.ncol();
757  T d = determinant(&m,&n,b.data(),&info);
758  if (info!=0) {
759  std::cerr<<"Cannot compute LU decomposition"<<std::endl;
760  return 0;
761  }
762  return d;
763  }
773  friend inline std::ostream& operator<<(std::ostream& out,const Matrix& A) {
774  out<<A.nrow()<<" "<<A.ncol()<< std::endl;
775  for (int i=0;i<A.nrow();i++) {
776  for (int j=0;j<A.ncol();j++) {
777  out<<A(i,j);
778  if (j<A.ncol()-1) out<<" ";
779  }
780  out<<std::endl;
781  }
782  return out;
783  }
793  friend inline std::istream& operator>>(std::istream& in,Matrix& A) {
794  int m,n;
795  in>>m>>n;
796  A.setSize(m,n);
797  for (int i=0;i<A.nrow();i++)
798  for (int j=0;j<A.ncol();j++)
799  in>>A(i,j);
800  return in;
801  }
802  };
803 
813  template <typename T>
814  inline Matrix<T> Diagonal(const Vector<T>& d){
815  int n=int(d.size());
816  Matrix<T> D=Matrix<T>::Zero(n,n);
817  for (int i=0;i<n;i++)
818  D(i,i)=d[i];
819  return D;
820  }
821 
829  template <typename T>
830  Matrix<T> outerProduct(const Vector<T>& v1,const Vector<T>& v2) {
831  assert(v1.size()==v2.size());
832  int n=(int)v1.size();
833  Matrix<T> A(n,n);
834  A.fill(T(0));
835  matrixProduct(n,n,1,T(1),v1.data(),n,v2.data(),1,T(0),A.data(),'N','N');
836  return A;
837  }
838 
852  template <typename T,int N>
854  if (N<=3)
855  return inverse(A); // Direct inversion
856  FMatrix<T,N,N> invA(A);
857  // LU
858  int *pivots=new int[N];
859  int info,n=N;
860  // Inverse
861  matrixInverse(&n,invA.data(),&n,pivots,&info);
862  if (info) {
863  std::cerr<<"Cannot invert matrix"<<std::endl;
864  invA.fill(T(0));
865  }
866  delete[] pivots;
867  return invA;
868  }
869 
885  template <typename T,int M,int N>
886  void svd(const FMatrix<T,M,N>& A, FMatrix<T,M,M> &U, FVector<T,M> &S, FMatrix<T,N,N> &Vt, bool all=false) {
887  FMatrix<T,M,N> cpy(A);
888  int nr=M,nc=N;
889  int info;
890  singularValuesDecomposition(all?"A":"S",&nr,&nc,cpy.data(),S.data(),U.data(),Vt.data(),&info);
891  if (info) {
892  std::cerr<<"Cannot compute SVD"<<std::endl;
893  return;
894  }
895  }
896 
906  template <typename T>
907  void eigenvalues(const Matrix<T>& A, Vector<T>& wr, Vector<T>& wi)
908  {
909  assert(A.ncol()==A.nrow());
910  Matrix<T> cpy = A.clone();
911  wr = Vector<T>(A.nrow());
912  wi = Vector<T>(A.nrow());
913 
914  int n = A.nrow();
915  int info;
916  eigenValues(&n, cpy.data(), wr.data(), wi.data(), 0, &info);
917  if (info)
918  std::cerr << "Cannot compute eigenvalues" << std::endl;
919  }
920 
934  template <typename T,int M,int N>
936  {
937  if(M==N)
938  {
939  // Linear system is not over- or under- determined
940  FMatrix<T,M,N> tmpA(A);
941  FVector<T,M> tmpb(b);
942  // LU factorization
943  int *pivots=new int[M];
944  int info,n=M;
945  // Solve linear system using factorization
946  LUSystemSolve(&n,tmpA.data(),pivots,tmpb.data(),&info);
947  if (info) {
948  std::cerr<<"Cannot solve linear system"<<std::endl;
949  delete [] pivots;
950  return FVector<T,N>(T(0));
951  tmpb.fill(T(0));
952  std::cerr<<"Cannot solve linear system"<<std::endl;
953  }
954  delete[] pivots;
955  return FVector<T,N>(tmpb.data()); //Safe since M==N but necessary for instantiation if M!=N
956  }
957  // Least-squares solution in the over-determined case, minimum norm solution in the under-determined one
958  int info,m=M,n=N;
959  FMatrix<T,M,N> tmpA(A);
960  if(M > N) { // Over-determined system: least squares solution
961  FVector<T,M> tmpb;
962  for (int i=0;i<m;i++) tmpb[i] = b[i];
963  QRSystemSolve(&m,&n,tmpA.data(),tmpb.data(),&info);
964  if (info) {
965  tmpb.fill(T(0));
966  std::cerr<<"Cannot solve linear system"<<std::endl;
967  }
968  return FVector<T,N>( tmpb.data() );
969  }
970  // Under-determined system (M<N): minimum norm solution
971  FVector<T,N> tmpb;
972  for (int i=0;i<m;i++) tmpb[i] = b[i];
973  QRSystemSolve(&m,&n,tmpA.data(),tmpb.data(),&info);
974  if (info) {
975  tmpb.fill(T(0));
976  std::cerr<<"Cannot solve linear system"<<std::endl;
977  }
978  return tmpb;
979  }
994  template <typename T,int M,int N>
996  if (N > M)
997  return transpose(pseudoInverse(transpose(A)));
999  FVector<T,M> S;
1000  svd(A,U,S,Vt);
1001  if (tolrel==0)
1002  tolrel=std::numeric_limits<T>::epsilon();
1003  T tol = S[0] * tolrel;
1004  for (int i=0;i<N;i++) {
1005  T s=(S[i]>tol)?T(1)/S[i]:T(0);
1006  for (int j=0;j<N;j++)
1007  Vt(i,j)*=s;
1008  }
1009  return transpose(FMatrix<T,M,N>(U.data())*Vt);
1010  }
1021  template <typename T,int M,int N>
1023  assert(M>=N);
1024  int info,m=M,n=N;
1025  Q = A;
1026  QRFactorization(&m,&n,&n,Q.data(),R.data(),&info);
1027  if(info!=1) {
1028  std::cerr << "Failed QR" << std::endl;
1029  return false;
1030  }
1031  return true;
1032  }
1043  template <typename T,int M,int N>
1045  assert(M>=N);
1046  int info,m=M,n=N;
1047  equalize(M*N,A.data(),1,Q.data(),1);
1048  QRFactorization(&m,&m,&n,Q.data(),R.data(),&info);
1049  if(info!=1) {
1050  std::cerr << "Failed QR" << std::endl;
1051  return false;
1052  }
1053  return true;
1054  }
1066  template <typename T,int N>
1067  FMatrix<T,N,N> cholesky(const FMatrix<T,N,N>& A,bool low=true) {
1068  FMatrix<T,N,N> b(A);
1069  int n=N;
1070  int info;
1071  choleskyDecomposition(low?"L":"U",&n,b.data(),&info);
1072  if (info) {
1073  std::cerr << "Bad cholesky" << std::endl;
1074  b.fill(T(0));
1075  return b;
1076  }
1077  if (low) {
1078  for (int j=1; j<n; j++)
1079  for (int i=0; i<j; i++)
1080  b(i,j)=0;
1081  } else {
1082  for (int j=1; j<n; j++)
1083  for (int i=0; i<j; i++)
1084  b(j,i)=0;
1085  }
1086  return b;
1087  }
1088 
1098  template <typename T,int N>
1100  int n=N;
1101  FMatrix<T,N,N> b(A);
1102  int info;
1103  // Try cholesky factorization -> only for positive matrices!
1104  choleskyDecomposition((char*)"L",&n,b.data(),&info);
1105  if (!info) {
1106  T d=T(1);
1107  for (int i=0; i<n; i++) d*=b(i,i);
1108  return d*d;
1109  }
1110  // Try SVD
1111  FMatrix<T,N,N> U,V;
1112  FVector<T,N> s;
1113  svd(A,U,s,V);
1114  T d=T(1);
1115  for (int i=0; i<N; i++) d*=s[i];
1116  return d;
1117  }
1118 
1120 }
FMatrix< T, N, N > inverseFMatrix(const FMatrix< T, N, N > &A)
Inverse.
Definition: Matrix.h:853
Vector & fill(T x)
Filling.
Definition: Vector.h:93
friend Matrix multt(const Matrix &A, const Matrix &B)
Product.
Definition: Matrix.h:379
Vector of fixed size.
Definition: FVector.h:17
Vector clone() const
Cloning.
Definition: Vector.h:80
Vector of variable size.
Definition: Vector.h:24
friend Matrix tmult(const Matrix &A, const Matrix &B)
Product.
Definition: Matrix.h:363
Matrix & setCol(int j, const Vector< T > &v)
Set column.
Definition: Matrix.h:189
friend void svd(const Matrix &A, Matrix &U, Vector< T > &S, Matrix &Vt, bool all=false)
SVD.
Definition: Matrix.h:589
FVector & fill(const T &v)
Filling.
Definition: FVector.h:75
friend T rcond(const Matrix &A)
Reciprocal condition number.
Definition: Matrix.h:558
friend Matrix pseudoInverse(const Matrix &A, T tolrel=0)
Pseudo inverse.
Definition: Matrix.h:660
void setSize(const Coords< dim > &sz)
Change sizes.
Definition: MultiArray.h:146
friend Matrix tmultt(const Matrix &A, const Matrix &B)
Product.
Definition: Matrix.h:395
Base::iterator iterator
Iterator type.
Definition: Matrix.h:27
Matrix< T > outerProduct(const Vector< T > &v1, const Vector< T > &v2)
Outer product matrix.
Definition: Matrix.h:830
Matrix & setDiagonal(const Vector< T > &v)
Set diagonal.
Definition: Matrix.h:214
Vector< T > operator*(const Vector< T > &v) const
Product with vector.
Definition: Matrix.h:283
Matrix operator+(const Matrix &B) const
Addition.
Definition: Matrix.h:418
~Matrix()
Destructor.
Definition: Matrix.h:74
Matrix & operator/=(T x)
Scalar in place division.
Definition: Matrix.h:508
FMatrix< T, M, M > Diagonal(const FVector< T, M > &d)
Diagonal.
Definition: FMatrix.h:743
size_t size() const
Size.
Definition: Array.h:194
MultiArray & operator=(const MultiArray &A)
Assignment.
Definition: MultiArray.h:175
friend Matrix inverse(const Matrix &A)
Inverse.
Definition: Matrix.h:324
friend Matrix cholesky(const Matrix &A, bool low=true)
Cholesky decomposition.
Definition: Matrix.h:720
int ncol() const
Number of columns.
Definition: Matrix.h:108
Vector< T > getSubColRef(int j, int i0, int m)
Get sub column by reference.
Definition: Matrix.h:243
Matrix & operator*=(T x)
Scalar in place multiplication.
Definition: Matrix.h:497
friend T det(const Matrix &A)
Determinant.
Definition: Matrix.h:752
T detFMatrix(const FMatrix< T, N, N > &A)
Determinant.
Definition: Matrix.h:1099
FMatrix & fill(const T &v)
Filling.
Definition: FMatrix.h:78
bool QRAll(const FMatrix< T, M, N > &A, FMatrix< T, M, M > &Q, FMatrix< T, M, N > &R)
QR decomposition.
Definition: Matrix.h:1044
int nrow() const
Number of rows.
Definition: Matrix.h:101
Matrix(const Base &A)
Copy constructor.
Definition: Matrix.h:64
Matrix getColsRef(int j0, int n)
Get part of columns by reference.
Definition: Matrix.h:272
size_t totalSize() const
Total Size.
Definition: MultiArray.h:218
Vector< T > getCol(int j) const
Get column.
Definition: Matrix.h:149
Vector< T > getRow(int i) const
Get row.
Definition: Matrix.h:162
Matrix()
Empty constructor.
Definition: Matrix.h:36
T * data()
Data pointer (read/write).
Definition: Array.h:217
Matrix operator-() const
Opposite.
Definition: Matrix.h:465
Matrix & operator=(const Matrix &A)
Assignment.
Definition: Matrix.h:82
Vector getSubVect(size_t offset, size_t size) const
Sub vector.
Definition: Vector.h:114
nD array of variable size.
Definition: MultiArray.h:21
friend std::ostream & operator<<(std::ostream &out, const Matrix &A)
ASCII write.
Definition: Matrix.h:773
friend std::istream & operator>>(std::istream &in, Matrix &A)
ASCII read.
Definition: Matrix.h:793
const T * data() const
Data pointer (read).
Definition: FMatrix.h:182
Matrix & setRow(int i, const Vector< T > &v)
Set row.
Definition: Matrix.h:202
Matrix getSubMat(int i0, int m, int j0, int n) const
Get sub matrix.
Definition: Matrix.h:256
MultiArray & fill(T x)
Filling.
Definition: MultiArray.h:204
friend bool QR(const Matrix &A, Matrix &Q, Matrix &R, bool all=false)
QR decomposition.
Definition: Matrix.h:688
Base::const_iterator const_iterator
Const iterator type.
Definition: Matrix.h:29
int height() const
Size alias 1.
Definition: MultiArray.h:239
friend Vector< T > tmult(const Matrix &A, const Vector< T > &v)
Product with vector.
Definition: Matrix.h:298
T * iterator
Iterator type.
Definition: MultiArray.h:45
friend T dot(const Matrix &a, const Matrix &b)
Scalar product.
Definition: Matrix.h:529
Vector< T > getDiagonal() const
Get diagonal.
Definition: Matrix.h:174
Matrix(T *t, int M, int N, bool handleDelete=false)
Constructor (pre-allocated).
Definition: Matrix.h:57
static Matrix Zero(int M, int N)
Zero matrix.
Definition: Matrix.h:124
void eigenvalues(const Matrix< T > &A, Vector< T > &wr, Vector< T > &wi)
Eigenvalues calculation Given square matrix, only eigenvalues are obtained.
Definition: Matrix.h:907
Matrix & operator+=(const Matrix &B)
In place Addition.
Definition: Matrix.h:440
Matrix(int M, int N)
Constructor (known size).
Definition: Matrix.h:44
friend Matrix transpose(const Matrix &A)
Transpose.
Definition: Matrix.h:312
static Matrix Identity(int N)
Identity.
Definition: Matrix.h:136
friend Vector< T > linSolve(const Matrix &A, const Vector< T > &b)
Linear system.
Definition: Matrix.h:615
Imagine++ namespace.
Definition: Array.h:7
Matrix of variable size.
Definition: Matrix.h:23
Matrix & operator-=(const Matrix &B)
In place Substraction.
Definition: Matrix.h:453
Matrix clone() const
Cloning.
Definition: Matrix.h:89
Matrix & fill(T x)
Filling.
Definition: Matrix.h:116
const T * data() const
Data pointer (read).
Definition: FArray.h:188
Vector< T > getColRef(int j)
Get column by reference.
Definition: Matrix.h:229
int width() const
Size alias 0.
Definition: MultiArray.h:232
Matrix operator/(T x) const
Scalar division.
Definition: Matrix.h:486
const T * const_iterator
Const iterator type.
Definition: MultiArray.h:47
Symmetric Matrix of variable size.
Definition: Matrix.h:12
friend T norm(const Matrix &A, char type='F')
Norm.
Definition: Matrix.h:546
Matrix of fixed dimension.
Definition: FMatrix.h:17