9 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
10 r
+= sqr(dat
->elem(i
,j
));
16 Matrix::operator()(int i
,int j
) const
18 assert(i
>= 0 && j
>= 0);
19 assert(i
< rows() && j
< cols());
20 return dat
->elem(i
,j
);
25 Matrix::operator()(int i
, int j
)
27 assert(i
>= 0 && j
>= 0);
28 assert(i
< rows() && j
< cols());
29 return dat
->elem(i
,j
);
35 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
40 Matrix::set_diag(Real r
)
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;
47 Matrix::set_diag(Vector d
)
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;
54 Matrix::operator+=(const Matrix
&m
)
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
);
63 Matrix::operator-=(const Matrix
&m
)
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
);
73 Matrix::operator*=(Real a
)
75 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
80 Matrix::operator=(const Matrix
&m
)
88 Matrix::Matrix(const Matrix
&m
)
96 Matrix::Matrix(int n
, int m
)
98 dat
= virtual_smat::get_full(n
,m
);
102 Matrix::Matrix(int n
)
104 dat
= virtual_smat::get_full(n
,n
);
108 Matrix::Matrix(Vector v
, Vector w
)
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
);
117 Matrix::row(int k
) const
123 for(int i
=0; i
< n
; i
++)
130 Matrix::col(int k
) const
134 for(int i
=0; i
< n
; i
++)
140 Matrix::left_multiply(const Vector
& v
) const
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
);
150 Matrix::operator *(const Vector
& v
) const
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
);
160 operator /(Matrix
const& m1
,Real a
)
168 Matrix::transpose() // delegate to storage?
170 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
)) {
173 Real r
=dat
->elem(i
,j
);
174 dat
->elem(i
,j
) = dat
->elem(j
,i
);
180 Matrix::operator-() const
189 Matrix::transposed() const
197 /* should do something smarter: bandmatrix * bandmatrix is also banded matrix. */
199 operator *(const Matrix
&m1
, const Matrix
&m2
)
201 Matrix
result(m1
.rows(), m2
.cols());
202 result
.set_product(m1
,m2
);
207 Matrix::set_product(const Matrix
&m1
, const Matrix
&m2
)
209 assert(m1
.cols()==m2
.rows());
210 assert(cols()==m2
.cols() && rows()==m1
.rows());
212 for (int i
=0, j
=0; dat
->mult_ok(i
,j
);
213 dat
->mult_next(i
,j
)) {
215 for (int k
= 0; k
< m1
.cols(); k
++)
216 r
+= m1(i
,k
)*m2(k
,j
);
222 Matrix::insert_row(Vector v
, int k
)
224 assert(v
.dim()==cols());
226 for (int j
=0; j
< cols(); j
++)
232 Matrix::swap_columns(int c1
, int c2
)
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
);
243 Matrix::swap_rows(int c1
, int c2
)
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
);
257 assert(cols() == rows());
261 Matrix::operator String() const
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 ");
278 Matrix::print() const