Back to index

texmacs  1.0.7.15
matrix.hpp
Go to the documentation of this file.
00001 
00002 /******************************************************************************
00003 * MODULE     : matrix.hpp
00004 * DESCRIPTION: matrixs with reference counting and pointer copying
00005 * COPYRIGHT  : (C) 2006  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license version 3 or later.
00008 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
00009 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
00010 ******************************************************************************/
00011 
00012 #ifndef MATRIX_H
00013 #define MATRIX_H
00014 #include "vector.hpp"
00015 #include "ntuple.hpp"
00016 #define TMPL template<typename T>
00017 #define BINARY_TMPL template<typename T, typename U>
00018 #define R typename unary_properties<T>::norm_type
00019 #define M typename binary_properties<T,U>::product_type
00020 
00021 TMPL class matrix;
00022 TMPL int NR (matrix<T> m);
00023 TMPL int NC (matrix<T> m);
00024 TMPL T*  A  (matrix<T> m);
00025 
00026 /******************************************************************************
00027 * The matrix class
00028 ******************************************************************************/
00029 
00030 TMPL
00031 class matrix_rep: concrete_struct {
00032   int rows;
00033   int cols;
00034   T* a;
00035 public:
00036   inline matrix_rep (T* a2, int rows2, int cols2):
00037     rows (rows2), cols (cols2), a (a2) {}
00038   inline ~matrix_rep () { if (a != NULL) tm_delete_array (a); }
00039   friend class matrix<T>;
00040   friend int NR LESSGTR (matrix<T> m);
00041   friend int NC LESSGTR (matrix<T> m);
00042   friend T*  A  LESSGTR (matrix<T> m);
00043 };
00044 
00045 TMPL
00046 class matrix {
00047 CONCRETE_TEMPLATE(matrix,T);
00048   inline matrix (T *a, int rows, int cols):
00049     rep (tm_new<matrix_rep<T> > (a, rows, cols)) {}
00050   inline matrix (T c, int rows, int cols) {
00051     int i, n= rows * cols;
00052     T* a= (n == 0? (T*) NULL: tm_new_array<T> (n));
00053     for (i=0; i<n; i++)
00054       a[i]= ((i%(cols+1)) == 0? c: T(0));
00055     rep= tm_new<matrix_rep<T> > (a, rows, cols); }
00056   inline matrix () {
00057     rep= tm_new<matrix_rep<T> > ((T*) NULL, 0, 0); }
00058   inline T& operator () (int i, int j) {
00059     return rep->a[i*rep->cols + j]; }
00060 };
00061 CONCRETE_TEMPLATE_CODE(matrix,typename,T);
00062 
00063 TMPL inline int NR (matrix<T> m) { return m->rows; }
00064 TMPL inline int NC (matrix<T> m) { return m->cols; }
00065 TMPL inline T*  A  (matrix<T> m) { return m->a; }
00066 
00067 TMPL
00068 class unary_properties<matrix<T> > {
00069 public:
00070   typedef T scalar_type;
00071   typedef typename unary_properties<T>::norm_type norm_type;
00072   typedef pair<int,int> index_type;
00073   static inline tree index_name (index_type i) {
00074     return tree (RSUB, tree (RSUB, "x", as_string (i.x1 + 1)),
00075                as_string (i.x2 + 1)); }
00076   static inline scalar_type access (matrix<T> m, index_type var) {
00077     return m (var.x1, var.x2); }
00078 };
00079 
00080 BINARY_TMPL
00081 class binary_properties<matrix<T>,matrix<U> > {
00082 public:
00083   typedef matrix<M > product_type;
00084 };
00085 
00086 BINARY_TMPL
00087 class binary_properties<T,matrix<U> > {
00088 public:
00089   typedef matrix<M > product_type;
00090 };
00091 
00092 BINARY_TMPL
00093 class binary_properties<matrix<T>,U > {
00094 public:
00095   typedef matrix<M > product_type;
00096 };
00097 
00098 /******************************************************************************
00099 * Conversions
00100 ******************************************************************************/
00101 
00102 TMPL tree
00103 as_tree (matrix<T> m) {
00104   int i, rows= NR (m);
00105   int j, cols= NC (m);
00106   tree t (TUPLE, rows);
00107   for (i=0; i<rows; i++) {
00108     t[i]= tree (TUPLE, cols);
00109     for (j=0; j<cols; j++)
00110       t[i][j]= as_tree (m(i,j));
00111   }
00112   return t;
00113 }
00114 
00115 TMPL inline tm_ostream&
00116 operator << (tm_ostream& out, matrix<T> m) {
00117   return out << as_math_string (as_tree (m));
00118 }
00119 
00120 TMPL void
00121 parse (tree t, matrix<T>& m) {
00122   ASSERT (is_tuple (t) && N(t)>0 && is_tuple (t[0]), "not a matrix");
00123   int i, j, rows= N(t), cols= N(t[0]);
00124   m= matrix<T> (T(0), rows, cols);
00125   for (i=0; i<rows; i++)
00126     for (j=0; j<cols; j++)
00127       m(i,j)= parse_as<T> (t[i][j]);
00128 }
00129 
00130 TMPL inline matrix<T>
00131 as_matrix (tree t) {
00132   matrix<T> result;
00133   parse (t, result);
00134   return result;
00135 }
00136 
00137 /******************************************************************************
00138 * Vectorial operations on matrices
00139 ******************************************************************************/
00140 
00141 template<typename T, typename Op> matrix<T>
00142 unary (matrix<T> m) {
00143   int i, n= NR(m) * NC(m);
00144   T* a= A(m);
00145   T* r= tm_new_array<T> (n);
00146   for (i=0; i<n; i++)
00147     r[i]= Op::eval (a[i]);
00148   return matrix<T> (r, NR(m), NC(m));
00149 }
00150 
00151 TMPL inline matrix<T>
00152 copy (matrix<T> m) {
00153   return unary<T,copy_op> (m); }
00154 
00155 TMPL inline matrix<T>
00156 operator - (matrix<T> m) {
00157   return unary<T,neg_op> (m); }
00158 
00159 template<typename T, typename Op> matrix<T>
00160 binary (matrix<T> m1, matrix<T> m2) {
00161   int i, n= NR(m1) * NC(m1);
00162   ASSERT (NR(m1) == NR(m2) && NC(m1) == NC(m2), "matrix sizes don't match");
00163   T* a= A(m1);
00164   T* b= A(m2);
00165   T* r= tm_new_array<T> (n);
00166   for (i=0; i<n; i++)
00167     r[i]= Op::eval (a[i], b[i]);
00168   return matrix<T> (r, NR(m1), NC(m1));
00169 }
00170 
00171 TMPL inline matrix<T>
00172 operator + (matrix<T> m1, matrix<T> m2) {
00173   return binary<T,add_op> (m1, m2); }
00174 
00175 TMPL inline matrix<T>
00176 operator - (matrix<T> m1, matrix<T> m2) {
00177   return binary<T,sub_op> (m1, m2); }
00178 
00179 /******************************************************************************
00180 * Multiplication
00181 ******************************************************************************/
00182 
00183 TMPL inline matrix<T>
00184 operator * (matrix<T> m1, matrix<T> m2) {
00185   int i, j, k, rows= NR (m1), aux= NC (m1), cols= NC (m2);
00186   ASSERT (NR (m2) == aux, "dimensions don't match");
00187   matrix<T> prod (T(0), rows, cols);
00188   for (i=0; i<rows; i++)
00189     for (j=0; j<cols; j++)
00190       for (k=0; k<aux; k++)
00191        prod (i, j) += m1 (i, k) * m2 (k, j);
00192   return prod;
00193 }
00194 
00195 TMPL inline vector<T>
00196 operator * (matrix<T> m, vector<T> v) {
00197   int i, j, rows= NR (m), cols= NC (m);
00198   ASSERT (N (v) == cols, "dimensions don't match");
00199   vector<T> prod (T(0), rows);
00200   for (i=0; i<rows; i++)
00201     for (j=0; j<cols; j++)
00202       prod [i] += m (i, j) * v[j];
00203   return prod;
00204 }
00205 
00206 TMPL inline array<T>
00207 operator * (matrix<T> m, array<T> v) {
00208   int i, j, rows= NR (m), cols= NC (m);
00209   ASSERT (N (v) == cols, "dimensions don't match");
00210   array<T> prod (T(0), rows);
00211   for (i=0; i<rows; i++)
00212     for (j=0; j<cols; j++)
00213       prod [i] += m (i, j) * v[j];
00214   return prod;
00215 }
00216 
00217 #undef TMPL
00218 #undef BINARY_TMPL
00219 #undef R
00220 #undef M
00221 #endif // defined MATRIX_H