Imagine++
Matrix.h
1 // ===========================================================================
2// Imagine++ Libraries
3// Copyright (C) Imagine
4// For detailed information: http://imagine.enpc.fr/software
5// ===========================================================================
6
7
8namespace 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) {}
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 }
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 }
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 }
441 assert(ncol()==B.ncol());
442 assert(nrow()==B.nrow());
443 combine((int)this->totalSize(), 1, B.data(), this->data());
444 return *this;
445 }
454 assert(ncol()==B.ncol());
455 assert(nrow()==B.nrow());
456 combine((int)this->totalSize(), -1, B.data(), this->data());
457 return *this;
458 }
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)
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());
817 for (int i=0;i<n;i++)
818 D(i,i)=d[i];
819 return D;
820 }
821
829 template <typename T>
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)
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}
size_t size() const
Size.
Definition: Array.h:194
T * data()
Data pointer (read/write).
Definition: Array.h:217
const T * data() const
Data pointer (read).
Definition: FArray.h:188
Matrix of fixed dimension.
Definition: FMatrix.h:17
const T * data() const
Data pointer (read).
Definition: FMatrix.h:182
FMatrix & fill(const T &v)
Filling.
Definition: FMatrix.h:78
Vector of fixed size.
Definition: FVector.h:17
FVector & fill(const T &v)
Filling.
Definition: FVector.h:75
Matrix of variable size.
Definition: Matrix.h:23
Vector< T > getCol(int j) const
Get column.
Definition: Matrix.h:149
Matrix & operator/=(T x)
Scalar in place division.
Definition: Matrix.h:508
Matrix & operator=(const Matrix &A)
Assignment.
Definition: Matrix.h:82
Vector< T > getRow(int i) const
Get row.
Definition: Matrix.h:162
int nrow() const
Number of rows.
Definition: Matrix.h:101
friend Matrix cholesky(const Matrix &A, bool low=true)
Cholesky decomposition.
Definition: Matrix.h:720
friend Matrix pseudoInverse(const Matrix &A, T tolrel=0)
Pseudo inverse.
Definition: Matrix.h:660
friend std::ostream & operator<<(std::ostream &out, const Matrix &A)
ASCII write.
Definition: Matrix.h:773
Matrix getSubMat(int i0, int m, int j0, int n) const
Get sub matrix.
Definition: Matrix.h:256
friend Vector< T > linSolve(const Matrix &A, const Vector< T > &b)
Linear system.
Definition: Matrix.h:615
friend T norm(const Matrix &A, char type='F')
Norm.
Definition: Matrix.h:546
Vector< T > operator*(const Vector< T > &v) const
Product with vector.
Definition: Matrix.h:283
Matrix & operator*=(T x)
Scalar in place multiplication.
Definition: Matrix.h:497
static Matrix Identity(int N)
Identity.
Definition: Matrix.h:136
friend Matrix tmult(const Matrix &A, const Matrix &B)
Product.
Definition: Matrix.h:363
friend Matrix tmultt(const Matrix &A, const Matrix &B)
Product.
Definition: Matrix.h:395
int ncol() const
Number of columns.
Definition: Matrix.h:108
friend void svd(const Matrix &A, Matrix &U, Vector< T > &S, Matrix &Vt, bool all=false)
SVD.
Definition: Matrix.h:589
Vector< T > getColRef(int j)
Get column by reference.
Definition: Matrix.h:229
Base::iterator iterator
Iterator type.
Definition: Matrix.h:27
Matrix clone() const
Cloning.
Definition: Matrix.h:89
Matrix & setDiagonal(const Vector< T > &v)
Set diagonal.
Definition: Matrix.h:214
friend Matrix inverse(const Matrix &A)
Inverse.
Definition: Matrix.h:324
~Matrix()
Destructor.
Definition: Matrix.h:74
friend Matrix multt(const Matrix &A, const Matrix &B)
Product.
Definition: Matrix.h:379
Matrix(T *t, int M, int N, bool handleDelete=false)
Constructor (pre-allocated).
Definition: Matrix.h:57
Matrix(const SymMatrix< T > &A)
Conversion from SymMatrix.
friend T dot(const Matrix &a, const Matrix &b)
Scalar product.
Definition: Matrix.h:529
Matrix & fill(T x)
Filling.
Definition: Matrix.h:116
Vector< T > getDiagonal() const
Get diagonal.
Definition: Matrix.h:174
Matrix operator+(const Matrix &B) const
Addition.
Definition: Matrix.h:418
static Matrix Zero(int M, int N)
Zero matrix.
Definition: Matrix.h:124
friend std::istream & operator>>(std::istream &in, Matrix &A)
ASCII read.
Definition: Matrix.h:793
Matrix(const Base &A)
Copy constructor.
Definition: Matrix.h:64
Matrix & operator-=(const Matrix &B)
In place Substraction.
Definition: Matrix.h:453
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 bool QR(const Matrix &A, Matrix &Q, Matrix &R, bool all=false)
QR decomposition.
Definition: Matrix.h:688
Matrix operator/(T x) const
Scalar division.
Definition: Matrix.h:486
Base::const_iterator const_iterator
Const iterator type.
Definition: Matrix.h:29
friend T det(const Matrix &A)
Determinant.
Definition: Matrix.h:752
Vector< T > getSubColRef(int j, int i0, int m)
Get sub column by reference.
Definition: Matrix.h:243
friend Matrix transpose(const Matrix &A)
Transpose.
Definition: Matrix.h:312
Matrix operator-() const
Opposite.
Definition: Matrix.h:465
Matrix & setCol(int j, const Vector< T > &v)
Set column.
Definition: Matrix.h:189
Matrix()
Empty constructor.
Definition: Matrix.h:36
friend Vector< T > tmult(const Matrix &A, const Vector< T > &v)
Product with vector.
Definition: Matrix.h:298
friend T rcond(const Matrix &A)
Reciprocal condition number.
Definition: Matrix.h:558
Matrix getColsRef(int j0, int n)
Get part of columns by reference.
Definition: Matrix.h:272
Matrix & setRow(int i, const Vector< T > &v)
Set row.
Definition: Matrix.h:202
nD array of variable size.
Definition: MultiArray.h:21
MultiArray & fill(T x)
Filling.
Definition: MultiArray.h:204
void setSize(const Coords< dim > &sz)
Change sizes.
Definition: MultiArray.h:146
size_t totalSize() const
Total Size.
Definition: MultiArray.h:218
const T * const_iterator
Const iterator type.
Definition: MultiArray.h:47
T * iterator
Iterator type.
Definition: MultiArray.h:45
int width() const
Size alias 0.
Definition: MultiArray.h:232
int height() const
Size alias 1.
Definition: MultiArray.h:239
MultiArray & operator=(const MultiArray &A)
Assignment.
Definition: MultiArray.h:175
Symmetric Matrix of variable size.
Definition: SymMatrix.h:22
Vector of variable size.
Definition: Vector.h:24
Vector & fill(T x)
Filling.
Definition: Vector.h:93
Vector clone() const
Cloning.
Definition: Vector.h:80
Vector getSubVect(size_t offset, size_t size) const
Sub vector.
Definition: Vector.h:114
FMatrix< T, M, M > Diagonal(const FVector< T, M > &d)
Diagonal.
Definition: FMatrix.h:743
FMatrix< T, M, M > inverse(const FMatrix< T, M, M > &A)
Inverse.
Definition: FMatrix.h:678
FMatrix< T, N, M > transpose(const FMatrix< T, M, N > &A)
Transpose.
Definition: FMatrix.h:564
T detFMatrix(const FMatrix< T, N, N > &A)
Determinant.
Definition: Matrix.h:1099
FMatrix< T, N, M > pseudoInverse(const FMatrix< T, M, N > &A, T tolrel=0)
Pseudo inverse.
Definition: Matrix.h:995
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
FVector< T, N > linSolve(const FMatrix< T, M, N > &A, const FVector< T, M > &b)
Linear system.
Definition: Matrix.h:935
FMatrix< T, N, N > cholesky(const FMatrix< T, N, N > &A, bool low=true)
Cholesky decomposition.
Definition: Matrix.h:1067
bool QRAll(const FMatrix< T, M, N > &A, FMatrix< T, M, M > &Q, FMatrix< T, M, N > &R)
QR decomposition.
Definition: Matrix.h:1044
bool QR(const FMatrix< T, M, N > &A, FMatrix< T, M, N > &Q, FMatrix< T, N, N > &R)
QR decomposition.
Definition: Matrix.h:1022
FMatrix< T, N, N > inverseFMatrix(const FMatrix< T, N, N > &A)
Inverse.
Definition: Matrix.h:853
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)
SVD.
Definition: Matrix.h:886
Matrix< T > outerProduct(const Vector< T > &v1, const Vector< T > &v2)
Outer product matrix.
Definition: Matrix.h:830
Imagine++ namespace.