GNU Octave  4.2.1
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-2017 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 #if defined (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  (*current_liboctave_error_handler) ("range error for fill");
84 
85  for (octave_idx_type i = beg; i <= end; i++)
86  elem (i, i) = val;
87 
88  return *this;
89 }
90 
94 {
95  if (beg < 0 || end >= length () || end < beg)
96  (*current_liboctave_error_handler) ("range error for fill");
97 
98  for (octave_idx_type i = beg; i <= end; i++)
99  elem (i, i) = val;
100 
101  return *this;
102 }
103 
106 {
107  octave_idx_type len = length ();
108  if (a.numel () != len)
109  (*current_liboctave_error_handler) ("range error for fill");
110 
111  for (octave_idx_type i = 0; i < len; i++)
112  elem (i, i) = a.elem (i);
113 
114  return *this;
115 }
116 
119 {
120  octave_idx_type len = length ();
121  if (a.numel () != len)
122  (*current_liboctave_error_handler) ("range error for fill");
123 
124  for (octave_idx_type i = 0; i < len; i++)
125  elem (i, i) = a.elem (i);
126 
127  return *this;
128 }
129 
132 {
133  octave_idx_type len = length ();
134  if (a.numel () != len)
135  (*current_liboctave_error_handler) ("range error for fill");
136 
137  for (octave_idx_type i = 0; i < len; i++)
138  elem (i, i) = a.elem (i);
139 
140  return *this;
141 }
142 
145 {
146  octave_idx_type len = length ();
147  if (a.numel () != len)
148  (*current_liboctave_error_handler) ("range error for fill");
149 
150  for (octave_idx_type i = 0; i < len; i++)
151  elem (i, i) = a.elem (i);
152 
153  return *this;
154 }
155 
158 {
159  octave_idx_type a_len = a.numel ();
160  if (beg < 0 || beg + a_len >= length ())
161  (*current_liboctave_error_handler) ("range error for fill");
162 
163  for (octave_idx_type i = 0; i < a_len; i++)
164  elem (i+beg, i+beg) = a.elem (i);
165 
166  return *this;
167 }
168 
171  octave_idx_type beg)
172 {
173  octave_idx_type a_len = a.numel ();
174  if (beg < 0 || beg + a_len >= length ())
175  (*current_liboctave_error_handler) ("range error for fill");
176 
177  for (octave_idx_type i = 0; i < a_len; i++)
178  elem (i+beg, i+beg) = a.elem (i);
179 
180  return *this;
181 }
182 
185 {
186  octave_idx_type a_len = a.numel ();
187  if (beg < 0 || beg + a_len >= length ())
188  (*current_liboctave_error_handler) ("range error for fill");
189 
190  for (octave_idx_type i = 0; i < a_len; i++)
191  elem (i+beg, i+beg) = a.elem (i);
192 
193  return *this;
194 }
195 
198  octave_idx_type beg)
199 {
200  octave_idx_type a_len = a.numel ();
201  if (beg < 0 || beg + a_len >= length ())
202  (*current_liboctave_error_handler) ("range error for fill");
203 
204  for (octave_idx_type i = 0; i < a_len; i++)
205  elem (i+beg, i+beg) = a.elem (i);
206 
207  return *this;
208 }
209 
212 {
213  return FloatDiagMatrix (extract_diag ().abs (), rows (), columns ());
214 }
215 
218 {
219  return FloatComplexDiagMatrix (conj (a.extract_diag ()), a.rows (),
220  a.columns ());
221 }
222 
223 // resize is the destructive analog for this one
224 
227  octave_idx_type r2, octave_idx_type c2) const
228 {
229  if (r1 > r2) { std::swap (r1, r2); }
230  if (c1 > c2) { std::swap (c1, c2); }
231 
232  octave_idx_type new_r = r2 - r1 + 1;
233  octave_idx_type new_c = c2 - c1 + 1;
234 
235  FloatComplexMatrix result (new_r, new_c);
236 
237  for (octave_idx_type j = 0; j < new_c; j++)
238  for (octave_idx_type i = 0; i < new_r; i++)
239  result.elem (i, j) = elem (r1+i, c1+j);
240 
241  return result;
242 }
243 
244 // extract row or column i.
245 
248 {
249  octave_idx_type r = rows ();
250  octave_idx_type c = cols ();
251  if (i < 0 || i >= r)
252  (*current_liboctave_error_handler) ("invalid row selection");
253 
254  FloatComplexRowVector retval (c, 0.0);
255  if (r <= c || (r > c && i < c))
256  retval.elem (i) = elem (i, i);
257 
258  return retval;
259 }
260 
263 {
264  if (! s)
265  (*current_liboctave_error_handler) ("invalid row selection");
266 
267  char c = *s;
268  if (c == 'f' || c == 'F')
269  return row (static_cast<octave_idx_type>(0));
270  else if (c == 'l' || c == 'L')
271  return row (rows () - 1);
272  else
273  (*current_liboctave_error_handler) ("invalid row selection");
274 }
275 
278 {
279  octave_idx_type r = rows ();
280  octave_idx_type c = cols ();
281  if (i < 0 || i >= c)
282  (*current_liboctave_error_handler) ("invalid column selection");
283 
285  if (r >= c || (r < c && i < r))
286  retval.elem (i) = elem (i, i);
287 
288  return retval;
289 }
290 
293 {
294  if (! s)
295  (*current_liboctave_error_handler) ("invalid column selection");
296 
297  char c = *s;
298  if (c == 'f' || c == 'F')
299  return column (static_cast<octave_idx_type>(0));
300  else if (c == 'l' || c == 'L')
301  return column (cols () - 1);
302  else
303  (*current_liboctave_error_handler) ("invalid column selection");
304 }
305 
308 {
309  octave_idx_type info;
310  return inverse (info);
311 }
312 
315 {
316  octave_idx_type r = rows ();
317  octave_idx_type c = cols ();
318  if (r != c)
319  (*current_liboctave_error_handler) ("inverse requires square matrix");
320 
322 
323  info = 0;
324  for (octave_idx_type i = 0; i < length (); i++)
325  {
326  if (elem (i, i) == 0.0f)
327  {
328  info = -1;
329  return *this;
330  }
331  else
332  retval.elem (i, i) = 1.0f / elem (i, i);
333  }
334 
335  return retval;
336 }
337 
340 {
341  octave_idx_type r = rows ();
342  octave_idx_type c = cols ();
343  octave_idx_type len = length ();
344 
346 
347  for (octave_idx_type i = 0; i < len; i++)
348  {
349  float val = std::abs (elem (i, i));
350  if (val < tol || val == 0.0f)
351  retval.elem (i, i) = 0.0f;
352  else
353  retval.elem (i, i) = 1.0f / elem (i, i);
354  }
355 
356  return retval;
357 }
358 
359 bool
361 {
362  return mx_inline_all_real (length (), data ());
363 }
364 
365 // diagonal matrix by diagonal matrix -> diagonal matrix operations
366 
369 {
370  octave_idx_type r = rows ();
371  octave_idx_type c = cols ();
372 
373  octave_idx_type a_nr = a.rows ();
374  octave_idx_type a_nc = a.cols ();
375 
376  if (r != a_nr || c != a_nc)
377  octave::err_nonconformant ("operator +=", r, c, a_nr, a_nc);
378 
379  if (r == 0 || c == 0)
380  return *this;
381 
382  FloatComplex *d = fortran_vec (); // Ensures only 1 reference to my privates!
383 
384  mx_inline_add2 (length (), d, a.data ());
385  return *this;
386 }
387 
390 {
391  octave_idx_type a_nr = a.rows ();
392  octave_idx_type a_nc = a.cols ();
393 
394  octave_idx_type b_nr = b.rows ();
395  octave_idx_type b_nc = b.cols ();
396 
397  if (a_nc != b_nr)
398  octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
399 
400  FloatComplexDiagMatrix c (a_nr, b_nc);
401 
402  octave_idx_type len = c.length ();
403  octave_idx_type lenm = len < a_nc ? len : a_nc;
404 
405  for (octave_idx_type i = 0; i < lenm; i++)
406  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
407  for (octave_idx_type i = lenm; i < len; i++)
408  c.dgxelem (i) = 0.0f;
409 
410  return c;
411 }
412 
415 {
416  octave_idx_type a_nr = a.rows ();
417  octave_idx_type a_nc = a.cols ();
418 
419  octave_idx_type b_nr = b.rows ();
420  octave_idx_type b_nc = b.cols ();
421 
422  if (a_nc != b_nr)
423  octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
424 
425  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
426  return FloatComplexDiagMatrix (a_nr, a_nc, 0.0);
427 
428  FloatComplexDiagMatrix c (a_nr, b_nc);
429 
430  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
431 
432  for (octave_idx_type i = 0; i < len; i++)
433  {
434  float a_element = a.elem (i, i);
435  FloatComplex b_element = b.elem (i, i);
436 
437  c.elem (i, i) = a_element * b_element;
438  }
439 
440  return c;
441 }
442 
445 {
446  octave_idx_type a_nr = a.rows ();
447  octave_idx_type a_nc = a.cols ();
448 
449  octave_idx_type b_nr = b.rows ();
450  octave_idx_type b_nc = b.cols ();
451 
452  if (a_nc != b_nr)
453  octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
454 
455  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
456  return FloatComplexDiagMatrix (a_nr, a_nc, 0.0);
457 
458  FloatComplexDiagMatrix c (a_nr, b_nc);
459 
460  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
461 
462  for (octave_idx_type i = 0; i < len; i++)
463  {
464  FloatComplex a_element = a.elem (i, i);
465  FloatComplex b_element = b.elem (i, i);
466 
467  c.elem (i, i) = a_element * b_element;
468  }
469 
470  return c;
471 }
472 
473 // other operations
474 
477 {
478  FloatComplexDET det (1.0f);
479  if (rows () != cols ())
480  (*current_liboctave_error_handler) ("determinant requires square matrix");
481 
482  octave_idx_type len = length ();
483  for (octave_idx_type i = 0; i < len; i++)
484  det *= elem (i, i);
485 
486  return det;
487 }
488 
489 float
491 {
492  FloatColumnVector av = extract_diag (0).map<float> (std::abs);
493  float amx = av.max ();
494  float amn = av.min ();
495  return amx == 0 ? 0.0f : amn / amx;
496 }
497 
498 // i/o
499 
500 std::ostream&
501 operator << (std::ostream& os, const FloatComplexDiagMatrix& a)
502 {
503  FloatComplex ZERO (0.0);
504 // int field_width = os.precision () + 7;
505  for (octave_idx_type i = 0; i < a.rows (); i++)
506  {
507  for (octave_idx_type j = 0; j < a.cols (); j++)
508  {
509  if (i == j)
510  os << " " /* setw (field_width) */ << a.elem (i, i);
511  else
512  os << " " /* setw (field_width) */ << ZERO;
513  }
514  os << "\n";
515  }
516  return os;
517 }
FloatComplexColumnVector column(octave_idx_type i) const
void mx_inline_add2(size_t n, R *r, const X *x)
Definition: mx-inlines.cc:127
FloatComplexDiagMatrix inverse(void) const
const FloatComplex * fortran_vec(void) const
Definition: DiagArray2.h:178
FloatComplex elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:114
FloatComplexRowVector row(octave_idx_type i) const
float rcond(void) const
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
identity matrix If supplied two scalar respectively For allows like xample val
Definition: data.cc:5068
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
FloatComplexDiagMatrix & fill(float val)
Definition: fCDiagMatrix.cc:63
octave_idx_type b_nr
Definition: sylvester.cc:74
octave_idx_type rows(void) const
Definition: DiagArray2.h:88
const FloatComplex * data(void) const
Definition: DiagArray2.h:176
T & elem(octave_idx_type n)
Definition: Array.h:482
s
Definition: file-io.cc:2682
octave_idx_type a_nc
Definition: sylvester.cc:72
bool operator==(const FloatComplexDiagMatrix &a) const
Definition: fCDiagMatrix.cc:48
octave_idx_type rows(void) const
Definition: Array.h:401
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
bool swap
Definition: load-save.cc:725
bool mx_inline_all_real(size_t n, const std::complex< T > *x)
Definition: mx-inlines.cc:295
T dgelem(octave_idx_type i) const
Definition: DiagArray2.h:121
float min(void) const
Definition: fColVector.cc:244
bool all_elements_are_real(void) const
FloatComplexDiagMatrix pseudo_inverse(float tol=0.0f) const
octave_idx_type a_nr
Definition: sylvester.cc:71
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:760
Definition: DET.h:33
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_value retval
Definition: data.cc:6294
Template for two dimensional diagonal array with math operators.
Definition: MDiagArray2.h:33
T & dgxelem(octave_idx_type i)
Definition: DiagArray2.h:159
FloatComplexDiagMatrix operator*(const FloatComplexDiagMatrix &a, const FloatDiagMatrix &b)
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
FloatComplexColumnVector extract_diag(octave_idx_type k=0) const
Definition: fCDiagMatrix.h:138
With real return the complex result
Definition: data.cc:3375
FloatComplexDiagMatrix & operator+=(const FloatDiagMatrix &a)
FloatComplexDET determinant(void) const
FloatComplexDiagMatrix conj(const FloatComplexDiagMatrix &a)
octave_idx_type cols(void) const
Definition: DiagArray2.h:89
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave_idx_type b_nc
Definition: sylvester.cc:75
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
float max(void) const
Definition: fColVector.cc:260
octave_idx_type columns(void) const
Definition: DiagArray2.h:90
b
Definition: cellfun.cc:398
FloatDiagMatrix abs(void) const
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::ostream & operator<<(std::ostream &os, const FloatComplexDiagMatrix &a)
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:532
octave_idx_type cols(void) const
Definition: Array.h:409
FloatComplexMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
octave_idx_type length(void) const
Definition: DiagArray2.h:94