2 matrix.cc -- implement Matrix
4 source file of the Flower Library
6 (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.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 ());