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