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