matrix.cc

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