Mercurial > octave-nkf
annotate liboctave/CDiagMatrix.cc @ 8840:c690e3772583
support diagonal matrices in pinv
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 23 Feb 2009 14:54:56 +0100 |
parents | 20dfb885f877 |
children | eb63fbe60fab |
rev | line source |
---|---|
1993 | 1 // DiagMatrix manipulations. |
458 | 2 /* |
3 | |
7017 | 4 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2004, |
5 2005, 2007 John W. Eaton | |
458 | 6 |
7 This file is part of Octave. | |
8 | |
9 Octave is free software; you can redistribute it and/or modify it | |
10 under the terms of the GNU General Public License as published by the | |
7016 | 11 Free Software Foundation; either version 3 of the License, or (at your |
12 option) any later version. | |
458 | 13 |
14 Octave is distributed in the hope that it will be useful, but WITHOUT | |
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
17 for more details. | |
18 | |
19 You should have received a copy of the GNU General Public License | |
7016 | 20 along with Octave; see the file COPYING. If not, see |
21 <http://www.gnu.org/licenses/>. | |
458 | 22 |
23 */ | |
24 | |
25 #ifdef HAVE_CONFIG_H | |
1192 | 26 #include <config.h> |
458 | 27 #endif |
28 | |
3503 | 29 #include <iostream> |
458 | 30 |
4669 | 31 #include "Array-util.h" |
1368 | 32 #include "lo-error.h" |
458 | 33 #include "mx-base.h" |
34 #include "mx-inlines.cc" | |
1650 | 35 #include "oct-cmplx.h" |
458 | 36 |
1360 | 37 // Complex Diagonal Matrix class |
458 | 38 |
39 ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a) | |
1989 | 40 : MDiagArray2<Complex> (a.rows (), a.cols ()) |
458 | 41 { |
5275 | 42 for (octave_idx_type i = 0; i < length (); i++) |
458 | 43 elem (i, i) = a.elem (i, i); |
44 } | |
45 | |
2386 | 46 bool |
458 | 47 ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const |
48 { | |
49 if (rows () != a.rows () || cols () != a.cols ()) | |
50 return 0; | |
51 | |
3769 | 52 return mx_inline_equal (data (), a.data (), length ()); |
458 | 53 } |
54 | |
2386 | 55 bool |
458 | 56 ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const |
57 { | |
58 return !(*this == a); | |
59 } | |
60 | |
61 ComplexDiagMatrix& | |
62 ComplexDiagMatrix::fill (double val) | |
63 { | |
5275 | 64 for (octave_idx_type i = 0; i < length (); i++) |
458 | 65 elem (i, i) = val; |
66 return *this; | |
67 } | |
68 | |
69 ComplexDiagMatrix& | |
70 ComplexDiagMatrix::fill (const Complex& val) | |
71 { | |
5275 | 72 for (octave_idx_type i = 0; i < length (); i++) |
458 | 73 elem (i, i) = val; |
74 return *this; | |
75 } | |
76 | |
77 ComplexDiagMatrix& | |
5275 | 78 ComplexDiagMatrix::fill (double val, octave_idx_type beg, octave_idx_type end) |
458 | 79 { |
80 if (beg < 0 || end >= length () || end < beg) | |
81 { | |
82 (*current_liboctave_error_handler) ("range error for fill"); | |
83 return *this; | |
84 } | |
85 | |
5275 | 86 for (octave_idx_type i = beg; i <= end; i++) |
458 | 87 elem (i, i) = val; |
88 | |
89 return *this; | |
90 } | |
91 | |
92 ComplexDiagMatrix& | |
5275 | 93 ComplexDiagMatrix::fill (const Complex& val, octave_idx_type beg, octave_idx_type end) |
458 | 94 { |
95 if (beg < 0 || end >= length () || end < beg) | |
96 { | |
97 (*current_liboctave_error_handler) ("range error for fill"); | |
98 return *this; | |
99 } | |
100 | |
5275 | 101 for (octave_idx_type i = beg; i <= end; i++) |
458 | 102 elem (i, i) = val; |
103 | |
104 return *this; | |
105 } | |
106 | |
107 ComplexDiagMatrix& | |
108 ComplexDiagMatrix::fill (const ColumnVector& a) | |
109 { | |
5275 | 110 octave_idx_type len = length (); |
458 | 111 if (a.length () != len) |
112 { | |
113 (*current_liboctave_error_handler) ("range error for fill"); | |
114 return *this; | |
115 } | |
116 | |
5275 | 117 for (octave_idx_type i = 0; i < len; i++) |
458 | 118 elem (i, i) = a.elem (i); |
119 | |
120 return *this; | |
121 } | |
122 | |
123 ComplexDiagMatrix& | |
124 ComplexDiagMatrix::fill (const ComplexColumnVector& a) | |
125 { | |
5275 | 126 octave_idx_type len = length (); |
458 | 127 if (a.length () != len) |
128 { | |
129 (*current_liboctave_error_handler) ("range error for fill"); | |
130 return *this; | |
131 } | |
132 | |
5275 | 133 for (octave_idx_type i = 0; i < len; i++) |
458 | 134 elem (i, i) = a.elem (i); |
135 | |
136 return *this; | |
137 } | |
138 | |
139 ComplexDiagMatrix& | |
140 ComplexDiagMatrix::fill (const RowVector& a) | |
141 { | |
5275 | 142 octave_idx_type len = length (); |
458 | 143 if (a.length () != len) |
144 { | |
145 (*current_liboctave_error_handler) ("range error for fill"); | |
146 return *this; | |
147 } | |
148 | |
5275 | 149 for (octave_idx_type i = 0; i < len; i++) |
458 | 150 elem (i, i) = a.elem (i); |
151 | |
152 return *this; | |
153 } | |
154 | |
155 ComplexDiagMatrix& | |
156 ComplexDiagMatrix::fill (const ComplexRowVector& a) | |
157 { | |
5275 | 158 octave_idx_type len = length (); |
458 | 159 if (a.length () != len) |
160 { | |
161 (*current_liboctave_error_handler) ("range error for fill"); | |
162 return *this; | |
163 } | |
164 | |
5275 | 165 for (octave_idx_type i = 0; i < len; i++) |
458 | 166 elem (i, i) = a.elem (i); |
167 | |
168 return *this; | |
169 } | |
170 | |
171 ComplexDiagMatrix& | |
5275 | 172 ComplexDiagMatrix::fill (const ColumnVector& a, octave_idx_type beg) |
458 | 173 { |
5275 | 174 octave_idx_type a_len = a.length (); |
458 | 175 if (beg < 0 || beg + a_len >= length ()) |
176 { | |
177 (*current_liboctave_error_handler) ("range error for fill"); | |
178 return *this; | |
179 } | |
180 | |
5275 | 181 for (octave_idx_type i = 0; i < a_len; i++) |
458 | 182 elem (i+beg, i+beg) = a.elem (i); |
183 | |
184 return *this; | |
185 } | |
186 | |
187 ComplexDiagMatrix& | |
5275 | 188 ComplexDiagMatrix::fill (const ComplexColumnVector& a, octave_idx_type beg) |
458 | 189 { |
5275 | 190 octave_idx_type a_len = a.length (); |
458 | 191 if (beg < 0 || beg + a_len >= length ()) |
192 { | |
193 (*current_liboctave_error_handler) ("range error for fill"); | |
194 return *this; | |
195 } | |
196 | |
5275 | 197 for (octave_idx_type i = 0; i < a_len; i++) |
458 | 198 elem (i+beg, i+beg) = a.elem (i); |
199 | |
200 return *this; | |
201 } | |
202 | |
203 ComplexDiagMatrix& | |
5275 | 204 ComplexDiagMatrix::fill (const RowVector& a, octave_idx_type beg) |
458 | 205 { |
5275 | 206 octave_idx_type a_len = a.length (); |
458 | 207 if (beg < 0 || beg + a_len >= length ()) |
208 { | |
209 (*current_liboctave_error_handler) ("range error for fill"); | |
210 return *this; | |
211 } | |
212 | |
5275 | 213 for (octave_idx_type i = 0; i < a_len; i++) |
458 | 214 elem (i+beg, i+beg) = a.elem (i); |
215 | |
216 return *this; | |
217 } | |
218 | |
219 ComplexDiagMatrix& | |
5275 | 220 ComplexDiagMatrix::fill (const ComplexRowVector& a, octave_idx_type beg) |
458 | 221 { |
5275 | 222 octave_idx_type a_len = a.length (); |
458 | 223 if (beg < 0 || beg + a_len >= length ()) |
224 { | |
225 (*current_liboctave_error_handler) ("range error for fill"); | |
226 return *this; | |
227 } | |
228 | |
5275 | 229 for (octave_idx_type i = 0; i < a_len; i++) |
458 | 230 elem (i+beg, i+beg) = a.elem (i); |
231 | |
232 return *this; | |
233 } | |
234 | |
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
235 DiagMatrix |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
236 ComplexDiagMatrix::abs (void) const |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
237 { |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
238 DiagMatrix retval (rows (), cols ()); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
239 for (octave_idx_type i = 0; i < rows (); i++) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
240 retval(i, i) = std::abs (elem (i, i)); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
241 return retval; |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
242 } |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
243 |
458 | 244 ComplexDiagMatrix |
245 conj (const ComplexDiagMatrix& a) | |
246 { | |
247 ComplexDiagMatrix retval; | |
5275 | 248 octave_idx_type a_len = a.length (); |
458 | 249 if (a_len > 0) |
3769 | 250 retval = ComplexDiagMatrix (mx_inline_conj_dup (a.data (), a_len), |
458 | 251 a.rows (), a.cols ()); |
252 return retval; | |
253 } | |
254 | |
255 // resize is the destructive analog for this one | |
256 | |
257 ComplexMatrix | |
5275 | 258 ComplexDiagMatrix::extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const |
458 | 259 { |
5275 | 260 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } |
261 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } | |
458 | 262 |
5275 | 263 octave_idx_type new_r = r2 - r1 + 1; |
264 octave_idx_type new_c = c2 - c1 + 1; | |
458 | 265 |
266 ComplexMatrix result (new_r, new_c); | |
267 | |
5275 | 268 for (octave_idx_type j = 0; j < new_c; j++) |
269 for (octave_idx_type i = 0; i < new_r; i++) | |
458 | 270 result.elem (i, j) = elem (r1+i, c1+j); |
271 | |
272 return result; | |
273 } | |
274 | |
275 // extract row or column i. | |
276 | |
277 ComplexRowVector | |
5275 | 278 ComplexDiagMatrix::row (octave_idx_type i) const |
458 | 279 { |
5275 | 280 octave_idx_type r = rows (); |
281 octave_idx_type c = cols (); | |
3504 | 282 if (i < 0 || i >= r) |
458 | 283 { |
284 (*current_liboctave_error_handler) ("invalid row selection"); | |
3585 | 285 return ComplexRowVector (); |
458 | 286 } |
287 | |
3504 | 288 ComplexRowVector retval (c, 0.0); |
289 if (r <= c || (r > c && i < c)) | |
458 | 290 retval.elem (i) = elem (i, i); |
291 | |
292 return retval; | |
293 } | |
294 | |
295 ComplexRowVector | |
296 ComplexDiagMatrix::row (char *s) const | |
297 { | |
533 | 298 if (! s) |
458 | 299 { |
300 (*current_liboctave_error_handler) ("invalid row selection"); | |
301 return ComplexRowVector (); | |
302 } | |
303 | |
304 char c = *s; | |
305 if (c == 'f' || c == 'F') | |
5275 | 306 return row (static_cast<octave_idx_type>(0)); |
458 | 307 else if (c == 'l' || c == 'L') |
308 return row (rows () - 1); | |
309 else | |
310 { | |
311 (*current_liboctave_error_handler) ("invalid row selection"); | |
312 return ComplexRowVector (); | |
313 } | |
314 } | |
315 | |
316 ComplexColumnVector | |
5275 | 317 ComplexDiagMatrix::column (octave_idx_type i) const |
458 | 318 { |
5275 | 319 octave_idx_type r = rows (); |
320 octave_idx_type c = cols (); | |
3504 | 321 if (i < 0 || i >= c) |
458 | 322 { |
323 (*current_liboctave_error_handler) ("invalid column selection"); | |
3585 | 324 return ComplexColumnVector (); |
458 | 325 } |
326 | |
3504 | 327 ComplexColumnVector retval (r, 0.0); |
328 if (r >= c || (r < c && i < r)) | |
458 | 329 retval.elem (i) = elem (i, i); |
330 | |
331 return retval; | |
332 } | |
333 | |
334 ComplexColumnVector | |
335 ComplexDiagMatrix::column (char *s) const | |
336 { | |
533 | 337 if (! s) |
458 | 338 { |
339 (*current_liboctave_error_handler) ("invalid column selection"); | |
3585 | 340 return ComplexColumnVector (); |
458 | 341 } |
342 | |
343 char c = *s; | |
344 if (c == 'f' || c == 'F') | |
5275 | 345 return column (static_cast<octave_idx_type>(0)); |
458 | 346 else if (c == 'l' || c == 'L') |
347 return column (cols () - 1); | |
348 else | |
349 { | |
350 (*current_liboctave_error_handler) ("invalid column selection"); | |
3585 | 351 return ComplexColumnVector (); |
458 | 352 } |
353 } | |
354 | |
355 ComplexDiagMatrix | |
356 ComplexDiagMatrix::inverse (void) const | |
357 { | |
8811 | 358 octave_idx_type info; |
458 | 359 return inverse (info); |
360 } | |
361 | |
362 ComplexDiagMatrix | |
8811 | 363 ComplexDiagMatrix::inverse (octave_idx_type& info) const |
458 | 364 { |
5275 | 365 octave_idx_type r = rows (); |
366 octave_idx_type c = cols (); | |
3504 | 367 if (r != c) |
458 | 368 { |
369 (*current_liboctave_error_handler) ("inverse requires square matrix"); | |
3585 | 370 return ComplexDiagMatrix (); |
458 | 371 } |
372 | |
3504 | 373 ComplexDiagMatrix retval (r, c); |
458 | 374 |
375 info = 0; | |
5275 | 376 for (octave_idx_type i = 0; i < length (); i++) |
458 | 377 { |
378 if (elem (i, i) == 0.0) | |
379 { | |
380 info = -1; | |
381 return *this; | |
382 } | |
383 else | |
384 retval.elem (i, i) = 1.0 / elem (i, i); | |
385 } | |
386 | |
1627 | 387 return retval; |
458 | 388 } |
389 | |
8840
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
390 ComplexDiagMatrix |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
391 ComplexDiagMatrix::pseudo_inverse (void) const |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
392 { |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
393 octave_idx_type r = rows (); |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
394 octave_idx_type c = cols (); |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
395 octave_idx_type len = length (); |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
396 |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
397 ComplexDiagMatrix retval (c, r); |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
398 |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
399 for (octave_idx_type i = 0; i < len; i++) |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
400 { |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
401 if (elem (i, i) != 0.0) |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
402 retval.elem (i, i) = 1.0 / elem (i, i); |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
403 else |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
404 retval.elem (i, i) = 0.0; |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
405 } |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
406 |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
407 return retval; |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
408 } |
c690e3772583
support diagonal matrices in pinv
Jaroslav Hajek <highegg@gmail.com>
parents:
8811
diff
changeset
|
409 |
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
410 bool |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
411 ComplexDiagMatrix::all_elements_are_real (void) const |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
412 { |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
413 octave_idx_type len = length (); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
414 for (octave_idx_type i = 0; i < len; i++) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
415 { |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
416 double ip = std::imag (elem (i, i)); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
417 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
418 if (ip != 0.0 || lo_ieee_signbit (ip)) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
419 return false; |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
420 } |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
421 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
422 return true; |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
423 } |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
424 |
458 | 425 // diagonal matrix by diagonal matrix -> diagonal matrix operations |
426 | |
427 ComplexDiagMatrix& | |
428 ComplexDiagMatrix::operator += (const DiagMatrix& a) | |
429 { | |
5275 | 430 octave_idx_type r = rows (); |
431 octave_idx_type c = cols (); | |
2386 | 432 |
5275 | 433 octave_idx_type a_nr = a.rows (); |
434 octave_idx_type a_nc = a.cols (); | |
2386 | 435 |
3504 | 436 if (r != a_nr || c != a_nc) |
458 | 437 { |
3504 | 438 gripe_nonconformant ("operator +=", r, c, a_nr, a_nc); |
458 | 439 return *this; |
440 } | |
441 | |
3504 | 442 if (r == 0 || c == 0) |
458 | 443 return *this; |
444 | |
445 Complex *d = fortran_vec (); // Ensures only one reference to my privates! | |
446 | |
3769 | 447 mx_inline_add2 (d, a.data (), length ()); |
458 | 448 return *this; |
449 } | |
450 | |
451 ComplexDiagMatrix | |
452 operator * (const ComplexDiagMatrix& a, const DiagMatrix& b) | |
453 { | |
5275 | 454 octave_idx_type a_nr = a.rows (); |
455 octave_idx_type a_nc = a.cols (); | |
2386 | 456 |
5275 | 457 octave_idx_type b_nr = b.rows (); |
458 octave_idx_type b_nc = b.cols (); | |
2386 | 459 |
3504 | 460 if (a_nc != b_nr) |
458 | 461 { |
3504 | 462 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
458 | 463 return ComplexDiagMatrix (); |
464 } | |
465 | |
3504 | 466 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
467 return ComplexDiagMatrix (a_nr, a_nc, 0.0); | |
458 | 468 |
3504 | 469 ComplexDiagMatrix c (a_nr, b_nc); |
458 | 470 |
5275 | 471 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
458 | 472 |
5275 | 473 for (octave_idx_type i = 0; i < len; i++) |
458 | 474 { |
475 Complex a_element = a.elem (i, i); | |
476 double b_element = b.elem (i, i); | |
477 | |
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
478 c.elem (i, i) = a_element * b_element; |
458 | 479 } |
480 | |
481 return c; | |
482 } | |
483 | |
484 ComplexDiagMatrix | |
1205 | 485 operator * (const DiagMatrix& a, const ComplexDiagMatrix& b) |
486 { | |
5275 | 487 octave_idx_type a_nr = a.rows (); |
488 octave_idx_type a_nc = a.cols (); | |
2386 | 489 |
5275 | 490 octave_idx_type b_nr = b.rows (); |
491 octave_idx_type b_nc = b.cols (); | |
2386 | 492 |
3504 | 493 if (a_nc != b_nr) |
1205 | 494 { |
3504 | 495 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
1205 | 496 return ComplexDiagMatrix (); |
497 } | |
498 | |
3504 | 499 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
500 return ComplexDiagMatrix (a_nr, a_nc, 0.0); | |
1205 | 501 |
3504 | 502 ComplexDiagMatrix c (a_nr, b_nc); |
1205 | 503 |
5275 | 504 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
1205 | 505 |
5275 | 506 for (octave_idx_type i = 0; i < len; i++) |
1205 | 507 { |
508 double a_element = a.elem (i, i); | |
509 Complex b_element = b.elem (i, i); | |
510 | |
8366
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
511 c.elem (i, i) = a_element * b_element; |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
512 } |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
513 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
514 return c; |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
515 } |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
516 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
517 ComplexDiagMatrix |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
518 operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
519 { |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
520 octave_idx_type a_nr = a.rows (); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
521 octave_idx_type a_nc = a.cols (); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
522 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
523 octave_idx_type b_nr = b.rows (); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
524 octave_idx_type b_nc = b.cols (); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
525 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
526 if (a_nc != b_nr) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
527 { |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
528 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
529 return ComplexDiagMatrix (); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
530 } |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
531 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
532 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
533 return ComplexDiagMatrix (a_nr, a_nc, 0.0); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
534 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
535 ComplexDiagMatrix c (a_nr, b_nc); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
536 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
537 octave_idx_type len = a_nr < b_nc ? a_nr : b_nc; |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
538 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
539 for (octave_idx_type i = 0; i < len; i++) |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
540 { |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
541 Complex a_element = a.elem (i, i); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
542 Complex b_element = b.elem (i, i); |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
543 |
8b1a2555c4e2
implement diagonal matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
544 c.elem (i, i) = a_element * b_element; |
1205 | 545 } |
546 | |
547 return c; | |
548 } | |
549 | |
458 | 550 // other operations |
551 | |
8371
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
552 ComplexDET |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
553 ComplexDiagMatrix::determinant (void) const |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
554 { |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
555 ComplexDET det (1.0); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
556 if (rows () != cols ()) |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
557 { |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
558 (*current_liboctave_error_handler) ("determinant requires square matrix"); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
559 det = ComplexDET (0.0); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
560 } |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
561 else |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
562 { |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
563 octave_idx_type len = length (); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
564 for (octave_idx_type i = 0; i < len; i++) |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
565 det *= elem (i, i); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
566 } |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
567 |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
568 return det; |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
569 } |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
570 |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
571 double |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
572 ComplexDiagMatrix::rcond (void) const |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
573 { |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
574 ColumnVector av = diag (0).map (std::abs); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
575 double amx = av.max (), amn = av.min (); |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
576 return amx == 0 ? 0.0 : amn / amx; |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
577 } |
c3f7e2549abb
make det & inv aware of diagonal & permutation matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
8366
diff
changeset
|
578 |
458 | 579 // i/o |
580 | |
3504 | 581 std::ostream& |
582 operator << (std::ostream& os, const ComplexDiagMatrix& a) | |
458 | 583 { |
584 Complex ZERO (0.0); | |
585 // int field_width = os.precision () + 7; | |
5275 | 586 for (octave_idx_type i = 0; i < a.rows (); i++) |
458 | 587 { |
5275 | 588 for (octave_idx_type j = 0; j < a.cols (); j++) |
458 | 589 { |
590 if (i == j) | |
591 os << " " /* setw (field_width) */ << a.elem (i, i); | |
592 else | |
593 os << " " /* setw (field_width) */ << ZERO; | |
594 } | |
595 os << "\n"; | |
596 } | |
597 return os; | |
598 } | |
599 | |
600 /* | |
601 ;;; Local Variables: *** | |
602 ;;; mode: C++ *** | |
603 ;;; End: *** | |
604 */ |