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"
11 #include "diagonal-storage.hh"
16 return dat->is_type_b( Diagonal_storage::static_name());
20 Matrix::set_full() const
22 if ( dat->name() != Full_storage::static_name()) {
23 Matrix_storage::set_full( ((Matrix*)this)->dat );
28 Matrix::try_set_band()const
36 // it only looks constant
37 Matrix*self = (Matrix*)this;
38 Matrix_storage::set_band(self->dat,b);
45 for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
46 r += sqr(dat->elem(i,j));
53 for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
58 Matrix::set_diag(Real r)
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;
65 Matrix::set_diag(Vector d)
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;
72 Matrix::operator+=(Matrix const &m)
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);
82 Matrix::operator-=(Matrix const &m)
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);
93 Matrix::operator*=(Real a)
95 for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
100 Matrix::operator=(Matrix const &m)
105 dat = m.dat->clone();
109 Matrix::band_i()const
112 while (starty >= 0 ) {
113 for ( int i = starty, j = 0; i < dim(); i++, j++ )
114 if (dat->elem( i,j ))
116 for ( int i=0, j = starty; j < dim(); i++,j++)
125 Matrix::Matrix(Matrix const &m)
129 dat = m.dat->clone();
133 Matrix::Matrix(int n, int m)
135 dat = Matrix_storage::get_full(n,m);
139 Matrix::Matrix(Matrix_storage*stor_p)
144 Matrix::Matrix(int n)
146 dat = Matrix_storage::get_full(n,n);
150 Matrix::Matrix(Vector v, Vector w)
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);
159 Matrix::row(int k) const
165 for(int i=0; i < n; i++)
172 Matrix::col(int k) const
176 for(int i=0; i < n; i++)
182 Matrix::left_multiply(Vector const & v) const
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);
192 Matrix::operator *(Vector const & v) const
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);
202 operator /(Matrix const& m1,Real a)
210 Matrix::transpose() // delegate to storage?
212 for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
215 Real r=dat->elem(i,j);
216 dat->elem(i,j) = dat->elem(j,i);
222 Matrix::operator-() const
231 Matrix::transposed() const
239 operator *(Matrix const &m1, Matrix const &m2)
241 Matrix result(Matrix_storage::get_product_result(m1.dat, m2.dat));
243 result.set_product(m1,m2);
248 Matrix::set_product(Matrix const &m1, Matrix const &m2)
250 assert(m1.cols()==m2.rows());
251 assert(cols()==m2.cols() && rows()==m1.rows());
253 if (m1.dat->try_right_multiply(dat, m2.dat))
255 for (int i=0, j=0; dat->mult_ok(i,j);
256 dat->mult_next(i,j)) {
258 for (int k = 0; k < m1.cols(); k++)
259 r += m1(i,k)*m2(k,j);
265 Matrix::insert_row(Vector v, int k)
268 assert(v.dim()==cols());
270 for (int j=0; j < c; j++)
276 Matrix::swap_columns(int c1, int c2)
278 assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
280 for (int i=0; i< r; i++) {
281 Real r=dat->elem(i,c1);
282 dat->elem(i,c1) = dat->elem(i,c2);
288 Matrix::swap_rows(int c1, int c2)
290 assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0);
292 for (int i=0; i< c; i++) {
293 Real r=dat->elem(c1,i);
294 dat->elem(c1,i) = dat->elem(c2,i);
303 assert(cols() == rows());