|
458
|
1 // DiagMatrix manipulations. -*- C++ -*- |
|
|
2 /* |
|
|
3 |
|
|
4 Copyright (C) 1992, 1993, 1994 John W. Eaton |
|
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
|
11 later version. |
|
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
|
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include "config.h" |
|
|
26 #endif |
|
|
27 |
|
|
28 #include <iostream.h> |
|
|
29 |
|
|
30 #include <Complex.h> |
|
|
31 |
|
|
32 #include "mx-base.h" |
|
|
33 #include "mx-inlines.cc" |
|
|
34 #include "lo-error.h" |
|
|
35 |
|
|
36 /* |
|
|
37 * Complex Diagonal Matrix class |
|
|
38 */ |
|
|
39 |
|
|
40 #define KLUDGE_DIAG_MATRICES |
|
|
41 #define TYPE Complex |
|
|
42 #define KL_DMAT_TYPE ComplexDiagMatrix |
|
|
43 #include "mx-kludge.cc" |
|
|
44 #undef KLUDGE_DIAG_MATRICES |
|
|
45 #undef TYPE |
|
|
46 #undef KL_DMAT_TYPE |
|
|
47 |
|
|
48 ComplexDiagMatrix::ComplexDiagMatrix (const RowVector& a) |
|
|
49 : DiagArray<Complex> (a.length ()) |
|
|
50 { |
|
|
51 for (int i = 0; i < length (); i++) |
|
|
52 elem (i, i) = a.elem (i); |
|
|
53 } |
|
|
54 |
|
|
55 ComplexDiagMatrix::ComplexDiagMatrix (const ColumnVector& a) |
|
|
56 : DiagArray<Complex> (a.length ()) |
|
|
57 { |
|
|
58 for (int i = 0; i < length (); i++) |
|
|
59 elem (i, i) = a.elem (i); |
|
|
60 } |
|
|
61 |
|
|
62 ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a) |
|
|
63 : DiagArray<Complex> (a.rows (), a.cols ()) |
|
|
64 { |
|
|
65 for (int i = 0; i < length (); i++) |
|
|
66 elem (i, i) = a.elem (i, i); |
|
|
67 } |
|
|
68 |
|
|
69 #if 0 |
|
|
70 ComplexDiagMatrix& |
|
|
71 ComplexDiagMatrix::resize (int r, int c) |
|
|
72 { |
|
|
73 if (r < 0 || c < 0) |
|
|
74 { |
|
|
75 (*current_liboctave_error_handler) |
|
|
76 ("can't resize to negative dimensions"); |
|
|
77 return *this; |
|
|
78 } |
|
|
79 |
|
|
80 int new_len = r < c ? r : c; |
|
533
|
81 Complex *new_data = 0; |
|
458
|
82 if (new_len > 0) |
|
|
83 { |
|
|
84 new_data = new Complex [new_len]; |
|
|
85 |
|
|
86 int min_len = new_len < len ? new_len : len; |
|
|
87 |
|
|
88 for (int i = 0; i < min_len; i++) |
|
|
89 new_data[i] = data[i]; |
|
|
90 } |
|
|
91 |
|
|
92 delete [] data; |
|
|
93 nr = r; |
|
|
94 nc = c; |
|
|
95 len = new_len; |
|
|
96 data = new_data; |
|
|
97 |
|
|
98 return *this; |
|
|
99 } |
|
|
100 |
|
|
101 ComplexDiagMatrix& |
|
|
102 ComplexDiagMatrix::resize (int r, int c, double val) |
|
|
103 { |
|
|
104 if (r < 0 || c < 0) |
|
|
105 { |
|
|
106 (*current_liboctave_error_handler) |
|
|
107 ("can't resize to negative dimensions"); |
|
|
108 return *this; |
|
|
109 } |
|
|
110 |
|
|
111 int new_len = r < c ? r : c; |
|
533
|
112 Complex *new_data = 0; |
|
458
|
113 if (new_len > 0) |
|
|
114 { |
|
|
115 new_data = new Complex [new_len]; |
|
|
116 |
|
|
117 int min_len = new_len < len ? new_len : len; |
|
|
118 |
|
|
119 for (int i = 0; i < min_len; i++) |
|
|
120 new_data[i] = data[i]; |
|
|
121 |
|
|
122 for (i = min_len; i < new_len; i++) |
|
|
123 new_data[i] = val; |
|
|
124 } |
|
|
125 |
|
|
126 delete [] data; |
|
|
127 nr = r; |
|
|
128 nc = c; |
|
|
129 len = new_len; |
|
|
130 data = new_data; |
|
|
131 |
|
|
132 return *this; |
|
|
133 } |
|
|
134 |
|
|
135 ComplexDiagMatrix& |
|
|
136 ComplexDiagMatrix::resize (int r, int c, const Complex& val) |
|
|
137 { |
|
|
138 if (r < 0 || c < 0) |
|
|
139 { |
|
|
140 (*current_liboctave_error_handler) |
|
|
141 ("can't resize to negative dimensions"); |
|
|
142 return *this; |
|
|
143 } |
|
|
144 |
|
|
145 int new_len = r < c ? r : c; |
|
533
|
146 Complex *new_data = 0; |
|
458
|
147 if (new_len > 0) |
|
|
148 { |
|
|
149 new_data = new Complex [new_len]; |
|
|
150 |
|
|
151 int min_len = new_len < len ? new_len : len; |
|
|
152 |
|
|
153 for (int i = 0; i < min_len; i++) |
|
|
154 new_data[i] = data[i]; |
|
|
155 |
|
|
156 for (i = min_len; i < new_len; i++) |
|
|
157 new_data[i] = val; |
|
|
158 } |
|
|
159 |
|
|
160 delete [] data; |
|
|
161 nr = r; |
|
|
162 nc = c; |
|
|
163 len = new_len; |
|
|
164 data = new_data; |
|
|
165 |
|
|
166 return *this; |
|
|
167 } |
|
|
168 #endif |
|
|
169 |
|
|
170 int |
|
|
171 ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const |
|
|
172 { |
|
|
173 if (rows () != a.rows () || cols () != a.cols ()) |
|
|
174 return 0; |
|
|
175 |
|
|
176 return equal (data (), a.data (), length ()); |
|
|
177 } |
|
|
178 |
|
|
179 int |
|
|
180 ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const |
|
|
181 { |
|
|
182 return !(*this == a); |
|
|
183 } |
|
|
184 |
|
|
185 ComplexDiagMatrix |
|
|
186 ComplexDiagMatrix::hermitian (void) const |
|
|
187 { |
|
|
188 return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ()); |
|
|
189 } |
|
|
190 |
|
|
191 ComplexDiagMatrix& |
|
|
192 ComplexDiagMatrix::fill (double val) |
|
|
193 { |
|
|
194 for (int i = 0; i < length (); i++) |
|
|
195 elem (i, i) = val; |
|
|
196 return *this; |
|
|
197 } |
|
|
198 |
|
|
199 ComplexDiagMatrix& |
|
|
200 ComplexDiagMatrix::fill (const Complex& val) |
|
|
201 { |
|
|
202 for (int i = 0; i < length (); i++) |
|
|
203 elem (i, i) = val; |
|
|
204 return *this; |
|
|
205 } |
|
|
206 |
|
|
207 ComplexDiagMatrix& |
|
|
208 ComplexDiagMatrix::fill (double val, int beg, int end) |
|
|
209 { |
|
|
210 if (beg < 0 || end >= length () || end < beg) |
|
|
211 { |
|
|
212 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
213 return *this; |
|
|
214 } |
|
|
215 |
|
|
216 for (int i = beg; i < end; i++) |
|
|
217 elem (i, i) = val; |
|
|
218 |
|
|
219 return *this; |
|
|
220 } |
|
|
221 |
|
|
222 ComplexDiagMatrix& |
|
|
223 ComplexDiagMatrix::fill (const Complex& val, int beg, int end) |
|
|
224 { |
|
|
225 if (beg < 0 || end >= length () || end < beg) |
|
|
226 { |
|
|
227 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
228 return *this; |
|
|
229 } |
|
|
230 |
|
|
231 for (int i = beg; i < end; i++) |
|
|
232 elem (i, i) = val; |
|
|
233 |
|
|
234 return *this; |
|
|
235 } |
|
|
236 |
|
|
237 ComplexDiagMatrix& |
|
|
238 ComplexDiagMatrix::fill (const ColumnVector& a) |
|
|
239 { |
|
|
240 int len = length (); |
|
|
241 if (a.length () != len) |
|
|
242 { |
|
|
243 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
244 return *this; |
|
|
245 } |
|
|
246 |
|
|
247 for (int i = 0; i < len; i++) |
|
|
248 elem (i, i) = a.elem (i); |
|
|
249 |
|
|
250 return *this; |
|
|
251 } |
|
|
252 |
|
|
253 ComplexDiagMatrix& |
|
|
254 ComplexDiagMatrix::fill (const ComplexColumnVector& a) |
|
|
255 { |
|
|
256 int len = length (); |
|
|
257 if (a.length () != len) |
|
|
258 { |
|
|
259 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
260 return *this; |
|
|
261 } |
|
|
262 |
|
|
263 for (int i = 0; i < len; i++) |
|
|
264 elem (i, i) = a.elem (i); |
|
|
265 |
|
|
266 return *this; |
|
|
267 } |
|
|
268 |
|
|
269 ComplexDiagMatrix& |
|
|
270 ComplexDiagMatrix::fill (const RowVector& a) |
|
|
271 { |
|
|
272 int len = length (); |
|
|
273 if (a.length () != len) |
|
|
274 { |
|
|
275 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
276 return *this; |
|
|
277 } |
|
|
278 |
|
|
279 for (int i = 0; i < len; i++) |
|
|
280 elem (i, i) = a.elem (i); |
|
|
281 |
|
|
282 return *this; |
|
|
283 } |
|
|
284 |
|
|
285 ComplexDiagMatrix& |
|
|
286 ComplexDiagMatrix::fill (const ComplexRowVector& a) |
|
|
287 { |
|
|
288 int len = length (); |
|
|
289 if (a.length () != len) |
|
|
290 { |
|
|
291 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
292 return *this; |
|
|
293 } |
|
|
294 |
|
|
295 for (int i = 0; i < len; i++) |
|
|
296 elem (i, i) = a.elem (i); |
|
|
297 |
|
|
298 return *this; |
|
|
299 } |
|
|
300 |
|
|
301 ComplexDiagMatrix& |
|
|
302 ComplexDiagMatrix::fill (const ColumnVector& a, int beg) |
|
|
303 { |
|
|
304 int a_len = a.length (); |
|
|
305 if (beg < 0 || beg + a_len >= length ()) |
|
|
306 { |
|
|
307 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
308 return *this; |
|
|
309 } |
|
|
310 |
|
|
311 for (int i = 0; i < a_len; i++) |
|
|
312 elem (i+beg, i+beg) = a.elem (i); |
|
|
313 |
|
|
314 return *this; |
|
|
315 } |
|
|
316 |
|
|
317 ComplexDiagMatrix& |
|
|
318 ComplexDiagMatrix::fill (const ComplexColumnVector& a, int beg) |
|
|
319 { |
|
|
320 int a_len = a.length (); |
|
|
321 if (beg < 0 || beg + a_len >= length ()) |
|
|
322 { |
|
|
323 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
324 return *this; |
|
|
325 } |
|
|
326 |
|
|
327 for (int i = 0; i < a_len; i++) |
|
|
328 elem (i+beg, i+beg) = a.elem (i); |
|
|
329 |
|
|
330 return *this; |
|
|
331 } |
|
|
332 |
|
|
333 ComplexDiagMatrix& |
|
|
334 ComplexDiagMatrix::fill (const RowVector& a, int beg) |
|
|
335 { |
|
|
336 int a_len = a.length (); |
|
|
337 if (beg < 0 || beg + a_len >= length ()) |
|
|
338 { |
|
|
339 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
340 return *this; |
|
|
341 } |
|
|
342 |
|
|
343 for (int i = 0; i < a_len; i++) |
|
|
344 elem (i+beg, i+beg) = a.elem (i); |
|
|
345 |
|
|
346 return *this; |
|
|
347 } |
|
|
348 |
|
|
349 ComplexDiagMatrix& |
|
|
350 ComplexDiagMatrix::fill (const ComplexRowVector& a, int beg) |
|
|
351 { |
|
|
352 int a_len = a.length (); |
|
|
353 if (beg < 0 || beg + a_len >= length ()) |
|
|
354 { |
|
|
355 (*current_liboctave_error_handler) ("range error for fill"); |
|
|
356 return *this; |
|
|
357 } |
|
|
358 |
|
|
359 for (int i = 0; i < a_len; i++) |
|
|
360 elem (i+beg, i+beg) = a.elem (i); |
|
|
361 |
|
|
362 return *this; |
|
|
363 } |
|
|
364 |
|
|
365 ComplexDiagMatrix |
|
|
366 ComplexDiagMatrix::transpose (void) const |
|
|
367 { |
|
|
368 return ComplexDiagMatrix (dup (data (), length ()), cols (), rows ()); |
|
|
369 } |
|
|
370 |
|
|
371 DiagMatrix |
|
|
372 real (const ComplexDiagMatrix& a) |
|
|
373 { |
|
|
374 DiagMatrix retval; |
|
|
375 int a_len = a.length (); |
|
|
376 if (a_len > 0) |
|
|
377 retval = DiagMatrix (real_dup (a.data (), a_len), a.rows (), |
|
|
378 a.cols ()); |
|
|
379 return retval; |
|
|
380 } |
|
|
381 |
|
|
382 DiagMatrix |
|
|
383 imag (const ComplexDiagMatrix& a) |
|
|
384 { |
|
|
385 DiagMatrix retval; |
|
|
386 int a_len = a.length (); |
|
|
387 if (a_len > 0) |
|
|
388 retval = DiagMatrix (imag_dup (a.data (), a_len), a.rows (), |
|
|
389 a.cols ()); |
|
|
390 return retval; |
|
|
391 } |
|
|
392 |
|
|
393 ComplexDiagMatrix |
|
|
394 conj (const ComplexDiagMatrix& a) |
|
|
395 { |
|
|
396 ComplexDiagMatrix retval; |
|
|
397 int a_len = a.length (); |
|
|
398 if (a_len > 0) |
|
|
399 retval = ComplexDiagMatrix (conj_dup (a.data (), a_len), |
|
|
400 a.rows (), a.cols ()); |
|
|
401 return retval; |
|
|
402 } |
|
|
403 |
|
|
404 // resize is the destructive analog for this one |
|
|
405 |
|
|
406 ComplexMatrix |
|
|
407 ComplexDiagMatrix::extract (int r1, int c1, int r2, int c2) const |
|
|
408 { |
|
|
409 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
|
410 if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } |
|
|
411 |
|
|
412 int new_r = r2 - r1 + 1; |
|
|
413 int new_c = c2 - c1 + 1; |
|
|
414 |
|
|
415 ComplexMatrix result (new_r, new_c); |
|
|
416 |
|
|
417 for (int j = 0; j < new_c; j++) |
|
|
418 for (int i = 0; i < new_r; i++) |
|
|
419 result.elem (i, j) = elem (r1+i, c1+j); |
|
|
420 |
|
|
421 return result; |
|
|
422 } |
|
|
423 |
|
|
424 // extract row or column i. |
|
|
425 |
|
|
426 ComplexRowVector |
|
|
427 ComplexDiagMatrix::row (int i) const |
|
|
428 { |
|
|
429 int nr = rows (); |
|
|
430 int nc = cols (); |
|
|
431 if (i < 0 || i >= nr) |
|
|
432 { |
|
|
433 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
434 return RowVector (); |
|
|
435 } |
|
|
436 |
|
|
437 ComplexRowVector retval (nc, 0.0); |
|
|
438 if (nr <= nc || (nr > nc && i < nc)) |
|
|
439 retval.elem (i) = elem (i, i); |
|
|
440 |
|
|
441 return retval; |
|
|
442 } |
|
|
443 |
|
|
444 ComplexRowVector |
|
|
445 ComplexDiagMatrix::row (char *s) const |
|
|
446 { |
|
533
|
447 if (! s) |
|
458
|
448 { |
|
|
449 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
450 return ComplexRowVector (); |
|
|
451 } |
|
|
452 |
|
|
453 char c = *s; |
|
|
454 if (c == 'f' || c == 'F') |
|
|
455 return row (0); |
|
|
456 else if (c == 'l' || c == 'L') |
|
|
457 return row (rows () - 1); |
|
|
458 else |
|
|
459 { |
|
|
460 (*current_liboctave_error_handler) ("invalid row selection"); |
|
|
461 return ComplexRowVector (); |
|
|
462 } |
|
|
463 } |
|
|
464 |
|
|
465 ComplexColumnVector |
|
|
466 ComplexDiagMatrix::column (int i) const |
|
|
467 { |
|
|
468 int nr = rows (); |
|
|
469 int nc = cols (); |
|
|
470 if (i < 0 || i >= nc) |
|
|
471 { |
|
|
472 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
473 return ColumnVector (); |
|
|
474 } |
|
|
475 |
|
|
476 ComplexColumnVector retval (nr, 0.0); |
|
|
477 if (nr >= nc || (nr < nc && i < nr)) |
|
|
478 retval.elem (i) = elem (i, i); |
|
|
479 |
|
|
480 return retval; |
|
|
481 } |
|
|
482 |
|
|
483 ComplexColumnVector |
|
|
484 ComplexDiagMatrix::column (char *s) const |
|
|
485 { |
|
533
|
486 if (! s) |
|
458
|
487 { |
|
|
488 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
489 return ColumnVector (); |
|
|
490 } |
|
|
491 |
|
|
492 char c = *s; |
|
|
493 if (c == 'f' || c == 'F') |
|
|
494 return column (0); |
|
|
495 else if (c == 'l' || c == 'L') |
|
|
496 return column (cols () - 1); |
|
|
497 else |
|
|
498 { |
|
|
499 (*current_liboctave_error_handler) ("invalid column selection"); |
|
|
500 return ColumnVector (); |
|
|
501 } |
|
|
502 } |
|
|
503 |
|
|
504 ComplexDiagMatrix |
|
|
505 ComplexDiagMatrix::inverse (void) const |
|
|
506 { |
|
|
507 int info; |
|
|
508 return inverse (info); |
|
|
509 } |
|
|
510 |
|
|
511 ComplexDiagMatrix |
|
|
512 ComplexDiagMatrix::inverse (int& info) const |
|
|
513 { |
|
|
514 int nr = rows (); |
|
|
515 int nc = cols (); |
|
|
516 if (nr != nc) |
|
|
517 { |
|
|
518 (*current_liboctave_error_handler) ("inverse requires square matrix"); |
|
|
519 return DiagMatrix (); |
|
|
520 } |
|
|
521 |
|
|
522 ComplexDiagMatrix retval (nr, nc); |
|
|
523 |
|
|
524 info = 0; |
|
|
525 for (int i = 0; i < length (); i++) |
|
|
526 { |
|
|
527 if (elem (i, i) == 0.0) |
|
|
528 { |
|
|
529 info = -1; |
|
|
530 return *this; |
|
|
531 } |
|
|
532 else |
|
|
533 retval.elem (i, i) = 1.0 / elem (i, i); |
|
|
534 } |
|
|
535 |
|
|
536 return *this; |
|
|
537 } |
|
|
538 |
|
|
539 // diagonal matrix by diagonal matrix -> diagonal matrix operations |
|
|
540 |
|
|
541 ComplexDiagMatrix& |
|
|
542 ComplexDiagMatrix::operator += (const DiagMatrix& a) |
|
|
543 { |
|
|
544 int nr = rows (); |
|
|
545 int nc = cols (); |
|
|
546 if (nr != a.rows () || nc != a.cols ()) |
|
|
547 { |
|
|
548 (*current_liboctave_error_handler) |
|
|
549 ("nonconformant matrix += operation attempted"); |
|
|
550 return *this; |
|
|
551 } |
|
|
552 |
|
|
553 if (nr == 0 || nc == 0) |
|
|
554 return *this; |
|
|
555 |
|
|
556 Complex *d = fortran_vec (); // Ensures only one reference to my privates! |
|
|
557 |
|
|
558 add2 (d, a.data (), length ()); |
|
|
559 return *this; |
|
|
560 } |
|
|
561 |
|
|
562 ComplexDiagMatrix& |
|
|
563 ComplexDiagMatrix::operator -= (const DiagMatrix& a) |
|
|
564 { |
|
|
565 int nr = rows (); |
|
|
566 int nc = cols (); |
|
|
567 if (nr != a.rows () || nc != a.cols ()) |
|
|
568 { |
|
|
569 (*current_liboctave_error_handler) |
|
|
570 ("nonconformant matrix -= operation attempted"); |
|
|
571 return *this; |
|
|
572 } |
|
|
573 |
|
|
574 if (nr == 0 || nc == 0) |
|
|
575 return *this; |
|
|
576 |
|
|
577 Complex *d = fortran_vec (); // Ensures only one reference to my privates! |
|
|
578 |
|
|
579 subtract2 (d, a.data (), length ()); |
|
|
580 return *this; |
|
|
581 } |
|
|
582 |
|
|
583 ComplexDiagMatrix& |
|
|
584 ComplexDiagMatrix::operator += (const ComplexDiagMatrix& a) |
|
|
585 { |
|
|
586 int nr = rows (); |
|
|
587 int nc = cols (); |
|
|
588 if (nr != a.rows () || nc != a.cols ()) |
|
|
589 { |
|
|
590 (*current_liboctave_error_handler) |
|
|
591 ("nonconformant matrix += operation attempted"); |
|
|
592 return *this; |
|
|
593 } |
|
|
594 |
|
|
595 if (nr == 0 || nc == 0) |
|
|
596 return *this; |
|
|
597 |
|
|
598 Complex *d = fortran_vec (); // Ensures only one reference to my privates! |
|
|
599 |
|
|
600 add2 (d, a.data (), length ()); |
|
|
601 return *this; |
|
|
602 } |
|
|
603 |
|
|
604 ComplexDiagMatrix& |
|
|
605 ComplexDiagMatrix::operator -= (const ComplexDiagMatrix& a) |
|
|
606 { |
|
|
607 int nr = rows (); |
|
|
608 int nc = cols (); |
|
|
609 if (nr != a.rows () || nc != a.cols ()) |
|
|
610 { |
|
|
611 (*current_liboctave_error_handler) |
|
|
612 ("nonconformant matrix -= operation attempted"); |
|
|
613 return *this; |
|
|
614 } |
|
|
615 |
|
|
616 if (nr == 0 || nc == 0) |
|
|
617 return *this; |
|
|
618 |
|
|
619 Complex *d = fortran_vec (); // Ensures only one reference to my privates! |
|
|
620 |
|
|
621 subtract2 (d, a.data (), length ()); |
|
|
622 return *this; |
|
|
623 } |
|
|
624 |
|
|
625 // diagonal matrix by scalar -> matrix operations |
|
|
626 |
|
|
627 ComplexMatrix |
|
|
628 operator + (const ComplexDiagMatrix& a, double s) |
|
|
629 { |
|
|
630 ComplexMatrix tmp (a.rows (), a.cols (), s); |
|
|
631 return a + tmp; |
|
|
632 } |
|
|
633 |
|
|
634 ComplexMatrix |
|
|
635 operator - (const ComplexDiagMatrix& a, double s) |
|
|
636 { |
|
|
637 ComplexMatrix tmp (a.rows (), a.cols (), -s); |
|
|
638 return a + tmp; |
|
|
639 } |
|
|
640 |
|
|
641 ComplexMatrix |
|
|
642 operator + (const ComplexDiagMatrix& a, const Complex& s) |
|
|
643 { |
|
|
644 ComplexMatrix tmp (a.rows (), a.cols (), s); |
|
|
645 return a + tmp; |
|
|
646 } |
|
|
647 |
|
|
648 ComplexMatrix |
|
|
649 operator - (const ComplexDiagMatrix& a, const Complex& s) |
|
|
650 { |
|
|
651 ComplexMatrix tmp (a.rows (), a.cols (), -s); |
|
|
652 return a + tmp; |
|
|
653 } |
|
|
654 |
|
|
655 // diagonal matrix by scalar -> diagonal matrix operations |
|
|
656 |
|
|
657 ComplexDiagMatrix |
|
|
658 operator * (const ComplexDiagMatrix& a, double s) |
|
|
659 { |
|
|
660 return ComplexDiagMatrix (multiply (a.data (), a.length (), s), |
|
|
661 a.rows (), a.cols ()); |
|
|
662 } |
|
|
663 |
|
|
664 ComplexDiagMatrix |
|
|
665 operator / (const ComplexDiagMatrix& a, double s) |
|
|
666 { |
|
|
667 return ComplexDiagMatrix (divide (a.data (), a.length (), s), |
|
|
668 a.rows (), a.cols ()); |
|
|
669 } |
|
|
670 |
|
|
671 // scalar by diagonal matrix -> matrix operations |
|
|
672 |
|
|
673 ComplexMatrix |
|
|
674 operator + (double s, const ComplexDiagMatrix& a) |
|
|
675 { |
|
|
676 ComplexMatrix tmp (a.rows (), a.cols (), s); |
|
|
677 return tmp + a; |
|
|
678 } |
|
|
679 |
|
|
680 ComplexMatrix |
|
|
681 operator - (double s, const ComplexDiagMatrix& a) |
|
|
682 { |
|
|
683 ComplexMatrix tmp (a.rows (), a.cols (), s); |
|
|
684 return tmp - a; |
|
|
685 } |
|
|
686 |
|
|
687 ComplexMatrix |
|
|
688 operator + (const Complex& s, const ComplexDiagMatrix& a) |
|
|
689 { |
|
|
690 ComplexMatrix tmp (a.rows (), a.cols (), s); |
|
|
691 return tmp + a; |
|
|
692 } |
|
|
693 |
|
|
694 ComplexMatrix |
|
|
695 operator - (const Complex& s, const ComplexDiagMatrix& a) |
|
|
696 { |
|
|
697 ComplexMatrix tmp (a.rows (), a.cols (), s); |
|
|
698 return tmp - a; |
|
|
699 } |
|
|
700 |
|
|
701 // scalar by diagonal matrix -> diagonal matrix operations |
|
|
702 |
|
|
703 ComplexDiagMatrix |
|
|
704 operator * (double s, const ComplexDiagMatrix& a) |
|
|
705 { |
|
|
706 return ComplexDiagMatrix (multiply (a.data (), a.length (), s), |
|
|
707 a.rows (), a.cols ()); |
|
|
708 } |
|
|
709 |
|
|
710 // diagonal matrix by column vector -> column vector operations |
|
|
711 |
|
|
712 ComplexColumnVector |
|
|
713 operator * (const ComplexDiagMatrix& m, const ColumnVector& a) |
|
|
714 { |
|
|
715 int nr = m.rows (); |
|
|
716 int nc = m.cols (); |
|
|
717 int a_len = a.length (); |
|
|
718 if (nc != a_len) |
|
|
719 { |
|
|
720 (*current_liboctave_error_handler) |
|
|
721 ("nonconformant matrix muliplication attempted"); |
|
|
722 return ComplexColumnVector (); |
|
|
723 } |
|
|
724 |
|
|
725 if (nc == 0 || nr == 0) |
|
|
726 return ComplexColumnVector (0); |
|
|
727 |
|
|
728 ComplexColumnVector result (nr); |
|
|
729 |
|
|
730 for (int i = 0; i < a_len; i++) |
|
|
731 result.elem (i) = a.elem (i) * m.elem (i, i); |
|
|
732 |
|
|
733 for (i = a_len; i < nr; i++) |
|
|
734 result.elem (i) = 0.0; |
|
|
735 |
|
|
736 return result; |
|
|
737 } |
|
|
738 |
|
|
739 ComplexColumnVector |
|
|
740 operator * (const ComplexDiagMatrix& m, const ComplexColumnVector& a) |
|
|
741 { |
|
|
742 int nr = m.rows (); |
|
|
743 int nc = m.cols (); |
|
|
744 int a_len = a.length (); |
|
|
745 if (nc != a_len) |
|
|
746 { |
|
|
747 (*current_liboctave_error_handler) |
|
|
748 ("nonconformant matrix muliplication attempted"); |
|
|
749 return ComplexColumnVector (); |
|
|
750 } |
|
|
751 |
|
|
752 if (nc == 0 || nr == 0) |
|
|
753 return ComplexColumnVector (0); |
|
|
754 |
|
|
755 ComplexColumnVector result (nr); |
|
|
756 |
|
|
757 for (int i = 0; i < a_len; i++) |
|
|
758 result.elem (i) = a.elem (i) * m.elem (i, i); |
|
|
759 |
|
|
760 for (i = a_len; i < nr; i++) |
|
|
761 result.elem (i) = 0.0; |
|
|
762 |
|
|
763 return result; |
|
|
764 } |
|
|
765 |
|
|
766 // diagonal matrix by diagonal matrix -> diagonal matrix operations |
|
|
767 |
|
|
768 ComplexDiagMatrix |
|
|
769 operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b) |
|
|
770 { |
|
|
771 int nr_a = a.rows (); |
|
|
772 int nc_a = a.cols (); |
|
|
773 int nr_b = b.rows (); |
|
|
774 int nc_b = b.cols (); |
|
|
775 if (nc_a != nr_b) |
|
|
776 { |
|
|
777 (*current_liboctave_error_handler) |
|
|
778 ("nonconformant matrix multiplication attempted"); |
|
|
779 return ComplexDiagMatrix (); |
|
|
780 } |
|
|
781 |
|
|
782 if (nr_a == 0 || nc_a == 0 || nc_b == 0) |
|
|
783 return ComplexDiagMatrix (nr_a, nc_a, 0.0); |
|
|
784 |
|
|
785 ComplexDiagMatrix c (nr_a, nc_b); |
|
|
786 |
|
|
787 int len = nr_a < nc_b ? nr_a : nc_b; |
|
|
788 |
|
|
789 for (int i = 0; i < len; i++) |
|
|
790 { |
|
|
791 Complex a_element = a.elem (i, i); |
|
|
792 Complex b_element = b.elem (i, i); |
|
|
793 |
|
|
794 if (a_element == 0.0 || b_element == 0.0) |
|
|
795 c.elem (i, i) = 0.0; |
|
|
796 else if (a_element == 1.0) |
|
|
797 c.elem (i, i) = b_element; |
|
|
798 else if (b_element == 1.0) |
|
|
799 c.elem (i, i) = a_element; |
|
|
800 else |
|
|
801 c.elem (i, i) = a_element * b_element; |
|
|
802 } |
|
|
803 |
|
|
804 return c; |
|
|
805 } |
|
|
806 |
|
|
807 ComplexDiagMatrix |
|
|
808 operator + (const ComplexDiagMatrix& m, const DiagMatrix& a) |
|
|
809 { |
|
|
810 int nr = m.rows (); |
|
|
811 int nc = m.cols (); |
|
|
812 if (nr != a.rows () || nc != a.cols ()) |
|
|
813 { |
|
|
814 (*current_liboctave_error_handler) |
|
|
815 ("nonconformant matrix addition attempted"); |
|
|
816 return ComplexDiagMatrix (); |
|
|
817 } |
|
|
818 |
|
|
819 if (nr == 0 || nc == 0) |
|
|
820 return ComplexDiagMatrix (nr, nc); |
|
|
821 |
|
|
822 return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc); |
|
|
823 } |
|
|
824 |
|
|
825 ComplexDiagMatrix |
|
|
826 operator - (const ComplexDiagMatrix& m, const DiagMatrix& a) |
|
|
827 { |
|
|
828 int nr = m.rows (); |
|
|
829 int nc = m.cols (); |
|
|
830 if (nr != a.rows () || nc != a.cols ()) |
|
|
831 { |
|
|
832 (*current_liboctave_error_handler) |
|
|
833 ("nonconformant matrix subtraction attempted"); |
|
|
834 return ComplexDiagMatrix (); |
|
|
835 } |
|
|
836 |
|
|
837 if (nr == 0 || nc == 0) |
|
|
838 return ComplexDiagMatrix (nr, nc); |
|
|
839 |
|
|
840 return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()), |
|
|
841 nr, nc); |
|
|
842 } |
|
|
843 |
|
|
844 ComplexDiagMatrix |
|
|
845 operator * (const ComplexDiagMatrix& a, const DiagMatrix& b) |
|
|
846 { |
|
|
847 int nr_a = a.rows (); |
|
|
848 int nc_a = a.cols (); |
|
|
849 int nr_b = b.rows (); |
|
|
850 int nc_b = b.cols (); |
|
|
851 if (nc_a != nr_b) |
|
|
852 { |
|
|
853 (*current_liboctave_error_handler) |
|
|
854 ("nonconformant matrix multiplication attempted"); |
|
|
855 return ComplexDiagMatrix (); |
|
|
856 } |
|
|
857 |
|
|
858 if (nr_a == 0 || nc_a == 0 || nc_b == 0) |
|
|
859 return ComplexDiagMatrix (nr_a, nc_a, 0.0); |
|
|
860 |
|
|
861 ComplexDiagMatrix c (nr_a, nc_b); |
|
|
862 |
|
|
863 int len = nr_a < nc_b ? nr_a : nc_b; |
|
|
864 |
|
|
865 for (int i = 0; i < len; i++) |
|
|
866 { |
|
|
867 Complex a_element = a.elem (i, i); |
|
|
868 double b_element = b.elem (i, i); |
|
|
869 |
|
|
870 if (a_element == 0.0 || b_element == 0.0) |
|
|
871 c.elem (i, i) = 0.0; |
|
|
872 else if (a_element == 1.0) |
|
|
873 c.elem (i, i) = b_element; |
|
|
874 else if (b_element == 1.0) |
|
|
875 c.elem (i, i) = a_element; |
|
|
876 else |
|
|
877 c.elem (i, i) = a_element * b_element; |
|
|
878 } |
|
|
879 |
|
|
880 return c; |
|
|
881 } |
|
|
882 |
|
|
883 ComplexDiagMatrix |
|
|
884 product (const ComplexDiagMatrix& m, const DiagMatrix& a) |
|
|
885 { |
|
|
886 int nr = m.rows (); |
|
|
887 int nc = m.cols (); |
|
|
888 if (nr != a.rows () || nc != a.cols ()) |
|
|
889 { |
|
|
890 (*current_liboctave_error_handler) |
|
|
891 ("nonconformant matrix product attempted"); |
|
|
892 return ComplexDiagMatrix (); |
|
|
893 } |
|
|
894 |
|
|
895 if (nr == 0 || nc == 0) |
|
|
896 return ComplexDiagMatrix (nr, nc); |
|
|
897 |
|
|
898 return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()), |
|
|
899 nr, nc); |
|
|
900 } |
|
|
901 |
|
|
902 // diagonal matrix by matrix -> matrix operations |
|
|
903 |
|
|
904 ComplexMatrix |
|
|
905 operator + (const ComplexDiagMatrix& m, const Matrix& a) |
|
|
906 { |
|
|
907 int nr = m.rows (); |
|
|
908 int nc = m.cols (); |
|
|
909 if (nr != a.rows () || nc != a.cols ()) |
|
|
910 { |
|
|
911 (*current_liboctave_error_handler) |
|
|
912 ("nonconformant matrix addition attempted"); |
|
|
913 return ComplexMatrix (); |
|
|
914 } |
|
|
915 |
|
|
916 if (nr == 0 || nc == 0) |
|
|
917 return ComplexMatrix (nr, nc); |
|
|
918 |
|
|
919 ComplexMatrix result (a); |
|
|
920 for (int i = 0; i < m.length (); i++) |
|
|
921 result.elem (i, i) += m.elem (i, i); |
|
|
922 |
|
|
923 return result; |
|
|
924 } |
|
|
925 |
|
|
926 ComplexMatrix |
|
|
927 operator - (const ComplexDiagMatrix& m, const Matrix& a) |
|
|
928 { |
|
|
929 int nr = m.rows (); |
|
|
930 int nc = m.cols (); |
|
|
931 if (nr != a.rows () || nc != a.cols ()) |
|
|
932 { |
|
|
933 (*current_liboctave_error_handler) |
|
|
934 ("nonconformant matrix subtraction attempted"); |
|
|
935 return ComplexMatrix (); |
|
|
936 } |
|
|
937 |
|
|
938 if (nr == 0 || nc == 0) |
|
|
939 return ComplexMatrix (nr, nc); |
|
|
940 |
|
|
941 ComplexMatrix result (-a); |
|
|
942 for (int i = 0; i < m.length (); i++) |
|
|
943 result.elem (i, i) += m.elem (i, i); |
|
|
944 |
|
|
945 return result; |
|
|
946 } |
|
|
947 |
|
|
948 ComplexMatrix |
|
|
949 operator * (const ComplexDiagMatrix& m, const Matrix& a) |
|
|
950 { |
|
|
951 int nr = m.rows (); |
|
|
952 int nc = m.cols (); |
|
|
953 int a_nr = a.rows (); |
|
|
954 int a_nc = a.cols (); |
|
|
955 if (nc != a_nr) |
|
|
956 { |
|
|
957 (*current_liboctave_error_handler) |
|
|
958 ("nonconformant matrix multiplication attempted"); |
|
|
959 return ComplexMatrix (); |
|
|
960 } |
|
|
961 |
|
|
962 if (nr == 0 || nc == 0 || a_nc == 0) |
|
|
963 return ComplexMatrix (nr, a_nc, 0.0); |
|
|
964 |
|
|
965 ComplexMatrix c (nr, a_nc); |
|
|
966 |
|
|
967 for (int i = 0; i < m.length (); i++) |
|
|
968 { |
|
|
969 if (m.elem (i, i) == 1.0) |
|
|
970 { |
|
|
971 for (int j = 0; j < a_nc; j++) |
|
|
972 c.elem (i, j) = a.elem (i, j); |
|
|
973 } |
|
|
974 else if (m.elem (i, i) == 0.0) |
|
|
975 { |
|
|
976 for (int j = 0; j < a_nc; j++) |
|
|
977 c.elem (i, j) = 0.0; |
|
|
978 } |
|
|
979 else |
|
|
980 { |
|
|
981 for (int j = 0; j < a_nc; j++) |
|
|
982 c.elem (i, j) = m.elem (i, i) * a.elem (i, j); |
|
|
983 } |
|
|
984 } |
|
|
985 |
|
|
986 if (nr > nc) |
|
|
987 { |
|
|
988 for (int j = 0; j < a_nc; j++) |
|
|
989 for (int i = a_nr; i < nr; i++) |
|
|
990 c.elem (i, j) = 0.0; |
|
|
991 } |
|
|
992 |
|
|
993 return c; |
|
|
994 } |
|
|
995 |
|
|
996 ComplexMatrix |
|
|
997 operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a) |
|
|
998 { |
|
|
999 int nr = m.rows (); |
|
|
1000 int nc = m.cols (); |
|
|
1001 if (nr != a.rows () || nc != a.cols ()) |
|
|
1002 { |
|
|
1003 (*current_liboctave_error_handler) |
|
|
1004 ("nonconformant matrix addition attempted"); |
|
|
1005 return ComplexMatrix (); |
|
|
1006 } |
|
|
1007 |
|
|
1008 if (nr == 0 || nc == 0) |
|
|
1009 return ComplexMatrix (nr, nc); |
|
|
1010 |
|
|
1011 ComplexMatrix result (a); |
|
|
1012 for (int i = 0; i < m.length (); i++) |
|
|
1013 result.elem (i, i) += m.elem (i, i); |
|
|
1014 |
|
|
1015 return result; |
|
|
1016 } |
|
|
1017 |
|
|
1018 ComplexMatrix |
|
|
1019 operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a) |
|
|
1020 { |
|
|
1021 int nr = m.rows (); |
|
|
1022 int nc = m.cols (); |
|
|
1023 if (nr != a.rows () || nc != a.cols ()) |
|
|
1024 { |
|
|
1025 (*current_liboctave_error_handler) |
|
|
1026 ("nonconformant matrix subtraction attempted"); |
|
|
1027 return ComplexMatrix (); |
|
|
1028 } |
|
|
1029 |
|
|
1030 if (nr == 0 || nc == 0) |
|
|
1031 return ComplexMatrix (nr, nc); |
|
|
1032 |
|
|
1033 ComplexMatrix result (-a); |
|
|
1034 for (int i = 0; i < m.length (); i++) |
|
|
1035 result.elem (i, i) += m.elem (i, i); |
|
|
1036 |
|
|
1037 return result; |
|
|
1038 } |
|
|
1039 |
|
|
1040 ComplexMatrix |
|
|
1041 operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a) |
|
|
1042 { |
|
|
1043 int nr = m.rows (); |
|
|
1044 int nc = m.cols (); |
|
|
1045 int a_nr = a.rows (); |
|
|
1046 int a_nc = a.cols (); |
|
|
1047 if (nc != a_nr) |
|
|
1048 { |
|
|
1049 (*current_liboctave_error_handler) |
|
|
1050 ("nonconformant matrix multiplication attempted"); |
|
|
1051 return ComplexMatrix (); |
|
|
1052 } |
|
|
1053 |
|
|
1054 if (nr == 0 || nc == 0 || a_nc == 0) |
|
|
1055 return ComplexMatrix (nr, a_nc, 0.0); |
|
|
1056 |
|
|
1057 ComplexMatrix c (nr, a_nc); |
|
|
1058 |
|
|
1059 for (int i = 0; i < m.length (); i++) |
|
|
1060 { |
|
|
1061 if (m.elem (i, i) == 1.0) |
|
|
1062 { |
|
|
1063 for (int j = 0; j < a_nc; j++) |
|
|
1064 c.elem (i, j) = a.elem (i, j); |
|
|
1065 } |
|
|
1066 else if (m.elem (i, i) == 0.0) |
|
|
1067 { |
|
|
1068 for (int j = 0; j < a_nc; j++) |
|
|
1069 c.elem (i, j) = 0.0; |
|
|
1070 } |
|
|
1071 else |
|
|
1072 { |
|
|
1073 for (int j = 0; j < a_nc; j++) |
|
|
1074 c.elem (i, j) = m.elem (i, i) * a.elem (i, j); |
|
|
1075 } |
|
|
1076 } |
|
|
1077 |
|
|
1078 if (nr > nc) |
|
|
1079 { |
|
|
1080 for (int j = 0; j < a_nc; j++) |
|
|
1081 for (int i = a_nr; i < nr; i++) |
|
|
1082 c.elem (i, j) = 0.0; |
|
|
1083 } |
|
|
1084 |
|
|
1085 return c; |
|
|
1086 } |
|
|
1087 |
|
|
1088 // other operations |
|
|
1089 |
|
|
1090 ComplexColumnVector |
|
|
1091 ComplexDiagMatrix::diag (void) const |
|
|
1092 { |
|
|
1093 return diag (0); |
|
|
1094 } |
|
|
1095 |
|
|
1096 // Could be optimized... |
|
|
1097 |
|
|
1098 ComplexColumnVector |
|
|
1099 ComplexDiagMatrix::diag (int k) const |
|
|
1100 { |
|
|
1101 int nnr = rows (); |
|
|
1102 int nnc = cols (); |
|
|
1103 if (k > 0) |
|
|
1104 nnc -= k; |
|
|
1105 else if (k < 0) |
|
|
1106 nnr += k; |
|
|
1107 |
|
|
1108 ComplexColumnVector d; |
|
|
1109 |
|
|
1110 if (nnr > 0 && nnc > 0) |
|
|
1111 { |
|
|
1112 int ndiag = (nnr < nnc) ? nnr : nnc; |
|
|
1113 |
|
|
1114 d.resize (ndiag); |
|
|
1115 |
|
|
1116 if (k > 0) |
|
|
1117 { |
|
|
1118 for (int i = 0; i < ndiag; i++) |
|
|
1119 d.elem (i) = elem (i, i+k); |
|
|
1120 } |
|
|
1121 else if ( k < 0) |
|
|
1122 { |
|
|
1123 for (int i = 0; i < ndiag; i++) |
|
|
1124 d.elem (i) = elem (i-k, i); |
|
|
1125 } |
|
|
1126 else |
|
|
1127 { |
|
|
1128 for (int i = 0; i < ndiag; i++) |
|
|
1129 d.elem (i) = elem (i, i); |
|
|
1130 } |
|
|
1131 } |
|
|
1132 else |
|
|
1133 cerr << "diag: requested diagonal out of range\n"; |
|
|
1134 |
|
|
1135 return d; |
|
|
1136 } |
|
|
1137 |
|
|
1138 // i/o |
|
|
1139 |
|
|
1140 ostream& |
|
|
1141 operator << (ostream& os, const ComplexDiagMatrix& a) |
|
|
1142 { |
|
|
1143 Complex ZERO (0.0); |
|
|
1144 // int field_width = os.precision () + 7; |
|
|
1145 for (int i = 0; i < a.rows (); i++) |
|
|
1146 { |
|
|
1147 for (int j = 0; j < a.cols (); j++) |
|
|
1148 { |
|
|
1149 if (i == j) |
|
|
1150 os << " " /* setw (field_width) */ << a.elem (i, i); |
|
|
1151 else |
|
|
1152 os << " " /* setw (field_width) */ << ZERO; |
|
|
1153 } |
|
|
1154 os << "\n"; |
|
|
1155 } |
|
|
1156 return os; |
|
|
1157 } |
|
|
1158 |
|
|
1159 /* |
|
|
1160 ;;; Local Variables: *** |
|
|
1161 ;;; mode: C++ *** |
|
|
1162 ;;; page-delimiter: "^/\\*" *** |
|
|
1163 ;;; End: *** |
|
|
1164 */ |
