2 matrix.cc -- implement Matrix
4 source file of the Flower Library
6 (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
10 #include "full-storage.hh"
16 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
17 r += sqr (dat_->elem (i,j));
24 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
29 Matrix::set_diag (Real r)
31 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
32 dat_->elem (i,j)=(i==j) ? r: 0.0;
36 Matrix::set_diag (Vector d)
38 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
39 dat_->elem (i,j)=(i==j) ? d (i): 0.0;
43 Matrix::operator+=(Matrix const &m)
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);
52 Matrix::operator-=(Matrix const &m)
54 assert (m.cols() == cols ());
55 assert (m.rows() == rows ());
56 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
57 dat_->elem (i,j) -= m (i,j);
62 Matrix::operator*=(Real a)
64 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
65 dat_->elem (i,j) *= a;
69 Matrix::operator=(Matrix const &m)
74 dat_ = new Full_storage (*m.dat_);
79 Matrix::band_i () const
87 dat_->band_i_ = calc_band_i ();
91 Matrix::calc_band_i() const
96 for (int i = starty, j = 0; i < dim(); i++, j++)
99 for (int i=0, j = starty; j < dim(); i++,j++)
100 if (dat_->elem (i,j))
108 Matrix::Matrix (Matrix const &m)
110 dat_ = new Full_storage (*m.dat_);
114 Matrix::Matrix (int n, int m)
116 dat_ = new Full_storage(n,m);
121 Matrix::Matrix (int n)
123 dat_= new Full_storage (n,n);
127 Matrix::Matrix (Vector v, Vector w)
131 dat_= new Full_storage (n,m);
132 for (int i=0; i < n; i++)
133 for (int j=0; j < m ; j++)
134 dat_->elem (i,j)=v (i)*w (j);
139 Matrix::row (int k) const
145 for (int i=0; i < n; i++)
146 v (i)=dat_->elem (k,i);
152 Matrix::col (int k) const
156 for (int i=0; i < n; i++)
157 v (i)=dat_->elem (i,k);
162 Matrix::left_multiply (Vector const & v) const
164 Vector dest (v.dim());
165 assert (dat_->cols()==v.dim ());
166 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
167 dest (i)+= dat_->elem (j,i)*v (j);
172 Matrix::operator *(Vector const & v) const
174 Vector dest (rows());
175 assert (dat_->cols()==v.dim ());
176 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
177 dest (i)+= dat_->elem (i,j)*v (j);
182 operator /(Matrix const& m1,Real a)
190 ugh. Only works for square matrices.
195 for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
199 Real r=dat_->elem (i,j);
200 dat_->elem (i,j) = dat_->elem (j,i);
206 Matrix::operator-() const
215 Matrix::transposed() const
223 operator *(Matrix const &m1, Matrix const &m2)
225 Matrix result (m1.rows (), m2.cols ());
226 result.set_product (m1,m2);
233 Matrix::set_product (Matrix const &m1, Matrix const &m2)
235 assert (m1.cols()==m2.rows ());
236 assert (cols()==m2.cols () && rows ()==m1.rows ());
238 for (int i=0; i < rows (); i++)
239 for (int j=0; j < cols (); j++)
242 for (int k = 0 >? i - m1.band_i () >? j - m2.band_i ();
245 r += m1(i,k)*m2(k,j);
252 Matrix::insert_row (Vector v, int k)
255 assert (v.dim()==cols ());
256 dat_->insert_row (k);
257 for (int j=0; j < c; j++)
258 dat_->elem (k,j)=v (j);
263 Matrix::swap_columns (int c1, int c2)
265 assert (c1>=0&& c1 < cols()&&c2 < cols () && c2 >=0);
267 for (int i=0; i< r; i++)
269 Real r=dat_->elem (i,c1);
270 dat_->elem (i,c1) = dat_->elem (i,c2);
276 Matrix::swap_rows (int c1, int c2)
278 assert (c1>=0&& c1 < rows()&&c2 < rows () && c2 >=0);
280 for (int i=0; i< c; i++)
282 Real r=dat_->elem (c1,i);
283 dat_->elem (c1,i) = dat_->elem (c2,i);
292 assert (cols() == rows ());