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
fDiagMatrix.cc
Go to the documentation of this file.
1 // FloatDiagMatrix 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 "mx-base.h"
34 #include "mx-inlines.cc"
35 #include "oct-cmplx.h"
36 
37 // Diagonal Matrix class.
38 
39 bool
41 {
42  if (rows () != a.rows () || cols () != a.cols ())
43  return 0;
44 
45  return mx_inline_equal (length (), data (), a.data ());
46 }
47 
48 bool
50 {
51  return !(*this == a);
52 }
53 
56 {
57  for (octave_idx_type i = 0; i < length (); i++)
58  elem (i, i) = val;
59  return *this;
60 }
61 
64 {
65  if (beg < 0 || end >= length () || end < beg)
66  (*current_liboctave_error_handler) ("range error for fill");
67 
68  for (octave_idx_type i = beg; i <= end; i++)
69  elem (i, i) = val;
70 
71  return *this;
72 }
73 
76 {
77  octave_idx_type len = length ();
78  if (a.numel () != len)
79  (*current_liboctave_error_handler) ("range error for fill");
80 
81  for (octave_idx_type i = 0; i < len; i++)
82  elem (i, i) = a.elem (i);
83 
84  return *this;
85 }
86 
89 {
90  octave_idx_type len = length ();
91  if (a.numel () != len)
92  (*current_liboctave_error_handler) ("range error for fill");
93 
94  for (octave_idx_type i = 0; i < len; i++)
95  elem (i, i) = a.elem (i);
96 
97  return *this;
98 }
99 
102 {
103  octave_idx_type a_len = a.numel ();
104  if (beg < 0 || beg + a_len >= length ())
105  (*current_liboctave_error_handler) ("range error for fill");
106 
107  for (octave_idx_type i = 0; i < a_len; i++)
108  elem (i+beg, i+beg) = a.elem (i);
109 
110  return *this;
111 }
112 
115 {
116  octave_idx_type a_len = a.numel ();
117  if (beg < 0 || beg + a_len >= length ())
118  (*current_liboctave_error_handler) ("range error for fill");
119 
120  for (octave_idx_type i = 0; i < a_len; i++)
121  elem (i+beg, i+beg) = a.elem (i);
122 
123  return *this;
124 }
125 
128 {
129  return FloatDiagMatrix (extract_diag ().abs (), rows (), columns ());
130 }
131 
134 {
135  return FloatDiagMatrix (real (a.extract_diag ()), a.rows (), a.columns ());
136 }
137 
140 {
141  return FloatDiagMatrix (imag (a.extract_diag ()), a.rows (), a.columns ());
142 }
143 
146  octave_idx_type r2, octave_idx_type c2) const
147 {
148  if (r1 > r2) { std::swap (r1, r2); }
149  if (c1 > c2) { std::swap (c1, c2); }
150 
151  octave_idx_type new_r = r2 - r1 + 1;
152  octave_idx_type new_c = c2 - c1 + 1;
153 
154  FloatMatrix result (new_r, new_c);
155 
156  for (octave_idx_type j = 0; j < new_c; j++)
157  for (octave_idx_type i = 0; i < new_r; i++)
158  result.elem (i, j) = elem (r1+i, c1+j);
159 
160  return result;
161 }
162 
163 // extract row or column i.
164 
167 {
168  octave_idx_type r = rows ();
169  octave_idx_type c = cols ();
170  if (i < 0 || i >= r)
171  (*current_liboctave_error_handler) ("invalid row selection");
172 
173  FloatRowVector retval (c, 0.0);
174  if (r <= c || (r > c && i < c))
175  retval.elem (i) = elem (i, i);
176 
177  return retval;
178 }
179 
181 FloatDiagMatrix::row (char *s) const
182 {
183  if (! s)
184  (*current_liboctave_error_handler) ("invalid row selection");
185 
186  char c = *s;
187  if (c == 'f' || c == 'F')
188  return row (static_cast<octave_idx_type>(0));
189  else if (c == 'l' || c == 'L')
190  return row (rows () - 1);
191  else
192  (*current_liboctave_error_handler) ("invalid row selection");
193 }
194 
197 {
198  octave_idx_type r = rows ();
199  octave_idx_type c = cols ();
200  if (i < 0 || i >= c)
201  (*current_liboctave_error_handler) ("invalid column selection");
202 
203  FloatColumnVector retval (r, 0.0);
204  if (r >= c || (r < c && i < r))
205  retval.elem (i) = elem (i, i);
206 
207  return retval;
208 }
209 
212 {
213  if (! s)
214  (*current_liboctave_error_handler) ("invalid column selection");
215 
216  char c = *s;
217  if (c == 'f' || c == 'F')
218  return column (static_cast<octave_idx_type>(0));
219  else if (c == 'l' || c == 'L')
220  return column (cols () - 1);
221  else
222  (*current_liboctave_error_handler) ("invalid column selection");
223 }
224 
227 {
228  octave_idx_type info;
229  return inverse (info);
230 }
231 
234 {
235  octave_idx_type r = rows ();
236  octave_idx_type c = cols ();
237  octave_idx_type len = length ();
238  if (r != c)
239  (*current_liboctave_error_handler) ("inverse requires square matrix");
240 
241  FloatDiagMatrix retval (r, c);
242 
243  info = 0;
244  for (octave_idx_type i = 0; i < len; i++)
245  {
246  if (elem (i, i) == 0.0)
247  retval.elem (i, i) = octave::numeric_limits<float>::Inf ();
248  else
249  retval.elem (i, i) = 1.0 / elem (i, i);
250  }
251 
252  return retval;
253 }
254 
257 {
258  octave_idx_type r = rows ();
259  octave_idx_type c = cols ();
260  octave_idx_type len = length ();
261 
262  FloatDiagMatrix retval (c, r);
263 
264  for (octave_idx_type i = 0; i < len; i++)
265  {
266  float val = std::abs (elem (i, i));
267  if (val < tol || val == 0.0f)
268  retval.elem (i, i) = 0.0f;
269  else
270  retval.elem (i, i) = 1.0f / elem (i, i);
271  }
272 
273  return retval;
274 }
275 
276 // diagonal matrix by diagonal matrix -> diagonal matrix operations
277 
278 // diagonal matrix by diagonal matrix -> diagonal matrix operations
279 
282 {
283  octave_idx_type a_nr = a.rows ();
284  octave_idx_type a_nc = a.cols ();
285 
286  octave_idx_type b_nr = b.rows ();
287  octave_idx_type b_nc = b.cols ();
288 
289  if (a_nc != b_nr)
290  octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
291 
292  FloatDiagMatrix c (a_nr, b_nc);
293 
294  octave_idx_type len = c.length ();
295  octave_idx_type lenm = len < a_nc ? len : a_nc;
296 
297  for (octave_idx_type i = 0; i < lenm; i++)
298  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
299  for (octave_idx_type i = lenm; i < len; i++)
300  c.dgxelem (i) = 0.0f;
301 
302  return c;
303 }
304 
305 // other operations
306 
307 FloatDET
309 {
310  FloatDET det (1.0f);
311  if (rows () != cols ())
312  (*current_liboctave_error_handler) ("determinant requires square matrix");
313 
314  octave_idx_type len = length ();
315  for (octave_idx_type i = 0; i < len; i++)
316  det *= elem (i, i);
317 
318  return det;
319 }
320 
321 float
323 {
324  FloatColumnVector av = extract_diag (0).map<float> (fabsf);
325  float amx = av.max ();
326  float amn = av.min ();
327  return amx == 0 ? 0.0f : amn / amx;
328 }
329 
330 std::ostream&
331 operator << (std::ostream& os, const FloatDiagMatrix& a)
332 {
333 // int field_width = os.precision () + 7;
334 
335  for (octave_idx_type i = 0; i < a.rows (); i++)
336  {
337  for (octave_idx_type j = 0; j < a.cols (); j++)
338  {
339  if (i == j)
340  os << " " /* setw (field_width) */ << a.elem (i, i);
341  else
342  os << " " /* setw (field_width) */ << 0.0;
343  }
344  os << "\n";
345  }
346  return os;
347 }
FloatDiagMatrix & fill(float val)
Definition: fDiagMatrix.cc:55
float elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:114
FloatDiagMatrix operator*(const FloatDiagMatrix &a, const FloatDiagMatrix &b)
Definition: fDiagMatrix.cc:281
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
FloatDiagMatrix abs(void) const
Definition: fDiagMatrix.cc:127
octave_idx_type b_nr
Definition: sylvester.cc:74
octave_idx_type rows(void) const
Definition: DiagArray2.h:88
const float * data(void) const
Definition: DiagArray2.h:176
FloatDET determinant(void) const
Definition: fDiagMatrix.cc:308
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
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
T dgelem(octave_idx_type i) const
Definition: DiagArray2.h:121
float min(void) const
Definition: fColVector.cc:244
octave_idx_type a_nr
Definition: sylvester.cc:71
Array< U > map(F fcn) const
Apply function fcn to each element of the Array.
Definition: Array.h:760
FloatMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition: fDiagMatrix.cc:145
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
FloatDiagMatrix real(const FloatComplexDiagMatrix &a)
Definition: fDiagMatrix.cc:133
float rcond(void) const
Definition: fDiagMatrix.cc:322
T & dgxelem(octave_idx_type i)
Definition: DiagArray2.h:159
FloatColumnVector extract_diag(octave_idx_type k=0) const
Definition: fDiagMatrix.h:105
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
bool operator!=(const FloatDiagMatrix &a) const
Definition: fDiagMatrix.cc:49
FloatComplexColumnVector extract_diag(octave_idx_type k=0) const
Definition: fCDiagMatrix.h:138
With real return the complex result
Definition: data.cc:3375
FloatDiagMatrix inverse(void) const
Definition: fDiagMatrix.cc:226
FloatDiagMatrix pseudo_inverse(float tol=0.0f) const
Definition: fDiagMatrix.cc:256
octave_idx_type cols(void) const
Definition: DiagArray2.h:89
#define Inf
Definition: Faddeeva.cc:247
=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
FloatColumnVector column(octave_idx_type i) const
Definition: fDiagMatrix.cc:196
octave_idx_type columns(void) const
Definition: DiagArray2.h:90
b
Definition: cellfun.cc:398
bool operator==(const FloatDiagMatrix &a) const
Definition: fDiagMatrix.cc:40
FloatDiagMatrix(void)
Definition: fDiagMatrix.h:42
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:532
FloatDiagMatrix imag(const FloatComplexDiagMatrix &a)
Definition: fDiagMatrix.cc:139
FloatRowVector row(octave_idx_type i) const
Definition: fDiagMatrix.cc:166
octave_idx_type length(void) const
Definition: DiagArray2.h:94
std::ostream & operator<<(std::ostream &os, const FloatDiagMatrix &a)
Definition: fDiagMatrix.cc:331