Imagine++
FMatrix.h
1// ===========================================================================
2// Imagine++ Libraries
3// Copyright (C) Imagine
4// For detailed information: http://imagine.enpc.fr/software
5// ===========================================================================
6
7namespace Imagine {
10
11
17 template <typename T,int M,int N> class FMatrix {
18 protected:
20 T _data[M*N];
21 public:
34 explicit FMatrix(const T& v) {
35 fill(v);
36 }
44 FMatrix(const T t[M*N]){
45 for (int k=0;k<M*N;k++)
46 _data[k] = t[k];
47 }
55 FMatrix(const T t[M][N]){
56 for (int j=0;j<N;j++) for (int i=0;i<M;i++)
57 (*this)(i,j) = t[i][j];
58 }
67 template <typename T2> FMatrix(const FMatrix<T2,M,N>& A) {
68 for (int k=0;k<M*N;k++)
69 _data[k] = T(A[k]);
70 }
78 FMatrix& fill(const T& v) {
79 for (int k=0;k<M*N;k++)
80 _data[k] = v;
81 return *this;
82 }
89 static FMatrix Zero() {
90 FMatrix Z;
91 std::memset(Z._data,0,M*N*sizeof(T));
92 return Z;
93 }
103 template <typename T2>
105 for (int k=0;k<M*N;k++)
106 _data[k] = T(B[k]);
107 return *this;
108 }
115 int size() const { return M*N; }
122 int nrow() const { return M; }
129 int ncol() const { return N; }
138 const T& operator () (int i, int j) const {
139 assert(i>=0 && i<M && j>=0 && j<N);
140 return (*this)[i+M*j];
141 }
150 T& operator () (int i, int j) {
151 assert(i>=0 && i<M && j>=0 && j<N);
152 return (*this)[i+M*j];
153 }
161 const T& operator [] (int k) const {
162 assert(k>=0 && k<M*N);
163 return _data[k];
164 }
172 T& operator [] (int k) {
173 assert(k>=0 && k<M*N);
174 return _data[k];
175 }
182 const T* data() const { return _data; }
189 T* data() { return _data; }
197 static FMatrix Identity() {
198 assert(M==N);
199 FMatrix I(T(0));
200 for (int i=0;i<M;i++)
201 I(i,i) = T(1);
202 return I;
203 }
212 static FMatrix CrossProd(const FVector<T,3>& v) {
213 assert(M==N && M==3);
214 FMatrix P(T(0));
215 P(1,0) = v.z();
216 P(0,1) = - v.z();
217 P(2,0) = - v.y();
218 P(0,2) = v.y();
219 P(2,1) = v.x();
220 P(1,2) = - v.x();
221 return P;
222 }
230 bool operator == (const FMatrix& B) const {
231 for (int k=0;k<M*N;k++)
232 if ((*this)[k] != B[k]) return false;
233 return true;
234 }
242 bool operator != (const FMatrix& B) const {
243 return !(*this == B);
244 }
252 FMatrix operator + (const FMatrix& B) const {
253 FMatrix C;
254 for (int k=0;k<M*N;k++) C[k] = _data[k] + B[k];
255 return C;
256 }
264 FMatrix operator - (const FMatrix& B) const {
265 FMatrix C;
266 for (int k=0;k<M*N;k++) C[k] = _data[k] - B[k];
267 return C;
268 }
277 for (int k=0;k<M*N;k++) _data[k] += B[k];
278 return *this;
279 }
288 for (int k=0;k<M*N;k++) _data[k] -= B[k];
289 return *this;
290 }
298 FMatrix operator + (T s) const {
299 FMatrix C;
300 for (int k=0;k<M*N;k++) C[k] = _data[k] + s;
301 return C;
302 }
310 FMatrix operator - (T s) const {
311 FMatrix C;
312 for (int k=0;k<M*N;k++) C[k] = _data[k] - s;
313 return C;
314 }
323 for (int k=0;k<M*N;k++) _data[k] += s;
324 return *this;
325 }
334 for (int k=0;k<M*N;k++) _data[k] -= s;
335 return *this;
336 }
344 FMatrix C;
345 for (int k=0;k<M*N;k++) C[k] = -_data[k];
346 return C;
347 }
355 FMatrix operator * (T s) const {
356 FMatrix C;
357 for (int k=0;k<M*N;k++) C[k] = _data[k] * s;
358 return C;
359 }
367 FMatrix operator / (T s) const {
368 FMatrix C;
369 for (int k=0;k<M*N;k++) C[k] = _data[k] / s;
370 return C;
371 }
380 for (int k=0;k<M*N;k++) _data[k] *= s;
381 return *this;
382 }
391 for (int k=0;k<M*N;k++) _data[k] /= s;
392 return *this;
393 }
401 template <int O>
403 FMatrix<T,M,O> C(T(0));
404 for (int j=0;j<O;j++) for (int i=0;i<M;i++) {
405 for (int k=0;k<N;k++) C(i,j) += (*this)(i,k) * B(k,j);
406 }
407 return C;
408 }
417 FVector<T,M> w(T(0));
418 for (int j=0;j<N;j++)
419 for (int i=0;i<M;i++) w[i] += (*this)(i,j) * v[j];
420 return w;
421 }
429 FVector<T,M> getCol(int j) const {
430 assert(j>=0 && j<N);
431 FVector<T,M> v;
432 for (int i=0;i<M;i++) v[i] = (*this)(i,j);
433 return v;
434 }
443 FMatrix& setCol(int j, const FVector<T,M>& v) {
444 assert(j>=0 && j<N);
445 for (int i=0;i<M;i++) (*this)(i,j) = v[i];
446 return *this;
447 }
455 FVector<T,N> getRow(int i) const {
456 assert(i>=0 && i<M);
457 FVector<T,N> v;
458 for (int j=0;j<N;j++) v[j] = (*this)(i,j);
459 return v;
460 }
469 FMatrix& setRow(int i, const FVector<T,N>& v) {
470 assert(i>=0 && i<M);
471 for (int j=0;j<N;j++) (*this)(i,j) = v[j];
472 return *this;
473 }
481 friend T norm2(const FMatrix& A) {
482 T n = 0;
483 for (int k=0;k<M*N;k++) n += A._data[k] * A._data[k];
484 return n;
485 }
493 friend inline T norm(const FMatrix& A) {
494 assert( !std::numeric_limits<T>::is_integer );
495 return ::sqrt(norm2(A));
496 }
505 friend void write(std::ostream& out, const FMatrix& A) {
506 out.write((const char*)A.data(),std::streamsize(M*N*sizeof(T)));
507 }
516 friend void read(std::istream& in, FMatrix& A) {
517 in.read((char*)A.data(),std::streamsize(M*N*sizeof(T)));
518 }
528 friend std::ostream& operator<<(std::ostream& out, const FMatrix& A) {
529 for (int i=0;i<A.nrow();i++) {
530 for (int j=0;j<A.ncol();j++) {
531 out<<A(i,j);
532 if (j<A.ncol()-1) out<<" ";
533 }
534 out<<std::endl;
535 }
536 return out;
537 }
547 friend std::istream& operator>>(std::istream& in, FMatrix& A) {
548 for (int i=0;i<M;i++)
549 for (int j=0;j<N;j++)
550 in>>A(i,j);
551 return in;
552 }
553
554 };
555
563 template <typename T,int M,int N>
566 for (int j=0;j<N;j++) for (int i=0;i<M;i++) C(j,i) = A(i,j);
567 return C;
568 }
569
578 template <typename T,int M,int N>
580 return A*s;
581 }
590 template <typename T,int M,int N>
592 return A+s;
593 }
602 template <typename T,int M,int N>
604 return (-A)+s;
605 }
614 template <typename T,int M,int N,int O>
616 FMatrix<T,N,O> C(T(0));
617 for (int j=0;j<O;j++) for (int i=0;i<N;i++)
618 for (int k=0;k<M;k++) C(i,j) += A(k,i) * B(k,j);
619 return C;
620 }
629 template <typename T,int M,int N,int O>
631 FMatrix<T,M,O> C(T(0));
632 for (int j=0;j<O;j++) for (int i=0;i<M;i++)
633 for (int k=0;k<N;k++) C(i,j) += A(i,k) * B(j,k);
634 return C;
635 }
644 template <typename T,int M,int N,int O>
646 FMatrix<T,N,O> C(T(0));
647 for (int j=0;j<O;j++) for (int i=0;i<N;i++)
648 for (int k=0;k<M;k++) C(i,j) += A(k,i) * B(j,k);
649 return C;
650 }
659 template <typename T,int M,int N>
661 FVector<T,N> w(T(0));
662 for (int j=0;j<N;j++)
663 for (int i=0;i<M;i++) w[j] += A(i,j) * v[i];
664 return w;
665 }
677 template <typename T,int M>
679 FMatrix<T,M,M> I(0.);
680 if (M<=3){
681 const T _det=det(A);
682 if (_det==0) {
683 std::cerr << "Non invertible matrix" << std::endl;
684 } else {
685 if (M==1) {
686 I(0,0)=1/_det;
687 } else if (M==2) {
688 I(0,0)=A(1,1)/_det;
689 I(1,1)=A(0,0)/_det;
690 I(0,1)=-A(0,1)/_det;
691 I(1,0)=-A(1,0)/_det;
692 } else if (M==3) {
693 // could be made faster with direct formula
694 for (int i=0;i<3;i++) for (int j=0;j<3;j++) {
695 const int ia = (i+1)%3;
696 const int ib = (i+2)%3;
697 const int ja = (j+1)%3;
698 const int jb = (j+2)%3;
699 I(j,i) = ( A(ia,ja) * A(ib,jb) - A(ib,ja) * A(ia,jb) ) / _det;
700 }
701 }
702 }
703 } else {
704 std::cerr << "Code for inverting matrices of size " << M << "*" << M << " is not available" << std::endl;
705 std::cerr << "Consider using inverseFMatrix" << std::endl;
706 }
707 return I;
708 }
716 template <typename T,int M>
717 T det(const FMatrix<T,M,M>& A) {
718 if (M==1)
719 return A(0,0);
720 if(M==2){
721 return A(0,0)*A(1,1)-A(0,1)*A(1,0);
722 }
723 if (M==3) {
724 return
725 A(0,0) * (A(1,1) * A(2,2) - A(1,2) * A(2,1))
726 - A(1,0) * (A(0,1) * A(2,2) - A(2,1) * A(0,2))
727 + A(2,0) * (A(0,1) * A(1,2) - A(1,1) * A(0,2));
728 }
729 std::cerr << "Code for determinant of matrices of size " << M << "*" << M << " is not available" << std::endl;
730 std::cerr << "Consider using Matrix class instead of FMatrix" << std::endl;
731 return T(0);
732 }
733
742 template <typename T,int M>
745 for (int i=0;i<M;i++)
746 D(i,i)=d[i];
747 return D;
748 }
749
751}
Matrix of fixed dimension.
Definition: FMatrix.h:17
friend void read(std::istream &in, FMatrix &A)
Binary read.
Definition: FMatrix.h:516
bool operator==(const FMatrix &B) const
Equality test.
Definition: FMatrix.h:230
FMatrix(const T &v)
Constructor with constant value.
Definition: FMatrix.h:34
static FMatrix CrossProd(const FVector< T, 3 > &v)
Cross product matrix.
Definition: FMatrix.h:212
bool operator!=(const FMatrix &B) const
Inequality test.
Definition: FMatrix.h:242
FMatrix & setRow(int i, const FVector< T, N > &v)
Set row.
Definition: FMatrix.h:469
friend T norm2(const FMatrix &A)
Squared Frobenius norm.
Definition: FMatrix.h:481
FMatrix & setCol(int j, const FVector< T, M > &v)
Set column.
Definition: FMatrix.h:443
const T * data() const
Data pointer (read).
Definition: FMatrix.h:182
FMatrix operator+(const FMatrix &B) const
Addition.
Definition: FMatrix.h:252
FMatrix & fill(const T &v)
Filling.
Definition: FMatrix.h:78
int ncol() const
Number of columns.
Definition: FMatrix.h:129
T _data[M *N]
internal storage.
Definition: FMatrix.h:20
FMatrix(const T t[M *N])
Constructor from C array.
Definition: FMatrix.h:44
friend T norm(const FMatrix &A)
Frobenius norm.
Definition: FMatrix.h:493
const T & operator()(int i, int j) const
Read access.
Definition: FMatrix.h:138
FVector< T, M > getCol(int j) const
Get column.
Definition: FMatrix.h:429
friend std::istream & operator>>(std::istream &in, FMatrix &A)
ASCII read.
Definition: FMatrix.h:547
FMatrix(const T t[M][N])
Constructor from bidim C array.
Definition: FMatrix.h:55
static FMatrix Identity()
Identity.
Definition: FMatrix.h:197
FMatrix & operator-=(const FMatrix &B)
In place Substraction.
Definition: FMatrix.h:287
int nrow() const
Number of rows.
Definition: FMatrix.h:122
const T & operator[](int k) const
1D read access.
Definition: FMatrix.h:161
FMatrix()
Empty constructor.
Definition: FMatrix.h:27
FMatrix & operator=(const FMatrix< T2, M, N > &B)
Assignment.
Definition: FMatrix.h:104
FMatrix & operator/=(T s)
Scalar in place division.
Definition: FMatrix.h:390
FMatrix(const FMatrix< T2, M, N > &A)
Copy constructor Constructs from another FMatrix (with a possibly different type)
Definition: FMatrix.h:67
FMatrix operator-() const
Opposite.
Definition: FMatrix.h:343
FVector< T, N > getRow(int i) const
Get row.
Definition: FMatrix.h:455
T * data()
Data pointer (write).
Definition: FMatrix.h:189
int size() const
Size.
Definition: FMatrix.h:115
static FMatrix Zero()
Zero matrix Matrix with constant 0 value.
Definition: FMatrix.h:89
FMatrix & operator*=(T s)
Scalar in place multiplication.
Definition: FMatrix.h:379
FMatrix operator*(T s) const
Scalar multiplication.
Definition: FMatrix.h:355
FMatrix & operator+=(const FMatrix &B)
In place Addition.
Definition: FMatrix.h:276
FMatrix operator/(T s) const
Scalar division.
Definition: FMatrix.h:367
friend void write(std::ostream &out, const FMatrix &A)
Binary write.
Definition: FMatrix.h:505
friend std::ostream & operator<<(std::ostream &out, const FMatrix &A)
ASCII write.
Definition: FMatrix.h:528
const T & z() const
Read alias.
Definition: FVector.h:132
const T & y() const
Read alias.
Definition: FVector.h:118
const T & x() const
Read alias.
Definition: FVector.h:104
FMatrix< T, M, O > multt(const FMatrix< T, M, N > &A, const FMatrix< T, O, N > &B)
Product.
Definition: FMatrix.h:630
FMatrix< T, M, M > Diagonal(const FVector< T, M > &d)
Diagonal.
Definition: FMatrix.h:743
T det(const FMatrix< T, M, M > &A)
Determinant.
Definition: FMatrix.h:717
FMatrix< T, M, N > operator+(T s, const FMatrix< T, M, N > &A)
Scalar addition.
Definition: FMatrix.h:591
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
FMatrix< T, M, N > operator-(T s, const FMatrix< T, M, N > &A)
Scalar substraction.
Definition: FMatrix.h:603
FMatrix< T, N, O > tmult(const FMatrix< T, M, N > &A, const FMatrix< T, M, O > &B)
Product.
Definition: FMatrix.h:615
FMatrix< T, N, O > tmultt(const FMatrix< T, M, N > &A, const FMatrix< T, O, M > &B)
Product.
Definition: FMatrix.h:645
FMatrix< T, M, N > operator*(T s, const FMatrix< T, M, N > &A)
Scalar multiplication.
Definition: FMatrix.h:579
Imagine++ namespace.