7 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
8 r
+= sqr(dat
->elem(i
,j
));
15 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
20 Matrix::set_diag(Real r
)
22 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
23 dat
->elem(i
,j
)=(i
==j
) ? r
: 0.0;
27 Matrix::set_diag(Vector d
)
29 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
30 dat
->elem(i
,j
)=(i
==j
) ? d(i
): 0.0;
34 Matrix::operator+=(Matrix
const &m
)
36 assert(m
.cols() == cols());
37 assert(m
.rows() == rows());
38 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
39 dat
->elem(i
,j
) += m(i
,j
);
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
);
53 Matrix::operator*=(Real a
)
55 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
60 Matrix::operator=(Matrix
const &m
)
68 Matrix::Matrix(Matrix
const &m
)
76 Matrix::Matrix(int n
, int m
)
78 dat
= Matrix_storage::get_full(n
,m
);
84 dat
= Matrix_storage::get_full(n
,n
);
88 Matrix::Matrix(Vector v
, Vector w
)
90 dat
= Matrix_storage::get_full(v
.dim(), w
.dim());
91 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
92 dat
->elem(i
,j
)=v(i
)*w(j
);
97 Matrix::row(int k
) const
103 for(int i
=0; i
< n
; i
++)
110 Matrix::col(int k
) const
114 for(int i
=0; i
< n
; i
++)
120 Matrix::left_multiply(Vector
const & v
) const
122 Vector
dest(v
.dim());
123 assert(dat
->cols()==v
.dim());
124 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
125 dest(i
)+= dat
->elem(j
,i
)*v(j
);
130 Matrix::operator *(Vector
const & v
) const
133 assert(dat
->cols()==v
.dim());
134 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
))
135 dest(i
)+= dat
->elem(i
,j
)*v(j
);
140 operator /(Matrix
const& m1
,Real a
)
148 Matrix::transpose() // delegate to storage?
150 for (int i
=0, j
=0; dat
->mult_ok(i
,j
); dat
->mult_next(i
,j
)) {
153 Real r
=dat
->elem(i
,j
);
154 dat
->elem(i
,j
) = dat
->elem(j
,i
);
160 Matrix::operator-() const
169 Matrix::transposed() const
177 /* should do something smarter: bandmatrix * bandmatrix is also banded matrix. */
179 operator *(Matrix
const &m1
, Matrix
const &m2
)
181 Matrix
result(m1
.rows(), m2
.cols());
182 result
.set_product(m1
,m2
);
187 Matrix::set_product(Matrix
const &m1
, Matrix
const &m2
)
189 assert(m1
.cols()==m2
.rows());
190 assert(cols()==m2
.cols() && rows()==m1
.rows());
192 if (m1
.dat
->try_right_multiply(dat
, m2
.dat
))
194 for (int i
=0, j
=0; dat
->mult_ok(i
,j
);
195 dat
->mult_next(i
,j
)) {
197 for (int k
= 0; k
< m1
.cols(); k
++)
198 r
+= m1(i
,k
)*m2(k
,j
);
204 Matrix::insert_row(Vector v
, int k
)
207 assert(v
.dim()==cols());
209 for (int j
=0; j
< c
; j
++)
215 Matrix::swap_columns(int c1
, int c2
)
217 assert(c1
>=0&& c1
< cols()&&c2
< cols() && c2
>=0);
219 for (int i
=0; i
< r
; i
++) {
220 Real r
=dat
->elem(i
,c1
);
221 dat
->elem(i
,c1
) = dat
->elem(i
,c2
);
227 Matrix::swap_rows(int c1
, int c2
)
229 assert(c1
>=0&& c1
< rows()&&c2
< rows() && c2
>=0);
231 for (int i
=0; i
< c
; i
++) {
232 Real r
=dat
->elem(c1
,i
);
233 dat
->elem(c1
,i
) = dat
->elem(c2
,i
);
242 assert(cols() == rows());