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_CMAT_H
00025 #define TNT_CMAT_H
00026
00027 #include "tnt_subscript.h"
00028 #include "tnt_vec.h"
00029 #include <cstdlib>
00030 #include <cassert>
00031 #include <iostream>
00032 #include <sstream>
00033
00034 namespace TNT
00035 {
00036
00037
00038 template <class T>
00039 class Matrix
00040 {
00041
00042
00043 public:
00044
00045 typedef Subscript size_type;
00046 typedef T value_type;
00047 typedef T element_type;
00048 typedef T* pointer;
00049 typedef T* iterator;
00050 typedef T& reference;
00051 typedef const T* const_iterator;
00052 typedef const T& const_reference;
00053
00054 Subscript lbound() const { return 1;}
00055
00056 protected:
00057 Subscript m_;
00058 Subscript n_;
00059 Subscript mn_;
00060 T* v_;
00061 T** row_;
00062 T* vm1_ ;
00063 T** rowm1_;
00064
00065
00066
00067
00068 void initialize(Subscript M, Subscript N)
00069 {
00070 mn_ = M*N;
00071 m_ = M;
00072 n_ = N;
00073
00074 v_ = new T[mn_];
00075 row_ = new T*[M];
00076 rowm1_ = new T*[M];
00077
00078 assert(v_ != NULL);
00079 assert(row_ != NULL);
00080 assert(rowm1_ != NULL);
00081
00082 T* p = v_;
00083 vm1_ = v_ - 1;
00084 for (Subscript i=0; i<M; i++)
00085 {
00086 row_[i] = p;
00087 rowm1_[i] = p-1;
00088 p += N ;
00089
00090 }
00091
00092 rowm1_ -- ;
00093 }
00094
00095 void copy(const T* v)
00096 {
00097 Subscript N = m_ * n_;
00098 Subscript i;
00099
00100 #ifdef TNT_UNROLL_LOOPS
00101 Subscript Nmod4 = N & 3;
00102 Subscript N4 = N - Nmod4;
00103
00104 for (i=0; i<N4; i+=4)
00105 {
00106 v_[i] = v[i];
00107 v_[i+1] = v[i+1];
00108 v_[i+2] = v[i+2];
00109 v_[i+3] = v[i+3];
00110 }
00111
00112 for (i=N4; i< N; i++)
00113 v_[i] = v[i];
00114 #else
00115
00116 for (i=0; i< N; i++)
00117 v_[i] = v[i];
00118 #endif
00119 }
00120
00121 void set(const T& val)
00122 {
00123 Subscript N = m_ * n_;
00124 Subscript i;
00125
00126 #ifdef TNT_UNROLL_LOOPS
00127 Subscript Nmod4 = N & 3;
00128 Subscript N4 = N - Nmod4;
00129
00130 for (i=0; i<N4; i+=4)
00131 {
00132 v_[i] = val;
00133 v_[i+1] = val;
00134 v_[i+2] = val;
00135 v_[i+3] = val;
00136 }
00137
00138 for (i=N4; i< N; i++)
00139 v_[i] = val;
00140 #else
00141
00142 for (i=0; i< N; i++)
00143 v_[i] = val;
00144
00145 #endif
00146 }
00147
00148
00149
00150 void destroy()
00151 {
00152
00153 if (v_ == NULL) return ;
00154
00155
00156 if (v_ != NULL) delete [] (v_);
00157 if (row_ != NULL) delete [] (row_);
00158
00159
00160 rowm1_ ++;
00161 if (rowm1_ != NULL ) delete [] (rowm1_);
00162 }
00163
00164
00165 public:
00166
00167 operator T**(){ return row_; }
00168 operator T**() const { return row_; }
00169
00170
00171 Subscript size() const { return mn_; }
00172
00173
00174
00175 Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
00176
00177 Matrix(const Matrix<T> &A)
00178 {
00179 initialize(A.m_, A.n_);
00180 copy(A.v_);
00181 }
00182
00183 Matrix(Subscript M, Subscript N, const T& value = T())
00184 {
00185 initialize(M,N);
00186 set(value);
00187 }
00188
00189 Matrix(Subscript M, Subscript N, const T* v)
00190 {
00191 initialize(M,N);
00192 copy(v);
00193 }
00194
00195 Matrix(Subscript M, Subscript N, const char *s)
00196 {
00197 initialize(M,N);
00198
00199 std::istringstream ins(s);
00200
00201 Subscript i, j;
00202
00203 for (i=0; i<M; i++)
00204 for (j=0; j<N; j++)
00205 ins >> row_[i][j];
00206 }
00207
00208
00209
00210 ~Matrix()
00211 {
00212 destroy();
00213 }
00214
00215
00216
00217
00218 Matrix<T>& newsize(Subscript M, Subscript N)
00219 {
00220 if (num_rows() == M && num_cols() == N)
00221 return *this;
00222
00223 destroy();
00224 initialize(M,N);
00225
00226 return *this;
00227 }
00228
00229
00230
00231
00232
00233
00234 Matrix<T>& operator=(const Matrix<T> &A)
00235 {
00236 if (v_ == A.v_)
00237 return *this;
00238
00239 if (m_ == A.m_ && n_ == A.n_)
00240 copy(A.v_);
00241
00242 else
00243 {
00244 destroy();
00245 initialize(A.m_, A.n_);
00246 copy(A.v_);
00247 }
00248
00249 return *this;
00250 }
00251
00252 Matrix<T>& operator=(const T& scalar)
00253 {
00254 set(scalar);
00255 return *this;
00256 }
00257
00258
00259 Subscript dim(Subscript d) const
00260 {
00261 #ifdef TNT_BOUNDS_CHECK
00262 assert( d >= 1);
00263 assert( d <= 2);
00264 #endif
00265 return (d==1) ? m_ : ((d==2) ? n_ : 0);
00266 }
00267
00268 Subscript num_rows() const { return m_; }
00269 Subscript num_cols() const { return n_; }
00270
00271
00272
00273
00274 inline T* operator[](Subscript i)
00275 {
00276 #ifdef TNT_BOUNDS_CHECK
00277 assert(0<=i);
00278 assert(i < m_) ;
00279 #endif
00280 return row_[i];
00281 }
00282
00283 inline const T* operator[](Subscript i) const
00284 {
00285 #ifdef TNT_BOUNDS_CHECK
00286 assert(0<=i);
00287 assert(i < m_) ;
00288 #endif
00289 return row_[i];
00290 }
00291
00292 inline reference operator()(Subscript i)
00293 {
00294 #ifdef TNT_BOUNDS_CHECK
00295 assert(1<=i);
00296 assert(i <= mn_) ;
00297 #endif
00298 return vm1_[i];
00299 }
00300
00301 inline const_reference operator()(Subscript i) const
00302 {
00303 #ifdef TNT_BOUNDS_CHECK
00304 assert(1<=i);
00305 assert(i <= mn_) ;
00306 #endif
00307 return vm1_[i];
00308 }
00309
00310
00311
00312 inline reference operator()(Subscript i, Subscript j)
00313 {
00314 #ifdef TNT_BOUNDS_CHECK
00315 assert(1<=i);
00316 assert(i <= m_) ;
00317 assert(1<=j);
00318 assert(j <= n_);
00319 #endif
00320 return rowm1_[i][j];
00321 }
00322
00323
00324
00325 inline const_reference operator() (Subscript i, Subscript j) const
00326 {
00327 #ifdef TNT_BOUNDS_CHECK
00328 assert(1<=i);
00329 assert(i <= m_) ;
00330 assert(1<=j);
00331 assert(j <= n_);
00332 #endif
00333 return rowm1_[i][j];
00334 }
00335
00336
00337
00338
00339 };
00340
00341
00342
00343
00344 template <class T>
00345 std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
00346 {
00347 Subscript M=A.num_rows();
00348 Subscript N=A.num_cols();
00349
00350 s << M << " " << N << "\n";
00351
00352 for (Subscript i=0; i<M; i++)
00353 {
00354 for (Subscript j=0; j<N; j++)
00355 {
00356 s << A[i][j] << " ";
00357 }
00358 s << "\n";
00359 }
00360
00361
00362 return s;
00363 }
00364
00365 template <class T>
00366 std::istream& operator>>(std::istream &s, Matrix<T> &A)
00367 {
00368
00369 Subscript M, N;
00370
00371 s >> M >> N;
00372
00373 if ( !(M == A.num_rows() && N == A.num_cols() ))
00374 {
00375 A.newsize(M,N);
00376 }
00377
00378
00379 for (Subscript i=0; i<M; i++)
00380 for (Subscript j=0; j<N; j++)
00381 {
00382 s >> A[i][j];
00383 }
00384
00385
00386 return s;
00387 }
00388
00389
00390
00391
00392 template <class T>
00393 Matrix<T> operator+(const Matrix<T> &A,
00394 const Matrix<T> &B)
00395 {
00396 Subscript M = A.num_rows();
00397 Subscript N = A.num_cols();
00398
00399 assert(M==B.num_rows());
00400 assert(N==B.num_cols());
00401
00402 Matrix<T> tmp(M,N);
00403 Subscript i,j;
00404
00405 for (i=0; i<M; i++)
00406 for (j=0; j<N; j++)
00407 tmp[i][j] = A[i][j] + B[i][j];
00408
00409 return tmp;
00410 }
00411
00412 template <class T>
00413 Matrix<T> operator-(const Matrix<T> &A,
00414 const Matrix<T> &B)
00415 {
00416 Subscript M = A.num_rows();
00417 Subscript N = A.num_cols();
00418
00419 assert(M==B.num_rows());
00420 assert(N==B.num_cols());
00421
00422 Matrix<T> tmp(M,N);
00423 Subscript i,j;
00424
00425 for (i=0; i<M; i++)
00426 for (j=0; j<N; j++)
00427 tmp[i][j] = A[i][j] - B[i][j];
00428
00429 return tmp;
00430 }
00431
00432 template <class T>
00433 Matrix<T> mult_element(const Matrix<T> &A,
00434 const Matrix<T> &B)
00435 {
00436 Subscript M = A.num_rows();
00437 Subscript N = A.num_cols();
00438
00439 assert(M==B.num_rows());
00440 assert(N==B.num_cols());
00441
00442 Matrix<T> tmp(M,N);
00443 Subscript i,j;
00444
00445 for (i=0; i<M; i++)
00446 for (j=0; j<N; j++)
00447 tmp[i][j] = A[i][j] * B[i][j];
00448
00449 return tmp;
00450 }
00451
00452
00453 template <class T>
00454 Matrix<T> transpose(const Matrix<T> &A)
00455 {
00456 Subscript M = A.num_rows();
00457 Subscript N = A.num_cols();
00458
00459 Matrix<T> S(N,M);
00460 Subscript i, j;
00461
00462 for (i=0; i<M; i++)
00463 for (j=0; j<N; j++)
00464 S[j][i] = A[i][j];
00465
00466 return S;
00467 }
00468
00469
00470
00471 template <class T>
00472 inline Matrix<T> matmult(const Matrix<T> &A,
00473 const Matrix<T> &B)
00474 {
00475
00476 #ifdef TNT_BOUNDS_CHECK
00477 assert(A.num_cols() == B.num_rows());
00478 #endif
00479
00480 Subscript M = A.num_rows();
00481 Subscript N = A.num_cols();
00482 Subscript K = B.num_cols();
00483
00484 Matrix<T> tmp(M,K);
00485 T sum;
00486
00487 for (Subscript i=0; i<M; i++)
00488 for (Subscript k=0; k<K; k++)
00489 {
00490 sum = 0;
00491 for (Subscript j=0; j<N; j++)
00492 sum = sum + A[i][j] * B[j][k];
00493
00494 tmp[i][k] = sum;
00495 }
00496
00497 return tmp;
00498 }
00499
00500 template <class T>
00501 inline Matrix<T> operator*(const Matrix<T> &A,
00502 const Matrix<T> &B)
00503 {
00504 return matmult(A,B);
00505 }
00506
00507 template <class T>
00508 inline int matmult(Matrix<T>& C, const Matrix<T> &A,
00509 const Matrix<T> &B)
00510 {
00511
00512 assert(A.num_cols() == B.num_rows());
00513
00514 Subscript M = A.num_rows();
00515 Subscript N = A.num_cols();
00516 Subscript K = B.num_cols();
00517
00518 C.newsize(M,K);
00519
00520 T sum;
00521
00522 const T* row_i;
00523 const T* col_k;
00524
00525 for (Subscript i=0; i<M; i++)
00526 for (Subscript k=0; k<K; k++)
00527 {
00528 row_i = &(A[i][0]);
00529 col_k = &(B[0][k]);
00530 sum = 0;
00531 for (Subscript j=0; j<N; j++)
00532 {
00533 sum += *row_i * *col_k;
00534 row_i++;
00535 col_k += K;
00536 }
00537 C[i][k] = sum;
00538 }
00539
00540 return 0;
00541 }
00542
00543
00544 template <class T>
00545 Vector<T> matmult(const Matrix<T> &A, const Vector<T> &x)
00546 {
00547
00548 #ifdef TNT_BOUNDS_CHECK
00549 assert(A.num_cols() == x.dim());
00550 #endif
00551
00552 Subscript M = A.num_rows();
00553 Subscript N = A.num_cols();
00554
00555 Vector<T> tmp(M);
00556 T sum;
00557
00558 for (Subscript i=0; i<M; i++)
00559 {
00560 sum = 0;
00561 const T* rowi = A[i];
00562 for (Subscript j=0; j<N; j++)
00563 sum = sum + rowi[j] * x[j];
00564
00565 tmp[i] = sum;
00566 }
00567
00568 return tmp;
00569 }
00570
00571 template <class T>
00572 inline Vector<T> operator*(const Matrix<T> &A, const Vector<T> &x)
00573 {
00574 return matmult(A,x);
00575 }
00576
00577 }
00578
00579 #endif
00580