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
fCDiagMatrix.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 // FloatComplex Diagonal Matrix class
39 
41  : MDiagArray2<FloatComplex> (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 
81 {
82  if (beg < 0 || end >= length () || end < beg)
83  {
84  (*current_liboctave_error_handler) ("range error for fill");
85  return *this;
86  }
87 
88  for (octave_idx_type i = beg; i <= end; i++)
89  elem (i, i) = val;
90 
91  return *this;
92 }
93 
97 {
98  if (beg < 0 || end >= length () || end < beg)
99  {
100  (*current_liboctave_error_handler) ("range error for fill");
101  return *this;
102  }
103 
104  for (octave_idx_type i = beg; i <= end; i++)
105  elem (i, i) = val;
106 
107  return *this;
108 }
109 
112 {
113  octave_idx_type len = length ();
114  if (a.length () != len)
115  {
116  (*current_liboctave_error_handler) ("range error for fill");
117  return *this;
118  }
119 
120  for (octave_idx_type i = 0; i < len; i++)
121  elem (i, i) = a.elem (i);
122 
123  return *this;
124 }
125 
128 {
129  octave_idx_type len = length ();
130  if (a.length () != len)
131  {
132  (*current_liboctave_error_handler) ("range error for fill");
133  return *this;
134  }
135 
136  for (octave_idx_type i = 0; i < len; i++)
137  elem (i, i) = a.elem (i);
138 
139  return *this;
140 }
141 
144 {
145  octave_idx_type len = length ();
146  if (a.length () != len)
147  {
148  (*current_liboctave_error_handler) ("range error for fill");
149  return *this;
150  }
151 
152  for (octave_idx_type i = 0; i < len; i++)
153  elem (i, i) = a.elem (i);
154 
155  return *this;
156 }
157 
160 {
161  octave_idx_type len = length ();
162  if (a.length () != len)
163  {
164  (*current_liboctave_error_handler) ("range error for fill");
165  return *this;
166  }
167 
168  for (octave_idx_type i = 0; i < len; i++)
169  elem (i, i) = a.elem (i);
170 
171  return *this;
172 }
173 
176 {
177  octave_idx_type a_len = a.length ();
178  if (beg < 0 || beg + a_len >= length ())
179  {
180  (*current_liboctave_error_handler) ("range error for fill");
181  return *this;
182  }
183 
184  for (octave_idx_type i = 0; i < a_len; i++)
185  elem (i+beg, i+beg) = a.elem (i);
186 
187  return *this;
188 }
189 
192  octave_idx_type beg)
193 {
194  octave_idx_type a_len = a.length ();
195  if (beg < 0 || beg + a_len >= length ())
196  {
197  (*current_liboctave_error_handler) ("range error for fill");
198  return *this;
199  }
200 
201  for (octave_idx_type i = 0; i < a_len; i++)
202  elem (i+beg, i+beg) = a.elem (i);
203 
204  return *this;
205 }
206 
209 {
210  octave_idx_type a_len = a.length ();
211  if (beg < 0 || beg + a_len >= length ())
212  {
213  (*current_liboctave_error_handler) ("range error for fill");
214  return *this;
215  }
216 
217  for (octave_idx_type i = 0; i < a_len; i++)
218  elem (i+beg, i+beg) = a.elem (i);
219 
220  return *this;
221 }
222 
225  octave_idx_type beg)
226 {
227  octave_idx_type a_len = a.length ();
228  if (beg < 0 || beg + a_len >= length ())
229  {
230  (*current_liboctave_error_handler) ("range error for fill");
231  return *this;
232  }
233 
234  for (octave_idx_type i = 0; i < a_len; i++)
235  elem (i+beg, i+beg) = a.elem (i);
236 
237  return *this;
238 }
239 
242 {
243  return FloatDiagMatrix (extract_diag ().abs (), rows (), columns ());
244 }
245 
248 {
249  return FloatComplexDiagMatrix (conj (a.extract_diag ()), a.rows (),
250  a.columns ());
251 }
252 
253 // resize is the destructive analog for this one
254 
258 {
259  if (r1 > r2) { std::swap (r1, r2); }
260  if (c1 > c2) { std::swap (c1, c2); }
261 
262  octave_idx_type new_r = r2 - r1 + 1;
263  octave_idx_type new_c = c2 - c1 + 1;
264 
265  FloatComplexMatrix result (new_r, new_c);
266 
267  for (octave_idx_type j = 0; j < new_c; j++)
268  for (octave_idx_type i = 0; i < new_r; i++)
269  result.elem (i, j) = elem (r1+i, c1+j);
270 
271  return result;
272 }
273 
274 // extract row or column i.
275 
278 {
279  octave_idx_type r = rows ();
280  octave_idx_type c = cols ();
281  if (i < 0 || i >= r)
282  {
283  (*current_liboctave_error_handler) ("invalid row selection");
284  return FloatComplexRowVector ();
285  }
286 
287  FloatComplexRowVector retval (c, 0.0);
288  if (r <= c || (r > c && i < c))
289  retval.elem (i) = elem (i, i);
290 
291  return retval;
292 }
293 
296 {
297  if (! s)
298  {
299  (*current_liboctave_error_handler) ("invalid row selection");
300  return FloatComplexRowVector ();
301  }
302 
303  char c = *s;
304  if (c == 'f' || c == 'F')
305  return row (static_cast<octave_idx_type>(0));
306  else if (c == 'l' || c == 'L')
307  return row (rows () - 1);
308  else
309  {
310  (*current_liboctave_error_handler) ("invalid row selection");
311  return FloatComplexRowVector ();
312  }
313 }
314 
317 {
318  octave_idx_type r = rows ();
319  octave_idx_type c = cols ();
320  if (i < 0 || i >= c)
321  {
322  (*current_liboctave_error_handler) ("invalid column selection");
323  return FloatComplexColumnVector ();
324  }
325 
326  FloatComplexColumnVector retval (r, 0.0);
327  if (r >= c || (r < c && i < r))
328  retval.elem (i) = elem (i, i);
329 
330  return retval;
331 }
332 
335 {
336  if (! s)
337  {
338  (*current_liboctave_error_handler) ("invalid column selection");
339  return FloatComplexColumnVector ();
340  }
341 
342  char c = *s;
343  if (c == 'f' || c == 'F')
344  return column (static_cast<octave_idx_type>(0));
345  else if (c == 'l' || c == 'L')
346  return column (cols () - 1);
347  else
348  {
349  (*current_liboctave_error_handler) ("invalid column selection");
350  return FloatComplexColumnVector ();
351  }
352 }
353 
356 {
357  octave_idx_type info;
358  return inverse (info);
359 }
360 
363 {
364  octave_idx_type r = rows ();
365  octave_idx_type c = cols ();
366  if (r != c)
367  {
368  (*current_liboctave_error_handler) ("inverse requires square matrix");
369  return FloatComplexDiagMatrix ();
370  }
371 
372  FloatComplexDiagMatrix retval (r, c);
373 
374  info = 0;
375  for (octave_idx_type i = 0; i < length (); i++)
376  {
377  if (elem (i, i) == 0.0f)
378  {
379  info = -1;
380  return *this;
381  }
382  else
383  retval.elem (i, i) = 1.0f / elem (i, i);
384  }
385 
386  return retval;
387 }
388 
391 {
392  octave_idx_type r = rows ();
393  octave_idx_type c = cols ();
394  octave_idx_type len = length ();
395 
396  FloatComplexDiagMatrix retval (c, r);
397 
398  for (octave_idx_type i = 0; i < len; i++)
399  {
400  float val = std::abs (elem (i, i));
401  if (val < tol || val == 0.0f)
402  retval.elem (i, i) = 0.0f;
403  else
404  retval.elem (i, i) = 1.0f / elem (i, i);
405  }
406 
407  return retval;
408 }
409 
410 bool
412 {
413  return mx_inline_all_real (length (), data ());
414 }
415 
416 // diagonal matrix by diagonal matrix -> diagonal matrix operations
417 
420 {
421  octave_idx_type r = rows ();
422  octave_idx_type c = cols ();
423 
424  octave_idx_type a_nr = a.rows ();
425  octave_idx_type a_nc = a.cols ();
426 
427  if (r != a_nr || c != a_nc)
428  {
429  gripe_nonconformant ("operator +=", r, c, a_nr, a_nc);
430  return *this;
431  }
432 
433  if (r == 0 || c == 0)
434  return *this;
435 
436  FloatComplex *d = fortran_vec (); // Ensures only 1 reference to my privates!
437 
438  mx_inline_add2 (length (), d, a.data ());
439  return *this;
440 }
441 
444 {
445  octave_idx_type a_nr = a.rows ();
446  octave_idx_type a_nc = a.cols ();
447 
448  octave_idx_type b_nr = b.rows ();
449  octave_idx_type b_nc = b.cols ();
450 
451  if (a_nc != b_nr)
452  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
453 
454  FloatComplexDiagMatrix c (a_nr, b_nc);
455 
456  octave_idx_type len = c.length ();
457  octave_idx_type lenm = len < a_nc ? len : a_nc;
458 
459  for (octave_idx_type i = 0; i < lenm; i++)
460  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
461  for (octave_idx_type i = lenm; i < len; i++)
462  c.dgxelem (i) = 0.0f;
463 
464  return c;
465 }
466 
469 {
470  octave_idx_type a_nr = a.rows ();
471  octave_idx_type a_nc = a.cols ();
472 
473  octave_idx_type b_nr = b.rows ();
474  octave_idx_type b_nc = b.cols ();
475 
476  if (a_nc != b_nr)
477  {
478  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
479  return FloatComplexDiagMatrix ();
480  }
481 
482  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
483  return FloatComplexDiagMatrix (a_nr, a_nc, 0.0);
484 
485  FloatComplexDiagMatrix c (a_nr, b_nc);
486 
487  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
488 
489  for (octave_idx_type i = 0; i < len; i++)
490  {
491  float a_element = a.elem (i, i);
492  FloatComplex b_element = b.elem (i, i);
493 
494  c.elem (i, i) = a_element * b_element;
495  }
496 
497  return c;
498 }
499 
502 {
503  octave_idx_type a_nr = a.rows ();
504  octave_idx_type a_nc = a.cols ();
505 
506  octave_idx_type b_nr = b.rows ();
507  octave_idx_type b_nc = b.cols ();
508 
509  if (a_nc != b_nr)
510  {
511  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
512  return FloatComplexDiagMatrix ();
513  }
514 
515  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
516  return FloatComplexDiagMatrix (a_nr, a_nc, 0.0);
517 
518  FloatComplexDiagMatrix c (a_nr, b_nc);
519 
520  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
521 
522  for (octave_idx_type i = 0; i < len; i++)
523  {
524  FloatComplex a_element = a.elem (i, i);
525  FloatComplex b_element = b.elem (i, i);
526 
527  c.elem (i, i) = a_element * b_element;
528  }
529 
530  return c;
531 }
532 
533 // other operations
534 
537 {
538  FloatComplexDET det (1.0f);
539  if (rows () != cols ())
540  {
541  (*current_liboctave_error_handler) ("determinant requires square matrix");
542  det = FloatComplexDET (0.0);
543  }
544  else
545  {
546  octave_idx_type len = length ();
547  for (octave_idx_type i = 0; i < len; i++)
548  det *= elem (i, i);
549  }
550 
551  return det;
552 }
553 
554 float
556 {
557  FloatColumnVector av = extract_diag (0).map<float> (std::abs);
558  float amx = av.max ();
559  float amn = av.min ();
560  return amx == 0 ? 0.0f : amn / amx;
561 }
562 
563 // i/o
564 
565 std::ostream&
566 operator << (std::ostream& os, const FloatComplexDiagMatrix& a)
567 {
568  FloatComplex ZERO (0.0);
569 // int field_width = os.precision () + 7;
570  for (octave_idx_type i = 0; i < a.rows (); i++)
571  {
572  for (octave_idx_type j = 0; j < a.cols (); j++)
573  {
574  if (i == j)
575  os << " " /* setw (field_width) */ << a.elem (i, i);
576  else
577  os << " " /* setw (field_width) */ << ZERO;
578  }
579  os << "\n";
580  }
581  return os;
582 }
FloatComplexColumnVector column(octave_idx_type i) const
void mx_inline_add2(size_t n, R *r, const X *x)
Definition: mx-inlines.cc:95
FloatComplexDiagMatrix inverse(void) const
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
base_det< FloatComplex > FloatComplexDET
Definition: DET.h:88
const FloatComplex * fortran_vec(void) const
Definition: DiagArray2.h:182
FloatComplex elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:110
FloatComplexRowVector row(octave_idx_type i) const
float rcond(void) const
FloatComplexDiagMatrix & fill(float val)
Definition: fCDiagMatrix.cc:63
octave_idx_type rows(void) const
Definition: DiagArray2.h:86
const FloatComplex * data(void) const
Definition: DiagArray2.h:180
T & elem(octave_idx_type n)
Definition: Array.h:380
bool operator==(const FloatComplexDiagMatrix &a) const
Definition: fCDiagMatrix.cc:48
octave_idx_type rows(void) const
Definition: Array.h:313
F77_RET_T const double const double double * d
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
float min(void) const
Definition: fColVector.cc:267
bool all_elements_are_real(void) const
FloatComplexDiagMatrix pseudo_inverse(float tol=0.0f) const
F77_RET_T const double const double * f
bool operator!=(const FloatComplexDiagMatrix &a) const
Definition: fCDiagMatrix.cc:57
Array< U > map(F fcn) const
Apply function fcn to each element of the Array.
Definition: Array.h:659
Definition: DET.h:31
T & dgxelem(octave_idx_type i)
Definition: DiagArray2.h:163
FloatComplexDiagMatrix operator*(const FloatComplexDiagMatrix &a, const FloatDiagMatrix &b)
FloatComplexColumnVector extract_diag(octave_idx_type k=0) const
Definition: fCDiagMatrix.h:136
FloatComplexDiagMatrix & operator+=(const FloatDiagMatrix &a)
FloatComplexDET determinant(void) const
friend class FloatComplexColumnVector
Definition: fCRowVector.h:35
FloatComplexDiagMatrix conj(const FloatComplexDiagMatrix &a)
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
float max(void) const
Definition: fColVector.cc:283
octave_idx_type columns(void) const
Definition: DiagArray2.h:88
FloatDiagMatrix abs(void) const
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
std::ostream & operator<<(std::ostream &os, const FloatComplexDiagMatrix &a)
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
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
FloatComplexMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
T abs(T x)
Definition: pr-output.cc:3062
octave_idx_type length(void) const
Definition: DiagArray2.h:92