GNU Octave  3.8.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-2013 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  if (elem (i, i) != 0.0f)
306  retval.elem (i, i) = 1.0f / elem (i, i);
307  else
308  retval.elem (i, i) = 0.0f;
309  }
310 
311  return retval;
312 }
313 
314 // diagonal matrix by diagonal matrix -> diagonal matrix operations
315 
316 // diagonal matrix by diagonal matrix -> diagonal matrix operations
317 
320 {
321  octave_idx_type a_nr = a.rows ();
322  octave_idx_type a_nc = a.cols ();
323 
324  octave_idx_type b_nr = b.rows ();
325  octave_idx_type b_nc = b.cols ();
326 
327  if (a_nc != b_nr)
328  gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
329 
330  FloatDiagMatrix c (a_nr, b_nc);
331 
332  octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
333 
334  for (octave_idx_type i = 0; i < lenm; i++)
335  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
336  for (octave_idx_type i = lenm; i < len; i++)
337  c.dgxelem (i) = 0.0f;
338 
339  return c;
340 }
341 
342 // other operations
343 
344 FloatDET
346 {
347  FloatDET det (1.0f);
348  if (rows () != cols ())
349  {
350  (*current_liboctave_error_handler) ("determinant requires square matrix");
351  det = 0.0f;
352  }
353  else
354  {
355  octave_idx_type len = length ();
356  for (octave_idx_type i = 0; i < len; i++)
357  det *= elem (i, i);
358  }
359 
360  return det;
361 }
362 
363 float
365 {
366  FloatColumnVector av = extract_diag (0).map<float> (fabsf);
367  float amx = av.max (), amn = av.min ();
368  return amx == 0 ? 0.0f : amn / amx;
369 }
370 
371 std::ostream&
372 operator << (std::ostream& os, const FloatDiagMatrix& a)
373 {
374 // int field_width = os.precision () + 7;
375 
376  for (octave_idx_type i = 0; i < a.rows (); i++)
377  {
378  for (octave_idx_type j = 0; j < a.cols (); j++)
379  {
380  if (i == j)
381  os << " " /* setw (field_width) */ << a.elem (i, i);
382  else
383  os << " " /* setw (field_width) */ << 0.0;
384  }
385  os << "\n";
386  }
387  return os;
388 }