GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dRowVector.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2024 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 #include <type_traits>
33 
34 #include "Array-util.h"
35 #include "lo-blas-proto.h"
36 #include "lo-error.h"
37 #include "mx-base.h"
38 #include "mx-inlines.cc"
39 #include "oct-cmplx.h"
40 
41 // Row Vector class.
42 
43 bool
45 {
47  if (len != a.numel ())
48  return 0;
49  return mx_inline_equal (len, data (), a.data ());
50 }
51 
52 bool
54 {
55  return !(*this == a);
56 }
57 
58 RowVector&
60 {
61  octave_idx_type a_len = a.numel ();
62 
63  if (c < 0 || c + a_len > numel ())
64  (*current_liboctave_error_handler) ("range error for insert");
65 
66  if (a_len > 0)
67  {
68  make_unique ();
69 
70  for (octave_idx_type i = 0; i < a_len; i++)
71  xelem (c+i) = a.elem (i);
72  }
73 
74  return *this;
75 }
76 
77 RowVector&
78 RowVector::fill (double val)
79 {
81 
82  if (len > 0)
83  {
84  make_unique ();
85 
86  for (octave_idx_type i = 0; i < len; i++)
87  xelem (i) = val;
88  }
89 
90  return *this;
91 }
92 
93 RowVector&
95 {
97 
98  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
99  (*current_liboctave_error_handler) ("range error for fill");
100 
101  if (c1 > c2) { std::swap (c1, c2); }
102 
103  if (c2 >= c1)
104  {
105  make_unique ();
106 
107  for (octave_idx_type i = c1; i <= c2; i++)
108  xelem (i) = val;
109  }
110 
111  return *this;
112 }
113 
114 RowVector
115 RowVector::append (const RowVector& a) const
116 {
118  octave_idx_type nc_insert = len;
119  RowVector retval (len + a.numel ());
120  retval.insert (*this, 0);
121  retval.insert (a, nc_insert);
122  return retval;
123 }
124 
127 {
128  return MArray<double>::transpose ();
129 }
130 
131 RowVector
133 {
134  return do_mx_unary_op<double, Complex> (a, mx_inline_real);
135 }
136 
137 RowVector
139 {
140  return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
141 }
142 
143 RowVector
145 {
146  if (c1 > c2) { std::swap (c1, c2); }
147 
148  octave_idx_type new_c = c2 - c1 + 1;
149 
150  RowVector result (new_c);
151 
152  for (octave_idx_type i = 0; i < new_c; i++)
153  result.xelem (i) = elem (c1+i);
154 
155  return result;
156 }
157 
158 RowVector
160 {
161  RowVector result (n);
162 
163  for (octave_idx_type i = 0; i < n; i++)
164  result.xelem (i) = elem (r1+i);
165 
166  return result;
167 }
168 
169 // row vector by matrix -> row vector
170 
171 RowVector
172 operator * (const RowVector& v, const Matrix& a)
173 {
174  RowVector retval;
175 
176  F77_INT len = octave::to_f77_int (v.numel ());
177 
178  F77_INT a_nr = octave::to_f77_int (a.rows ());
179  F77_INT a_nc = octave::to_f77_int (a.cols ());
180 
181  if (a_nr != len)
182  octave::err_nonconformant ("operator *", 1, len, a_nr, a_nc);
183 
184  if (len == 0)
185  retval.resize (a_nc, 0.0);
186  else
187  {
188  // Transpose A to form A'*x == (x'*A)'
189 
190  F77_INT ld = a_nr;
191 
192  retval.resize (a_nc);
193  double *y = retval.fortran_vec ();
194 
195  F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
196  a_nr, a_nc, 1.0, a.data (),
197  ld, v.data (), 1, 0.0, y, 1
198  F77_CHAR_ARG_LEN (1)));
199  }
200 
201  return retval;
202 }
203 
204 // other operations
205 
206 double
208 {
210  if (len == 0)
211  return 0;
212 
213  double res = elem (0);
214 
215  for (octave_idx_type i = 1; i < len; i++)
216  if (elem (i) < res)
217  res = elem (i);
218 
219  return res;
220 }
221 
222 double
224 {
226  if (len == 0)
227  return 0;
228 
229  double res = elem (0);
230 
231  for (octave_idx_type i = 1; i < len; i++)
232  if (elem (i) > res)
233  res = elem (i);
234 
235  return res;
236 }
237 
238 std::ostream&
239 operator << (std::ostream& os, const RowVector& a)
240 {
241 // int field_width = os.precision () + 7;
242 
243  for (octave_idx_type i = 0; i < a.numel (); i++)
244  os << ' ' /* setw (field_width) */ << a.elem (i);
245  return os;
246 }
247 
248 std::istream&
249 operator >> (std::istream& is, RowVector& a)
250 {
251  octave_idx_type len = a.numel ();
252 
253  if (len > 0)
254  {
255  double tmp;
256  for (octave_idx_type i = 0; i < len; i++)
257  {
258  is >> tmp;
259  if (is)
260  a.elem (i) = tmp;
261  else
262  break;
263  }
264  }
265  return is;
266 }
267 
268 // other operations
269 
270 RowVector
271 linspace (double x1, double x2, octave_idx_type n_in)
272 {
273  RowVector retval;
274 
275  if (n_in < 1)
276  return retval;
277  else if (n_in == 1)
278  {
279  retval.resize (1, x2);
280  return retval;
281  }
282  else if (x1 == x2)
283  {
284  retval.resize (n_in, x2);
285  return retval;
286  }
287 
288  // Use unsigned type (guaranteed n_in > 1 at this point) so that divisions
289  // by 2 can be replaced by compiler with shift right instructions.
290  typedef std::make_unsigned<octave_idx_type>::type unsigned_octave_idx_type;
291 
292  unsigned_octave_idx_type n = n_in;
293 
294  // Set endpoints, rather than calculate, for maximum accuracy.
295  retval.clear (n);
296  retval.xelem (0) = x1;
297  retval.xelem (n-1) = x2;
298 
299  // Construct linspace symmetrically from both ends.
300  bool isnan_delta = false;
301  double delta = (x2 - x1) / (n - 1);
302  if (octave::math::isinf (delta))
303  {
305  isnan_delta = true;
306  }
307 
308  unsigned_octave_idx_type n2 = n/2;
309  for (unsigned_octave_idx_type i = 1; i < n2; i++)
310  {
311  retval.xelem (i) = x1 + i*delta;
312  retval.xelem (n-1-i) = x2 - i*delta;
313  }
314  if (n % 2 == 1) // Middle element if number of elements is odd.
315  {
316  if (x1 == -x2)
317  retval.xelem (n2) = 0;
318  else if (isnan_delta)
320  else
321  retval.xelem (n2) = (x1 + x2) / 2;
322  }
323 
324  return retval;
325 }
326 
327 // row vector by column vector -> scalar
328 
329 double
330 operator * (const RowVector& v, const ColumnVector& a)
331 {
332  double retval = 0.0;
333 
334  F77_INT len = octave::to_f77_int (v.numel ());
335 
336  F77_INT a_len = octave::to_f77_int (a.numel ());
337 
338  if (len != a_len)
339  octave::err_nonconformant ("operator *", len, a_len);
340 
341  if (len != 0)
342  F77_FUNC (xddot, XDDOT) (len, v.data (), 1, a.data (), 1, retval);
343 
344  return retval;
345 }
346 
347 Complex
349 {
350  ComplexRowVector tmp (v);
351  return tmp * a;
352 }
#define NaN
Definition: Faddeeva.cc:261
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
octave_idx_type rows() const
Definition: Array.h:459
const double * data() const
Definition: Array.h:663
octave_idx_type cols() const
Definition: Array.h:469
void make_unique()
Definition: Array.h:216
double & xelem(octave_idx_type n)
Definition: Array.h:524
octave_idx_type numel() const
Number of elements in the array.
Definition: Array.h:414
MArray< T > transpose() const
Definition: MArray.h:99
Definition: dMatrix.h:42
bool operator==(const RowVector &a) const
Definition: dRowVector.cc:44
double min() const
Definition: dRowVector.cc:207
ColumnVector transpose() const
Definition: dRowVector.cc:126
bool operator!=(const RowVector &a) const
Definition: dRowVector.cc:53
double max() const
Definition: dRowVector.cc:223
RowVector & insert(const RowVector &a, octave_idx_type c)
Definition: dRowVector.cc:59
RowVector extract_n(octave_idx_type c1, octave_idx_type n) const
Definition: dRowVector.cc:159
void clear(octave_idx_type n)
Definition: dRowVector.h:107
RowVector extract(octave_idx_type c1, octave_idx_type c2) const
Definition: dRowVector.cc:144
RowVector & fill(double val)
Definition: dRowVector.cc:78
RowVector append(const RowVector &a) const
Definition: dRowVector.cc:115
void resize(octave_idx_type n, const double &rfv=0)
Definition: dRowVector.h:102
std::istream & operator>>(std::istream &is, RowVector &a)
Definition: dRowVector.cc:249
RowVector imag(const ComplexRowVector &a)
Definition: dRowVector.cc:138
std::ostream & operator<<(std::ostream &os, const RowVector &a)
Definition: dRowVector.cc:239
RowVector linspace(double x1, double x2, octave_idx_type n_in)
Definition: dRowVector.cc:271
RowVector real(const ComplexRowVector &a)
Definition: dRowVector.cc:132
RowVector operator*(const RowVector &v, const Matrix &a)
Definition: dRowVector.cc:172
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
bool isinf(double x)
Definition: lo-mappers.h:203
bool mx_inline_equal(std::size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:577
void mx_inline_real(std::size_t n, T *r, const std::complex< T > *x)
Definition: mx-inlines.cc:325
void mx_inline_imag(std::size_t n, T *r, const std::complex< T > *x)
Definition: mx-inlines.cc:333
octave_idx_type n
Definition: mx-inlines.cc:761
std::complex< double > Complex
Definition: oct-cmplx.h:33
subroutine xddot(n, dx, incx, dy, incy, retval)
Definition: xddot.f:2
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg
F77_RET_T len
Definition: xerbla.cc:61