Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members   Examples  

matrices.h

Go to the documentation of this file.
00001 //N.Yu.Zolotykh 1999, 2000
00002 //University of Nizhni Novgorod, Russia
00003 
00004 #include<iostream.h>
00005 #include<stdlib.h>
00006 #include<stdio.h>
00007 #include<strstrea.h>
00008 #include<string.h>
00009 #include"arageli.h"
00010 
00011 #ifndef MATRICES_H_
00012 #define MATRICES_H_
00013 
00014 
00035 
00038 const int st_ok     =  1;                    
00039 const int st_out_of_memory  = -1;            
00040 const int st_out_of_range = -2;              
00041 const int st_inconsistent_sizes = -3;        
00042 const int st_input_error  = -4;              
00043 const int st_output_error   = -5;            
00044 const int st_incorrect_index  = -6;          
00045 const int st_matrix_must_be_square = -7;     
00046 const int st_matrix_is_singular = -8;        
00047 
00048 
00053 int matrix_status = st_ok;
00054 
00055 
00057 typedef int index; // index >= -1
00058 
00063 index const n_max_size = 10000;
00064 
00066 index const m_max_size = 10000;
00067 
00068 
00069 #if defined(__BORLANDC__) && defined(__WIN32__) //bc 5.5
00070    template <class T> class vector;
00071    template <class T> class matrix;
00072    template <class T> ostream & operator << (ostream & s, const vector<T> & x);
00073    template <class T> istream & operator >> (istream & s, vector<T> & x);
00074    template <class T> vector<T> operator +  (const vector<T> & a, const vector<T> & b);
00075    template <class T> vector<T> operator -  (const vector<T> & a, const vector<T> & b);
00076    template <class T> vector<T> operator * (const T alpha, const vector<T> & a);
00077    template <class T> T operator * (const vector<T> & a, const vector<T> & b);
00078    template <class T> int operator == (const vector<T> & a, const vector<T> & b);
00079    template <class T> int operator != (const vector<T> & a, const vector<T> & b);
00080    template <class T> void pretty(ostr & s, const vector<T>& A);
00081    template <class T> ostream & operator << (ostream & s, const matrix<T> & x);
00082    template <class T> istream & operator >> (istream & s, matrix<T> & x); // ?
00083    template <class T> matrix<T> operator + (const matrix<T> &, const matrix<T> &);
00084    template <class T> matrix<T> operator * (const matrix<T> &, const matrix<T> &);
00085    template <class T> matrix<T> operator * (T alpha, const matrix<T> & a);
00086    template <class T> matrix<T> operator * (const matrix<T> & a, const vector<T> & v);
00087    template <class T> matrix<T> operator * ( const vector<T> & v, const matrix<T> & a);
00088    template <class T> matrix<T> transpose (const matrix<T> & a);
00089    template <class T> int operator == (const matrix<T> & a, const matrix<T> & b);
00090    template <class T> int operator != (const matrix<T> & a, const matrix<T> & b);
00091    template <class T> matrix<T> diag(const vector<T>& d);
00092    template <class T> matrix<T> diag(T a, index m, index n);
00093    template <class T> matrix<T> diag(T a, index m);
00094    template <class T> matrix<T> eye (T a, index m, index n);
00095    template <class T> matrix<T> eye (T a, index m);
00096    template <class T> matrix<T> fill(T a, index m, index n);
00097    template <class T> matrix<T> fill(T a, index m);
00098    template <class T> matrix<T> zero(T a, index m, index n);
00099    template <class T> matrix<T> zero(T a, index m);
00100    template <class T> void pretty(ostr & s, const matrix<T>& A);
00101    template <class T> void print2(ostr & s, const matrix<T>& A, const matrix<T>& B);
00102    template <class T> void print3(ostr & s, const matrix<T>& A, const matrix<T>& B, const matrix<T>& C);
00103 #endif
00104    
00105    
00106 /***************************/
00107 /*                         */
00108 /*     class vector<>      */
00109 /*                         */
00110 /***************************/
00111 
00115 template <class T> class vector
00116 {
00117 public:
00118 
00120 
00124   vector (index n_init = 0, index n_delta_init = 2);
00125 
00127 
00131   vector (index n_init, T* a, index n_delta_init = 2);
00132 
00134   vector (const vector<T> &x);
00135 
00137   vector<T> & operator = (const vector<T>  &x);
00138 
00140   ~vector(void);
00141 
00143   index get_n (void) const {return n;};
00144 
00146   inline T & operator[](index j);
00147 
00149   inline T at (index j) const;
00150 
00152   inline T access(index j) const;
00153 
00155 
00159   void resize (index new_n = 0);
00160 
00162   void ins_item (index j, T item = 0);
00163 
00167   void join (T item);
00168 
00170   int member(T item) const;
00171 
00173   index index_of_item(T item) const;
00174 
00176   void del_item (index j);
00177 
00179   T take_item (index j);
00180 
00182   void swap_items (index j1, index j2);
00183 
00185   vector<T> & operator += (const vector<T> & x);
00186 
00188   vector<T> & operator -= (const vector<T> & x);
00189 
00190 #if !defined(__BORLANDC__) || !defined(__WIN32__) 
00191 
00194   friend ostream & operator << (ostream & s, const vector<T> & x);
00195 
00199   friend istream & operator >> (istream & s, vector<T> & x);
00200     // почему-то, если сделать эту функцию-друга обычной функцией,
00201     // команды вида cin >> x не воспринимаются
00202 
00204   friend vector<T> operator +  (const vector<T> & a, const vector<T> & b);
00205 
00207   friend vector<T> operator -  (const vector<T> & a, const vector<T> & b);
00208 
00210   friend vector<T> operator * (const T alpha, const vector<T> & a);
00211 
00213   friend T operator * (const vector<T> & a, const vector<T> & b);
00214 
00216   friend int operator == (const vector<T> & a, const vector<T> & b);
00217 
00219   friend int operator != (const vector<T> & a, const vector<T> & b);
00220     // почему-то, если сделать эти операции-друзья обычными операциями,
00221     // команды вида a == b * c не воспринимаются
00222 
00224   friend void pretty(ostr & s, const vector<T>& A);
00225 #else // bc 5.5
00226   friend ostream & operator << <T> (ostream & s, const vector<T> & x);
00227   friend istream & operator >> <T> (istream & s, vector<T> & x);
00228   friend vector<T> operator +  <T> (const vector<T> & a, const vector<T> & b);
00229   friend vector<T> operator -  <T> (const vector<T> & a, const vector<T> & b);
00230   friend vector<T> operator *  <T>(const T alpha, const vector<T> & a);
00231   friend T operator * <T> (const vector<T> & a, const vector<T> & b);
00232   friend int operator == <T> (const vector<T> & a, const vector<T> & b);
00233   friend int operator != <T> (const vector<T> & a, const vector<T> & b);
00234   friend void pretty<T>(ostr & s, const vector<T>& A);
00235 #endif
00236 
00237 private:
00238 
00239   index n, n_limit, n_delta;
00240   T *data;
00241   void dispose_all (void);
00242   inline int index_correct (index j) const;
00243 };
00244 
00245 
00246 /***************************/
00247 /*                         */
00248 /*      class matrix<>     */
00249 /*                         */
00250 /***************************/
00251 
00253 template <class T>
00254   class matrix
00255 {
00256 public:
00257 
00264   matrix (index m_init = 0, index n_init = 0, index m_delta_init = 2,
00265           index n_delta_init = 2);
00266 
00268 
00272   matrix (index m_init, index n_init, T* A, index m_delta_init = 2,
00273           index n_delta_init = 2);
00274 
00275 
00277   matrix (const matrix<T>  &x);
00278 
00280   matrix<T>  &operator = (const matrix<T>  &x);
00281 
00283   ~matrix (void);
00284 
00286   inline index get_m (void) const {return m;};
00288   inline index get_n (void) const {return n;};
00289 
00291 
00293   inline T & operator ()(index i, index j);
00295 
00297   inline T at (index i, index j) const;
00298 
00300 
00303   T access(index i, index j) const;
00304 
00306 
00310   void resize(index new_m = 0, index new_n = 0);
00311 
00313   vector<T> row (index j) const;
00315   matrix<T> rows (vector<index> indices) const;
00317   void ins_row (index i);
00319   void ins_row (index i, vector<T>  x);
00321 
00322   void join_row (vector<T>  x);
00324   void del_row (index i);
00326   void del_rows (vector<index> indices);
00328   vector<T> take_row (index i);
00330 
00332   matrix<T> take_rows (vector<index> indices);
00334   void add_row (index i1, index i2, T alpha);
00336 
00339   void add_rows (index i1, T alpha1, index i2, T alpha2);
00341   void mult_row (index i, T alpha);
00343   void div_row (index i, T alpha);
00345   void swap_rows (index i1, index i2);
00346 
00348   vector<T> col (index j) const;
00350   matrix<T> cols (vector<index> indices) const;
00352   void ins_col (index j);
00354   void ins_col (index j, vector<T>  x);
00356 
00357   void join_col (vector<T>  x);
00359   void del_col (index j);
00361   void del_cols (vector<index> indices);
00363   vector<T> take_col (index j);
00365   matrix<T> take_cols (vector<index> indices);
00367   void add_col (index j1, index j2, T alpha);
00369 
00373   void add_cols (index j1, T alpha1, index j2, T alpha2);
00375   void mult_col (index j, T alpha);
00377   void div_col (index j, T alpha);
00379   void swap_cols (index j1, index j2);
00380 
00382   matrix<T> sumbatrix(vector<index> r, vector<index> c);
00383 
00385   matrix operator + () {return *this;}
00387   matrix operator - ();
00388 
00389 #if !defined(__BORLANDC__) || !defined(__WIN32__) 
00390 
00391 
00393   friend ostream & operator << (ostream & s, const matrix & x);
00395 
00397   friend istream & operator >> (istream & s, matrix & x); // ?
00398     // почему-то, если сделать эту функцию-друга обычной функцией,
00399     // команды вида cin >> x не воспринимаются
00400 
00402   friend matrix<T> operator +(const matrix<T> &, const matrix<T> &);
00404   friend matrix<T> operator *(const matrix<T> &, const matrix<T> &);
00406   friend matrix<T> operator *(T alpha, const matrix<T> & a);
00408   friend matrix<T> operator *(const matrix<T> & a, const vector<T> & v);
00410   friend matrix<T> operator *( const vector<T> & v, const matrix<T> & a);
00412   friend matrix<T> transpose (const matrix<T> & a);
00414   friend int operator == (const matrix<T> & a, const matrix<T> & b);
00416   friend int operator != (const matrix<T> & a, const matrix<T> & b);
00417     // почему-то, если сделать эти операции-друзья обычными операциями,
00418     // команды вида a == b * c не воспринимаются
00419 
00421   friend matrix<T> diag(const vector<T>& d);
00423   friend matrix<T> diag(T a, index m, index n);
00425   friend matrix<T> diag(T a, index m);
00427   friend matrix<T> eye (T a, index m, index n);
00429   friend matrix<T> eye (T a, index m);
00431   friend matrix<T> fill(T a, index m, index n);
00433   friend matrix<T> fill(T a, index m);
00435   friend matrix<T> zero(T a, index m, index n);
00437   friend matrix<T> zero(T a, index m);
00439   friend void pretty(ostr & s, const matrix<T>& A);
00440 
00441   friend
00442     void print2(ostr & s,
00443       const matrix<T>& A, const matrix<T>& B);
00444   friend
00445     void print3(ostr & s,
00446       const matrix<T>& A, const matrix<T>& B, const matrix<T>& C);
00447 #else //bc5.5
00448    friend ostream & operator <<<T>(ostream & s, const matrix & x);
00449    friend istream & operator >><T> (istream & s, matrix & x); // ?
00450    friend matrix<T> operator + <T>(const matrix<T> &, const matrix<T> &);
00451    friend matrix<T> operator * <T>(const matrix<T> &, const matrix<T> &);
00452    friend matrix<T> operator * <T>(T alpha, const matrix<T> & a);
00453    friend matrix<T> operator * <T>(const matrix<T> & a, const vector<T> & v);
00454    friend matrix<T> operator * <T>( const vector<T> & v, const matrix<T> & a);
00455    friend matrix<T> transpose<T> (const matrix<T> & a);
00456    friend int operator == <T>(const matrix<T> & a, const matrix<T> & b);
00457    friend int operator != <T>(const matrix<T> & a, const matrix<T> & b);
00458    friend matrix<T> diag<T>(const vector<T>& d);
00459    friend matrix<T> diag<T>(T a, index m, index n);
00460    friend matrix<T> diag<T>(T a, index m);
00461    friend matrix<T> eye <T>(T a, index m, index n);
00462    friend matrix<T> eye <T>(T a, index m);
00463    friend matrix<T> fill<T>(T a, index m, index n);
00464    friend matrix<T> fill<T>(T a, index m);
00465    friend matrix<T> zero<T>(T a, index m, index n);
00466    friend matrix<T> zero<T>(T a, index m);
00467    friend void pretty<T>(ostr & s, const matrix<T>& A);
00468    friend void print2<T>(ostr & s,
00469                 const matrix<T>& A, const matrix<T>& B);
00470    friend void print3<T>(ostr & s,
00471                 const matrix<T>& A, const matrix<T>& B, const matrix<T>& C);
00472 #endif
00473 
00474 private:
00475 
00476   index m, n, m_limit, n_limit, m_delta, n_delta;
00477   T **data;
00478   void dispose_all (void);
00479   void ins_this_row (index i, T * p);
00480   inline int indices_correct (index i, index j) const;
00481   int max_format() const;
00482 };
00483 
00488 template <class T>
00489   void print2(ostr & s,
00490   const matrix<T>& A, const matrix<T>& B);
00491 
00492 
00498 template <class T>
00499   void print3(ostr & s,
00500   const matrix<T>& A, const matrix<T>& B, const matrix<T>& C);
00501 
00502 
00504 
00508 template <class T> void mesh_grid(vector<T>& v, T begin, T end, T step);
00509 
00510 
00511 
00512 /***************************/
00513 /*                         */
00514 /*     implementation      */
00515 /*                         */
00516 /***************************/
00517 
00518 
00519 void matrix_error (int new_status, char *message)
00520 {
00521   if (matrix_status == st_ok)
00522     matrix_status = new_status;
00523   cerr << "matrix error: " << message << "\n";
00524 }
00525 
00526 void matrix_fatal_error (int new_status, char *message)
00527 {
00528   if (matrix_status == st_ok)
00529     matrix_status = new_status;
00530   cerr << "matrix fatal error: " << message << "\n";
00531   exit(1);
00532 }
00533 
00534 inline int sizes_correct (index n_size)
00535 {
00536   return (n_size >= 0) && (n_size <= n_max_size);
00537 }
00538 
00539 inline int sizes_correct (index m_size, index n_size)
00540 {
00541   return (m_size >= 0) && (n_size >= 0) && (m_size <= m_max_size) &&
00542   (n_size <= n_max_size);
00543 }
00544 
00545 /***************************/
00546 /*                         */
00547 /*     class vector<>      */
00548 /*                         */
00549 /***************************/
00550 
00551 template <class T>
00552   vector<T>::vector(index n_init, index n_delta_init)
00553 {
00554   n = n_limit = 0;
00555   n_delta = n_delta_init;
00556   data = NULL;
00557 
00558   resize (n_init);
00559 }
00560 
00561 template <class T>
00562   vector<T>::vector (index n_init, T* a, index n_delta_init)
00563 
00564 {
00565   n = n_limit = 0;
00566   n_delta = n_delta_init;
00567   data = NULL;
00568 
00569   resize (n_init);
00570   if (matrix_status != st_ok)
00571     return;
00572 
00573   for (index i = 0; i < get_n(); i++)
00574     data[i] = a[i];
00575 }
00576 
00577 
00578 template <class T>
00579   vector<T>::vector(const vector<T> &x)
00580 {
00581   n = n_limit = 0;
00582   n_delta = x.n_delta;
00583   data = NULL;
00584   resize (x.get_n());
00585   if (matrix_status == st_ok)
00586     for (index j = 0; j < n; j++)
00587       data[j] = x.data[j];
00588 }
00589 
00590 template <class T>
00591   vector<T>  &vector<T>::operator = (const vector<T>  &x)
00592 {
00593   if (this == &x)
00594     return *this;
00595   dispose_all ();
00596   resize (x.get_n());
00597   if (matrix_status != st_ok)
00598   {
00599     matrix_error(st_out_of_memory, "vector::operator = : can't assign matrix");
00600     return *this;
00601   }
00602   for (index j = 0; j < n; j++)
00603     data[j] = x.data[j];
00604   return *this;
00605 }
00606 
00607 template <class T>
00608   void vector<T> :: dispose_all ()
00609 {
00610   delete [] data;
00611   data = NULL;
00612   n = n_limit = 0;
00613 }
00614 
00615 template <class T>
00616   vector<T> :: ~vector ()
00617 {
00618   dispose_all ();
00619 }
00620 
00621 template <class T>
00622   T & vector<T> :: operator[](index j)
00623 {
00624   if (!index_correct (j))
00625     matrix_fatal_error (st_incorrect_index,
00626       "vector::operator[]: incorrect index");
00627   return data[j];
00628 }
00629 
00630 template <class T>
00631   T vector<T>::at(index j) const
00632 {
00633   if (!index_correct (j))
00634   {
00635     matrix_error (st_incorrect_index,
00636       "vector::at: incorrect index");
00637     return T(0);
00638   }
00639   return data[j];
00640 }
00641 
00642 template <class T>
00643   void vector<T>::resize (index new_n)
00644 {
00645   if (!sizes_correct (new_n + n_delta))
00646   {
00647     matrix_error (st_out_of_range,
00648       "vector::can't resize vector:\n incorrect sizes");
00649     return;
00650   }
00651   T *new_data = new T[new_n + n_delta];
00652   if (!new_data)
00653   {
00654     matrix_error (st_out_of_memory, "vector::can't resize vector: out of memory");
00655     return;
00656   }
00657   dispose_all ();
00658   data = new_data;
00659   n = new_n;
00660   n_limit = new_n + n_delta;
00661   return;
00662 }
00663 
00664 template <class T>
00665 ostream & operator << (ostream & s, const vector<T>  &x)
00666 {
00667     s << x.get_n() << '\n';
00668   for (index j = 0; j < x.get_n(); j++)
00669     s << x.data[j] << '\t';
00670   s << '\n';
00671   if (s.fail ())
00672     matrix_error (st_output_error, "vector::output error");
00673   return s;
00674 }
00675 
00676 template <class T>
00677   istream & operator >> (istream & s, vector<T>  &x)
00678 {
00679   index new_n=0;
00680   s >> new_n;
00681   x.resize (new_n);
00682   if (matrix_status == st_ok)
00683     for (index j = 0; j < x.get_n(); j++)
00684       s >> x[j];
00685   if (s.fail ())
00686     matrix_error (st_input_error, "vector::operator >> : input error");
00687   return s;
00688 }
00689 
00690 template <class T>
00691   void vector<T> :: ins_item (index j, T item)
00692 {
00693   if (!index_correct(j) && j != n)
00694   {
00695     matrix_error(st_out_of_range, "vector::ins_item: incorrect index");
00696     return;
00697   }
00698   if (n < n_limit)
00699   {
00700     for (index jj = n; jj > j; jj--)
00701       data[jj] = data[jj - 1];
00702     n++;
00703     data[j] = item;
00704     return;
00705   }
00706   if (!sizes_correct (n + 1))
00707   {
00708     matrix_error (st_out_of_range, "vector::ins_item: too many items");
00709     return;
00710   }
00711   T *new_data = new T[n_limit + n_delta];
00712   if (!new_data)
00713   {
00714     matrix_error (st_out_of_memory,
00715       "vector::can't insert new item: out of memory \n vector is not correct");
00716     return;
00717   }
00718   index jj;
00719   for (jj = 0; jj < j; jj++)
00720     new_data[jj] = data[jj];
00721   new_data[jj] = item;
00722   for (jj = j; jj < n; jj++)
00723     new_data[jj + 1] = data[jj];
00724   delete [] data;
00725   data = new_data;
00726   n++;
00727   n_limit += n_delta;
00728   return;
00729 }
00730 
00731 template <class T>
00732   void vector<T>::join (T item)
00733 {
00734   ins_item(get_n(), item);
00735 }
00736 
00737 template <class T>
00738   int vector<T>::member (T item) const
00739 {
00740   if (index_of_item(item) == -1)
00741     return 0;
00742   else
00743     return 1;
00744 }
00745 
00746 template <class T>
00747   index vector<T>::index_of_item(T item) const
00748 {
00749   for (index j = 0; j < n; j++)
00750     if (data[j] == item)
00751       return j;
00752   return -1;
00753 }
00754 
00755 template <class T>
00756   void vector<T> :: del_item (index j)
00757 {
00758   if (!index_correct(j))
00759   {
00760     matrix_error(st_out_of_range, "vector::del_item: incorrect index");
00761     return;
00762   }
00763   for (index jj = j + 1; jj < n; jj++)
00764     data[jj - 1] = data[jj];
00765   n--;
00766 }
00767 
00768 template <class T>
00769   T vector<T> :: take_item (index j)
00770 {
00771   T buf;
00772     if (!index_correct(j))
00773   {
00774     matrix_fatal_error(st_out_of_range, "vector::take_item: out of range");
00775     return buf;
00776   }
00777   buf = data[j];
00778   for (index jj = j + 1; jj < n; jj++)
00779     data[jj - 1] = data[jj];
00780   n--;
00781   return buf;
00782 }
00783 
00784 template <class T>
00785   void vector<T> :: swap_items (index j1, index j2)
00786 {
00787   if (!index_correct(j1) || !index_correct(j2))
00788   {
00789     matrix_fatal_error(st_out_of_range, "vector::swap_items: out of range");
00790     return;
00791   }
00792   T buf;
00793   buf = data[j1];
00794   data[j1] = data[j2];
00795   data[j2] = buf;
00796 }
00797 
00798 template <class T>
00799   vector<T> operator + (const vector<T>  &a, const vector<T>  &b)
00800 {
00801     if (a.get_n() != b.get_n())
00802   {
00803     matrix_error (st_inconsistent_sizes,
00804                   "vector::operator+ : different sizes of summand vectors");
00805     return a;
00806   }
00807   vector<T>  sum;
00808   sum.resize (a.get_n());
00809   if (matrix_status != st_ok)
00810     return a;
00811   for (index j = 0; j < a.get_n(); j++)
00812     sum[j] = a.data[j] + b.data[j];
00813   return sum;
00814 }
00815 
00816 template <class T>
00817   vector<T> operator - (const vector<T>  &a, const vector<T>  &b)
00818 {
00819     if (a.get_n() != b.get_n())
00820   {
00821     matrix_error (st_inconsistent_sizes,
00822                   "vector::operator- :different sizes of vectors");
00823     return a;
00824   }
00825   vector<T>  dif;
00826   dif.resize (a.get_n());
00827   if (matrix_status != st_ok)
00828     return a;
00829   for (index j = 0; j < a.get_n(); j++)
00830     dif[j] = a.data[j] - b.data[j];
00831   return dif;
00832 }
00833 
00834 template <class T>
00835 vector<T> & vector<T>::operator += (const vector<T> & x)
00836 {
00837   if (n != x.get_n())
00838   {
00839     matrix_error (st_inconsistent_sizes,
00840                   "vector::operator+= :different sizes of vectors");
00841     return *this;
00842   }
00843 
00844   for (index j = 0; j < n; j++)
00845     data[j] = data[j] + x.data[j];
00846 
00847   return *this;
00848 }
00849 
00850 template <class T>
00851   vector<T> & vector<T>::operator -= (const vector<T> & x)
00852 {
00853   if (n != x.get_n())
00854   {
00855     matrix_error (st_inconsistent_sizes,
00856                   "vector::operator-= :different sizes of vectors");
00857     return *this;
00858   }
00859 
00860   for (index j = 0; j < n; j++)
00861     data[j] = data[j] - x.data[j];
00862 
00863   return *this;
00864 }
00865 
00866 template <class T>
00867   vector<T> operator * (T alpha, const vector<T> & a)
00868 {
00869   vector<T>  result;
00870   result.resize (a.get_n());
00871   if (matrix_status != st_ok)
00872     return a;
00873   for (index j = 0; j < a.get_n(); j++)
00874     result[j] = alpha * a.data[j];
00875   return result;
00876 }
00877 
00878 template <class T>
00879   T operator * (const vector<T> & a, const vector<T> & b)
00880 {
00881   T prod = 0;
00882   for (index j = 0; j < a.get_n(); j++)
00883     prod = prod + a.data[j] * b.data[j];
00884   return prod;
00885 }
00886 
00887 template <class T>
00888   int operator == (const vector<T>  &a, const vector<T>  &b)
00889 {
00890 if (a.get_n() != b.get_n())
00891   return 0;
00892 for (index j = 0; j < a.get_n(); j++)
00893   if (a.data[j] != b.data[j])
00894     return 0;
00895 return 1;
00896 }
00897 
00898 template <class T>
00899   int  operator != (const vector<T>  &a, const vector<T>  &b)
00900 {
00901   return !(a == b);
00902 }
00903 
00904 template <class T>
00905   int vector<T>::index_correct (index j) const
00906 {
00907   return ((j >= 0) && (j < n));
00908 }
00909 
00910 template <class T>
00911   T vector<T>::access (index j) const
00912 {
00913   return data[j];
00914 }
00915 
00916 /***************************/
00917 /*                         */
00918 /*     class matrix<>      */
00919 /*                         */
00920 /***************************/
00921 
00922 template <class T>
00923   matrix<T>::matrix (index m_init,      index n_init,
00924   index m_delta_init, index n_delta_init)
00925 {
00926   m = n = m_limit = n_limit = 0;
00927   m_delta = m_delta_init;
00928   n_delta = n_delta_init;
00929   data = NULL;
00930 
00931   resize (m_init, n_init);
00932 }
00933 
00934 template <class T>
00935   matrix<T>::matrix (index m_init, index n_init, T* A,
00936   index m_delta_init, index n_delta_init)
00937 {
00938   m = n = m_limit = n_limit = 0;
00939   m_delta = m_delta_init;
00940   n_delta = n_delta_init;
00941   data = NULL;
00942 
00943   resize (m_init, n_init);
00944   if (matrix_status != st_ok)
00945     return;
00946 
00947   for (index i = 0; i < m; i++)
00948     for (index j = 0; j < n; j++)
00949       data[i][j] = A[i*n + j];
00950 }
00951 
00952 
00953 template <class T>
00954   matrix<T> ::matrix (const matrix<T>  &x)
00955 {
00956   n = n_limit = m = m_limit = 0;
00957   n_delta = x.n_delta;
00958   m_delta = x.m_delta;
00959   data = NULL;
00960   resize (x.get_m(), x.get_n());
00961   if (matrix_status != st_ok)
00962     return;
00963   for (index i = 0; i < m; i++)
00964     for (index j = 0; j < n; j++)
00965       data[i][j] = x.access(i, j);
00966 }
00967 
00968 template <class T>
00969   matrix<T> & matrix<T> :: operator = (const matrix<T>  &x)
00970 {
00971   if (this == &x)
00972     return *this;
00973   dispose_all ();
00974   resize (x.get_m(), x.get_n());
00975   if (matrix_status != st_ok)
00976   {
00977     matrix_error(st_out_of_range, "matrix::operator= :can't assign matrix");
00978     return *this;
00979   }
00980   for (index i = 0; i < m; i++)
00981     for (index j = 0; j < n; j++)
00982       data[i][j] = x.access(i, j);
00983   return *this;
00984 }
00985 
00986 template <class T>
00987   void matrix<T> ::dispose_all ()
00988 {
00989   if (!data)
00990     return;
00991   for (index i = 0; i < m; i++)
00992     delete [] data[i];
00993   delete [] data;
00994   data = NULL;
00995   m = n = m_limit = n_limit = 0;
00996 }
00997 
00998 template <class T>
00999   matrix<T> ::~matrix ()
01000 {
01001   dispose_all ();
01002 }
01003 
01004 template <class T>
01005   T & matrix<T>::operator ()(index i, index j)
01006 {
01007   if (!indices_correct (i, j))
01008     matrix_fatal_error (st_incorrect_index, "matrix::operator() :incorrect indices");
01009   return data[i][j];
01010 }
01011 
01012 template <class T>
01013   T matrix<T>::at(index i, index j) const
01014 {
01015   if (!indices_correct (i, j))
01016   {
01017     matrix_error (st_incorrect_index, "matrix::at: incorrect indices");
01018     return T(0);
01019   }
01020   return data[i][j];
01021 }
01022 
01023 template <class T>
01024   void matrix<T>::resize (index new_m, index new_n)
01025 {
01026 
01027   if (!sizes_correct(new_m + m_delta, new_n + n_delta))
01028   {
01029     matrix_error (st_out_of_range, "can't resize matrix:\n incorrect sizes");
01030     return;
01031   }
01032   T **new_data = new T *[new_m + m_delta];
01033   if (!new_data)
01034   {
01035     matrix_error (st_out_of_memory, "can't resize matrix: out of memory");
01036     return;
01037   }
01038   dispose_all ();
01039   data = new_data;
01040   /* now we are inserting rows */
01041   index m_count;
01042   for (m_count = 0; m_count < new_m; m_count++)
01043   {
01044     T *p = new T[new_n + n_delta];
01045     if (!p)
01046     {
01047       matrix_error (st_out_of_memory,
01048                     "can't correctly resize matrix: out of memory");
01049       break;
01050     }
01051     data[m_count] = p;
01052   }
01053   n = new_n;
01054   n_limit = new_n + n_delta;
01055   m = m_count;
01056   m_limit = new_m + m_delta;
01057 }
01058 
01059 template <class T> int matrix<T>::max_format() const
01060 {
01061   int format = 0;
01062   for (int i = 0; i < m; i++)
01063     for (int j = 0; j < n; j++)
01064     {
01065       ostrstream s;
01066 #ifdef __WATCOMC__
01067 #  pragma warn -665
01068       (ostream)s << at(i, j) << ends; // I guess Watcom C bug
01069 #  pragma warn +665
01070 #else
01071       s << at(i, j) << ends;
01072 #endif
01073       char *str = s.str();
01074       if (strlen(str) > format)
01075         format = strlen(str);
01076       delete [] str;
01077     }
01078   return format;
01079 }
01080 
01081 
01082 template <class T>
01083 ostream & operator << (ostream & s, const matrix<T> &x)
01084 {
01085   s << x.get_m() << ' ' << x.get_n() << '\n';
01086   for (index i = 0; i < x.get_m(); i++)
01087   {
01088     for (index j = 0; j < x.get_n(); j++)
01089       s << x.access(i, j) << '\t';
01090     s << '\n';
01091   }
01092   if (s.fail ())
01093     matrix_error (st_output_error, "matrix::output error");
01094   return s;
01095 }
01096 
01097 template <class T>
01098 istream & operator >> (istream & s, matrix<T>& x)
01099 {
01100   index new_m = 0, new_n = 0;
01101   s >> new_m >> new_n;
01102   x.resize (new_m, new_n);
01103   if (matrix_status == st_ok)
01104     for (index i = 0; i < x.get_m(); i++)
01105       for (index j = 0; j < x.get_n(); j++)
01106         s >> x (i, j);
01107   if (s.fail ())
01108     matrix_error (st_input_error, "matrix::operator >> : input error");
01109   return s;
01110 }
01111 
01112 template <class T>
01113   vector<T> matrix<T>::row (index i) const
01114 {
01115   vector<T>  x;
01116   if (!indices_correct (i, -1))
01117   {
01118     matrix_error (st_incorrect_index, "matrix::row: incorrect index");
01119     return x;
01120   }
01121   x.resize(n);
01122   if (matrix_status == st_ok)
01123     for (index j = 0; j < n; j++)
01124       x[j] = data[i][j];
01125   return x;
01126 }
01127 
01128 template <class T>
01129   matrix<T> matrix<T>::rows (vector<index> indices) const
01130 {
01131   matrix<T> result(indices.get_n(), n);
01132   if (matrix_status == st_ok)
01133     for (index i = 0; i < indices.get_n(); i++)
01134       for (index j = 0; j < n; j++)
01135         result(i, j) = at(indices[i], j);
01136   return result;
01137 }
01138 
01139 template <class T>
01140   void matrix<T> ::ins_row (index i)
01141 {
01142   if (!sizes_correct (m + 1, n))
01143   {
01144     matrix_error (st_out_of_range, "matrix::ins_row: matrix is too big");
01145     return;
01146   }
01147   T *p = new T[n_limit];
01148   //storage for a  new row
01149     if (!p)
01150     {
01151       matrix_error (st_out_of_memory, "can't insert new row: out of memory");
01152       return;
01153     }
01154   ins_this_row (i, p);
01155   if (matrix_status != st_ok)
01156   {
01157     delete [] p;
01158     return;
01159   }
01160   return;
01161 }
01162 
01163 template <class T>
01164   void matrix<T> :: ins_row (index i, vector<T>  x)
01165 {
01166   if (n != x.get_n ())
01167   {
01168     matrix_error (st_inconsistent_sizes,
01169                   "matrix::ins_row: can't insert row: sizes are inconsistent");
01170     return;
01171   }
01172   ins_row(i);
01173   if (matrix_status != st_ok)
01174     return;
01175   else
01176   {
01177     for (index j = 0; j < n; j++)
01178       data[i][j] = x[j];
01179     return;
01180   }
01181 }
01182 
01183 template <class T>
01184   void matrix<T> :: join_row (vector<T>  x)
01185 {
01186   ins_row(m, x);
01187 }
01188 
01189 template <class T>
01190   void matrix<T>::ins_this_row (index i, T * p)
01191 {
01192   if (!indices_correct (i, -1) && i != m)
01193   {
01194     matrix_error (st_incorrect_index, "matrix::ins_this_row: incorrect index");
01195     return;
01196   }
01197   if (m < m_limit)
01198   {
01199     memmove (&(data[i + 1]), &(data[i]), sizeof (T *) * (m - i));
01200     m++;
01201     data[i] = p;
01202   }
01203   else
01204   {
01205     index new_m = m_limit + m_delta;
01206     T **new_data = new T *[new_m];
01207     if (!new_data)
01208     {
01209       matrix_error (st_out_of_memory, "can't insert new row: out of memory");
01210       return;
01211     }
01212     memmove (new_data, data, sizeof (T *) * i);
01213     memmove (&(new_data[i + 1]), &(data[i]), sizeof (T *) * (m - i));
01214     delete [] data;
01215     data = new_data;
01216     m++;
01217     m_limit += m_delta;
01218     data[i] = p;
01219   }
01220   return;
01221 }
01222 
01223 template <class T>
01224   void matrix<T> :: del_row (index i)
01225 {
01226   take_row(i);
01227 }
01228 
01229 template <class T>
01230   vector<T> matrix<T> :: take_row (index i)
01231 {
01232   vector<T>  x;
01233     x = row (i);
01234   if (matrix_status != st_ok)
01235     return x;
01236   if (data)
01237     delete [] data[i];
01238   memmove (&(data[i]), &(data[i + 1]), sizeof (T *) * (m - i - 1));
01239   m--;
01240   return x;
01241 }
01242 
01243 template <class T>
01244   void matrix<T>::del_rows (vector<index> indices)
01245 {
01246   // first we dispose rows
01247   for (index i = 0; i < indices.get_n(); i++)
01248   {
01249     if (!indices_correct (indices[i], -1))
01250     {
01251       matrix_error (st_incorrect_index, "matrix::del_rows: incorrect index");
01252       continue;
01253     }
01254     delete [] data[indices[i]];
01255     data[indices[i]] = NULL;
01256   }
01257 
01258   // second we del null rows
01259   index new_i = 0;
01260   for (index i = 0; i < m; i++)
01261   {
01262     if (!indices_correct (i, -1))
01263       continue;
01264     if (data[i])
01265     {
01266       data[new_i] = data[i];
01267       new_i++;
01268     }
01269   }
01270   m = new_i;
01271 
01272 }
01273 
01274 template <class T>
01275   matrix<T> matrix<T>::take_rows (vector<index> indices)
01276 {
01277   matrix<T> result(0, n);
01278   if (matrix_status != st_ok)
01279     return result;
01280 
01281   for (index i = 0; i < indices.get_n(); i++)
01282   {
01283     if (!indices_correct (indices[i], -1))
01284     {
01285       matrix_error (st_incorrect_index, "matrix::take_rows: incorrect index");
01286       continue;
01287     }
01288     if (data[indices[i]])
01289     {
01290       T* current_row = data[indices[i]];
01291       result.ins_this_row(i, current_row);
01292       data[indices[i]] = NULL;
01293     }
01294     else
01295     {
01296       // the row indices[i] has been just taken !
01297       index old_copy = indices.index_of_item(indices[i]);
01298       result.join_row(result.row(old_copy));
01299     }
01300   }
01301 
01302   // delete null rows
01303   index new_i = 0;
01304   for (index i = 0; i < m; i++)
01305   {
01306     if (data[i])
01307     {
01308       data[new_i] = data[i];
01309       new_i++;
01310     }
01311   }
01312   m = new_i;
01313 
01314   return result;
01315 }
01316 
01317 
01318 
01319 template <class T>
01320   void matrix<T> :: add_row (index i1, index i2, T alpha)
01321 {
01322   if (indices_correct (i1, -1) && indices_correct (i2, -1))
01323     for (index j = 0; j < n; j++)
01324       data[i1][j] = data[i1][j] + data[i2][j] * alpha;
01325   else
01326     matrix_error (st_incorrect_index, "matrix::can't add rows: incorrect indices");
01327 }
01328 
01329 template <class T>
01330  void matrix<T> ::add_rows (index i1, T alpha1, index i2, T alpha2)
01331 {
01332   if (!indices_correct (i1, -1) || !indices_correct (i2, -1))
01333   {
01334     matrix_error (st_incorrect_index, "matrix::can't add rows: incorrect indices");
01335     return;
01336   }
01337   ins_row(m);
01338   if (matrix_status == st_ok)
01339   {
01340     for (index j = 0; j < n; j++)
01341       data[m - 1][j] =
01342         data[i1][j] * alpha1 + data[i2][j] * alpha2;
01343     return;
01344   }
01345   else
01346     return;
01347 }
01348 
01349 template <class T>
01350   void matrix<T> ::mult_row (index i, T alpha)
01351 {
01352   if (indices_correct (i, -1))
01353     for (index j = 0; j < n; j++)
01354       data[i][j] = data[i][j] * alpha;
01355   else
01356     matrix_error (st_incorrect_index, "matrix::can't mult row: incorrect index");
01357 }
01358 
01359 template <class T>
01360   void matrix<T> ::div_row (index i, T alpha)
01361 {
01362   if (indices_correct (i, -1))
01363     for (index j = 0; j < n; j++)
01364       data[i][j] = data[i][j] / alpha;
01365   else
01366     matrix_error (st_incorrect_index, "matrix::can't div row: incorrect index");
01367 }
01368 
01369 template <class T>
01370   void matrix<T> ::swap_rows (index i1, index i2)
01371 {
01372   if (!indices_correct (i1, -1) || !indices_correct (i2, -1))
01373   {
01374     matrix_error (st_incorrect_index, "matrix::can't swap rows: incorrect indices");
01375     return;
01376   }
01377   T * buf = data[i1];
01378   data[i1] = data[i2];
01379   data[i2] = buf;
01380 }
01381 
01382 template <class T>
01383   vector<T> matrix<T>::col(index j) const
01384 {
01385   vector<T>  x;
01386   if (!indices_correct (-1, j))
01387   {
01388     matrix_error (st_incorrect_index, "matrix::col: incorrect index");
01389     return x;
01390   }
01391   x.resize (m);
01392   if (matrix_status == st_ok)
01393     for (index i = 0; i < m; i++)
01394       x[i] = data[i][j];
01395   return x;
01396 }
01397 
01398 template <class T>
01399   matrix<T> matrix<T>::cols (vector<index> indices) const
01400 {
01401   matrix<T> result(m, indices.get_n());
01402   if (matrix_status == st_ok)
01403     for (index j = 0; j < indices.get_n(); j++)
01404       for (index i = 0; i < m; i++)
01405         result.data[i][j] = at(i, indices[j]);
01406   return result;
01407 }
01408 
01409 template <class T>
01410  void matrix<T> ::ins_col (index j)
01411 {
01412   if (!indices_correct(-1, j) && j !=n)
01413   {
01414     matrix_error (st_incorrect_index, "matrix::can't ins col: incorrect index");
01415     return;
01416   }
01417   if (n < n_limit)
01418   {
01419     for (index ii = 0; ii < m; ii++)
01420       for (index jj = n; jj > j; jj--)
01421         data[ii][jj] = data[ii][jj - 1];
01422     n++;
01423     return;
01424   }
01425 
01426   if (!sizes_correct (m, n + 1))
01427   {
01428     matrix_error (st_out_of_range, "matrix::ins_col: too many columns");
01429     return;
01430   }
01431 
01432   for (index i = 0; i < m; i++)
01433   {
01434     T *p = new T[n_limit + n_delta];
01435     if (!p)
01436     {
01437       matrix_error (st_out_of_memory,
01438          "matrix::can't insert new column: out of memory \n matrix is not correct");
01439       return;
01440     }
01441     index jj;
01442     for (jj = 0; jj < j; jj++)
01443       p[jj] = data[i][jj];
01444     for (jj = j; jj < n; jj++)
01445       p[jj + 1] = data[i][jj];
01446     delete [] data[i];
01447     data[i] = p;
01448   }
01449   n++;
01450   n_limit += n_delta;
01451 
01452   return;
01453 }
01454 
01455 template <class T>
01456  void matrix<T> ::ins_col (index j, vector<T>  x)
01457 {
01458   if (m != x.get_n ())
01459   {
01460     matrix_error (st_inconsistent_sizes,
01461       "matrix::can't insert column: sizes are inconsistent");
01462     return;
01463   }
01464   ins_col (j);
01465   if (matrix_status != st_ok)
01466     return;
01467   else
01468   {
01469     for (index i = 0; i < m; i++)
01470       data[i][j] = x[i];
01471     return;
01472   }
01473 }
01474 
01475 template <class T>
01476  void matrix<T> :: join_col (vector<T>  x)
01477 {
01478   ins_col(n, x);
01479 }
01480 
01481 template <class T>
01482   void matrix<T>::del_col (index j)
01483 {
01484   take_col(j);
01485 }
01486 
01487 template <class T>
01488   vector<T> matrix<T>::take_col (index j)
01489 {
01490   vector<T> x;
01491   x = col (j);
01492   if (matrix_status != st_ok)
01493     return x;
01494   for (index i = 0; i < m; i++)
01495     for (index jj = j + 1; jj < n; jj++)
01496       data[i][jj - 1] = data[i][jj];
01497   n--;
01498   return x;
01499 }
01500 
01501 template <class T>
01502   void matrix<T>::del_cols (vector<index> indices)
01503 {
01504   vector<int> disposed_cols(n);
01505 
01506   if (matrix_status != st_ok)
01507     return;
01508 
01509 
01510   // first we remember cols to delete
01511   for (index j = 0; j < n; j++)
01512     disposed_cols[j] = 0;
01513   for (index j = 0; j < indices.get_n(); j++)
01514     disposed_cols[indices[j]] = 1;
01515 
01516   // second we delete cols
01517   index new_j = 0;
01518   for (index j = 0; j < n; j++)
01519   {
01520     if (disposed_cols[j])
01521     {
01522       for (index i = 0; i < m; i++)
01523         data[i][new_j] = data[i][j];
01524       new_j++;
01525     }
01526   }
01527   n = new_j;
01528 }
01529 
01530 template <class T>
01531   matrix<T> matrix<T>::take_cols (vector<index> indices)
01532 {
01533   matrix<T> result = cols(indices);
01534   del_cols(indices);
01535   return result;
01536 }
01537 
01538 
01539 template <class T>
01540   void matrix<T> ::add_col (index j1, index j2, T alpha)
01541 {
01542   if (!indices_correct (-1, j1) || !indices_correct (-1, j2))
01543   {
01544     matrix_error (st_incorrect_index, "matrix::can't add cols: incorrect indices");
01545     return;
01546   }
01547   for (index i = 0; i < m; i++)
01548     data[i][j1] = data[i][j1] + data[i][j2] * alpha;
01549 }
01550 
01551 template <class T>
01552  void matrix<T> ::add_cols (index j1, T alpha1, index j2, T alpha2)
01553 {
01554   if (!indices_correct (-1, j1) || !indices_correct (-1, j2))
01555   {
01556     matrix_error (st_incorrect_index,
01557                   "matrix::can't add columns: incorrect indices");
01558     return;
01559   }
01560   ins_col (n);
01561   if (matrix_status == st_ok)
01562   {
01563     for (index i = 0; i < m; i++)
01564       data[i][n - 1] =
01565         data[i][j1] * alpha1 + data[i][j2] * alpha2;
01566     return;
01567   }
01568   else
01569     return;
01570 }
01571 
01572 template <class T>
01573   void matrix<T> ::mult_col (index j, T alpha)
01574 {
01575   if (!indices_correct (-1, j))
01576   {
01577     matrix_error (st_incorrect_index, "matrix::can't mult col: incorrect indices");
01578     return;
01579   }
01580   for (index i = 0; i < m; i++)
01581     data[i][j] = data[i][j] * alpha;
01582 }
01583 
01584 template <class T>
01585   void matrix<T> ::div_col (index j, T alpha)
01586 {
01587   if (!indices_correct (-1, j))
01588   {
01589     matrix_error (st_incorrect_index, "matrix::can't div col: incorrect indices");
01590     return;
01591   }
01592     for (index i = 0; i < m; i++)
01593     data[i][j] = data[i][j] / alpha;
01594 }
01595 
01596 template <class T>
01597   void matrix<T> ::swap_cols (index j1, index j2)
01598 {
01599   if (!indices_correct (-1, j1) || !indices_correct (-1, j2))
01600   {
01601     matrix_error (st_incorrect_index,
01602                   "matrix::can't swap columns: incorrect indices");
01603     return;
01604   }
01605   T buf;
01606   for (index i = 0; i < m; i++)
01607   {
01608     buf = data[i][j1];
01609     data[i][j1] = data[i][j2];
01610     data[i][j2] = buf;
01611   }
01612 }
01613 
01614 template <class T>
01615   matrix<T> matrix<T>::sumbatrix(vector<index> r, vector<index> c)
01616 {
01617   matrix<T> result(r.get_n(), c.get_n());
01618   if (matrix_status != st_ok)
01619     return result;
01620   for (index i = 0; i < r.get_n(); i++)
01621     if (!indices_correct(r[i], -1))
01622     {
01623       matrix_fatal_error (st_incorrect_index, "matrix::submatrix: incorrect index");
01624       return result;
01625     }
01626   for (index j = 0; j < c.get_n(); j++)
01627     if (!indices_correct(c[j], -1))
01628     {
01629       matrix_fatal_error (st_incorrect_index, "matrix::submatrix: incorrect index");
01630       return result;
01631     }
01632   for (index i = 0; i < r.get_n(); i++)
01633     for (index j = 0; j < c.get_n(); j++)
01634       result.data[i][j] = data[r[i]][c[j]];
01635   return result;
01636 }
01637 
01638 template <class T>
01639   matrix<T>  matrix<T> ::operator - ()
01640 {
01641   matrix<T>  x (m, n);
01642   for (index i = 0; i < m; i++)
01643     for (index j = 0; j < n; j++)
01644 #ifdef __BORLANDC__
01645       x (i, j) = 0-data[i][j];
01646 #else
01647       x (i, j) = -data[i][j];
01648 #endif
01649   return x;
01650 }
01651 
01652 template <class T>
01653   matrix<T> operator + (const matrix<T>  &a, const matrix<T>  &b)
01654 {
01655     if ((a.get_n() != b.get_n()) || (a.get_m() != b.get_m()))
01656   {
01657     matrix_error (st_inconsistent_sizes,
01658                   "matrix::operator+ : different sizes of summand matrices");
01659     return a;
01660   }
01661   matrix<T>  sum;
01662   sum.resize (a.get_m(), a.get_n());
01663   if (matrix_status != st_ok)
01664     return a;
01665   for (index i = 0; i < a.get_m(); i++)
01666     for (index j = 0; j < a.get_n(); j++)
01667       sum (i, j) = a.access(i, j) + b.access(i, j);
01668   return sum;
01669 }
01670 
01671 template <class T>
01672   matrix<T> operator * (const matrix<T>  &a, const matrix<T>  &b)
01673 {
01674     if (a.get_n() != b.get_m())
01675   {
01676     matrix_error (st_inconsistent_sizes,
01677       "operator * : inconsistent sizes of prodotted matrices");
01678     return a;
01679   }
01680   matrix<T>  result;
01681   result.resize (a.get_m(), b.get_n());
01682   if (matrix_status != st_ok)
01683     return a;
01684   for (index i = 0; i < result.get_m(); i++)
01685     for (index j = 0; j < result.get_n(); j++)
01686     {
01687       result (i, j) = 0;
01688       for (index k = 0; k < a.get_n(); k++)
01689         result (i, j) = result(i, j) + a.access(i, k) * b.access(k, j);
01690     }
01691   return result;
01692 }
01693 
01694 
01695 /*
01696 template <class T>
01697   matrix<T> operator *(const matrix<T> & a, const vector<T> & v)
01698 {
01699     if (a.get_n() != v.get_n())
01700   {
01701     matrix_error (st_inconsistent_sizes,
01702       "operator * : inconsistent sizes of prodotted matrices");
01703     return a;
01704   }
01705   matrix<T> result;
01706   result.resize (a.get_m(), v.get_n());
01707   if (matrix_status != st_ok)
01708     return a;
01709   for (index i = 0; i < result.get_m(); i++)
01710     for (index j = 0; j < result.get_n(); j++)
01711     {
01712       result (i, j) = 0;
01713       for (index k = 0; k < a.get_n(); k++)
01714         result (i, j) = result(i, j) + a.data[i][k] * b.data[k][j];
01715     }
01716   return result;
01717 }
01718 */
01719 
01720 //friend matrix<T> operator *(const vector<T> & v, const matrix<T> & a)
01721 
01722 
01723 template <class T>
01724   matrix<T> operator * (T alpha, const matrix<T>  &a)
01725 {
01726   matrix<T>  result;
01727   result.resize (a.get_m(), a.get_n());
01728   if (matrix_status != st_ok)
01729     return a;
01730   for (index i = 0; i < a.get_m(); i++)
01731     for (index j = 0; j < a.get_n(); j++)
01732       result (i, j) = alpha * a.data[i][j];
01733   return result;
01734 }
01735 
01736 template <class T>
01737   matrix<T> transpose (const matrix<T>  &a)
01738 {
01739   matrix<T>  result;
01740   result.resize (a.get_n(), a.get_m());
01741   if (matrix_status != st_ok)
01742     return a;
01743   for (index i = 0; i < a.get_n(); i++)
01744     for (index j = 0; j < a.get_m(); j++)
01745       result (i, j) = a.data[j][i];
01746   return result;
01747 }
01748 
01749 template <class T>
01750   int  operator == (const matrix<T>  &a, const matrix<T>  &b)
01751 {
01752     if ((a.get_n() != b.get_n()) || (a.get_m() != b.get_m()))
01753     return 0;
01754   for (index i = 0; i < a.get_m(); i++)
01755     for (index j = 0; j < a.get_n(); j++)
01756       if (a.data[i][j] != b.data[i][j])
01757         return 0;
01758   return 1;
01759 }
01760 
01761 template <class T>
01762   int  operator != (const matrix<T>  &a, const matrix<T>  &b)
01763 {
01764   return !(a == b);
01765 }
01766 
01767 template <class T>
01768   int matrix<T> :: indices_correct (index i, index j) const
01769 {
01770   return (
01771        i == -1 && j >= 0 && j < n
01772     || j == -1 && i >= 0 && i < m
01773     || i >= 0 && i < m && j >= 0 && j < n);
01774 }
01775 
01776 template <class T>
01777   T matrix<T> :: access (index i, index j) const
01778 {
01779   return data[i][j];
01780 }
01781 
01782 template <class T>
01783   matrix<T> diag(const vector<T>& d)
01784 {
01785   index n = d.get_n();
01786   matrix<T> D(n, n);
01787   if (matrix_status != st_ok)
01788     return D;
01789   for (index i = 0; i < n; i++)
01790     for (index j = 0; j < n; j++)
01791       if (i != j)
01792         D.data[i][j] = 0;
01793       else
01794         D.data[i][j] = d.at(i);
01795   return D;
01796 }
01797 
01798 template <class T>
01799   matrix<T> diag(T a, index m, index n)
01800 {
01801   matrix<T> D(m, n);
01802   if (matrix_status != st_ok)
01803     return D;
01804   for (index i = 0; i < m; i++)
01805     for (index j = 0; j < n; j++)
01806       if (i != j)
01807         D.data[i][j] = T(0);
01808       else
01809         D.data[i][j] = a;
01810   return D;
01811 }
01812 
01813 template <class T>
01814   matrix<T> diag(T a, index m)
01815 {
01816   return diag(a, m, m);
01817 }
01818 
01819 template <class T>
01820   matrix<T> eye (T a, index m, index n)
01821 {
01822   return diag(T(1), m, n);
01823 }
01824 
01825 template <class T>
01826   matrix<T> eye(T a, index m)
01827 {
01828   return diag(T(1), m);
01829 }
01830 
01831 template <class T>
01832   matrix<T> fill(T a, index m, index n)
01833 {
01834   matrix<T> D(m, n);
01835   if (matrix_status != st_ok)
01836     return D;
01837   for (index i = 0; i < m; i++)
01838     for (index j = 0; j < n; j++)
01839       D.data[i][j] = a;
01840   return D;
01841 }
01842 
01843 template <class T>
01844   matrix<T> fill(T a, index m)
01845 {
01846   return fill(a, m, m);
01847 }
01848 
01849 template <class T>
01850   matrix<T> zero(T a, index m, index n)
01851 {
01852   return fill(T(0), m, n);
01853 }
01854 
01855 template <class T>
01856   matrix<T> zero(T a, index m)
01857 {
01858   return fill(T(0), m);
01859 }
01860 
01861 template <class T>
01862   void mesh_grid(vector<T>& v, T begin, T end, T step)
01863 {
01864   index size = (end - begin)/step + 1;
01865   v.resize(size);
01866   for (index j = 0; j < size; j++)
01867   {
01868     v[j] = begin;
01869     begin = begin + step;
01870   }
01871 }
01872 
01873 template <class T>
01874   void pretty(ostr & s, const vector<T>& v)
01875 {
01876   s << "\\left(";
01877   for (index j = 0; j < v.get_n(); j++)
01878   {
01879     if (j > 0)
01880       s << ", & ";
01881     s << v.at(j);
01882   }
01883   s << "\\right) ^{\rm T} \n";
01884 }
01885 
01886 template <class T>
01887   void pretty(ostr & s, const matrix<T>& A)
01888 {
01889   s << "\\left(\n\\begin{array}{";
01890   for (index j = 0; j < A.get_n(); j++)
01891     s << "r";
01892   s << "} \n";
01893   for (index i = 0; i < A.get_m(); i++)
01894   {
01895     for (index j = 0; j < A.get_n(); j++)
01896     {
01897       if (j > 0)
01898         s << "\t & ";
01899       s << A.at(i, j);
01900     }
01901     s << " \\\\ \n";
01902   }
01903   s << "\\end{array}\n\\right) \n";
01904 }
01905 
01906 
01911 template <class T>
01912   void print2(ostr & s,
01913   const matrix<T>& A, const matrix<T>& B)
01914 {
01915   int format = A.max_format();
01916   if (format < B.max_format())
01917     format = B.max_format();
01918 
01919   s << " ┌";
01920   for (int j = 0; j < (format + 1) * A.get_n(); j++)
01921     s << "─";
01922   s << "─┬";
01923   for (int j = 0; j < (format + 1) * B.get_n(); j++)
01924     s << "─";
01925   s <<  "─┐\n";
01926 
01927   for (int i = 0; i < A.get_m(); i++)
01928   {
01929     s << " │";
01930     for (int j = 0; j < A.get_n(); j++)
01931     {
01932       s << " ";
01933       s.width(format);
01934       s << A.at(i, j);
01935     }
01936     s << " │";
01937     for (int j = 0; j < B.get_n(); j++)
01938     {
01939       s << " ";
01940       s.width(format);
01941       s << B.at(i, j);
01942     }
01943     s << " │\n";
01944   };
01945 
01946   s << " └";
01947   for (int j = 0; j < (format + 1) * A.get_n(); j++)
01948     s << "─";
01949   s << "─┴";
01950   for (int j = 0; j < (format + 1) * B.get_n(); j++)
01951     s << "─";
01952   s << "─┘\n";
01953 
01954 };
01955 
01961 template <class T>
01962   void print3(ostr & s,
01963   const matrix<T>& A, const matrix<T>& B, const matrix<T>& C)
01964 {
01965   int format = A.max_format();
01966   if (format < B.max_format())
01967     format = B.max_format();
01968   if (format < C.max_format())
01969     format = C.max_format();
01970 
01971   s << "  ";
01972   for (int j = 0; j < (format + 1) * B.get_n(); j++)
01973     s << " ";
01974   s << " ┌";
01975   for (int j = 0; j < (format + 1) * A.get_n(); j++)
01976     s << "─";
01977   s << "─┐\n";
01978 
01979   for (int i = 0; i < A.get_m(); i++)
01980   {
01981     s << "  ";
01982     for (int j = 0; j < (format + 1) * B.get_n(); j++)
01983       s << " ";
01984     s << " │";
01985     for (int j = 0; j < A.get_n(); j++)
01986     {
01987       s.width(format);
01988       s << A.at(i, j) << " ";
01989     }
01990     s << " │\n";
01991   };
01992 
01993   s << " ┌";
01994   for (int j = 0; j < (format + 1) * B.get_n(); j++)
01995     s << "─";
01996   s << "─┼";
01997   for (int j = 0; j < (format + 1) * A.get_n(); j++)
01998     s << "─";
01999   s <<  "─┤\n";
02000 
02001   for (int i = 0; i < B.get_m(); i++)
02002   {
02003     s << " │";
02004     for (int j = 0; j < B.get_n(); j++)
02005     {
02006       s.width(format);
02007       s << B.at(i, j) << " ";
02008     }
02009     s << " │";
02010     for (int j = 0; j < C.get_n(); j++)
02011     {
02012       s.width(format);
02013       s << C.at(i, j) << " ";
02014     }
02015     s << " │\n";
02016   };
02017 
02018   s << " └";
02019   for (int j = 0; j < (format + 1) * B.get_n(); j++)
02020     s << "─";
02021   s << "─┴";
02022   for (int j = 0; j < (format + 1) * A.get_n(); j++)
02023     s << "─";
02024   s << "─┘\n";
02025 
02026 };
02027 
02028 
02029 #endif
02030 

Generated at Tue Jan 22 20:37:04 2002 for Arageli by doxygen1.2.9.1 written by Dimitri van Heesch, © 1997-2001