lib-src/numeric/matrix.cc

   1 //////////////////////////////////////////////////////////////////////////////
   2 // NOTICE:
   3 //
   4 // ADLib, Prop and their related set of tools and documentation are in the
   5 // public domain.   The author(s) of this software reserve no copyrights on
   6 // the source code and any code generated using the tools.  You are encouraged
   7 // to use ADLib and Prop to develop software, in both academic and commercial
   8 // settings, and are free to incorporate any part of ADLib and Prop into
   9 // your programs.
  10 //
  11 // Although you are under no obligation to do so, we strongly recommend that
  12 // you give away all software developed using our tools.
  13 //
  14 // We also ask that credit be given to us when ADLib and/or Prop are used in
  15 // your programs, and that this notice be preserved intact in all the source
  16 // code.
  17 //
  18 // This software is still under development and we welcome any suggestions
  19 // and help from the users.
  20 //
  21 // Allen Leung
  22 // 1994
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 #include <stdio.h>
  26 #include <stdarg.h>
  27 #include <stdlib.h>
  28 #include <iostream.h>
  29 #include <AD/numeric/matrix.h>
  30
  31 ////////////////////////////////////////////////////////////////////////////
  32 // Default handler
  33 ////////////////////////////////////////////////////////////////////////////
  34 static void error_handler(const char format[], ...)
  35 {  va_list args;
  36    va_start(args,format);
  37    fprintf  (stderr, "[ Matrix error: ");
  38    vfprintf (stderr, format, args);
  39    fprintf  (stderr, "]\n");
  40    va_end(args);
  41    exit(1);
  42 }
  43
  44 void (*Matrix::error)(const char [], ...) = error_handler;
  45
  46 ////////////////////////////////////////////////////////////////////////////
  47 // Bounds checking
  48 ////////////////////////////////////////////////////////////////////////////
  49 void Matrix::can_add(const Matrix& m) const
  50 {  if (m.rows() != rows() || m.columns() != columns())
  51       Matrix::error(
  52          "Can't add/subtract/compare matrices with different dimensions [%d,%d], [%d,%d]",
  53          rows(), columns(), m.rows(), m.columns());
  54
  55 }
  56
  57 void Matrix::can_mult(const Matrix& m) const
  58 {  if (columns() != m.rows())
  59       Matrix::error(
  60          "Can't multiply matrices with incompatible dimensions [%d,%d], [%d,%d]",
  61          rows(), columns(), m.rows(), m.columns());
  62 }
  63
  64 ////////////////////////////////////////////////////////////////////////////
  65 // Negate
  66 ////////////////////////////////////////////////////////////////////////////
  67 Matrix Matrix::operator - ()
  68 {  Matrix r(M,N);
  69    for (int i = len - 1; i >= 0; i--) r.array[i] = -array[i]; return r;
  70 }
  71
  72 ////////////////////////////////////////////////////////////////////////////
  73 // Matrix operations with scalars
  74 ////////////////////////////////////////////////////////////////////////////
  75 Matrix operator + (double a, const Matrix& m)
  76 {  Matrix r(m.M,m.N);
  77    for (int i = m.len - 1; i >= 0; i--) r.array[i] = a + m.array[i];
  78    return r;
  79 }
  80 Matrix operator + (const Matrix& m, double a)
  81 {  Matrix r(m.M,m.N);
  82    for (int i = m.len - 1; i >= 0; i--) r.array[i] = m.array[i] + a;
  83    return r;
  84 }
  85 Matrix operator - (double a, const Matrix& m)
  86 {  Matrix r(m.M,m.N);
  87    for (int i = m.len - 1; i >= 0; i--) r.array[i] = a - m.array[i];
  88    return r;
  89 }
  90 Matrix operator - (const Matrix& m, double a)
  91 {  Matrix r(m.M,m.N);
  92    for (int i = m.len - 1; i >= 0; i--) r.array[i] = m.array[i] - a;
  93    return r;
  94 }
  95 Matrix operator * (double a, const Matrix& m)
  96 {  Matrix r(m.M,m.N);
  97    for (int i = m.len - 1; i >= 0; i--) r.array[i] = a * m.array[i];
  98    return r;
  99 }
 100 Matrix operator * (const Matrix& m, double a)
 101 {  Matrix r(m.M,m.N);
 102    for (int i = m.len - 1; i >= 0; i--) r.array[i] = m.array[i] * a;
 103    return r;
 104 }
 105 Matrix operator / (const Matrix& m, double a)
 106 {  Matrix r(m.M,m.N);
 107    for (int i = m.len - 1; i >= 0; i--) r.array[i] = m.array[i] / a;
 108    return r;
 109 }
 110
 111 //////////////////////////////////////////////////////////////////////////
 112 //  In place operations with scalars
 113 //////////////////////////////////////////////////////////////////////////
 114 Matrix& Matrix::negate()
 115 {  for (int i = len - 1; i >= 0; i--) array[i] = -array[i]; return *this; }
 116 Matrix& Matrix::operator += (double a)
 117 {  for (int i = len - 1; i >= 0; i--) array[i] += a; return *this; }
 118 Matrix& Matrix::operator -= (double a)
 119 {  for (int i = len - 1; i >= 0; i--) array[i] -= a; return *this; }
 120 Matrix& Matrix::operator *= (double a)
 121 {  for (int i = len - 1; i >= 0; i--) array[i] *= a; return *this; }
 122 Matrix& Matrix::operator /= (double a)
 123 {  for (int i = len - 1; i >= 0; i--) array[i] /= a; return *this; }
 124
 125 //////////////////////////////////////////////////////////////////////////
 126 //  Matrix/Matrix operations
 127 //////////////////////////////////////////////////////////////////////////
 128 Matrix operator + (const Matrix& a, const Matrix& b)
 129 {  a.can_add(b);
 130    Matrix r(a.M,a.N);
 131    for (int i = a.len - 1; i >= 0; i--) r.array[i] = a.array[i] + b.array[i];
 132    return r;
 133 }
 134
 135 Matrix operator - (const Matrix& a, const Matrix& b)
 136 {  a.can_add(b);
 137    Matrix r(a.M,a.N);
 138    for (int i = a.len - 1; i >= 0; i--) r.array[i] = a.array[i] - b.array[i];
 139    return r;
 140 }
 141
 142 Matrix operator * (const Matrix& a, const Matrix& b)
 143 {  a.can_mult(b);
 144    Matrix r(a.N,b.M);
 145    for (int i = 0; i < a.N; i++) {
 146       int j;
 147       double sum = 0;
 148       for (j = 0; j < b.M; j++)
 149          for (int k = 0; k < a.M; k++)
 150             sum += a[i][k] * b[k][j];
 151       r[i][j] = sum;
 152    }
 153    return r;
 154 }
 155
 156 //////////////////////////////////////////////////////////////////////////
 157 //  In place matrix operations
 158 //////////////////////////////////////////////////////////////////////////
 159 double Matrix::det() const
 160 {  can_mult(*this);
 161    return 0;
 162 }
 163
 164 Matrix Matrix::inverse() const
 165 {  can_mult(*this);
 166    return *this;
 167 }
 168
 169 Matrix& Matrix::operator += (const Matrix& m)
 170 {  can_add(m);
 171    for (int i = len - 1; i >= 0; i--) array[i] += m.array[i];
 172    return *this;
 173 }
 174
 175 Matrix& Matrix::operator -= (const Matrix& m)
 176 {  can_add(m);
 177    for (int i = len - 1; i >= 0; i--) array[i] -= m.array[i]; return *this;
 178 }
 179
 180 Matrix& Matrix::operator *= (const Matrix& m)
 181 { *this = *this * m; return *this; }
 182
 183 ////////////////////////////////////////////////////////////////////////////
 184 //  Comparison
 185 ////////////////////////////////////////////////////////////////////////////
 186 Bool operator == (const Matrix& a, const Matrix& b)
 187 {  a.can_add(b);
 188    for (int i = a.len - 1; i >= 0; i--)
 189       if (a.array[i] != b.array[i]) return false;
 190    return true;
 191 }