GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
fCColVector.cc
Go to the documentation of this file.
1 // ColumnVector manipulations.
2 /*
3 
4 Copyright (C) 1994-2015 John W. Eaton
5 Copyright (C) 2010 VZLU Prague
6 
7 This file is part of Octave.
8 
9 Octave is free software; you can redistribute it and/or modify it
10 under the terms of the GNU General Public License as published by the
11 Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13 
14 Octave is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18 
19 You should have received a copy of the GNU General Public License
20 along with Octave; see the file COPYING. If not, see
21 <http://www.gnu.org/licenses/>.
22 
23 */
24 
25 #ifdef HAVE_CONFIG_H
26 #include <config.h>
27 #endif
28 
29 #include <iostream>
30 
31 #include "Array-util.h"
32 #include "f77-fcn.h"
33 #include "functor.h"
34 #include "lo-error.h"
35 #include "mx-base.h"
36 #include "mx-inlines.cc"
37 #include "oct-cmplx.h"
38 
39 // Fortran functions we call.
40 
41 extern "C"
42 {
43  F77_RET_T
44  F77_FUNC (cgemv, CGEMV) (F77_CONST_CHAR_ARG_DECL,
46  const FloatComplex&, const FloatComplex*,
47  const octave_idx_type&, const FloatComplex*,
48  const octave_idx_type&, const FloatComplex&,
49  FloatComplex*, const octave_idx_type&
51 }
52 
53 // FloatComplex Column Vector class
54 
56  : MArray<FloatComplex> (a)
57 {
58 }
59 
60 bool
62 {
63  octave_idx_type len = length ();
64  if (len != a.length ())
65  return 0;
66  return mx_inline_equal (len, data (), a.data ());
67 }
68 
69 bool
71 {
72  return !(*this == a);
73 }
74 
75 // destructive insert/delete/reorder operations
76 
79 {
80  octave_idx_type a_len = a.length ();
81 
82  if (r < 0 || r + a_len > length ())
83  {
84  (*current_liboctave_error_handler) ("range error for insert");
85  return *this;
86  }
87 
88  if (a_len > 0)
89  {
90  make_unique ();
91 
92  for (octave_idx_type i = 0; i < a_len; i++)
93  xelem (r+i) = a.elem (i);
94  }
95 
96  return *this;
97 }
98 
101  octave_idx_type r)
102 {
103  octave_idx_type a_len = a.length ();
104 
105  if (r < 0 || r + a_len > length ())
106  {
107  (*current_liboctave_error_handler) ("range error for insert");
108  return *this;
109  }
110 
111  if (a_len > 0)
112  {
113  make_unique ();
114 
115  for (octave_idx_type i = 0; i < a_len; i++)
116  xelem (r+i) = a.elem (i);
117  }
118 
119  return *this;
120 }
121 
124 {
125  octave_idx_type len = length ();
126 
127  if (len > 0)
128  {
129  make_unique ();
130 
131  for (octave_idx_type i = 0; i < len; i++)
132  xelem (i) = val;
133  }
134 
135  return *this;
136 }
137 
140 {
141  octave_idx_type len = length ();
142 
143  if (len > 0)
144  {
145  make_unique ();
146 
147  for (octave_idx_type i = 0; i < len; i++)
148  xelem (i) = val;
149  }
150 
151 
152  return *this;
153 }
154 
158 {
159  octave_idx_type len = length ();
160 
161  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
162  {
163  (*current_liboctave_error_handler) ("range error for fill");
164  return *this;
165  }
166 
167  if (r1 > r2) { std::swap (r1, r2); }
168 
169  if (r2 >= r1)
170  {
171  make_unique ();
172 
173  for (octave_idx_type i = r1; i <= r2; i++)
174  xelem (i) = val;
175  }
176 
177  return *this;
178 }
179 
183 {
184  octave_idx_type len = length ();
185 
186  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
187  {
188  (*current_liboctave_error_handler) ("range error for fill");
189  return *this;
190  }
191 
192  if (r1 > r2) { std::swap (r1, r2); }
193 
194  if (r2 >= r1)
195  {
196  make_unique ();
197 
198  for (octave_idx_type i = r1; i <= r2; i++)
199  xelem (i) = val;
200  }
201 
202  return *this;
203 }
204 
207 {
208  octave_idx_type len = length ();
209  octave_idx_type nr_insert = len;
210  FloatComplexColumnVector retval (len + a.length ());
211  retval.insert (*this, 0);
212  retval.insert (a, nr_insert);
213  return retval;
214 }
215 
218 {
219  octave_idx_type len = length ();
220  octave_idx_type nr_insert = len;
221  FloatComplexColumnVector retval (len + a.length ());
222  retval.insert (*this, 0);
223  retval.insert (a, nr_insert);
224  return retval;
225 }
226 
229 {
231 }
232 
235 {
237 }
238 
241 {
242  return do_mx_unary_map<float, FloatComplex, std::abs> (*this);
243 }
244 
247 {
248  return do_mx_unary_map<FloatComplex, FloatComplex, std::conj<float> > (a);
249 }
250 
251 // resize is the destructive equivalent for this one
252 
255 {
256  if (r1 > r2) { std::swap (r1, r2); }
257 
258  octave_idx_type new_r = r2 - r1 + 1;
259 
260  FloatComplexColumnVector result (new_r);
261 
262  for (octave_idx_type i = 0; i < new_r; i++)
263  result.elem (i) = elem (r1+i);
264 
265  return result;
266 }
267 
270  octave_idx_type n) const
271 {
272  FloatComplexColumnVector result (n);
273 
274  for (octave_idx_type i = 0; i < n; i++)
275  result.elem (i) = elem (r1+i);
276 
277  return result;
278 }
279 
280 // column vector by column vector -> column vector operations
281 
284 {
285  octave_idx_type len = length ();
286 
287  octave_idx_type a_len = a.length ();
288 
289  if (len != a_len)
290  {
291  gripe_nonconformant ("operator +=", len, a_len);
292  return *this;
293  }
294 
295  if (len == 0)
296  return *this;
297 
298  FloatComplex *d = fortran_vec (); // Ensures only 1 reference to my privates!
299 
300  mx_inline_add2 (len, d, a.data ());
301  return *this;
302 }
303 
306 {
307  octave_idx_type len = length ();
308 
309  octave_idx_type a_len = a.length ();
310 
311  if (len != a_len)
312  {
313  gripe_nonconformant ("operator -=", len, a_len);
314  return *this;
315  }
316 
317  if (len == 0)
318  return *this;
319 
320  FloatComplex *d = fortran_vec (); // Ensures only 1 reference to my privates!
321 
322  mx_inline_sub2 (len, d, a.data ());
323  return *this;
324 }
325 
326 // matrix by column vector -> column vector operations
327 
330 {
331  FloatComplexColumnVector tmp (a);
332  return m * tmp;
333 }
334 
337 {
339 
340  octave_idx_type nr = m.rows ();
341  octave_idx_type nc = m.cols ();
342 
343  octave_idx_type a_len = a.length ();
344 
345  if (nc != a_len)
346  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
347  else
348  {
349  retval.clear (nr);
350 
351  if (nr != 0)
352  {
353  if (nc == 0)
354  retval.fill (0.0);
355  else
356  {
357  FloatComplex *y = retval.fortran_vec ();
358 
359  F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
360  nr, nc, 1.0f, m.data (), nr,
361  a.data (), 1, 0.0f, y, 1
362  F77_CHAR_ARG_LEN (1)));
363  }
364  }
365  }
366 
367  return retval;
368 }
369 
370 // matrix by column vector -> column vector operations
371 
374 {
375  FloatComplexMatrix tmp (m);
376  return tmp * a;
377 }
378 
379 // diagonal matrix by column vector -> column vector operations
380 
383 {
384  octave_idx_type nr = m.rows ();
385  octave_idx_type nc = m.cols ();
386 
387  octave_idx_type a_len = a.length ();
388 
389  if (nc != a_len)
390  {
391  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
392  return FloatComplexColumnVector ();
393  }
394 
395  if (nc == 0 || nr == 0)
396  return FloatComplexColumnVector (0);
397 
398  FloatComplexColumnVector result (nr);
399 
400  for (octave_idx_type i = 0; i < a_len; i++)
401  result.elem (i) = a.elem (i) * m.elem (i, i);
402 
403  for (octave_idx_type i = a_len; i < nr; i++)
404  result.elem (i) = 0.0;
405 
406  return result;
407 }
408 
411 {
412  octave_idx_type nr = m.rows ();
413  octave_idx_type nc = m.cols ();
414 
415  octave_idx_type a_len = a.length ();
416 
417  if (nc != a_len)
418  {
419  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
420  return FloatComplexColumnVector ();
421  }
422 
423  if (nc == 0 || nr == 0)
424  return FloatComplexColumnVector (0);
425 
426  FloatComplexColumnVector result (nr);
427 
428  for (octave_idx_type i = 0; i < a_len; i++)
429  result.elem (i) = a.elem (i) * m.elem (i, i);
430 
431  for (octave_idx_type i = a_len; i < nr; i++)
432  result.elem (i) = 0.0;
433 
434  return result;
435 }
436 
439 {
440  octave_idx_type nr = m.rows ();
441  octave_idx_type nc = m.cols ();
442 
443  octave_idx_type a_len = a.length ();
444 
445  if (nc != a_len)
446  {
447  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
448  return FloatComplexColumnVector ();
449  }
450 
451  if (nc == 0 || nr == 0)
452  return FloatComplexColumnVector (0);
453 
454  FloatComplexColumnVector result (nr);
455 
456  for (octave_idx_type i = 0; i < a_len; i++)
457  result.elem (i) = a.elem (i) * m.elem (i, i);
458 
459  for (octave_idx_type i = a_len; i < nr; i++)
460  result.elem (i) = 0.0;
461 
462  return result;
463 }
464 
465 // other operations
466 
469 {
470  octave_idx_type len = length ();
471  if (len == 0)
472  return 0.0;
473 
474  FloatComplex res = elem (0);
475  float absres = std::abs (res);
476 
477  for (octave_idx_type i = 1; i < len; i++)
478  if (std::abs (elem (i)) < absres)
479  {
480  res = elem (i);
481  absres = std::abs (res);
482  }
483 
484  return res;
485 }
486 
489 {
490  octave_idx_type len = length ();
491  if (len == 0)
492  return 0.0;
493 
494  FloatComplex res = elem (0);
495  float absres = std::abs (res);
496 
497  for (octave_idx_type i = 1; i < len; i++)
498  if (std::abs (elem (i)) > absres)
499  {
500  res = elem (i);
501  absres = std::abs (res);
502  }
503 
504  return res;
505 }
506 
507 // i/o
508 
509 std::ostream&
510 operator << (std::ostream& os, const FloatComplexColumnVector& a)
511 {
512 // int field_width = os.precision () + 7;
513  for (octave_idx_type i = 0; i < a.length (); i++)
514  os << /* setw (field_width) << */ a.elem (i) << "\n";
515  return os;
516 }
517 
518 std::istream&
520 {
521  octave_idx_type len = a.length ();
522 
523  if (len > 0)
524  {
525  float tmp;
526  for (octave_idx_type i = 0; i < len; i++)
527  {
528  is >> tmp;
529  if (is)
530  a.elem (i) = tmp;
531  else
532  break;
533  }
534  }
535  return is;
536 }
void mx_inline_add2(size_t n, R *r, const X *x)
Definition: mx-inlines.cc:95
#define F77_CHAR_ARG_LEN(l)
Definition: f77-fcn.h:253
FloatComplexColumnVector operator*(const FloatComplexMatrix &m, const FloatColumnVector &a)
Definition: fCColVector.cc:329
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
void mx_inline_sub2(size_t n, R *r, const X *x)
Definition: mx-inlines.cc:96
FloatComplexRowVector transpose(void) const
Definition: fCColVector.cc:234
T elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:110
bool operator!=(const FloatComplexColumnVector &a) const
Definition: fCColVector.cc:70
FloatComplexColumnVector & operator-=(const FloatColumnVector &a)
Definition: fCColVector.cc:305
octave_idx_type rows(void) const
Definition: DiagArray2.h:86
MArray< T > transpose(void) const
Definition: MArray.h:83
FloatColumnVector abs(void) const
Definition: fCColVector.cc:240
std::istream & operator>>(std::istream &is, FloatComplexColumnVector &a)
Definition: fCColVector.cc:519
FloatComplex min(void) const
Definition: fCColVector.cc:468
T & elem(octave_idx_type n)
Definition: Array.h:380
Definition: MArray.h:36
FloatComplexColumnVector conj(const FloatComplexColumnVector &a)
Definition: fCColVector.cc:246
bool operator==(const FloatComplexColumnVector &a) const
Definition: fCColVector.cc:61
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
FloatComplex max(void) const
Definition: fCColVector.cc:488
octave_idx_type rows(void) const
Definition: Array.h:313
F77_RET_T const double const double double * d
#define F77_CONST_CHAR_ARG2(x, l)
Definition: f77-fcn.h:251
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type r2
FloatComplexColumnVector & operator+=(const FloatColumnVector &a)
Definition: fCColVector.cc:283
F77_RET_T const double const double * f
void make_unique(void)
Definition: Array.h:104
const FloatComplex * data(void) const
Definition: Array.h:479
#define F77_RET_T
Definition: f77-fcn.h:264
MArray< T > hermitian(T(*fcn)(const T &)=0) const
Definition: MArray.h:86
FloatComplexColumnVector stack(const FloatColumnVector &a) const
Definition: fCColVector.cc:206
FloatComplexRowVector hermitian(void) const
Definition: fCColVector.cc:228
FloatComplex & xelem(octave_idx_type n)
Definition: Array.h:353
F77_RET_T F77_FUNC(cgemv, CGEMV)(F77_CONST_CHAR_ARG_DECL
void clear(octave_idx_type n)
Definition: fCColVector.h:151
octave_idx_type cols(void) const
Definition: DiagArray2.h:87
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
#define F77_CONST_CHAR_ARG_DECL
Definition: f77-fcn.h:255
FloatComplexColumnVector extract_n(octave_idx_type r1, octave_idx_type n) const
Definition: fCColVector.cc:269
F77_RET_T const octave_idx_type const octave_idx_type const FloatComplex const FloatComplex const octave_idx_type const FloatComplex const octave_idx_type const FloatComplex FloatComplex const octave_idx_type & F77_CHAR_ARG_LEN_DECL
Definition: fCColVector.cc:45
FloatComplexColumnVector & fill(float val)
Definition: fCColVector.cc:123
std::ostream & operator<<(std::ostream &os, const FloatComplexColumnVector &a)
Definition: fCColVector.cc:510
FloatComplexColumnVector & insert(const FloatColumnVector &a, octave_idx_type r)
Definition: fCColVector.cc:78
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type r1
const FloatComplex * fortran_vec(void) const
Definition: Array.h:481
FloatComplexColumnVector extract(octave_idx_type r1, octave_idx_type r2) const
Definition: fCColVector.cc:254
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:441
octave_idx_type cols(void) const
Definition: Array.h:321
T abs(T x)
Definition: pr-output.cc:3062