flower/matrix.cc

   1 #include "matrix.hh"
   2
   3 Real
   4 Matrix::norm() const
   5 {
   6     Real r =0.0;
   7     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
   8         r += sqr(dat->elem(i,j));
   9     return sqrt(r);
  10 }
  11
  12 void
  13 Matrix::fill(Real r)
  14 {
  15     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
  16         dat->elem(i,j)=r;
  17 }
  18
  19 void
  20 Matrix::set_diag(Real r)
  21 {
  22     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
  23         dat->elem(i,j)=(i==j) ? r: 0.0;
  24 }
  25
  26 void
  27 Matrix::set_diag(Vector d)
  28 {
  29     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
  30         dat->elem(i,j)=(i==j) ? d(i): 0.0;
  31 }
  32
  33 void
  34 Matrix::operator+=(Matrix const &m)
  35 {
  36     assert(m.cols() == cols());
  37     assert(m.rows() == rows());
  38     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
  39         dat->elem(i,j) += m(i,j);
  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
  52 void
  53 Matrix::operator*=(Real a)
  54 {
  55     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
  56         dat->elem(i,j) *= a;
  57 }
  58
  59 void
  60 Matrix::operator=(Matrix const &m)
  61 {
  62     if (&m == this)
  63         return ;
  64     delete dat;
  65     dat = m.dat->clone();
  66 }
  67
  68 Matrix::Matrix(Matrix const &m)
  69 {
  70     m.OK();
  71
  72     dat = m.dat->clone();
  73 }
  74
  75
  76 Matrix::Matrix(int n, int m)
  77 {
  78     dat = Matrix_storage::get_full(n,m);
  79     fill(0);
  80 }
  81
  82 Matrix::Matrix(int n)
  83 {
  84     dat = Matrix_storage::get_full(n,n);
  85     fill(0);
  86 }
  87
  88 Matrix::Matrix(Vector v, Vector w)
  89 {
  90     dat = Matrix_storage::get_full(v.dim(), w.dim());
  91     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
  92         dat->elem(i,j)=v(i)*w(j);
  93 }
  94
  95
  96 Vector
  97 Matrix::row(int k) const
  98 {
  99     int n=cols();
 100
 101
 102     Vector v(n);
 103     for(int i=0; i < n; i++)
 104         v(i)=dat->elem(k,i);
 105
 106     return v;
 107 }
 108
 109 Vector
 110 Matrix::col(int k) const
 111 {
 112     int n=rows();
 113     Vector v(n);
 114     for(int i=0; i < n; i++)
 115         v(i)=dat->elem(i,k);
 116     return v;
 117 }
 118
 119 Vector
 120 Matrix::left_multiply(Vector const & v) const
 121 {
 122      Vector dest(v.dim());
 123     assert(dat->cols()==v.dim());
 124     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
 125         dest(i)+= dat->elem(j,i)*v(j);
 126     return dest;
 127 }
 128
 129 Vector
 130 Matrix::operator *(Vector const & v) const
 131 {
 132     Vector dest(rows());
 133     assert(dat->cols()==v.dim());
 134     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
 135         dest(i)+= dat->elem(i,j)*v(j);
 136     return dest;
 137 }
 138
 139 Matrix
 140 operator /(Matrix const& m1,Real a)
 141 {
 142     Matrix m(m1);
 143     m /= a;
 144     return  m;
 145 }
 146
 147 void
 148 Matrix::transpose()             // delegate to storage?
 149 {
 150     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
 151         if (i >= j)
 152             continue;
 153         Real r=dat->elem(i,j);
 154         dat->elem(i,j) = dat->elem(j,i);
 155         dat->elem(j,i)=r;
 156     }
 157 }
 158
 159 Matrix
 160 Matrix::operator-() const
 161 {
 162     OK();
 163     Matrix m(*this);
 164     m*=-1.0;
 165     return m;
 166 }
 167
 168 Matrix
 169 Matrix::transposed() const
 170 {
 171     Matrix m(*this);
 172     m.transpose();
 173     return m;
 174 }
 175
 176
 177 /* should do something smarter: bandmatrix * bandmatrix is also banded matrix.  */
 178 Matrix
 179 operator *(Matrix const &m1, Matrix const &m2)
 180 {
 181     Matrix result(m1.rows(), m2.cols());
 182     result.set_product(m1,m2);
 183     return result;
 184 }
 185
 186 void
 187 Matrix::set_product(Matrix const &m1, Matrix const &m2)
 188 {
 189     assert(m1.cols()==m2.rows());
 190     assert(cols()==m2.cols() && rows()==m1.rows());
 191
 192     if (m1.dat->try_right_multiply(dat, m2.dat))
 193         return;
 194     for (int i=0, j=0; dat->mult_ok(i,j);
 195          dat->mult_next(i,j)) {
 196         Real r=0.0;
 197         for (int k = 0; k < m1.cols(); k++)
 198             r += m1(i,k)*m2(k,j);
 199         dat->elem(i,j)=r;
 200     }
 201 }
 202
 203 void
 204 Matrix::insert_row(Vector v, int k)
 205 {
 206     int c = cols();
 207     assert(v.dim()==cols());
 208     dat->insert_row(k);
 209     for (int j=0; j < c; j++)
 210         dat->elem(k,j)=v(j);
 211 }
 212
 213
 214 void
 215 Matrix::swap_columns(int c1, int c2)
 216 {
 217     assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
 218     int r = rows();
 219     for (int i=0; i< r; i++) {
 220         Real r=dat->elem(i,c1);
 221         dat->elem(i,c1) = dat->elem(i,c2);
 222         dat->elem(i,c2)=r;
 223     }
 224 }
 225
 226 void
 227 Matrix::swap_rows(int c1, int c2)
 228 {
 229     assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0);
 230     int c = cols();
 231     for (int i=0; i< c; i++) {
 232         Real r=dat->elem(c1,i);
 233         dat->elem(c1,i) = dat->elem(c2,i);
 234         dat->elem(c2,i)=r;
 235     }
 236 }
 237
 238
 239 int
 240 Matrix::dim() const
 241 {
 242     assert(cols() == rows());
 243     return rows();
 244 }
 245