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 = virtual_smat::get_full(n,m);
84 dat = virtual_smat::get_full(n,n);
88 Matrix::Matrix(Vector v, Vector w)
90 dat = virtual_smat::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 for (int i=0, j=0; dat->mult_ok(i,j);
193 dat->mult_next(i,j)) {
195 for (int k = 0; k < m1.cols(); k++)
196 r += m1(i,k)*m2(k,j);
202 Matrix::insert_row(Vector v, int k)
205 assert(v.dim()==cols());
207 for (int j=0; j < c; j++)
213 Matrix::swap_columns(int c1, int c2)
215 assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
217 for (int i=0; i< r; i++) {
218 Real r=dat->elem(i,c1);
219 dat->elem(i,c1) = dat->elem(i,c2);
225 Matrix::swap_rows(int c1, int c2)
227 assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0);
229 for (int i=0; i< c; i++) {
230 Real r=dat->elem(c1,i);
231 dat->elem(c1,i) = dat->elem(c2,i);
240 assert(cols() == rows());