GNU Octave  3.8.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
CColVector.cc
Go to the documentation of this file.
1 // ColumnVector manipulations.
2 /*
3 
4 Copyright (C) 1994-2013 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 (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
46  const Complex&, const Complex*,
47  const octave_idx_type&, const Complex*,
48  const octave_idx_type&, const Complex&,
49  Complex*, const octave_idx_type&
51 }
52 
53 // Complex Column Vector class
54 
56  : MArray<Complex> (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 {
102  octave_idx_type a_len = a.length ();
103 
104  if (r < 0 || r + a_len > length ())
105  {
106  (*current_liboctave_error_handler) ("range error for insert");
107  return *this;
108  }
109 
110  if (a_len > 0)
111  {
112  make_unique ();
113 
114  for (octave_idx_type i = 0; i < a_len; i++)
115  xelem (r+i) = a.elem (i);
116  }
117 
118  return *this;
119 }
120 
123 {
124  octave_idx_type len = length ();
125 
126  if (len > 0)
127  {
128  make_unique ();
129 
130  for (octave_idx_type i = 0; i < len; i++)
131  xelem (i) = val;
132  }
133 
134  return *this;
135 }
136 
139 {
140  octave_idx_type len = length ();
141 
142  if (len > 0)
143  {
144  make_unique ();
145 
146  for (octave_idx_type i = 0; i < len; i++)
147  xelem (i) = val;
148  }
149 
150 
151  return *this;
152 }
153 
156 {
157  octave_idx_type len = length ();
158 
159  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
160  {
161  (*current_liboctave_error_handler) ("range error for fill");
162  return *this;
163  }
164 
165  if (r1 > r2) { std::swap (r1, r2); }
166 
167  if (r2 >= r1)
168  {
169  make_unique ();
170 
171  for (octave_idx_type i = r1; i <= r2; i++)
172  xelem (i) = val;
173  }
174 
175  return *this;
176 }
177 
181 {
182  octave_idx_type len = length ();
183 
184  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
185  {
186  (*current_liboctave_error_handler) ("range error for fill");
187  return *this;
188  }
189 
190  if (r1 > r2) { std::swap (r1, r2); }
191 
192  if (r2 >= r1)
193  {
194  make_unique ();
195 
196  for (octave_idx_type i = r1; i <= r2; i++)
197  xelem (i) = val;
198  }
199 
200  return *this;
201 }
202 
205 {
206  octave_idx_type len = length ();
207  octave_idx_type nr_insert = len;
208  ComplexColumnVector retval (len + a.length ());
209  retval.insert (*this, 0);
210  retval.insert (a, nr_insert);
211  return retval;
212 }
213 
216 {
217  octave_idx_type len = length ();
218  octave_idx_type nr_insert = len;
219  ComplexColumnVector retval (len + a.length ());
220  retval.insert (*this, 0);
221  retval.insert (a, nr_insert);
222  return retval;
223 }
224 
227 {
229 }
230 
233 {
234  return MArray<Complex>::transpose ();
235 }
236 
239 {
240  return do_mx_unary_map<double, Complex, std::abs> (*this);
241 }
242 
245 {
246  return do_mx_unary_map<Complex, Complex, std::conj<double> > (a);
247 }
248 
249 // resize is the destructive equivalent for this one
250 
253 {
254  if (r1 > r2) { std::swap (r1, r2); }
255 
256  octave_idx_type new_r = r2 - r1 + 1;
257 
258  ComplexColumnVector result (new_r);
259 
260  for (octave_idx_type i = 0; i < new_r; i++)
261  result.elem (i) = elem (r1+i);
262 
263  return result;
264 }
265 
268 {
269  ComplexColumnVector result (n);
270 
271  for (octave_idx_type i = 0; i < n; i++)
272  result.elem (i) = elem (r1+i);
273 
274  return result;
275 }
276 
277 // column vector by column vector -> column vector operations
278 
281 {
282  octave_idx_type len = length ();
283 
284  octave_idx_type a_len = a.length ();
285 
286  if (len != a_len)
287  {
288  gripe_nonconformant ("operator +=", len, a_len);
289  return *this;
290  }
291 
292  if (len == 0)
293  return *this;
294 
295  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
296 
297  mx_inline_add2 (len, d, a.data ());
298  return *this;
299 }
300 
303 {
304  octave_idx_type len = length ();
305 
306  octave_idx_type a_len = a.length ();
307 
308  if (len != a_len)
309  {
310  gripe_nonconformant ("operator -=", len, a_len);
311  return *this;
312  }
313 
314  if (len == 0)
315  return *this;
316 
317  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
318 
319  mx_inline_sub2 (len, d, a.data ());
320  return *this;
321 }
322 
323 // matrix by column vector -> column vector operations
324 
327 {
328  ComplexColumnVector tmp (a);
329  return m * tmp;
330 }
331 
334 {
335  ComplexColumnVector retval;
336 
337  octave_idx_type nr = m.rows ();
338  octave_idx_type nc = m.cols ();
339 
340  octave_idx_type a_len = a.length ();
341 
342  if (nc != a_len)
343  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
344  else
345  {
346  retval.clear (nr);
347 
348  if (nr != 0)
349  {
350  if (nc == 0)
351  retval.fill (0.0);
352  else
353  {
354  Complex *y = retval.fortran_vec ();
355 
356  F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
357  nr, nc, 1.0, m.data (), nr,
358  a.data (), 1, 0.0, y, 1
359  F77_CHAR_ARG_LEN (1)));
360  }
361  }
362 
363  }
364 
365  return retval;
366 }
367 
368 // matrix by column vector -> column vector operations
369 
372 {
373  ComplexMatrix tmp (m);
374  return tmp * a;
375 }
376 
377 // diagonal matrix by column vector -> column vector operations
378 
381 {
382  octave_idx_type nr = m.rows ();
383  octave_idx_type nc = m.cols ();
384 
385  octave_idx_type a_len = a.length ();
386 
387  if (nc != a_len)
388  {
389  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
390  return ComplexColumnVector ();
391  }
392 
393  if (nc == 0 || nr == 0)
394  return ComplexColumnVector (0);
395 
396  ComplexColumnVector result (nr);
397 
398  for (octave_idx_type i = 0; i < a_len; i++)
399  result.elem (i) = a.elem (i) * m.elem (i, i);
400 
401  for (octave_idx_type i = a_len; i < nr; i++)
402  result.elem (i) = 0.0;
403 
404  return result;
405 }
406 
409 {
410  octave_idx_type nr = m.rows ();
411  octave_idx_type nc = m.cols ();
412 
413  octave_idx_type a_len = a.length ();
414 
415  if (nc != a_len)
416  {
417  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
418  return ComplexColumnVector ();
419  }
420 
421  if (nc == 0 || nr == 0)
422  return ComplexColumnVector (0);
423 
424  ComplexColumnVector result (nr);
425 
426  for (octave_idx_type i = 0; i < a_len; i++)
427  result.elem (i) = a.elem (i) * m.elem (i, i);
428 
429  for (octave_idx_type i = a_len; i < nr; i++)
430  result.elem (i) = 0.0;
431 
432  return result;
433 }
434 
437 {
438  octave_idx_type nr = m.rows ();
439  octave_idx_type nc = m.cols ();
440 
441  octave_idx_type a_len = a.length ();
442 
443  if (nc != a_len)
444  {
445  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
446  return ComplexColumnVector ();
447  }
448 
449  if (nc == 0 || nr == 0)
450  return ComplexColumnVector (0);
451 
452  ComplexColumnVector result (nr);
453 
454  for (octave_idx_type i = 0; i < a_len; i++)
455  result.elem (i) = a.elem (i) * m.elem (i, i);
456 
457  for (octave_idx_type i = a_len; i < nr; i++)
458  result.elem (i) = 0.0;
459 
460  return result;
461 }
462 
463 // other operations
464 
465 Complex
467 {
468  octave_idx_type len = length ();
469  if (len == 0)
470  return 0.0;
471 
472  Complex res = elem (0);
473  double absres = std::abs (res);
474 
475  for (octave_idx_type i = 1; i < len; i++)
476  if (std::abs (elem (i)) < absres)
477  {
478  res = elem (i);
479  absres = std::abs (res);
480  }
481 
482  return res;
483 }
484 
485 Complex
487 {
488  octave_idx_type len = length ();
489  if (len == 0)
490  return 0.0;
491 
492  Complex res = elem (0);
493  double absres = std::abs (res);
494 
495  for (octave_idx_type i = 1; i < len; i++)
496  if (std::abs (elem (i)) > absres)
497  {
498  res = elem (i);
499  absres = std::abs (res);
500  }
501 
502  return res;
503 }
504 
505 // i/o
506 
507 std::ostream&
508 operator << (std::ostream& os, const ComplexColumnVector& a)
509 {
510 // int field_width = os.precision () + 7;
511  for (octave_idx_type i = 0; i < a.length (); i++)
512  os << /* setw (field_width) << */ a.elem (i) << "\n";
513  return os;
514 }
515 
516 std::istream&
517 operator >> (std::istream& is, ComplexColumnVector& a)
518 {
519  octave_idx_type len = a.length ();
520 
521  if (len > 0)
522  {
523  double tmp;
524  for (octave_idx_type i = 0; i < len; i++)
525  {
526  is >> tmp;
527  if (is)
528  a.elem (i) = tmp;
529  else
530  break;
531  }
532  }
533  return is;
534 }