GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CRowVector.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 // Complex 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 // destructive insert/delete/reorder operations
59 
62 {
63  octave_idx_type a_len = a.numel ();
64 
65  if (c < 0 || c + a_len > numel ())
66  (*current_liboctave_error_handler) ("range error for insert");
67 
68  if (a_len > 0)
69  {
70  make_unique ();
71 
72  for (octave_idx_type i = 0; i < a_len; i++)
73  xelem (c+i) = a.elem (i);
74  }
75 
76  return *this;
77 }
78 
81 {
82  octave_idx_type a_len = a.numel ();
83 
84  if (c < 0 || c + a_len > numel ())
85  (*current_liboctave_error_handler) ("range error for insert");
86 
87  if (a_len > 0)
88  {
89  make_unique ();
90 
91  for (octave_idx_type i = 0; i < a_len; i++)
92  xelem (c+i) = a.elem (i);
93  }
94 
95  return *this;
96 }
97 
100 {
102 
103  if (len > 0)
104  {
105  make_unique ();
106 
107  for (octave_idx_type i = 0; i < len; i++)
108  xelem (i) = val;
109  }
110 
111  return *this;
112 }
113 
116 {
118 
119  if (len > 0)
120  {
121  make_unique ();
122 
123  for (octave_idx_type i = 0; i < len; i++)
124  xelem (i) = val;
125  }
126 
127  return *this;
128 }
129 
132 {
134 
135  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
136  (*current_liboctave_error_handler) ("range error for fill");
137 
138  if (c1 > c2) { std::swap (c1, c2); }
139 
140  if (c2 >= c1)
141  {
142  make_unique ();
143 
144  for (octave_idx_type i = c1; i <= c2; i++)
145  xelem (i) = val;
146  }
147 
148  return *this;
149 }
150 
154 {
156 
157  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
158  (*current_liboctave_error_handler) ("range error for fill");
159 
160  if (c1 > c2) { std::swap (c1, c2); }
161 
162  if (c2 >= c1)
163  {
164  make_unique ();
165 
166  for (octave_idx_type i = c1; i <= c2; i++)
167  xelem (i) = val;
168  }
169 
170  return *this;
171 }
172 
175 {
177  octave_idx_type nc_insert = len;
178  ComplexRowVector retval (len + a.numel ());
179  retval.insert (*this, 0);
180  retval.insert (a, nc_insert);
181  return retval;
182 }
183 
186 {
188  octave_idx_type nc_insert = len;
189  ComplexRowVector retval (len + a.numel ());
190  retval.insert (*this, 0);
191  retval.insert (a, nc_insert);
192  return retval;
193 }
194 
197 {
199 }
200 
203 {
204  return MArray<Complex>::transpose ();
205 }
206 
209 {
210  return do_mx_unary_map<Complex, Complex, std::conj<double>> (a);
211 }
212 
213 // resize is the destructive equivalent for this one
214 
217 {
218  if (c1 > c2) { std::swap (c1, c2); }
219 
220  octave_idx_type new_c = c2 - c1 + 1;
221 
222  ComplexRowVector result (new_c);
223 
224  for (octave_idx_type i = 0; i < new_c; i++)
225  result.elem (i) = elem (c1+i);
226 
227  return result;
228 }
229 
232 {
233  ComplexRowVector result (n);
234 
235  for (octave_idx_type i = 0; i < n; i++)
236  result.elem (i) = elem (r1+i);
237 
238  return result;
239 }
240 
241 // row vector by row vector -> row vector operations
242 
245 {
247 
248  octave_idx_type a_len = a.numel ();
249 
250  if (len != a_len)
251  octave::err_nonconformant ("operator +=", len, a_len);
252 
253  if (len == 0)
254  return *this;
255 
256  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
257 
258  mx_inline_add2 (len, d, a.data ());
259  return *this;
260 }
261 
264 {
266 
267  octave_idx_type a_len = a.numel ();
268 
269  if (len != a_len)
270  octave::err_nonconformant ("operator -=", len, a_len);
271 
272  if (len == 0)
273  return *this;
274 
275  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
276 
277  mx_inline_sub2 (len, d, a.data ());
278  return *this;
279 }
280 
281 // row vector by matrix -> row vector
282 
285 {
286  ComplexRowVector retval;
287 
288  F77_INT len = octave::to_f77_int (v.numel ());
289 
290  F77_INT a_nr = octave::to_f77_int (a.rows ());
291  F77_INT a_nc = octave::to_f77_int (a.cols ());
292 
293  if (a_nr != len)
294  octave::err_nonconformant ("operator *", 1, len, a_nr, a_nc);
295 
296  if (len == 0)
297  retval.resize (a_nc, 0.0);
298  else
299  {
300  // Transpose A to form A'*x == (x'*A)'
301 
302  F77_INT ld = a_nr;
303 
304  retval.resize (a_nc);
305  Complex *y = retval.fortran_vec ();
306 
307  F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
308  a_nr, a_nc, 1.0, F77_CONST_DBLE_CMPLX_ARG (a.data ()),
309  ld, F77_CONST_DBLE_CMPLX_ARG (v.data ()), 1, 0.0, F77_DBLE_CMPLX_ARG (y), 1
310  F77_CHAR_ARG_LEN (1)));
311  }
312 
313  return retval;
314 }
315 
318 {
319  ComplexRowVector tmp (v);
320  return tmp * a;
321 }
322 
323 // other operations
324 
325 Complex
327 {
329  if (len == 0)
330  return Complex (0.0);
331 
332  Complex res = elem (0);
333  double absres = std::abs (res);
334 
335  for (octave_idx_type i = 1; i < len; i++)
336  if (std::abs (elem (i)) < absres)
337  {
338  res = elem (i);
339  absres = std::abs (res);
340  }
341 
342  return res;
343 }
344 
345 Complex
347 {
349  if (len == 0)
350  return Complex (0.0);
351 
352  Complex res = elem (0);
353  double absres = std::abs (res);
354 
355  for (octave_idx_type i = 1; i < len; i++)
356  if (std::abs (elem (i)) > absres)
357  {
358  res = elem (i);
359  absres = std::abs (res);
360  }
361 
362  return res;
363 }
364 
365 // i/o
366 
367 std::ostream&
368 operator << (std::ostream& os, const ComplexRowVector& a)
369 {
370 // int field_width = os.precision () + 7;
371  for (octave_idx_type i = 0; i < a.numel (); i++)
372  os << ' ' /* setw (field_width) */ << a.elem (i);
373  return os;
374 }
375 
376 std::istream&
377 operator >> (std::istream& is, ComplexRowVector& a)
378 {
379  octave_idx_type len = a.numel ();
380 
381  if (len > 0)
382  {
383  Complex tmp;
384  for (octave_idx_type i = 0; i < len; i++)
385  {
386  is >> tmp;
387  if (is)
388  a.elem (i) = tmp;
389  else
390  break;
391  }
392  }
393  return is;
394 }
395 
396 // row vector by column vector -> scalar
397 
398 // row vector by column vector -> scalar
399 
400 Complex
402 {
403  ComplexColumnVector tmp (a);
404  return v * tmp;
405 }
406 
407 Complex
409 {
410  Complex retval (0.0, 0.0);
411 
412  F77_INT len = octave::to_f77_int (v.numel ());
413 
414  F77_INT a_len = octave::to_f77_int (a.numel ());
415 
416  if (len != a_len)
417  octave::err_nonconformant ("operator *", len, a_len);
418  if (len != 0)
419  F77_FUNC (xzdotu, XZDOTU) (len, F77_CONST_DBLE_CMPLX_ARG (v.data ()), 1,
420  F77_CONST_DBLE_CMPLX_ARG (a.data ()), 1, F77_DBLE_CMPLX_ARG (&retval));
421 
422  return retval;
423 }
424 
425 // other operations
426 
428 linspace (const Complex& x1, const Complex& x2, octave_idx_type n_in)
429 {
430  ComplexRowVector retval;
431 
432  if (n_in < 1)
433  return retval;
434  else if (n_in == 1)
435  {
436  retval.resize (1, x2);
437  return retval;
438  }
439  else if (x1 == x2)
440  {
441  retval.resize (n_in, x2);
442  return retval;
443  }
444 
445  // Use unsigned type (guaranteed n_in > 1 at this point) so that divisions
446  // by 2 can be replaced by compiler with shift right instructions.
447  typedef std::make_unsigned<octave_idx_type>::type unsigned_octave_idx_type;
448 
449  unsigned_octave_idx_type n = n_in;
450 
451  // Set endpoints, rather than calculate, for maximum accuracy.
452  retval.clear (n);
453  retval.xelem (0) = x1;
454  retval.xelem (n-1) = x2;
455 
456  // Construct linspace symmetrically from both ends.
457  bool isnan_delta = false;
458  Complex delta = (x2 - x1) / (n - 1.0);
459  if (octave::math::isinf (delta))
460  {
461  if (octave::math::isinf (delta.real ()))
462  delta.real (octave::numeric_limits<double>::NaN ());
463  if (octave::math::isinf (delta.imag ()))
464  delta.imag (octave::numeric_limits<double>::NaN ());
465  isnan_delta = true;
466  }
467 
468  unsigned_octave_idx_type n2 = n/2;
469  for (unsigned_octave_idx_type i = 1; i < n2; i++)
470  {
471  retval.xelem (i) = x1 + static_cast<double> (i)*delta;
472  retval.xelem (n-1-i) = x2 - static_cast<double> (i)*delta;
473  }
474  if (n % 2 == 1) // Middle element if number of elements is odd.
475  {
476  if (x1 == -x2)
477  retval.xelem (n2) = 0;
478  else
479  {
480  Complex c = (x1 + x2) / 2.0;
481  if (isnan_delta)
482  {
483  if (octave::math::isnan (delta.real ()))
485  if (octave::math::isnan (delta.imag ()))
487  }
488  retval.xelem (n2) = c;
489  }
490  }
491 
492  return retval;
493 }
std::istream & operator>>(std::istream &is, ComplexRowVector &a)
Definition: CRowVector.cc:377
ComplexRowVector conj(const ComplexRowVector &a)
Definition: CRowVector.cc:208
ComplexRowVector operator*(const ComplexRowVector &v, const ComplexMatrix &a)
Definition: CRowVector.cc:284
std::ostream & operator<<(std::ostream &os, const ComplexRowVector &a)
Definition: CRowVector.cc:368
ComplexRowVector linspace(const Complex &x1, const Complex &x2, octave_idx_type n_in)
Definition: CRowVector.cc:428
#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 T * data() const
Size of the specified dimension.
Definition: Array.h:663
octave_idx_type cols() const
Definition: Array.h:469
void make_unique()
Definition: Array.h:216
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:524
octave_idx_type numel() const
Number of elements in the array.
Definition: Array.h:414
ComplexRowVector & operator-=(const RowVector &a)
Definition: CRowVector.cc:263
ComplexRowVector & operator+=(const RowVector &a)
Definition: CRowVector.cc:244
ComplexRowVector extract(octave_idx_type c1, octave_idx_type c2) const
Definition: CRowVector.cc:216
Complex max() const
Definition: CRowVector.cc:346
ComplexRowVector append(const RowVector &a) const
Definition: CRowVector.cc:174
void resize(octave_idx_type n, const Complex &rfv=Complex(0))
Definition: CRowVector.h:129
void clear(octave_idx_type n)
Definition: CRowVector.h:134
ComplexColumnVector hermitian() const
Definition: CRowVector.cc:196
bool operator==(const ComplexRowVector &a) const
Definition: CRowVector.cc:44
ComplexRowVector extract_n(octave_idx_type c1, octave_idx_type n) const
Definition: CRowVector.cc:231
ComplexRowVector & fill(double val)
Definition: CRowVector.cc:99
ComplexRowVector & insert(const RowVector &a, octave_idx_type c)
Definition: CRowVector.cc:61
bool operator!=(const ComplexRowVector &a) const
Definition: CRowVector.cc:53
Complex min() const
Definition: CRowVector.cc:326
ComplexColumnVector transpose() const
Definition: CRowVector.cc:202
MArray< T > hermitian(T(*fcn)(const T &)=nullptr) const
Definition: MArray.h:102
MArray< T > transpose() const
Definition: MArray.h:99
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:319
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 isnan(bool)
Definition: lo-mappers.h:178
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
void mx_inline_sub2(std::size_t n, R *r, const X *x)
Definition: mx-inlines.cc:128
void mx_inline_add2(std::size_t n, R *r, const X *x)
Definition: mx-inlines.cc:127
bool mx_inline_equal(std::size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:577
octave_idx_type n
Definition: mx-inlines.cc:761
std::complex< double > Complex
Definition: oct-cmplx.h:33
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg
F77_RET_T len
Definition: xerbla.cc:61
subroutine xzdotu(n, zx, incx, zy, incy, retval)
Definition: xzdotu.f:2