00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef TNT_MATRIX_H
00025 #define TNT_MATRIX_H
00026
00027 #include "tnt_subscript.h"
00028 #include "tnt_vector.h"
00029 #include <cstdlib>
00030 #include <cassert>
00031 #include <iostream>
00032 #include <sstream>
00033 #include <fstream>
00034
00035 namespace TNT
00036 {
00037
00038
00039
00040
00055 template <class T>
00056 class Matrix
00057 {
00058
00059 private:
00060 Subscript m_;
00061 Subscript n_;
00062 Subscript mn_;
00063 T* v_;
00064 T** row_;
00065 T* vm1_ ;
00066 T** rowm1_;
00067
00068
00069
00070
00071 void initialize(Subscript M, Subscript N)
00072 {
00073 mn_ = M*N;
00074 m_ = M;
00075 n_ = N;
00076
00077 v_ = new T[mn_];
00078 row_ = new T*[M];
00079 rowm1_ = new T*[M];
00080
00081 assert(v_ != NULL);
00082 assert(row_ != NULL);
00083 assert(rowm1_ != NULL);
00084
00085 T* p = v_;
00086 vm1_ = v_ - 1;
00087 for (Subscript i=0; i<M; i++)
00088 {
00089 row_[i] = p;
00090 rowm1_[i] = p-1;
00091 p += N ;
00092
00093 }
00094
00095 rowm1_ -- ;
00096 }
00097
00098 void copy(const T* v)
00099 {
00100 Subscript N = m_ * n_;
00101 Subscript i;
00102
00103 #ifdef TNT_UNROLL_LOOPS
00104 Subscript Nmod4 = N & 3;
00105 Subscript N4 = N - Nmod4;
00106
00107 for (i=0; i<N4; i+=4)
00108 {
00109 v_[i] = v[i];
00110 v_[i+1] = v[i+1];
00111 v_[i+2] = v[i+2];
00112 v_[i+3] = v[i+3];
00113 }
00114
00115 for (i=N4; i< N; i++)
00116 v_[i] = v[i];
00117 #else
00118
00119 for (i=0; i< N; i++)
00120 v_[i] = v[i];
00121 #endif
00122 }
00123
00124 void set(const T& val)
00125 {
00126 Subscript N = m_ * n_;
00127 Subscript i;
00128
00129 #ifdef TNT_UNROLL_LOOPS
00130 Subscript Nmod4 = N & 3;
00131 Subscript N4 = N - Nmod4;
00132
00133 for (i=0; i<N4; i+=4)
00134 {
00135 v_[i] = val;
00136 v_[i+1] = val;
00137 v_[i+2] = val;
00138 v_[i+3] = val;
00139 }
00140
00141 for (i=N4; i< N; i++)
00142 v_[i] = val;
00143 #else
00144
00145 for (i=0; i< N; i++)
00146 v_[i] = val;
00147
00148 #endif
00149 }
00150
00151
00152
00153 void destroy()
00154 {
00155
00156 if (v_ == NULL) return ;
00157
00158
00159 if (v_ != NULL) delete [] (v_);
00160 if (row_ != NULL) delete [] (row_);
00161
00162
00163 rowm1_ ++;
00164 if (rowm1_ != NULL ) delete [] (rowm1_);
00165 }
00166
00167 public:
00168
00169 typedef Subscript size_type;
00170 typedef T value_type;
00171 typedef T element_type;
00172 typedef T* pointer;
00173 typedef T* iterator;
00174 typedef T& reference;
00175 typedef const T* const_iterator;
00176 typedef const T& const_reference;
00177
00178 Subscript lbound() const { return 1;}
00179
00180
00181
00182
00183 operator T**(){ return row_; }
00184 operator const T**() const { return row_; }
00185
00186
00190 Subscript size() const { return mn_; }
00191
00192
00193
00194 Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
00195
00196 Matrix(const Matrix<T> &A)
00197 {
00198 initialize(A.m_, A.n_);
00199 copy(A.v_);
00200 }
00201
00210 Matrix(Subscript M, Subscript N, const T& value = T(0))
00211 {
00212 initialize(M,N);
00213 set(value);
00214 }
00215
00224 Matrix(Subscript M, Subscript N, const T* v)
00225 {
00226 initialize(M,N);
00227 copy(v);
00228 }
00229
00238 Matrix(Subscript M, Subscript N, const char *s)
00239 {
00240 initialize(M,N);
00241
00242 std::istringstream ins(s);
00243
00244 Subscript i, j;
00245
00246 for (i=0; i<M; i++)
00247 for (j=0; j<N; j++)
00248 ins >> row_[i][j];
00249 }
00250
00251
00252
00253 ~Matrix()
00254 {
00255 destroy();
00256 }
00257
00258
00282 Matrix<T>& newsize(Subscript M, Subscript N)
00283 {
00284 if (num_rows() == M && num_cols() == N)
00285 return *this;
00286
00287 destroy();
00288 initialize(M,N);
00289
00290 return *this;
00291 }
00292
00293
00294
00295
00303 Matrix<T>& operator=(const Matrix<T> &B)
00304 {
00305 if (v_ == B.v_)
00306 return *this;
00307
00308 if (m_ == B.m_ && n_ == B.n_)
00309 copy(B.v_);
00310
00311 else
00312 {
00313 destroy();
00314 initialize(B.m_, B.n_);
00315 copy(B.v_);
00316 }
00317
00318 return *this;
00319 }
00320
00321 Matrix<T>& operator=(const T& scalar)
00322 {
00323 set(scalar);
00324 return *this;
00325 }
00326
00327
00328 Subscript dim(Subscript d) const
00329 {
00330 #ifdef TNT_BOUNDS_CHECK
00331 assert( d >= 1);
00332 assert( d <= 2);
00333 #endif
00334 return (d==1) ? m_ : ((d==2) ? n_ : 0);
00335 }
00336
00337 Subscript num_rows() const { return m_; }
00338 Subscript num_cols() const { return n_; }
00339
00340
00341
00342
00343 inline T* operator[](Subscript i)
00344 {
00345 #ifdef TNT_BOUNDS_CHECK
00346 assert(0<=i);
00347 assert(i < m_) ;
00348 #endif
00349 return row_[i];
00350 }
00351
00352 inline const T* operator[](Subscript i) const
00353 {
00354 #ifdef TNT_BOUNDS_CHECK
00355 assert(0<=i);
00356 assert(i < m_) ;
00357 #endif
00358 return row_[i];
00359 }
00360
00361 inline reference operator()(Subscript i)
00362 {
00363 #ifdef TNT_BOUNDS_CHECK
00364 assert(1<=i);
00365 assert(i <= mn_) ;
00366 #endif
00367 return vm1_[i];
00368 }
00369
00370 inline const_reference operator()(Subscript i) const
00371 {
00372 #ifdef TNT_BOUNDS_CHECK
00373 assert(1<=i);
00374 assert(i <= mn_) ;
00375 #endif
00376 return vm1_[i];
00377 }
00378
00379
00380
00381 inline reference operator()(Subscript i, Subscript j)
00382 {
00383 #ifdef TNT_BOUNDS_CHECK
00384 assert(1<=i);
00385 assert(i <= m_) ;
00386 assert(1<=j);
00387 assert(j <= n_);
00388 #endif
00389 return rowm1_[i][j];
00390 }
00391
00392
00393
00394 inline const_reference operator() (Subscript i, Subscript j) const
00395 {
00396 #ifdef TNT_BOUNDS_CHECK
00397 assert(1<=i);
00398 assert(i <= m_) ;
00399 assert(1<=j);
00400 assert(j <= n_);
00401 #endif
00402 return rowm1_[i][j];
00403 }
00404
00405
00406
00407 Vector<T> diag() const
00408 {
00409 Subscript N = n_ < m_ ? n_ : m_;
00410 Vector<T> d(N);
00411
00412 for (int i=0; i<N; i++)
00413 d[i] = row_[i][i];
00414
00415 return d;
00416 }
00417
00418 Matrix<T> upper_triangular() const;
00419 Matrix<T> lower_triangular() const;
00420
00421
00422
00423 };
00424
00425
00426
00427
00428
00429 template <class T>
00430 std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
00431 {
00432 Subscript M=A.num_rows();
00433 Subscript N=A.num_cols();
00434
00435 s << M << " " << N << "\n";
00436 for (Subscript i=0; i<M; i++)
00437 {
00438 for (Subscript j=0; j<N; j++)
00439 {
00440 s << A[i][j] << " ";
00441 }
00442 s << "\n";
00443 }
00444
00445
00446 return s;
00447 }
00448
00449
00450 template <class T>
00451 std::istream& operator>>(std::istream &s, Matrix<T> &A)
00452 {
00453
00454 Subscript M, N;
00455
00456 s >> M >> N;
00457
00458 if ( !(M == A.num_rows() && N == A.num_cols() ))
00459 {
00460 A.newsize(M,N);
00461 }
00462
00463
00464 for (Subscript i=0; i<M; i++)
00465 for (Subscript j=0; j<N; j++)
00466 {
00467 s >> A[i][j];
00468 }
00469
00470
00471 return s;
00472 }
00473
00474
00475
00476
00477
00491 template <class T>
00492 Matrix<T> & mult(Matrix<T>& C, const Matrix<T> &A, const Matrix<T> &B)
00493 {
00494
00495 #ifdef TNT_BOUNDS_CHECK
00496 assert(A.num_cols() == B.num_rows());
00497 #endif
00498
00499 Subscript M = A.num_rows();
00500 Subscript N = A.num_cols();
00501 Subscript K = B.num_cols();
00502
00503 #ifdef TNT_BOUNDS_CHECK
00504 assert(C.num_rows() == M);
00505 assert(C.num_cols() == K);
00506 #endif
00507
00508 T sum;
00509
00510 const T* row_i;
00511 const T* col_k;
00512
00513 for (Subscript i=0; i<M; i++)
00514 for (Subscript k=0; k<K; k++)
00515 {
00516 row_i = &(A[i][0]);
00517 col_k = &(B[0][k]);
00518 sum = 0;
00519 for (Subscript j=0; j<N; j++)
00520 {
00521 sum += *row_i * *col_k;
00522 row_i++;
00523 col_k += K;
00524 }
00525 C[i][k] = sum;
00526 }
00527
00528 return C;
00529 }
00530
00531
00540 template <class T>
00541 inline Matrix<T> mult(const Matrix<T> &A, const Matrix<T> &B)
00542 {
00543
00544 #ifdef TNT_BOUNDS_CHECK
00545 assert(A.num_cols() == B.num_rows());
00546 #endif
00547
00548 Subscript M = A.num_rows();
00549 Subscript N = A.num_cols();
00550 Subscript K = B.num_cols();
00551
00552 Matrix<T> tmp(M,K);
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566 mult(tmp, A, B);
00567
00568 return tmp;
00569 }
00570
00579 template <class T>
00580 inline Matrix<T> operator*(const Matrix<T> &A, const Matrix<T> &B)
00581 {
00582 return mult(A,B);
00583 }
00584
00594 template <class T>
00595 inline Vector<T> mult(const Matrix<T> &A, const Vector<T> &b)
00596 {
00597
00598 #ifdef TNT_BOUNDS_CHECK
00599 assert(A.num_cols() == b.dim());
00600 #endif
00601
00602 Subscript M = A.num_rows();
00603 Subscript N = A.num_cols();
00604
00605 Vector<T> tmp(M);
00606 T sum;
00607
00608 for (Subscript i=0; i<M; i++)
00609 {
00610 sum = 0;
00611 for (Subscript j=0; j<N; j++)
00612 sum = sum + A[i][j] * b[j];
00613
00614 tmp[i] = sum;
00615 }
00616
00617 return tmp;
00618 }
00619
00629 template <class T>
00630 inline Vector<T> operator*(const Matrix<T> &A, const Vector<T> &b)
00631 {
00632 return mult(A,b);
00633 }
00634
00635
00647 template <class T>
00648 inline Matrix<T> mult(const T& s, const Matrix<T> &A)
00649 {
00650 Subscript M = A.num_rows();
00651 Subscript N = A.num_cols();
00652
00653 Matrix<T> R(M,N);
00654 for (int i=0; i<M; i++)
00655 for (int j=0; j<N; j++)
00656 R[i][j] = s * A[i][j];
00657
00658 return R;
00659 }
00660
00675 template <class T>
00676 inline Matrix<T> mult(const Matrix<T> &A, const T& s)
00677 {
00678 return mult(s, A);
00679 }
00680
00681
00694 template <class T>
00695 inline Matrix<T> mult_eq(const T& s, Matrix<T> &A)
00696 {
00697 Subscript M = A.num_rows();
00698 Subscript N = A.num_cols();
00699
00700 for (int i=0; i<M; i++)
00701 for (int j=0; j<N; j++)
00702 A[i][j] *= s;
00703 }
00704
00705 template <class T>
00706 inline Matrix<T> mult_eq(Matrix<T> &A, const T&a)
00707 {
00708 return mult_eq(a, A);
00709 }
00710
00711
00723 template <class T>
00724 inline Matrix<T> transpose_mult(const Matrix<T> &A, const Matrix<T> &B)
00725 {
00726
00727 #ifdef TNT_BOUNDS_CHECK
00728 assert(A.num_rows() == B.num_rows());
00729 #endif
00730
00731 Subscript M = A.num_cols();
00732 Subscript N = A.num_rows();
00733 Subscript K = B.num_cols();
00734
00735 Matrix<T> tmp(M,K);
00736 T sum;
00737
00738 for (Subscript i=0; i<N; i++)
00739 for (Subscript k=0; k<K; k++)
00740 {
00741 sum = 0;
00742 for (Subscript j=0; j<M; j++)
00743 sum = sum + A[j][i] * B[j][k];
00744
00745 tmp[i][k] = sum;
00746 }
00747
00748 return tmp;
00749 }
00750
00762 template <class T>
00763 inline Vector<T> transpose_mult(const Matrix<T> &A, const Vector<T> &b)
00764 {
00765
00766 #ifdef TNT_BOUNDS_CHECK
00767 assert(A.num_rows() == b.dim());
00768 #endif
00769
00770 Subscript M = A.num_cols();
00771 Subscript N = A.num_rows();
00772
00773 Vector<T> tmp(M);
00774
00775 for (Subscript i=0; i<M; i++)
00776 {
00777 T sum = 0;
00778 for (Subscript j=0; j<N; j++)
00779 sum = sum + A[j][i] * b[j];
00780
00781 tmp[i] = sum;
00782 }
00783
00784 return tmp;
00785 }
00786
00787
00796 template <class T>
00797 Matrix<T> add(const Matrix<T> &A, const Matrix<T> &B)
00798 {
00799 Subscript M = A.num_rows();
00800 Subscript N = A.num_cols();
00801
00802 assert(M==B.num_rows());
00803 assert(N==B.num_cols());
00804
00805 Matrix<T> tmp(M,N);
00806 Subscript i,j;
00807
00808 for (i=0; i<M; i++)
00809 for (j=0; j<N; j++)
00810 tmp[i][j] = A[i][j] + B[i][j];
00811
00812 return tmp;
00813 }
00814
00826 template <class T>
00827 inline Matrix<T> operator+(const Matrix<T> &A, const Matrix<T> &B)
00828 {
00829 return add(A,B);
00830 }
00831
00832
00833
00843 template <class T>
00844 Matrix<T>& add_eq(Matrix<T> &A, const Matrix<T> &B)
00845 {
00846 Subscript M = A.num_rows();
00847 Subscript N = A.num_cols();
00848
00849 assert(M==B.num_rows());
00850 assert(N==B.num_cols());
00851
00852 Matrix<T> tmp(M,N);
00853 Subscript i,j;
00854
00855 for (i=0; i<M; i++)
00856 for (j=0; j<N; j++)
00857 tmp[i][j] = A[i][j] + B[i][j];
00858
00859 return A += tmp;
00860 }
00861
00873 template <class T>
00874 inline Matrix<T> operator+=(Matrix<T> &A, const Matrix<T> &B)
00875 {
00876 return add_eq(A,B);
00877 }
00878
00888 template <class T>
00889 Matrix<T> minus(const Matrix <T>& A, const Matrix<T> &B)
00890 {
00891 Subscript M = A.num_rows();
00892 Subscript N = A.num_cols();
00893
00894 assert(M==B.num_rows());
00895 assert(N==B.num_cols());
00896
00897 Matrix<T> tmp(M,N);
00898 Subscript i,j;
00899
00900 for (i=0; i<M; i++)
00901 for (j=0; j<N; j++)
00902 tmp[i][j] = A[i][j] - B[i][j];
00903
00904 return tmp;
00905
00906 }
00907
00919 template <class T>
00920 inline Matrix<T> operator-(const Matrix<T> &A,
00921 const Matrix<T> &B)
00922 {
00923 return minus(A,B);
00924 }
00925
00926
00937 template <class T>
00938 Matrix<T> mult_element(const Matrix<T> &A, const Matrix<T> &B)
00939 {
00940 Subscript M = A.num_rows();
00941 Subscript N = A.num_cols();
00942
00943 assert(M==B.num_rows());
00944 assert(N==B.num_cols());
00945
00946 Matrix<T> tmp(M,N);
00947 Subscript i,j;
00948
00949 for (i=0; i<M; i++)
00950 for (j=0; j<N; j++)
00951 tmp[i][j] = A[i][j] * B[i][j];
00952
00953 return tmp;
00954 }
00955
00966 template <class T>
00967 Matrix<T>& mult_element_eq(Matrix<T> &A, const Matrix<T> &B)
00968 {
00969 Subscript M = A.num_rows();
00970 Subscript N = A.num_cols();
00971
00972 assert(M==B.num_rows());
00973 assert(N==B.num_cols());
00974
00975 Subscript i,j;
00976
00977 for (i=0; i<M; i++)
00978 for (j=0; j<N; j++)
00979 A[i][j] *= B[i][j];
00980
00981 return A;
00982 }
00983
00984
00995 template <class T>
00996 double norm(const Matrix<T> &A)
00997 {
00998 Subscript M = A.num_rows();
00999 Subscript N = A.num_cols();
01000
01001 T sum = 0.0;
01002 for (int i=1; i<=M; i++)
01003 for (int j=1; j<=N; j++)
01004 sum += A(i,j) * A(i,j);
01005 return sqrt(sum);
01006
01007 }
01008
01018 template <class T>
01019 Matrix<T> transpose(const Matrix<T> &A)
01020 {
01021 Subscript M = A.num_rows();
01022 Subscript N = A.num_cols();
01023
01024 Matrix<T> S(N,M);
01025 Subscript i, j;
01026
01027 for (i=0; i<M; i++)
01028 for (j=0; j<N; j++)
01029 S[j][i] = A[i][j];
01030
01031 return S;
01032 }
01033
01034
01035
01036
01037
01038
01039
01051 template <class T>
01052 Vector<T> upper_triangular_solve(const Matrix<T> &A, const Vector<T> &b)
01053 {
01054 int n = A.num_rows() < A.num_cols() ? A.num_rows() : A.num_cols();
01055 Vector<T> x(b);
01056 for (int k=n; k >= 1; k--)
01057 {
01058 x(k) /= A(k,k);
01059 for (int i=1; i< k; i++ )
01060 x(i) -= x(k) * A(i,k);
01061 }
01062
01063 return x;
01064 }
01065
01077 template <class T>
01078 Vector<T> lower_triangular_solve(const Matrix<T> &A, const Vector<T> &b)
01079 {
01080 int n = A.num_rows() < A.num_cols() ? A.num_rows() : A.num_cols();
01081 Vector<T> x(b);
01082 for (int k=1; k <= n; k++)
01083 {
01084 x(k) /= A(k,k);
01085 for (int i=k+1; i<= n; i++ )
01086 x(i) -= x(k) * A(i,k);
01087 }
01088
01089 return x;
01090 }
01091
01092
01093
01094 }
01095
01096 #endif
01097