GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CDiagMatrix.cc
Go to the documentation of this file.
1 // DiagMatrix manipulations.
2 /*
3 
4 Copyright (C) 1994-2015 John W. Eaton
5 Copyright (C) 2009 VZLU Prague
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
11 Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
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
20 along with Octave; see the file COPYING. If not, see
21 <http://www.gnu.org/licenses/>.
22 
23 */
24 
25 #ifdef HAVE_CONFIG_H
26 #include <config.h>
27 #endif
28 
29 #include <iostream>
30 
31 #include "Array-util.h"
32 #include "lo-error.h"
33 #include "lo-ieee.h"
34 #include "mx-base.h"
35 #include "mx-inlines.cc"
36 #include "oct-cmplx.h"
37 
38 // Complex Diagonal Matrix class
39 
41  : MDiagArray2<Complex> (a.rows (), a.cols ())
42 {
43  for (octave_idx_type i = 0; i < length (); i++)
44  elem (i, i) = a.elem (i, i);
45 }
46 
47 bool
49 {
50  if (rows () != a.rows () || cols () != a.cols ())
51  return 0;
52 
53  return mx_inline_equal (length (), data (), a.data ());
54 }
55 
56 bool
58 {
59  return !(*this == a);
60 }
61 
64 {
65  for (octave_idx_type i = 0; i < length (); i++)
66  elem (i, i) = val;
67  return *this;
68 }
69 
72 {
73  for (octave_idx_type i = 0; i < length (); i++)
74  elem (i, i) = val;
75  return *this;
76 }
77 
80 {
81  if (beg < 0 || end >= length () || end < beg)
82  {
83  (*current_liboctave_error_handler) ("range error for fill");
84  return *this;
85  }
86 
87  for (octave_idx_type i = beg; i <= end; i++)
88  elem (i, i) = val;
89 
90  return *this;
91 }
92 
96 {
97  if (beg < 0 || end >= length () || end < beg)
98  {
99  (*current_liboctave_error_handler) ("range error for fill");
100  return *this;
101  }
102 
103  for (octave_idx_type i = beg; i <= end; i++)
104  elem (i, i) = val;
105 
106  return *this;
107 }
108 
111 {
112  octave_idx_type len = length ();
113  if (a.length () != len)
114  {
115  (*current_liboctave_error_handler) ("range error for fill");
116  return *this;
117  }
118 
119  for (octave_idx_type i = 0; i < len; i++)
120  elem (i, i) = a.elem (i);
121 
122  return *this;
123 }
124 
127 {
128  octave_idx_type len = length ();
129  if (a.length () != len)
130  {
131  (*current_liboctave_error_handler) ("range error for fill");
132  return *this;
133  }
134 
135  for (octave_idx_type i = 0; i < len; i++)
136  elem (i, i) = a.elem (i);
137 
138  return *this;
139 }
140 
143 {
144  octave_idx_type len = length ();
145  if (a.length () != len)
146  {
147  (*current_liboctave_error_handler) ("range error for fill");
148  return *this;
149  }
150 
151  for (octave_idx_type i = 0; i < len; i++)
152  elem (i, i) = a.elem (i);
153 
154  return *this;
155 }
156 
159 {
160  octave_idx_type len = length ();
161  if (a.length () != len)
162  {
163  (*current_liboctave_error_handler) ("range error for fill");
164  return *this;
165  }
166 
167  for (octave_idx_type i = 0; i < len; i++)
168  elem (i, i) = a.elem (i);
169 
170  return *this;
171 }
172 
175 {
176  octave_idx_type a_len = a.length ();
177  if (beg < 0 || beg + a_len >= length ())
178  {
179  (*current_liboctave_error_handler) ("range error for fill");
180  return *this;
181  }
182 
183  for (octave_idx_type i = 0; i < a_len; i++)
184  elem (i+beg, i+beg) = a.elem (i);
185 
186  return *this;
187 }
188 
191 {
192  octave_idx_type a_len = a.length ();
193  if (beg < 0 || beg + a_len >= length ())
194  {
195  (*current_liboctave_error_handler) ("range error for fill");
196  return *this;
197  }
198 
199  for (octave_idx_type i = 0; i < a_len; i++)
200  elem (i+beg, i+beg) = a.elem (i);
201 
202  return *this;
203 }
204 
207 {
208  octave_idx_type a_len = a.length ();
209  if (beg < 0 || beg + a_len >= length ())
210  {
211  (*current_liboctave_error_handler) ("range error for fill");
212  return *this;
213  }
214 
215  for (octave_idx_type i = 0; i < a_len; i++)
216  elem (i+beg, i+beg) = a.elem (i);
217 
218  return *this;
219 }
220 
223 {
224  octave_idx_type a_len = a.length ();
225  if (beg < 0 || beg + a_len >= length ())
226  {
227  (*current_liboctave_error_handler) ("range error for fill");
228  return *this;
229  }
230 
231  for (octave_idx_type i = 0; i < a_len; i++)
232  elem (i+beg, i+beg) = a.elem (i);
233 
234  return *this;
235 }
236 
239 {
240  return DiagMatrix (extract_diag ().abs (), rows (), columns ());
241 }
242 
245 {
246  return ComplexDiagMatrix (conj (a.extract_diag ()), a.rows (), a.columns ());
247 }
248 
249 // resize is the destructive analog for this one
250 
254 {
255  if (r1 > r2) { std::swap (r1, r2); }
256  if (c1 > c2) { std::swap (c1, c2); }
257 
258  octave_idx_type new_r = r2 - r1 + 1;
259  octave_idx_type new_c = c2 - c1 + 1;
260 
261  ComplexMatrix result (new_r, new_c);
262 
263  for (octave_idx_type j = 0; j < new_c; j++)
264  for (octave_idx_type i = 0; i < new_r; i++)
265  result.elem (i, j) = elem (r1+i, c1+j);
266 
267  return result;
268 }
269 
270 // extract row or column i.
271 
274 {
275  octave_idx_type r = rows ();
276  octave_idx_type c = cols ();
277  if (i < 0 || i >= r)
278  {
279  (*current_liboctave_error_handler) ("invalid row selection");
280  return ComplexRowVector ();
281  }
282 
283  ComplexRowVector retval (c, 0.0);
284  if (r <= c || (r > c && i < c))
285  retval.elem (i) = elem (i, i);
286 
287  return retval;
288 }
289 
291 ComplexDiagMatrix::row (char *s) const
292 {
293  if (! s)
294  {
295  (*current_liboctave_error_handler) ("invalid row selection");
296  return ComplexRowVector ();
297  }
298 
299  char c = *s;
300  if (c == 'f' || c == 'F')
301  return row (static_cast<octave_idx_type>(0));
302  else if (c == 'l' || c == 'L')
303  return row (rows () - 1);
304  else
305  {
306  (*current_liboctave_error_handler) ("invalid row selection");
307  return ComplexRowVector ();
308  }
309 }
310 
313 {
314  octave_idx_type r = rows ();
315  octave_idx_type c = cols ();
316  if (i < 0 || i >= c)
317  {
318  (*current_liboctave_error_handler) ("invalid column selection");
319  return ComplexColumnVector ();
320  }
321 
322  ComplexColumnVector retval (r, 0.0);
323  if (r >= c || (r < c && i < r))
324  retval.elem (i) = elem (i, i);
325 
326  return retval;
327 }
328 
331 {
332  if (! s)
333  {
334  (*current_liboctave_error_handler) ("invalid column selection");
335  return ComplexColumnVector ();
336  }
337 
338  char c = *s;
339  if (c == 'f' || c == 'F')
340  return column (static_cast<octave_idx_type>(0));
341  else if (c == 'l' || c == 'L')
342  return column (cols () - 1);
343  else
344  {
345  (*current_liboctave_error_handler) ("invalid column selection");
346  return ComplexColumnVector ();
347  }
348 }
349 
352 {
353  octave_idx_type info;
354  return inverse (info);
355 }
356 
359 {
360  octave_idx_type r = rows ();
361  octave_idx_type c = cols ();
362  if (r != c)
363  {
364  (*current_liboctave_error_handler) ("inverse requires square matrix");
365  return ComplexDiagMatrix ();
366  }
367 
368  ComplexDiagMatrix retval (r, c);
369 
370  info = 0;
371  for (octave_idx_type i = 0; i < length (); i++)
372  {
373  if (elem (i, i) == 0.0)
374  {
375  info = -1;
376  return *this;
377  }
378  else
379  retval.elem (i, i) = 1.0 / elem (i, i);
380  }
381 
382  return retval;
383 }
384 
387 {
388  octave_idx_type r = rows ();
389  octave_idx_type c = cols ();
390  octave_idx_type len = length ();
391 
392  ComplexDiagMatrix retval (c, r);
393 
394  for (octave_idx_type i = 0; i < len; i++)
395  {
396  double val = std::abs (elem (i, i));
397  if (val < tol || val == 0.0)
398  retval.elem (i, i) = 0.0;
399  else
400  retval.elem (i, i) = 1.0 / elem (i, i);
401  }
402 
403  return retval;
404 }
405 
406 bool
408 {
409  return mx_inline_all_real (length (), data ());
410 }
411 
412 // diagonal matrix by diagonal matrix -> diagonal matrix operations
413 
416 {
417  octave_idx_type r = rows ();
418  octave_idx_type c = cols ();
419 
420  octave_idx_type a_nr = a.rows ();
421  octave_idx_type a_nc = a.cols ();
422 
423  if (r != a_nr || c != a_nc)
424  {
425  gripe_nonconformant ("operator +=", r, c, a_nr, a_nc);
426  return *this;
427  }
428 
429  if (r == 0 || c == 0)
430  return *this;
431 
432  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
433 
434  mx_inline_add2 (length (), d, a.data ());
435  return *this;
436 }
437 
440 {
441  octave_idx_type a_nr = a.rows ();
442  octave_idx_type a_nc = a.cols ();
443 
444  octave_idx_type b_nr = b.rows ();
445  octave_idx_type b_nc = b.cols ();
446 
447  if (a_nc != b_nr)
448  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
449 
450  ComplexDiagMatrix c (a_nr, b_nc);
451 
452  octave_idx_type len = c.length ();
453  octave_idx_type lenm = len < a_nc ? len : a_nc;
454 
455  for (octave_idx_type i = 0; i < lenm; i++)
456  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
457  for (octave_idx_type i = lenm; i < len; i++)
458  c.dgxelem (i) = 0.0;
459 
460  return c;
461 }
462 
465 {
466  octave_idx_type a_nr = a.rows ();
467  octave_idx_type a_nc = a.cols ();
468 
469  octave_idx_type b_nr = b.rows ();
470  octave_idx_type b_nc = b.cols ();
471 
472  if (a_nc != b_nr)
473  {
474  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
475  return ComplexDiagMatrix ();
476  }
477 
478  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
479  return ComplexDiagMatrix (a_nr, a_nc, 0.0);
480 
481  ComplexDiagMatrix c (a_nr, b_nc);
482 
483  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
484 
485  for (octave_idx_type i = 0; i < len; i++)
486  {
487  double a_element = a.elem (i, i);
488  Complex b_element = b.elem (i, i);
489 
490  c.elem (i, i) = a_element * b_element;
491  }
492 
493  return c;
494 }
495 
498 {
499  octave_idx_type a_nr = a.rows ();
500  octave_idx_type a_nc = a.cols ();
501 
502  octave_idx_type b_nr = b.rows ();
503  octave_idx_type b_nc = b.cols ();
504 
505  if (a_nc != b_nr)
506  {
507  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
508  return ComplexDiagMatrix ();
509  }
510 
511  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
512  return ComplexDiagMatrix (a_nr, a_nc, 0.0);
513 
514  ComplexDiagMatrix c (a_nr, b_nc);
515 
516  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
517 
518  for (octave_idx_type i = 0; i < len; i++)
519  {
520  Complex a_element = a.elem (i, i);
521  Complex b_element = b.elem (i, i);
522 
523  c.elem (i, i) = a_element * b_element;
524  }
525 
526  return c;
527 }
528 
529 // other operations
530 
533 {
534  ComplexDET det (1.0);
535  if (rows () != cols ())
536  {
537  (*current_liboctave_error_handler) ("determinant requires square matrix");
538  det = ComplexDET (0.0);
539  }
540  else
541  {
542  octave_idx_type len = length ();
543  for (octave_idx_type i = 0; i < len; i++)
544  det *= elem (i, i);
545  }
546 
547  return det;
548 }
549 
550 double
552 {
553  ColumnVector av = extract_diag (0).map<double> (std::abs);
554  double amx = av.max ();
555  double amn = av.min ();
556  return amx == 0 ? 0.0 : amn / amx;
557 }
558 
559 // i/o
560 
561 std::ostream&
562 operator << (std::ostream& os, const ComplexDiagMatrix& a)
563 {
564  Complex ZERO (0.0);
565 // int field_width = os.precision () + 7;
566  for (octave_idx_type i = 0; i < a.rows (); i++)
567  {
568  for (octave_idx_type j = 0; j < a.cols (); j++)
569  {
570  if (i == j)
571  os << " " /* setw (field_width) */ << a.elem (i, i);
572  else
573  os << " " /* setw (field_width) */ << ZERO;
574  }
575  os << "\n";
576  }
577  return os;
578 }
void mx_inline_add2(size_t n, R *r, const X *x)
Definition: mx-inlines.cc:95
bool operator!=(const ComplexDiagMatrix &a) const
Definition: CDiagMatrix.cc:57
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
const Complex * fortran_vec(void) const
Definition: DiagArray2.h:182
Complex elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:110
std::ostream & operator<<(std::ostream &os, const ComplexDiagMatrix &a)
Definition: CDiagMatrix.cc:562
ComplexDiagMatrix conj(const ComplexDiagMatrix &a)
Definition: CDiagMatrix.cc:244
ComplexColumnVector extract_diag(octave_idx_type k=0) const
Definition: CDiagMatrix.h:130
ComplexDiagMatrix inverse(void) const
Definition: CDiagMatrix.cc:351
ComplexMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition: CDiagMatrix.cc:252
octave_idx_type rows(void) const
Definition: DiagArray2.h:86
const Complex * data(void) const
Definition: DiagArray2.h:180
ComplexRowVector row(octave_idx_type i) const
Definition: CDiagMatrix.cc:273
T & elem(octave_idx_type n)
Definition: Array.h:380
ComplexDiagMatrix pseudo_inverse(double tol=0.0) const
Definition: CDiagMatrix.cc:386
bool operator==(const ComplexDiagMatrix &a) const
Definition: CDiagMatrix.cc:48
octave_idx_type rows(void) const
Definition: Array.h:313
F77_RET_T const double const double double * d
ComplexDiagMatrix operator*(const ComplexDiagMatrix &a, const DiagMatrix &b)
Definition: CDiagMatrix.cc:439
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type r2
bool mx_inline_all_real(size_t n, const std::complex< T > *x)
Definition: mx-inlines.cc:234
T dgelem(octave_idx_type i) const
Definition: DiagArray2.h:121
ComplexRowVector(void)
Definition: CRowVector.h:39
Array< U > map(F fcn) const
Apply function fcn to each element of the Array.
Definition: Array.h:659
ComplexDiagMatrix & operator+=(const DiagMatrix &a)
Definition: CDiagMatrix.cc:415
Definition: DET.h:31
friend class ComplexColumnVector
Definition: CRowVector.h:35
ComplexDiagMatrix(void)
Definition: CDiagMatrix.h:42
T & dgxelem(octave_idx_type i)
Definition: DiagArray2.h:163
double rcond(void) const
Definition: CDiagMatrix.cc:551
ComplexDET determinant(void) const
Definition: CDiagMatrix.cc:532
octave_idx_type cols(void) const
Definition: DiagArray2.h:87
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
double min(void) const
Definition: dColVector.cc:268
double max(void) const
Definition: dColVector.cc:284
ComplexDiagMatrix & fill(double val)
Definition: CDiagMatrix.cc:63
octave_idx_type columns(void) const
Definition: DiagArray2.h:88
base_det< Complex > ComplexDET
Definition: DET.h:87
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type octave_idx_type c1
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type r1
std::complex< double > Complex
Definition: oct-cmplx.h:29
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:441
octave_idx_type cols(void) const
Definition: Array.h:321
bool all_elements_are_real(void) const
Definition: CDiagMatrix.cc:407
ComplexColumnVector column(octave_idx_type i) const
Definition: CDiagMatrix.cc:312
T abs(T x)
Definition: pr-output.cc:3062
octave_idx_type length(void) const
Definition: DiagArray2.h:92
DiagMatrix abs(void) const
Definition: CDiagMatrix.cc:238