GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dDiagMatrix.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <ostream>
31 
32 #include "Array-util.h"
33 #include "lo-error.h"
34 #include "mx-base.h"
35 #include "mx-inlines.cc"
36 #include "oct-cmplx.h"
37 
38 // Diagonal Matrix class.
39 
40 bool
42 {
43  if (rows () != a.rows () || cols () != a.cols ())
44  return 0;
45 
46  return mx_inline_equal (length (), data (), a.data ());
47 }
48 
49 bool
51 {
52  return !(*this == a);
53 }
54 
56 DiagMatrix::fill (double val)
57 {
58  for (octave_idx_type i = 0; i < length (); i++)
59  elem (i, i) = val;
60  return *this;
61 }
62 
65 {
66  if (beg < 0 || end >= length () || end < beg)
67  (*current_liboctave_error_handler) ("range error for fill");
68 
69  for (octave_idx_type i = beg; i <= end; i++)
70  elem (i, i) = val;
71 
72  return *this;
73 }
74 
77 {
79  if (a.numel () != len)
80  (*current_liboctave_error_handler) ("range error for fill");
81 
82  for (octave_idx_type i = 0; i < len; i++)
83  elem (i, i) = a.elem (i);
84 
85  return *this;
86 }
87 
90 {
92  if (a.numel () != len)
93  (*current_liboctave_error_handler) ("range error for fill");
94 
95  for (octave_idx_type i = 0; i < len; i++)
96  elem (i, i) = a.elem (i);
97 
98  return *this;
99 }
100 
101 DiagMatrix&
103 {
104  octave_idx_type a_len = a.numel ();
105  if (beg < 0 || beg + a_len >= length ())
106  (*current_liboctave_error_handler) ("range error for fill");
107 
108  for (octave_idx_type i = 0; i < a_len; i++)
109  elem (i+beg, i+beg) = a.elem (i);
110 
111  return *this;
112 }
113 
114 DiagMatrix&
116 {
117  octave_idx_type a_len = a.numel ();
118  if (beg < 0 || beg + a_len >= length ())
119  (*current_liboctave_error_handler) ("range error for fill");
120 
121  for (octave_idx_type i = 0; i < a_len; i++)
122  elem (i+beg, i+beg) = a.elem (i);
123 
124  return *this;
125 }
126 
128 DiagMatrix::abs (void) const
129 {
130  return DiagMatrix (extract_diag ().abs (), rows (), columns ());
131 }
132 
135 {
136  return DiagMatrix (real (a.extract_diag ()), a.rows (), a.cols ());
137 }
138 
141 {
142  return DiagMatrix (imag (a.extract_diag ()), a.rows (), a.cols ());
143 }
144 
145 Matrix
147  octave_idx_type r2, octave_idx_type c2) const
148 {
149  if (r1 > r2) { std::swap (r1, r2); }
150  if (c1 > c2) { std::swap (c1, c2); }
151 
152  octave_idx_type new_r = r2 - r1 + 1;
153  octave_idx_type new_c = c2 - c1 + 1;
154 
155  Matrix result (new_r, new_c);
156 
157  for (octave_idx_type j = 0; j < new_c; j++)
158  for (octave_idx_type i = 0; i < new_r; i++)
159  result.elem (i, j) = elem (r1+i, c1+j);
160 
161  return result;
162 }
163 
164 // extract row or column i.
165 
166 RowVector
168 {
169  octave_idx_type r = rows ();
170  octave_idx_type c = cols ();
171  if (i < 0 || i >= r)
172  (*current_liboctave_error_handler) ("invalid row selection");
173 
174  RowVector retval (c, 0.0);
175  if (r <= c || i < c)
176  retval.elem (i) = elem (i, i);
177 
178  return retval;
179 }
180 
181 RowVector
182 DiagMatrix::row (char *s) const
183 {
184  if (! s)
185  (*current_liboctave_error_handler) ("invalid row selection");
186 
187  char c = s[0];
188  if (c == 'f' || c == 'F')
189  return row (static_cast<octave_idx_type> (0));
190  else if (c == 'l' || c == 'L')
191  return row (rows () - 1);
192  else
193  (*current_liboctave_error_handler) ("invalid row selection");
194 }
195 
198 {
199  octave_idx_type r = rows ();
200  octave_idx_type c = cols ();
201  if (i < 0 || i >= c)
202  (*current_liboctave_error_handler) ("invalid column selection");
203 
204  ColumnVector retval (r, 0.0);
205  if (r >= c || i < r)
206  retval.elem (i) = elem (i, i);
207 
208  return retval;
209 }
210 
212 DiagMatrix::column (char *s) const
213 {
214  if (! s)
215  (*current_liboctave_error_handler) ("invalid column selection");
216 
217  char c = s[0];
218  if (c == 'f' || c == 'F')
219  return column (static_cast<octave_idx_type> (0));
220  else if (c == 'l' || c == 'L')
221  return column (cols () - 1);
222  else
223  (*current_liboctave_error_handler) ("invalid column selection");
224 }
225 
228 {
229  octave_idx_type info;
230  return inverse (info);
231 }
232 
235 {
236  octave_idx_type r = rows ();
237  octave_idx_type c = cols ();
238  if (r != c)
239  (*current_liboctave_error_handler) ("inverse requires square matrix");
240 
241  DiagMatrix retval (r, c);
242 
243  info = 0;
244  octave_idx_type len = r; // alias for readability
245  octave_idx_type z_count = 0; // zeros
246  octave_idx_type nz_count = 0; // non-zeros
247  for (octave_idx_type i = 0; i < len; i++)
248  {
249  if (xelem (i, i) == 0.0)
250  {
251  z_count++;
252  if (nz_count > 0)
253  break;
254  }
255  else
256  {
257  nz_count++;
258  if (z_count > 0)
259  break;
260  retval.elem (i, i) = 1.0 / xelem (i, i);
261  }
262  }
263  if (nz_count == 0)
264  {
265  (*current_liboctave_error_handler)
266  ("inverse of the null matrix not defined");
267  }
268  else if (z_count > 0)
269  {
270  info = -1;
273  }
274 
275  return retval;
276 }
277 
279 DiagMatrix::pseudo_inverse (double tol) const
280 {
281  octave_idx_type r = rows ();
282  octave_idx_type c = cols ();
284 
285  DiagMatrix retval (c, r);
286 
287  for (octave_idx_type i = 0; i < len; i++)
288  {
289  double val = std::abs (elem (i, i));
290  if (val < tol || val == 0.0)
291  retval.elem (i, i) = 0.0;
292  else
293  retval.elem (i, i) = 1.0 / elem (i, i);
294  }
295 
296  return retval;
297 }
298 
299 // diagonal matrix by diagonal matrix -> diagonal matrix operations
300 
301 // diagonal matrix by diagonal matrix -> diagonal matrix operations
302 
304 operator * (const DiagMatrix& a, const DiagMatrix& b)
305 {
306  octave_idx_type a_nr = a.rows ();
307  octave_idx_type a_nc = a.cols ();
308 
309  octave_idx_type b_nr = b.rows ();
310  octave_idx_type b_nc = b.cols ();
311 
312  if (a_nc != b_nr)
313  octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
314 
315  DiagMatrix c (a_nr, b_nc);
316 
317  octave_idx_type len = c.length ();
318  octave_idx_type lenm = (len < a_nc ? len : a_nc);
319 
320  for (octave_idx_type i = 0; i < lenm; i++)
321  c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
322  for (octave_idx_type i = lenm; i < len; i++)
323  c.dgxelem (i) = 0.0;
324 
325  return c;
326 }
327 
328 // other operations
329 
330 DET
332 {
333  DET det (1.0);
334  if (rows () != cols ())
335  (*current_liboctave_error_handler) ("determinant requires square matrix");
336 
338  for (octave_idx_type i = 0; i < len; i++)
339  det *= elem (i, i);
340 
341  return det;
342 }
343 
344 double
345 DiagMatrix::rcond (void) const
346 {
347  ColumnVector av = extract_diag (0).map<double> (fabs);
348  double amx = av.max ();
349  double amn = av.min ();
350  return amx == 0 ? 0.0 : amn / amx;
351 }
352 
353 std::ostream&
354 operator << (std::ostream& os, const DiagMatrix& a)
355 {
356 // int field_width = os.precision () + 7;
357 
358  for (octave_idx_type i = 0; i < a.rows (); i++)
359  {
360  for (octave_idx_type j = 0; j < a.cols (); j++)
361  {
362  if (i == j)
363  os << ' ' /* setw (field_width) */ << a.elem (i, i);
364  else
365  os << ' ' /* setw (field_width) */ << 0.0;
366  }
367  os << "\n";
368  }
369  return os;
370 }
#define Inf
Definition: Faddeeva.cc:247
Array< U > map(F fcn) const
Apply function fcn to each element of the Array<T>.
Definition: Array.h:759
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:377
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:499
T element_type
Definition: Array.h:202
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
double max(void) const
Definition: dColVector.cc:261
double min(void) const
Definition: dColVector.cc:245
ComplexColumnVector extract_diag(octave_idx_type k=0) const
Definition: CDiagMatrix.h:136
T dgelem(octave_idx_type i) const
Definition: DiagArray2.h:123
T elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:116
const T * data(void) const
Definition: DiagArray2.h:168
T & dgxelem(octave_idx_type i)
Definition: DiagArray2.h:151
T xelem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:146
octave_idx_type length(void) const
Definition: DiagArray2.h:94
octave_idx_type cols(void) const
Definition: DiagArray2.h:89
octave_idx_type columns(void) const
Definition: DiagArray2.h:90
octave_idx_type rows(void) const
Definition: DiagArray2.h:88
DiagMatrix(void)=default
DET determinant(void) const
Definition: dDiagMatrix.cc:331
RowVector row(octave_idx_type i) const
Definition: dDiagMatrix.cc:167
Matrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition: dDiagMatrix.cc:146
DiagMatrix abs(void) const
Definition: dDiagMatrix.cc:128
DiagMatrix & fill(double val)
Definition: dDiagMatrix.cc:56
DiagMatrix pseudo_inverse(double tol=0.0) const
Definition: dDiagMatrix.cc:279
DiagMatrix inverse(void) const
Definition: dDiagMatrix.cc:227
double rcond(void) const
Definition: dDiagMatrix.cc:345
ColumnVector extract_diag(octave_idx_type k=0) const
Definition: dDiagMatrix.h:105
bool operator!=(const DiagMatrix &a) const
Definition: dDiagMatrix.cc:50
bool operator==(const DiagMatrix &a) const
Definition: dDiagMatrix.cc:41
ColumnVector column(octave_idx_type i) const
Definition: dDiagMatrix.cc:197
Definition: dMatrix.h:42
Definition: DET.h:39
DiagMatrix real(const ComplexDiagMatrix &a)
Definition: dDiagMatrix.cc:134
DiagMatrix imag(const ComplexDiagMatrix &a)
Definition: dDiagMatrix.cc:140
std::ostream & operator<<(std::ostream &os, const DiagMatrix &a)
Definition: dDiagMatrix.cc:354
DiagMatrix operator*(const DiagMatrix &a, const DiagMatrix &b)
Definition: dDiagMatrix.cc:304
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:41
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:570
T * r
Definition: mx-inlines.cc:773
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_value::octave_value(const Array< char > &chm, char type) return retval
Definition: ov.cc:811
static T abs(T x)
Definition: pr-output.cc:1678
F77_RET_T len
Definition: xerbla.cc:61