flower/matrix.cc

   1 /*
   2   matrix.cc -- implement Matrix
   3
   4   source file of the Flower Library
   5
   6   (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
   7 */
   8
   9 #include "matrix.hh"
  10 #include "full-storage.hh"
  11
  12 Real
  13 Matrix::norm() const
  14 {
  15   Real r =0.0;
  16   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  17     r += sqr (dat_->elem (i,j));
  18   return sqrt (r);
  19 }
  20
  21 void
  22 Matrix::fill (Real r)
  23 {
  24   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  25     dat_->elem (i,j)=r;
  26 }
  27
  28 void
  29 Matrix::set_diag (Real r)
  30 {
  31   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  32     dat_->elem (i,j)=(i==j) ? r: 0.0;
  33 }
  34
  35 void
  36 Matrix::set_diag (Vector d)
  37 {
  38   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  39     dat_->elem (i,j)=(i==j) ? d (i): 0.0;
  40 }
  41
  42 void
  43 Matrix::operator+=(Matrix const &m)
  44 {
  45   assert (m.cols() == cols ());
  46   assert (m.rows() == rows ());
  47   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  48     dat_->elem (i,j) += m (i,j);
  49 }
  50
  51 void
  52 Matrix::operator-=(Matrix const &m)
  53 {
  54   assert (m.cols() == cols ());
  55   assert (m.rows() == rows ());
  56   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  57     dat_->elem (i,j) -= m (i,j);
  58 }
  59
  60
  61 void
  62 Matrix::operator*=(Real a)
  63 {
  64   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
  65     dat_->elem (i,j) *= a;
  66 }
  67
  68 void
  69 Matrix::operator=(Matrix const &m)
  70 {
  71   if (&m == this)
  72     return ;
  73   delete dat_;
  74   dat_ = new Full_storage (*m.dat_);
  75 }
  76
  77
  78 int
  79 Matrix::band_i () const
  80 {
  81   return dat_->band_i_;
  82 }
  83
  84 void
  85 Matrix::set_band ()
  86 {
  87   dat_->band_i_ = calc_band_i ();
  88 }
  89
  90 int
  91 Matrix::calc_band_i() const
  92 {
  93   int starty = dim();
  94   while (starty >= 0)
  95     {
  96       for (int i = starty, j = 0; i < dim(); i++, j++)
  97         if (dat_->elem (i,j))
  98           goto gotcha;
  99       for (int i=0, j = starty; j < dim(); i++,j++)
 100         if (dat_->elem (i,j))
 101           goto gotcha;
 102       starty --;
 103     }
 104  gotcha:
 105   return  starty;
 106 }
 107
 108 Matrix::Matrix (Matrix const &m)
 109 {
 110   dat_ = new Full_storage (*m.dat_);
 111 }
 112
 113
 114 Matrix::Matrix (int n, int m)
 115 {
 116   dat_ = new Full_storage(n,m);
 117   fill (0);
 118 }
 119
 120
 121 Matrix::Matrix (int n)
 122 {
 123   dat_= new Full_storage (n,n);
 124   fill (0);
 125 }
 126
 127 Matrix::Matrix (Vector v, Vector w)
 128 {
 129   int n = v.dim();
 130   int m = w.dim ();
 131   dat_= new Full_storage (n,m);
 132   for (int i=0; i < n; i++)
 133     for (int j=0; j < m ; j++)
 134       dat_->elem (i,j)=v (i)*w (j);
 135 }
 136
 137
 138 Vector
 139 Matrix::row (int k) const
 140 {
 141   int n=cols();
 142
 143
 144   Vector v (n);
 145   for (int i=0; i < n; i++)
 146     v (i)=dat_->elem (k,i);
 147
 148   return v;
 149 }
 150
 151 Vector
 152 Matrix::col (int k) const
 153 {
 154   int n=rows();
 155   Vector v (n);
 156   for (int i=0; i < n; i++)
 157     v (i)=dat_->elem (i,k);
 158   return v;
 159 }
 160
 161 Vector
 162 Matrix::left_multiply (Vector const & v) const
 163 {
 164   Vector dest (v.dim());
 165   assert (dat_->cols()==v.dim ());
 166   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
 167     dest (i)+= dat_->elem (j,i)*v (j);
 168   return dest;
 169 }
 170
 171 Vector
 172 Matrix::operator *(Vector const & v) const
 173 {
 174   Vector dest (rows());
 175   assert (dat_->cols()==v.dim ());
 176   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
 177     dest (i)+= dat_->elem (i,j)*v (j);
 178   return dest;
 179 }
 180
 181 Matrix
 182 operator /(Matrix const& m1,Real a)
 183 {
 184   Matrix m (m1);
 185   m /= a;
 186   return  m;
 187 }
 188
 189 /*
 190   ugh. Only works for square matrices.
 191  */
 192 void
 193 Matrix::transpose()
 194 {
 195   for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
 196     {
 197       if (i >= j)
 198         continue;
 199       Real r=dat_->elem (i,j);
 200       dat_->elem (i,j) = dat_->elem (j,i);
 201       dat_->elem (j,i)=r;
 202     }
 203 }
 204
 205 Matrix
 206 Matrix::operator-() const
 207 {
 208   OK();
 209   Matrix m (*this);
 210   m*=-1.0;
 211   return m;
 212 }
 213
 214 Matrix
 215 Matrix::transposed() const
 216 {
 217   Matrix m (*this);
 218   m.transpose();
 219   return m;
 220 }
 221
 222 Matrix
 223 operator *(Matrix const &m1, Matrix const &m2)
 224 {
 225   Matrix result (m1.rows (), m2.cols ());
 226   result.set_product (m1,m2);
 227   return result;
 228 }
 229
 230
 231
 232 void
 233 Matrix::set_product (Matrix const &m1, Matrix const &m2)
 234 {
 235   assert (m1.cols()==m2.rows ());
 236   assert (cols()==m2.cols () && rows ()==m1.rows ());
 237
 238   for (int i=0; i < rows (); i++)
 239     for (int j=0; j <  cols (); j++)
 240       {
 241         Real r=0.0;
 242         for (int k = 0 >? i - m1.band_i () >? j - m2.band_i ();
 243              k < m1.cols(); k++)
 244           {
 245             r += m1(i,k)*m2(k,j);
 246           }
 247         dat_->elem (i,j)=r;
 248       }
 249 }
 250
 251 void
 252 Matrix::insert_row (Vector v, int k)
 253 {
 254   int c = cols();
 255   assert (v.dim()==cols ());
 256   dat_->insert_row (k);
 257   for (int j=0; j < c; j++)
 258     dat_->elem (k,j)=v (j);
 259 }
 260
 261
 262 void
 263 Matrix::swap_columns (int c1, int c2)
 264 {
 265   assert (c1>=0&& c1 < cols()&&c2 < cols () && c2 >=0);
 266   int r = rows();
 267   for (int i=0; i< r; i++)
 268     {
 269       Real r=dat_->elem (i,c1);
 270       dat_->elem (i,c1) = dat_->elem (i,c2);
 271       dat_->elem (i,c2)=r;
 272     }
 273 }
 274
 275 void
 276 Matrix::swap_rows (int c1, int c2)
 277 {
 278   assert (c1>=0&& c1 < rows()&&c2 < rows () && c2 >=0);
 279   int c = cols();
 280   for (int i=0; i< c; i++)
 281     {
 282       Real r=dat_->elem (c1,i);
 283       dat_->elem (c1,i) = dat_->elem (c2,i);
 284       dat_->elem (c2,i)=r;
 285     }
 286 }
 287
 288
 289 int
 290 Matrix::dim() const
 291 {
 292   assert (cols() == rows ());
 293   return rows();
 294 }
 295
 296 Matrix::~Matrix ()
 297 {
 298   delete dat_;
 299 }