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
fDiagMatrix.cc
Go to the documentation of this file.
1 // FloatDiagMatrix 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 "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  {
67  (*current_liboctave_error_handler) ("range error for fill");
68  return *this;
69  }
70 
71  for (octave_idx_type i = beg; i <= end; i++)
72  elem (i, i) = val;
73 
74  return *this;
75 }
76 
79 {
80  octave_idx_type len = length ();
81  if (a.length () != len)
82  {
83  (*current_liboctave_error_handler) ("range error for fill");
84  return *this;
85  }
86 
87  for (octave_idx_type i = 0; i < len; i++)
88  elem (i, i) = a.elem (i);
89 
90  return *this;
91 }
92 
95 {
96  octave_idx_type len = length ();
97  if (a.length () != len)
98  {
99  (*current_liboctave_error_handler) ("range error for fill");
100  return *this;
101  }
102 
103  for (octave_idx_type i = 0; i < len; i++)
104  elem (i, i) = a.elem (i);
105 
106  return *this;
107 }
108 
111 {
112  octave_idx_type a_len = a.length ();
113  if (beg < 0 || beg + a_len >= length ())
114  {
115  (*current_liboctave_error_handler) ("range error for fill");
116  return *this;
117  }
118 
119  for (octave_idx_type i = 0; i < a_len; i++)
120  elem (i+beg, i+beg) = a.elem (i);
121 
122  return *this;
123 }
124 
127 {
128  octave_idx_type a_len = a.length ();
129  if (beg < 0 || beg + a_len >= length ())
130  {
131  (*current_liboctave_error_handler) ("range error for fill");
132  return *this;
133  }
134 
135  for (octave_idx_type i = 0; i < a_len; i++)
136  elem (i+beg, i+beg) = a.elem (i);
137 
138  return *this;
139 }
140 
143 {
144  return FloatDiagMatrix (extract_diag ().abs (), rows (), columns ());
145 }
146 
149 {
150  return FloatDiagMatrix (real (a.extract_diag ()), a.rows (), a.columns ());
151 }
152 
155 {
156  return FloatDiagMatrix (imag (a.extract_diag ()), a.rows (), a.columns ());
157 }
158 
162 {
163  if (r1 > r2) { std::swap (r1, r2); }
164  if (c1 > c2) { std::swap (c1, c2); }
165 
166  octave_idx_type new_r = r2 - r1 + 1;
167  octave_idx_type new_c = c2 - c1 + 1;
168 
169  FloatMatrix result (new_r, new_c);
170 
171  for (octave_idx_type j = 0; j < new_c; j++)
172  for (octave_idx_type i = 0; i < new_r; i++)
173  result.elem (i, j) = elem (r1+i, c1+j);
174 
175  return result;
176 }
177 
178 // extract row or column i.
179 
182 {
183  octave_idx_type r = rows ();
184  octave_idx_type c = cols ();
185  if (i < 0 || i >= r)
186  {
187  (*current_liboctave_error_handler) ("invalid row selection");
188  return FloatRowVector ();
189  }
190 
191  FloatRowVector retval (c, 0.0);
192  if (r <= c || (r > c && i < c))
193  retval.elem (i) = elem (i, i);
194 
195  return retval;
196 }
197 
199 FloatDiagMatrix::row (char *s) const
200 {
201  if (! s)
202  {
203  (*current_liboctave_error_handler) ("invalid row selection");
204  return FloatRowVector ();
205  }
206 
207  char c = *s;
208  if (c == 'f' || c == 'F')
209  return row (static_cast<octave_idx_type>(0));
210  else if (c == 'l' || c == 'L')
211  return row (rows () - 1);
212  else
213  {
214  (*current_liboctave_error_handler) ("invalid row selection");
215  return FloatRowVector ();
216  }
217 }
218 
221 {
222  octave_idx_type r = rows ();
223  octave_idx_type c = cols ();
224  if (i < 0 || i >= c)
225  {
226  (*current_liboctave_error_handler) ("invalid column selection");
227  return FloatColumnVector ();
228  }
229 
230  FloatColumnVector retval (r, 0.0);
231  if (r >= c || (r < c && i < r))
232  retval.elem (i) = elem (i, i);
233 
234  return retval;
235 }
236 
238 FloatDiagMatrix::column (char *s) const
239 {
240  if (! s)
241  {
242  (*current_liboctave_error_handler) ("invalid column selection");
243  return FloatColumnVector ();
244  }
245 
246  char c = *s;
247  if (c == 'f' || c == 'F')
248  return column (static_cast<octave_idx_type>(0));
249  else if (c == 'l' || c == 'L')
250  return column (cols () - 1);
251  else
252  {
253  (*current_liboctave_error_handler) ("invalid column selection");
254  return FloatColumnVector ();
255  }
256 }
257 
260 {
261  octave_idx_type info;
262  return inverse (info);
263 }
264 
267 {
268  octave_idx_type r = rows ();
269  octave_idx_type c = cols ();
270  octave_idx_type len = length ();
271  if (r != c)
272  {
273  (*current_liboctave_error_handler) ("inverse requires square matrix");
274  return FloatDiagMatrix ();
275  }
276 
277  FloatDiagMatrix retval (r, c);
278 
279  info = 0;
280  for (octave_idx_type i = 0; i < len; i++)
281  {
282  if (elem (i, i) == 0.0)
283  {
284  info = -1;
285  return *this;
286  }
287  else
288  retval.elem (i, i) = 1.0 / elem (i, i);
289  }
290 
291  return retval;
292 }
293 
296 {
297  octave_idx_type r = rows ();
298  octave_idx_type c = cols ();
299  octave_idx_type len = length ();
300 
301  FloatDiagMatrix retval (c, r);
302 
303  for (octave_idx_type i = 0; i < len; i++)
304  {
305  float val = std::abs (elem (i, i));
306  if (val < tol || val == 0.0f)
307  retval.elem (i, i) = 0.0f;
308  else
309  retval.elem (i, i) = 1.0f / elem (i, i);
310  }
311 
312  return retval;
313 }
314 
315 // diagonal matrix by diagonal matrix -> diagonal matrix operations
316 
317 // diagonal matrix by diagonal matrix -> diagonal matrix operations
318 
321 {
322  octave_idx_type a_nr = a.rows ();
323  octave_idx_type a_nc = a.cols ();
324 
325  octave_idx_type b_nr = b.rows ();
326  octave_idx_type b_nc = b.cols ();
327 
328  if (a_nc != b_nr)
329  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
330 
331  FloatDiagMatrix c (a_nr, b_nc);
332 
333  octave_idx_type len = c.length ();
334  octave_idx_type lenm = len < a_nc ? len : a_nc;
335 
336  for (octave_idx_type i = 0; i < lenm; i++)
337  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
338  for (octave_idx_type i = lenm; i < len; i++)
339  c.dgxelem (i) = 0.0f;
340 
341  return c;
342 }
343 
344 // other operations
345 
346 FloatDET
348 {
349  FloatDET det (1.0f);
350  if (rows () != cols ())
351  {
352  (*current_liboctave_error_handler) ("determinant requires square matrix");
353  det = 0.0f;
354  }
355  else
356  {
357  octave_idx_type len = length ();
358  for (octave_idx_type i = 0; i < len; i++)
359  det *= elem (i, i);
360  }
361 
362  return det;
363 }
364 
365 float
367 {
368  FloatColumnVector av = extract_diag (0).map<float> (fabsf);
369  float amx = av.max ();
370  float amn = av.min ();
371  return amx == 0 ? 0.0f : amn / amx;
372 }
373 
374 std::ostream&
375 operator << (std::ostream& os, const FloatDiagMatrix& a)
376 {
377 // int field_width = os.precision () + 7;
378 
379  for (octave_idx_type i = 0; i < a.rows (); i++)
380  {
381  for (octave_idx_type j = 0; j < a.cols (); j++)
382  {
383  if (i == j)
384  os << " " /* setw (field_width) */ << a.elem (i, i);
385  else
386  os << " " /* setw (field_width) */ << 0.0;
387  }
388  os << "\n";
389  }
390  return os;
391 }
FloatDiagMatrix & fill(float val)
Definition: fDiagMatrix.cc:55
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
float elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:110
FloatDiagMatrix operator*(const FloatDiagMatrix &a, const FloatDiagMatrix &b)
Definition: fDiagMatrix.cc:320
FloatDiagMatrix abs(void) const
Definition: fDiagMatrix.cc:142
octave_idx_type rows(void) const
Definition: DiagArray2.h:86
const float * data(void) const
Definition: DiagArray2.h:180
FloatDET determinant(void) const
Definition: fDiagMatrix.cc:347
T & elem(octave_idx_type n)
Definition: Array.h:380
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type r2
T dgelem(octave_idx_type i) const
Definition: DiagArray2.h:121
float min(void) const
Definition: fColVector.cc:267
F77_RET_T const double const double * f
Array< U > map(F fcn) const
Apply function fcn to each element of the Array.
Definition: Array.h:659
FloatMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition: fDiagMatrix.cc:160
Definition: DET.h:31
FloatDiagMatrix real(const FloatComplexDiagMatrix &a)
Definition: fDiagMatrix.cc:148
float rcond(void) const
Definition: fDiagMatrix.cc:366
T & dgxelem(octave_idx_type i)
Definition: DiagArray2.h:163
FloatColumnVector extract_diag(octave_idx_type k=0) const
Definition: fDiagMatrix.h:106
bool operator!=(const FloatDiagMatrix &a) const
Definition: fDiagMatrix.cc:49
FloatComplexColumnVector extract_diag(octave_idx_type k=0) const
Definition: fCDiagMatrix.h:136
FloatDiagMatrix inverse(void) const
Definition: fDiagMatrix.cc:259
FloatDiagMatrix pseudo_inverse(float tol=0.0f) const
Definition: fDiagMatrix.cc:295
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
FloatColumnVector column(octave_idx_type i) const
Definition: fDiagMatrix.cc:220
octave_idx_type columns(void) const
Definition: DiagArray2.h:88
bool operator==(const FloatDiagMatrix &a) const
Definition: fDiagMatrix.cc:40
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
FloatDiagMatrix(void)
Definition: fDiagMatrix.h:43
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:441
FloatDiagMatrix imag(const FloatComplexDiagMatrix &a)
Definition: fDiagMatrix.cc:154
T abs(T x)
Definition: pr-output.cc:3062
FloatRowVector row(octave_idx_type i) const
Definition: fDiagMatrix.cc:181
octave_idx_type length(void) const
Definition: DiagArray2.h:92
std::ostream & operator<<(std::ostream &os, const FloatDiagMatrix &a)
Definition: fDiagMatrix.cc:375