1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-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 // This file should not include config.h. It is only included in other
27 // C++ source files that should have included config.h before including
28 // this file.
29
30 #include <cassert>
31
32 #include <algorithm>
33
34 #include "DiagArray2.h"
35
36 #include "lo-error.h"
37
38 template <typename T>
DiagArray2(const Array<T> & a,octave_idx_type r,octave_idx_type c)39 DiagArray2<T>::DiagArray2 (const Array<T>& a, octave_idx_type r,
40 octave_idx_type c)
41 : Array<T> (a.as_column ()), d1 (r), d2 (c)
42 {
43 octave_idx_type rcmin = std::min (r, c);
44 if (rcmin != a.numel ())
45 Array<T>::resize (dim_vector (rcmin, 1));
46 }
47
48 template <typename T>
49 Array<T>
extract_diag(octave_idx_type k) const50 DiagArray2<T>::extract_diag (octave_idx_type k) const
51 {
52 Array<T> d;
53
54 if (k == 0)
55 // The main diagonal is shallow-copied.
56 d = *this;
57 else if (k > 0 && k < cols ())
58 d = Array<T> (dim_vector (std::min (cols () - k, rows ()), 1), T ());
59 else if (k < 0 && -k < rows ())
60 d = Array<T> (dim_vector (std::min (rows () + k, cols ()), 1), T ());
61 else // Matlab returns [] 0x1 for out-of-range diagonal
62 d.resize (dim_vector (0, 1));
63
64 return d;
65 }
66
67 template <typename T>
68 DiagArray2<T>
transpose(void) const69 DiagArray2<T>::transpose (void) const
70 {
71 return DiagArray2<T> (*this, d2, d1);
72 }
73
74 template <typename T>
75 DiagArray2<T>
hermitian(T (* fcn)(const T &)) const76 DiagArray2<T>::hermitian (T (* fcn) (const T&)) const
77 {
78 return DiagArray2<T> (Array<T>::template map<T> (fcn), d2, d1);
79 }
80
81 // A two-dimensional array with diagonal elements only.
82
83 template <typename T>
84 T&
elem(octave_idx_type r,octave_idx_type c)85 DiagArray2<T>::elem (octave_idx_type r, octave_idx_type c)
86 {
87 static T zero (0);
88 return (r == c) ? Array<T>::elem (r) : zero;
89 }
90
91 template <typename T>
92 T&
checkelem(octave_idx_type r,octave_idx_type c)93 DiagArray2<T>::checkelem (octave_idx_type r, octave_idx_type c)
94 {
95 static T zero (0);
96 return check_idx (r, c) ? elem (r, c) : zero;
97 }
98
99 template <typename T>
100 void
resize(octave_idx_type r,octave_idx_type c,const T & rfv)101 DiagArray2<T>::resize (octave_idx_type r, octave_idx_type c,
102 const T& rfv)
103 {
104 if (r < 0 || c < 0)
105 (*current_liboctave_error_handler) ("can't resize to negative dimensions");
106
107 if (r != dim1 () || c != dim2 ())
108 {
109 Array<T>::resize (dim_vector (std::min (r, c), 1), rfv);
110 d1 = r; d2 = c;
111 }
112 }
113
114 template <typename T>
115 Array<T>
array_value(void) const116 DiagArray2<T>::array_value (void) const
117 {
118 Array<T> result (dims (), T (0));
119
120 for (octave_idx_type i = 0, len = length (); i < len; i++)
121 result.xelem (i, i) = dgelem (i);
122
123 return result;
124 }
125
126 template <typename T>
127 bool
check_idx(octave_idx_type r,octave_idx_type c) const128 DiagArray2<T>::check_idx (octave_idx_type r, octave_idx_type c) const
129 {
130 bool ok = true;
131
132 if (r < 0 || r >= dim1 ())
133 octave::err_index_out_of_range (2, 1, r+1, dim1 (), dims ());
134
135 if (c < 0 || c >= dim2 ())
136 octave::err_index_out_of_range (2, 2, c+1, dim2 (), dims ());
137
138 return ok;
139 }
140