home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 3 / 3257 / Matrix.C next >
Encoding:
C/C++ Source or Header  |  1991-04-30  |  6.6 KB  |  358 lines
  1. /*
  2.  * C++ class for MxN matrices.
  3.  *
  4.  * $Id: Matrix.C,v 1.6 91/04/14 18:27:48 bjaspan Exp $
  5.  */
  6.  
  7.  
  8. /*
  9.  * (1) An m x n matrix has m rows and n columns -- ie height x width.
  10.  * 
  11.  * (2) nam == 1 --> nothing else is guaranteed (some code here might
  12.  * violate that.. sigh)
  13.  *
  14.  * (3) Two NAM's are equal.
  15.  */
  16.  
  17. #include <stdio.h>
  18. #include <stdlib.h>
  19. #include <stdarg.h>
  20. #include "Matrix.H"
  21.  
  22. Matrix *Matrix::M_NAM = Matrix::MakeNAM();
  23.  
  24. /* Constructors */
  25.  
  26. void Matrix::Init(const int m, const int n)
  27. {
  28.      M_m = m;
  29.      M_n = n;
  30.      nam = 0;
  31.      d = new double[m*n];
  32. }
  33.  
  34. Matrix::Matrix()
  35. {
  36.      nam = 1;
  37.      M_m = M_n = 0;
  38.      d = 0;
  39. }
  40.  
  41. Matrix::Matrix(int m, int n)
  42. {
  43.      Init(m, n);
  44. }
  45.  
  46. Matrix::Matrix(int m, int n, double first ...)
  47. {
  48.      int i;
  49.      va_list ap;
  50.      
  51.      Init(m, n);
  52.  
  53.      d[0] = first;
  54.      va_start(ap, first);
  55.      for (i=1; i < m*n; i += 1)
  56.       d[i] = va_arg(ap, double);
  57. }
  58.  
  59. Matrix::Matrix(int m, int n, int first ...)
  60. {
  61.      int i;
  62.      va_list ap;
  63.      
  64.      Init(m, n);
  65.  
  66.      d[0] = double(first);
  67.      va_start(ap, first);
  68.      for (i=1; i < m*n; i += 1)
  69.       d[i] = double(va_arg(ap, int));
  70. }
  71.  
  72. Matrix::Matrix(int m, int n, int *array)
  73. {
  74.      int i;
  75.  
  76.      Init(m, n);
  77.      for (i=0; i < (m*n); i += 1)
  78.       d[i] = double(array[i]);
  79. }
  80.  
  81. Matrix::Matrix(int m, int n, double *array)
  82. {
  83.      int i;
  84.  
  85.      Init(m, n);
  86.      for (i=0; i < (m*n); i += 1)
  87.       d[i] = array[i];
  88. }
  89.  
  90. Matrix::Matrix(const Matrix& m1)
  91. {
  92.      int i;
  93.  
  94.      Init(m1.M_m, m1.M_n);
  95.      for (i=0; i < (M_m*M_n); i += 1)
  96.       d[i] = m1.d[i];
  97. }
  98.  
  99. Matrix::Matrix(const Matrix& m1, int r, int c, int m, int n)
  100. {
  101.      int i,j;
  102.  
  103.      if (r<1 || c<1 || m1.Rows() < r-1+m || m1.Cols() < c-1+n)
  104.       nam = 1;
  105.      else {
  106.       Init(m,n);
  107.  
  108.       for (i=1;i<=m;i++)
  109.            for (j=1;j<=n;j++)
  110.             (*this)(i,j) = m1(r+i-1,c+j-1);
  111.      }
  112. }
  113.  
  114. Matrix::~Matrix()
  115. {
  116.      delete [M_m*M_n] d;
  117. }
  118.  
  119. /* Overloaded operators */
  120.  
  121. int operator==(const Matrix& m1, const Matrix& m2)
  122. {
  123.      int i;
  124.  
  125.      if (m1.nam != m2.nam)
  126.       return 0;
  127.      if (m1.nam == 1)
  128.       return 1;
  129.      if (! m1.SameSize(m2))
  130.       return 0;
  131.  
  132.      for (i = 0; i < m1.M_m*m1.M_n; i++)
  133.       if (m1.d[i] != m2.d[i])
  134.            return 0;
  135.  
  136.      return 1;
  137. }
  138.  
  139. Matrix operator+(const Matrix& m1, const Matrix& m2)
  140. {
  141.      if (m1.M_m != m2.M_m || m1.M_n != m2.M_n ||
  142.      m1 == Matrix::NAM() || m2 == Matrix::NAM())
  143.       return Matrix::NAM();
  144.  
  145.      Matrix m(m1);
  146.  
  147.      for (int i=0; i < (m1.M_m*m1.M_n); i++)
  148.       m.d[i] += m2.d[i];
  149.  
  150.      return m;
  151. }
  152.  
  153. Matrix operator-(const Matrix& m1, const Matrix& m2)
  154. {
  155.      if (m1.M_m != m2.M_m || m1.M_n != m2.M_n ||
  156.      m1 == Matrix::NAM() || m2 == Matrix::NAM())
  157.       return Matrix::NAM();
  158.  
  159.      Matrix m(m1);
  160.  
  161.      for (int i=0; i < (m1.M_m*m1.M_n); i++)
  162.       m.d[i] -= m2.d[i];
  163.  
  164.      return m;
  165. }
  166.  
  167. Matrix operator*(const Matrix& m1, const Matrix& m2)
  168. {
  169.      if (m1.M_n != m2.M_m || m1 == Matrix::NAM() || m2 == Matrix::NAM())
  170.       return Matrix::NAM();
  171.  
  172.      Matrix m(m1.M_m, m2.M_n);
  173.      double val = 0.0;
  174.  
  175.      for (int i=1; i <= m1.M_m; i += 1)
  176.      for (int j=1; j <= m2.M_n; j += 1) {
  177.       m(i, j) = 0.0;
  178.       for (int k=1; k <= m1.M_n; k += 1) 
  179.            m(i, j) += m1(i, k) * m2(k, j);
  180.      }
  181.      
  182.      return m;
  183. }
  184.  
  185. Matrix& Matrix::operator+=(const Matrix& m)
  186. {
  187.      if (M_m != m.M_m || M_n != m.M_n || *this == Matrix::NAM() ||
  188.      m == Matrix::NAM())
  189.       nam = 1;
  190.      else
  191.       for (int i=0; i < (M_m*M_n); i++)
  192.            d[i] += m.d[i];
  193.  
  194.      return *this;
  195. }
  196.  
  197. Matrix& Matrix::operator-=(const Matrix& m)
  198. {
  199.      if (M_m != m.M_m || M_n != m.M_n || *this == Matrix::NAM() ||
  200.      m == Matrix::NAM())
  201.       nam = 1;
  202.      else
  203.       for (int i=0; i < (M_m*M_n); i++)
  204.            d[i] -= m.d[i];
  205.  
  206.      return *this;
  207. }
  208.  
  209. Matrix operator*(double k, const Matrix& m)
  210. {
  211.    Matrix n(m);
  212.    int i;
  213.  
  214.    for (i=0;i<n.M_m*n.M_n;i++) n.d[i]*=k;
  215.    return(n);
  216. }
  217.  
  218. Matrix& Matrix::operator=(const Matrix& m)
  219. {
  220.      nam = m.nam;
  221.      
  222.      if (nam)
  223.       ;
  224.      else if (M_m == m.M_m && M_n == m.M_n)
  225.       bcopy(m.d, d, M_m*M_n*sizeof(double));
  226.      else if (M_m*M_n == m.M_m*m.M_n) {
  227.       bcopy(m.d, d, M_m*M_n*sizeof(double));
  228.       M_m = m.M_m;
  229.       M_n = m.M_n;
  230.      } else {
  231.       delete [M_m*M_n] d;
  232.       
  233.       M_m = m.M_m;
  234.       M_n = m.M_n;
  235.       
  236.       d = new double[M_m*M_n];
  237.       bcopy(m.d, d, M_m*M_n*sizeof(double));
  238.      }
  239.  
  240.      return(*this);
  241. }
  242.  
  243. Matrix operator|(const Matrix& m1, const Matrix& m2)
  244.    // Concatenates two matrices horizontally.  Requires that m1 and m2
  245.    // have the same number of rows.
  246.    // [ A ] | [ B ] --> [ A | B ]
  247. {
  248.      if (m1.Rows() != m2.Rows() || m1 == Matrix::NAM() || m2 == Matrix::NAM())
  249.       return Matrix::NAM();
  250.  
  251.      Matrix m(m1.Rows(), m1.Cols() + m2.Cols());
  252.      int r, c;
  253.  
  254.      for (r = 1; r <= m1.Rows(); r += 1)
  255.       for (c = 1; c <= m1.Cols(); c += 1)
  256.            m(r, c) = m1(r, c);
  257.  
  258.      for (r = 1; r <= m1.Rows(); r += 1)
  259.       for (c = 1; c <= m2.Cols(); c += 1)
  260.            m(r, c + m1.Cols()) = m2(r, c);
  261.  
  262.      return m;
  263. }
  264.  
  265. Matrix operator&(const Matrix& m1, const Matrix& m2)
  266.    // Concatenates two matrices vertically.  Requires that m1 and m2
  267.    // have the same number of columns.
  268.    //                   [ A ]
  269.    // [ A ] & [ B ] --> [ - ]
  270.    //                   [ B ]
  271. {
  272.      if (m1.Cols() != m2.Cols() || m1 == Matrix::NAM() || m2 == Matrix::NAM())
  273.       return Matrix::NAM();
  274.  
  275.      Matrix m(m1.Rows()+m2.Rows(), m1.Cols());
  276.      int r, c;
  277.  
  278.      for (r = 1; r <= m1.Rows(); r += 1)
  279.       for (c = 1; c <= m1.Cols(); c += 1)
  280.            m(r, c) = m1(r, c);
  281.      for (r = 1; r <= m2.Rows(); r += 1)
  282.       for (c = 1; c <= m2.Cols(); c += 1)
  283.            m(r + m1.Rows(), c) = m2(r, c);
  284.      
  285. #if 0
  286.      // This doesn't work, and I don't know why.
  287.      bcopy(m1.d, m.d, m1.Rows()*m1.Cols()*sizeof(double));
  288.      bcopy(m2.d, m.d + m1.Rows()*m1.Cols()*sizeof(double),
  289.        m2.Rows()*m2.Cols()*sizeof(double));
  290. #endif
  291.  
  292.      return m;
  293. }
  294.  
  295. ostream& operator<<(ostream& s, const Matrix& m)
  296. {
  297.      int i,j;
  298.      
  299.      if (m.nam) {
  300.       s << "NAM\n";
  301.       return(s);
  302.      }
  303.      
  304.      for (i = 1; i <= m.M_m; i++) {
  305.       for (j = 1; j <= m.M_n; j++)
  306. #ifdef __GNUG__
  307.            s.form("%6.2f ", m(i, j));
  308. #else
  309.            s << form("%6.2f ", m(i, j));
  310. #endif
  311.       s << "\n";
  312.      }
  313.      
  314.      return(s);
  315. }
  316.  
  317. /* Member functions */
  318.  
  319. Matrix Matrix::T() const
  320. {
  321.      Matrix t(M_n, M_m);
  322.  
  323.      for (int i = 1; i <= M_m; i++)
  324.       for (int j = 1; j <= M_n; j++)
  325.            t(j, i) = (*this)(i, j);
  326.  
  327.      return t;
  328. }
  329.  
  330. void Matrix::Print() const
  331. {
  332.    cout << *this;
  333. }
  334.  
  335. /* Special matrix constructors */
  336. Matrix *Matrix::MakeNAM()
  337. {
  338.      Matrix *m = new Matrix(1,1);
  339.      m->nam = 1;
  340.      return m;
  341. }
  342.  
  343. /* Matrix mutators */
  344. void Matrix::Identify()    
  345. {
  346.      if (M_m != M_n) {
  347.       nam = 1;
  348.       return;
  349.      }
  350.  
  351.      for (int i = 1; i <= M_m; i++)
  352.       for (int j = 1; j <= M_n; j++)
  353.            (*this)(i, j) = 0.0;
  354.  
  355.      for (i = 1; i <= M_m; i++)
  356.       (*this)(i, i) = 1.0;
  357. }
  358.