00001
00002
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;
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
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
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
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
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
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
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
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
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
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
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;
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
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
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
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
01297 index old_copy = indices.index_of_item(indices[i]);
01298 result.join_row(result.row(old_copy));
01299 }
01300 }
01301
01302
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
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
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
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
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