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     int starty = dim();
 112     while (starty >= 0 ) {
 113         for ( int i = starty, j = 0; i < dim(); i++, j++ )
 114             if (dat->elem( i,j ))
 115                 goto gotcha;
 116         for ( int i=0, j = starty; j < dim(); i++,j++)
 117             if (dat->elem(i,j))
 118                 goto gotcha;
 119         starty --;
 120     }
 121 gotcha:
 122     return  starty;
 123 }
 124
 125 Matrix::Matrix(Matrix const &m)
 126 {
 127     m.OK();
 128
 129     dat = m.dat->clone();
 130 }
 131
 132
 133 Matrix::Matrix(int n, int m)
 134 {
 135     dat = Matrix_storage::get_full(n,m);
 136     fill(0);
 137 }
 138
 139 Matrix::Matrix(Matrix_storage*stor_p)
 140 {
 141     dat = stor_p;
 142 }
 143
 144 Matrix::Matrix(int n)
 145 {
 146     dat = Matrix_storage::get_full(n,n);
 147     fill(0);
 148 }
 149
 150 Matrix::Matrix(Vector v, Vector w)
 151 {
 152     dat = Matrix_storage::get_full(v.dim(), w.dim());
 153     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
 154         dat->elem(i,j)=v(i)*w(j);
 155 }
 156
 157
 158 Vector
 159 Matrix::row(int k) const
 160 {
 161     int n=cols();
 162
 163
 164     Vector v(n);
 165     for(int i=0; i < n; i++)
 166         v(i)=dat->elem(k,i);
 167
 168     return v;
 169 }
 170
 171 Vector
 172 Matrix::col(int k) const
 173 {
 174     int n=rows();
 175     Vector v(n);
 176     for(int i=0; i < n; i++)
 177         v(i)=dat->elem(i,k);
 178     return v;
 179 }
 180
 181 Vector
 182 Matrix::left_multiply(Vector const & v) const
 183 {
 184      Vector dest(v.dim());
 185     assert(dat->cols()==v.dim());
 186     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
 187         dest(i)+= dat->elem(j,i)*v(j);
 188     return dest;
 189 }
 190
 191 Vector
 192 Matrix::operator *(Vector const & v) const
 193 {
 194     Vector dest(rows());
 195     assert(dat->cols()==v.dim());
 196     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
 197         dest(i)+= dat->elem(i,j)*v(j);
 198     return dest;
 199 }
 200
 201 Matrix
 202 operator /(Matrix const& m1,Real a)
 203 {
 204     Matrix m(m1);
 205     m /= a;
 206     return  m;
 207 }
 208
 209 void
 210 Matrix::transpose()             // delegate to storage?
 211 {
 212     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
 213         if (i >= j)
 214             continue;
 215         Real r=dat->elem(i,j);
 216         dat->elem(i,j) = dat->elem(j,i);
 217         dat->elem(j,i)=r;
 218     }
 219 }
 220
 221 Matrix
 222 Matrix::operator-() const
 223 {
 224     OK();
 225     Matrix m(*this);
 226     m*=-1.0;
 227     return m;
 228 }
 229
 230 Matrix
 231 Matrix::transposed() const
 232 {
 233     Matrix m(*this);
 234     m.transpose();
 235     return m;
 236 }
 237
 238 Matrix
 239 operator *(Matrix const &m1, Matrix const &m2)
 240 {
 241     Matrix result(Matrix_storage::get_product_result(m1.dat, m2.dat));
 242
 243     result.set_product(m1,m2);
 244     return result;
 245 }
 246
 247 void
 248 Matrix::set_product(Matrix const &m1, Matrix const &m2)
 249 {
 250     assert(m1.cols()==m2.rows());
 251     assert(cols()==m2.cols() && rows()==m1.rows());
 252
 253     if (m1.dat->try_right_multiply(dat, m2.dat))
 254         return;
 255     for (int i=0, j=0; dat->mult_ok(i,j);
 256          dat->mult_next(i,j)) {
 257         Real r=0.0;
 258         for (int k = 0; k < m1.cols(); k++)
 259             r += m1(i,k)*m2(k,j);
 260         dat->elem(i,j)=r;
 261     }
 262 }
 263
 264 void
 265 Matrix::insert_row(Vector v, int k)
 266 {
 267     int c = cols();
 268     assert(v.dim()==cols());
 269     dat->insert_row(k);
 270     for (int j=0; j < c; j++)
 271         dat->elem(k,j)=v(j);
 272 }
 273
 274
 275 void
 276 Matrix::swap_columns(int c1, int c2)
 277 {
 278     assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
 279     int r = rows();
 280     for (int i=0; i< r; i++) {
 281         Real r=dat->elem(i,c1);
 282         dat->elem(i,c1) = dat->elem(i,c2);
 283         dat->elem(i,c2)=r;
 284     }
 285 }
 286
 287 void
 288 Matrix::swap_rows(int c1, int c2)
 289 {
 290     assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0);
 291     int c = cols();
 292     for (int i=0; i< c; i++) {
 293         Real r=dat->elem(c1,i);
 294         dat->elem(c1,i) = dat->elem(c2,i);
 295         dat->elem(c2,i)=r;
 296     }
 297 }
 298
 299
 300 int
 301 Matrix::dim() const
 302 {
 303     assert(cols() == rows());
 304     return rows();
 305 }
 306
 307
 308