]> git.donarmstrong.com Git - lilypond.git/blob - flower/matrix.cc
release: 0.0.42.pre3
[lilypond.git] / flower / matrix.cc
1 #include "matrix.hh"
2
3 Real
4 Matrix::norm() const
5 {
6     Real r =0.0;
7     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
8         r += sqr(dat->elem(i,j));
9     return sqrt(r);
10 }
11
12 void
13 Matrix::fill(Real r)
14 {
15     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
16         dat->elem(i,j)=r;
17 }
18
19 void
20 Matrix::set_diag(Real r)
21 {
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;
24 }
25
26 void
27 Matrix::set_diag(Vector d)
28 {
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;
31 }
32
33 void
34 Matrix::operator+=(Matrix const &m)
35 {
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);
40 }
41  
42 void
43 Matrix::operator-=(Matrix const &m)
44 {
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);
49 }
50
51
52 void
53 Matrix::operator*=(Real a)
54 {
55     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
56         dat->elem(i,j) *= a;
57 }
58
59 void
60 Matrix::operator=(Matrix const &m)
61 {
62     if (&m == this)
63         return ;
64     delete dat;
65     dat = m.dat->clone();
66 }
67     
68 Matrix::Matrix(Matrix const &m)
69 {
70     m.OK();
71     
72     dat = m.dat->clone();
73 }
74
75
76 Matrix::Matrix(int n, int m)
77 {
78     dat = virtual_smat::get_full(n,m);
79     fill(0);
80 }
81
82 Matrix::Matrix(int n)
83 {
84     dat = virtual_smat::get_full(n,n);
85     fill(0);
86 }
87
88 Matrix::Matrix(Vector v, Vector w)
89 {   
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);
93 }
94
95
96 Vector
97 Matrix::row(int k) const
98 {
99     int n=cols();
100
101     
102     Vector v(n);
103     for(int i=0; i < n; i++)
104         v(i)=dat->elem(k,i);
105
106     return v;
107 }
108
109 Vector
110 Matrix::col(int k) const
111 {
112     int n=rows();
113     Vector v(n);
114     for(int i=0; i < n; i++)
115         v(i)=dat->elem(i,k);
116     return v;
117 }
118
119 Vector
120 Matrix::left_multiply(Vector const & v) const
121 {
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);
126     return dest;
127 }
128
129 Vector
130 Matrix::operator *(Vector const & v) const
131 {
132     Vector dest(rows());
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);
136     return dest;
137 }
138
139 Matrix
140 operator /(Matrix const& m1,Real a)
141 {
142     Matrix m(m1);
143     m /= a;
144     return  m;
145 }
146
147 void
148 Matrix::transpose()             // delegate to storage?
149 {
150     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
151         if (i >= j)
152             continue;
153         Real r=dat->elem(i,j);
154         dat->elem(i,j) = dat->elem(j,i);
155         dat->elem(j,i)=r;
156     }
157 }
158
159 Matrix
160 Matrix::operator-() const
161 {
162     OK();
163     Matrix m(*this);
164     m*=-1.0;    
165     return m;
166 }
167
168 Matrix
169 Matrix::transposed() const
170 {
171     Matrix m(*this);
172     m.transpose();
173     return m;
174 }
175
176
177 /* should do something smarter: bandmatrix * bandmatrix is also banded matrix.  */
178 Matrix
179 operator *(Matrix const &m1, Matrix const &m2)
180 {
181     Matrix result(m1.rows(), m2.cols());
182     result.set_product(m1,m2);
183     return result;
184 }
185
186 void
187 Matrix::set_product(Matrix const &m1, Matrix const &m2)
188 {
189     assert(m1.cols()==m2.rows());
190     assert(cols()==m2.cols() && rows()==m1.rows());
191     
192     for (int i=0, j=0; dat->mult_ok(i,j);
193          dat->mult_next(i,j)) {
194         Real r=0.0;
195         for (int k = 0; k < m1.cols(); k++)
196             r += m1(i,k)*m2(k,j);
197         dat->elem(i,j)=r;
198     }
199 }
200
201 void
202 Matrix::insert_row(Vector v, int k)
203 {
204     assert(v.dim()==cols());
205     dat->insert_row(k);
206     for (int j=0; j < cols(); j++)
207         dat->elem(k,j)=v(j);
208 }
209
210
211 void
212 Matrix::swap_columns(int c1, int c2)
213 {
214     assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
215     for (int i=0; i< rows(); i++) {
216         Real r=dat->elem(i,c1);
217         dat->elem(i,c1) = dat->elem(i,c2);
218         dat->elem(i,c2)=r;
219     }
220 }
221
222 void
223 Matrix::swap_rows(int c1, int c2)
224 {
225     assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0);
226     for (int i=0; i< cols(); i++) {
227         Real r=dat->elem(c1,i);
228         dat->elem(c1,i) = dat->elem(c2,i);
229         dat->elem(c2,i)=r;
230     }
231 }
232
233
234 int
235 Matrix::dim() const
236 {
237     assert(cols() == rows());
238     return rows();
239 }
240