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 <istream>
31 #include <ostream>
32
33 #include "Array-util.h"
34 #include "lo-blas-proto.h"
35 #include "lo-error.h"
36 #include "mx-base.h"
37 #include "mx-inlines.cc"
38 #include "oct-cmplx.h"
39
40 // Column Vector class.
41
42 bool
operator ==(const ColumnVector & a) const43 ColumnVector::operator == (const ColumnVector& a) const
44 {
45 octave_idx_type len = numel ();
46 if (len != a.numel ())
47 return 0;
48 return mx_inline_equal (len, data (), a.data ());
49 }
50
51 bool
operator !=(const ColumnVector & a) const52 ColumnVector::operator != (const ColumnVector& a) const
53 {
54 return !(*this == a);
55 }
56
57 ColumnVector&
insert(const ColumnVector & a,octave_idx_type r)58 ColumnVector::insert (const ColumnVector& a, octave_idx_type r)
59 {
60 octave_idx_type a_len = a.numel ();
61
62 if (r < 0 || r + a_len > numel ())
63 (*current_liboctave_error_handler) ("range error for insert");
64
65 if (a_len > 0)
66 {
67 make_unique ();
68
69 for (octave_idx_type i = 0; i < a_len; i++)
70 xelem (r+i) = a.elem (i);
71 }
72
73 return *this;
74 }
75
76 ColumnVector&
fill(double val)77 ColumnVector::fill (double val)
78 {
79 octave_idx_type len = numel ();
80
81 if (len > 0)
82 {
83 make_unique ();
84
85 for (octave_idx_type i = 0; i < len; i++)
86 xelem (i) = val;
87 }
88
89 return *this;
90 }
91
92 ColumnVector&
fill(double val,octave_idx_type r1,octave_idx_type r2)93 ColumnVector::fill (double val, octave_idx_type r1, octave_idx_type r2)
94 {
95 octave_idx_type len = numel ();
96
97 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
98 (*current_liboctave_error_handler) ("range error for fill");
99
100 if (r1 > r2) { std::swap (r1, r2); }
101
102 if (r2 >= r1)
103 {
104 make_unique ();
105
106 for (octave_idx_type i = r1; i <= r2; i++)
107 xelem (i) = val;
108 }
109
110 return *this;
111 }
112
113 ColumnVector
stack(const ColumnVector & a) const114 ColumnVector::stack (const ColumnVector& a) const
115 {
116 octave_idx_type len = numel ();
117 octave_idx_type nr_insert = len;
118 ColumnVector retval (len + a.numel ());
119 retval.insert (*this, 0);
120 retval.insert (a, nr_insert);
121 return retval;
122 }
123
124 RowVector
transpose(void) const125 ColumnVector::transpose (void) const
126 {
127 return MArray<double>::transpose ();
128 }
129
130 ColumnVector
abs(void) const131 ColumnVector::abs (void) const
132 {
133 return do_mx_unary_map<double, double, std::abs> (*this);
134 }
135
136 ColumnVector
real(const ComplexColumnVector & a)137 real (const ComplexColumnVector& a)
138 {
139 return do_mx_unary_op<double, Complex> (a, mx_inline_real);
140 }
141
142 ColumnVector
imag(const ComplexColumnVector & a)143 imag (const ComplexColumnVector& a)
144 {
145 return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
146 }
147
148 // resize is the destructive equivalent for this one
149
150 ColumnVector
extract(octave_idx_type r1,octave_idx_type r2) const151 ColumnVector::extract (octave_idx_type r1, octave_idx_type r2) const
152 {
153 if (r1 > r2) { std::swap (r1, r2); }
154
155 octave_idx_type new_r = r2 - r1 + 1;
156
157 ColumnVector result (new_r);
158
159 for (octave_idx_type i = 0; i < new_r; i++)
160 result.xelem (i) = elem (r1+i);
161
162 return result;
163 }
164
165 ColumnVector
extract_n(octave_idx_type r1,octave_idx_type n) const166 ColumnVector::extract_n (octave_idx_type r1, octave_idx_type n) const
167 {
168 ColumnVector result (n);
169
170 for (octave_idx_type i = 0; i < n; i++)
171 result.xelem (i) = elem (r1+i);
172
173 return result;
174 }
175
176 // matrix by column vector -> column vector operations
177
178 ColumnVector
operator *(const Matrix & m,const ColumnVector & a)179 operator * (const Matrix& m, const ColumnVector& a)
180 {
181 ColumnVector retval;
182
183 F77_INT nr = octave::to_f77_int (m.rows ());
184 F77_INT nc = octave::to_f77_int (m.cols ());
185
186 F77_INT a_len = octave::to_f77_int (a.numel ());
187
188 if (nc != a_len)
189 octave::err_nonconformant ("operator *", nr, nc, a_len, 1);
190
191 retval.clear (nr);
192
193 if (nr != 0)
194 {
195 if (nc == 0)
196 retval.fill (0.0);
197 else
198 {
199 double *y = retval.fortran_vec ();
200
201 F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
202 nr, nc, 1.0, m.data (), nr,
203 a.data (), 1, 0.0, y, 1
204 F77_CHAR_ARG_LEN (1)));
205 }
206 }
207
208 return retval;
209 }
210
211 // diagonal matrix by column vector -> column vector operations
212
213 ColumnVector
operator *(const DiagMatrix & m,const ColumnVector & a)214 operator * (const DiagMatrix& m, const ColumnVector& a)
215 {
216 ColumnVector retval;
217
218 F77_INT nr = octave::to_f77_int (m.rows ());
219 F77_INT nc = octave::to_f77_int (m.cols ());
220
221 F77_INT a_len = octave::to_f77_int (a.numel ());
222
223 if (nc != a_len)
224 octave::err_nonconformant ("operator *", nr, nc, a_len, 1);
225
226 if (nr == 0 || nc == 0)
227 retval.resize (nr, 0.0);
228 else
229 {
230 retval.resize (nr);
231
232 for (octave_idx_type i = 0; i < a_len; i++)
233 retval.elem (i) = a.elem (i) * m.elem (i, i);
234
235 for (octave_idx_type i = a_len; i < nr; i++)
236 retval.elem (i) = 0.0;
237 }
238
239 return retval;
240 }
241
242 // other operations
243
244 double
min(void) const245 ColumnVector::min (void) const
246 {
247 octave_idx_type len = numel ();
248 if (len == 0)
249 return 0.0;
250
251 double res = elem (0);
252
253 for (octave_idx_type i = 1; i < len; i++)
254 if (elem (i) < res)
255 res = elem (i);
256
257 return res;
258 }
259
260 double
max(void) const261 ColumnVector::max (void) const
262 {
263 octave_idx_type len = numel ();
264 if (len == 0)
265 return 0.0;
266
267 double res = elem (0);
268
269 for (octave_idx_type i = 1; i < len; i++)
270 if (elem (i) > res)
271 res = elem (i);
272
273 return res;
274 }
275
276 std::ostream&
operator <<(std::ostream & os,const ColumnVector & a)277 operator << (std::ostream& os, const ColumnVector& a)
278 {
279 // int field_width = os.precision () + 7;
280 for (octave_idx_type i = 0; i < a.numel (); i++)
281 os << /* setw (field_width) << */ a.elem (i) << "\n";
282 return os;
283 }
284
285 std::istream&
operator >>(std::istream & is,ColumnVector & a)286 operator >> (std::istream& is, ColumnVector& a)
287 {
288 octave_idx_type len = a.numel ();
289
290 if (len > 0)
291 {
292 double tmp;
293 for (octave_idx_type i = 0; i < len; i++)
294 {
295 is >> tmp;
296 if (is)
297 a.elem (i) = tmp;
298 else
299 break;
300 }
301 }
302 return is;
303 }
304