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
CRowVector.cc
Go to the documentation of this file.
1 // RowVector manipulations.
2 /*
3 
4 Copyright (C) 1994-2013 John W. Eaton
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <iostream>
29 
30 #include "Array-util.h"
31 #include "f77-fcn.h"
32 #include "functor.h"
33 #include "lo-error.h"
34 #include "mx-base.h"
35 #include "mx-inlines.cc"
36 #include "oct-cmplx.h"
37 
38 // Fortran functions we call.
39 
40 extern "C"
41 {
42  F77_RET_T
43  F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
45  const Complex&, const Complex*,
46  const octave_idx_type&, const Complex*,
47  const octave_idx_type&, const Complex&, Complex*,
48  const octave_idx_type&
50 
51  F77_RET_T
52  F77_FUNC (xzdotu, XZDOTU) (const octave_idx_type&, const Complex*,
53  const octave_idx_type&, const Complex*,
54  const octave_idx_type&, Complex&);
55 }
56 
57 // Complex Row Vector class
58 
59 bool
61 {
62  octave_idx_type len = length ();
63  if (len != a.length ())
64  return 0;
65  return mx_inline_equal (len, data (), a.data ());
66 }
67 
68 bool
70 {
71  return !(*this == a);
72 }
73 
74 // destructive insert/delete/reorder operations
75 
78 {
79  octave_idx_type a_len = a.length ();
80 
81  if (c < 0 || c + a_len > length ())
82  {
83  (*current_liboctave_error_handler) ("range error for insert");
84  return *this;
85  }
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 {
101  octave_idx_type a_len = a.length ();
102 
103  if (c < 0 || c + a_len > length ())
104  {
105  (*current_liboctave_error_handler) ("range error for insert");
106  return *this;
107  }
108 
109  if (a_len > 0)
110  {
111  make_unique ();
112 
113  for (octave_idx_type i = 0; i < a_len; i++)
114  xelem (c+i) = a.elem (i);
115  }
116 
117  return *this;
118 }
119 
122 {
123  octave_idx_type len = length ();
124 
125  if (len > 0)
126  {
127  make_unique ();
128 
129  for (octave_idx_type i = 0; i < len; i++)
130  xelem (i) = val;
131  }
132 
133  return *this;
134 }
135 
138 {
139  octave_idx_type len = length ();
140 
141  if (len > 0)
142  {
143  make_unique ();
144 
145  for (octave_idx_type i = 0; i < len; i++)
146  xelem (i) = val;
147  }
148 
149  return *this;
150 }
151 
154 {
155  octave_idx_type len = length ();
156 
157  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
158  {
159  (*current_liboctave_error_handler) ("range error for fill");
160  return *this;
161  }
162 
163  if (c1 > c2) { std::swap (c1, c2); }
164 
165  if (c2 >= c1)
166  {
167  make_unique ();
168 
169  for (octave_idx_type i = c1; i <= c2; i++)
170  xelem (i) = val;
171  }
172 
173  return *this;
174 }
175 
179 {
180  octave_idx_type len = length ();
181 
182  if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
183  {
184  (*current_liboctave_error_handler) ("range error for fill");
185  return *this;
186  }
187 
188  if (c1 > c2) { std::swap (c1, c2); }
189 
190  if (c2 >= c1)
191  {
192  make_unique ();
193 
194  for (octave_idx_type i = c1; i <= c2; i++)
195  xelem (i) = val;
196  }
197 
198  return *this;
199 }
200 
203 {
204  octave_idx_type len = length ();
205  octave_idx_type nc_insert = len;
206  ComplexRowVector retval (len + a.length ());
207  retval.insert (*this, 0);
208  retval.insert (a, nc_insert);
209  return retval;
210 }
211 
214 {
215  octave_idx_type len = length ();
216  octave_idx_type nc_insert = len;
217  ComplexRowVector retval (len + a.length ());
218  retval.insert (*this, 0);
219  retval.insert (a, nc_insert);
220  return retval;
221 }
222 
225 {
227 }
228 
231 {
232  return MArray<Complex>::transpose ();
233 }
234 
237 {
238  return do_mx_unary_map<Complex, Complex, std::conj<double> > (a);
239 }
240 
241 // resize is the destructive equivalent for this one
242 
245 {
246  if (c1 > c2) { std::swap (c1, c2); }
247 
248  octave_idx_type new_c = c2 - c1 + 1;
249 
250  ComplexRowVector result (new_c);
251 
252  for (octave_idx_type i = 0; i < new_c; i++)
253  result.elem (i) = elem (c1+i);
254 
255  return result;
256 }
257 
260 {
261  ComplexRowVector result (n);
262 
263  for (octave_idx_type i = 0; i < n; i++)
264  result.elem (i) = elem (r1+i);
265 
266  return result;
267 }
268 
269 // row vector by row vector -> row vector operations
270 
273 {
274  octave_idx_type len = length ();
275 
276  octave_idx_type a_len = a.length ();
277 
278  if (len != a_len)
279  {
280  gripe_nonconformant ("operator +=", len, a_len);
281  return *this;
282  }
283 
284  if (len == 0)
285  return *this;
286 
287  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
288 
289  mx_inline_add2 (len, d, a.data ());
290  return *this;
291 }
292 
295 {
296  octave_idx_type len = length ();
297 
298  octave_idx_type a_len = a.length ();
299 
300  if (len != a_len)
301  {
302  gripe_nonconformant ("operator -=", len, a_len);
303  return *this;
304  }
305 
306  if (len == 0)
307  return *this;
308 
309  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
310 
311  mx_inline_sub2 (len, d, a.data ());
312  return *this;
313 }
314 
315 // row vector by matrix -> row vector
316 
319 {
320  ComplexRowVector retval;
321 
322  octave_idx_type len = v.length ();
323 
324  octave_idx_type a_nr = a.rows ();
325  octave_idx_type a_nc = a.cols ();
326 
327  if (a_nr != len)
328  gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
329  else
330  {
331  if (len == 0)
332  retval.resize (a_nc, 0.0);
333  else
334  {
335  // Transpose A to form A'*x == (x'*A)'
336 
337  octave_idx_type ld = a_nr;
338 
339  retval.resize (a_nc);
340  Complex *y = retval.fortran_vec ();
341 
342  F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
343  a_nr, a_nc, 1.0, a.data (),
344  ld, v.data (), 1, 0.0, y, 1
345  F77_CHAR_ARG_LEN (1)));
346  }
347  }
348 
349  return retval;
350 }
351 
354 {
355  ComplexRowVector tmp (v);
356  return tmp * a;
357 }
358 
359 // other operations
360 
361 Complex
363 {
364  octave_idx_type len = length ();
365  if (len == 0)
366  return Complex (0.0);
367 
368  Complex res = elem (0);
369  double absres = std::abs (res);
370 
371  for (octave_idx_type i = 1; i < len; i++)
372  if (std::abs (elem (i)) < absres)
373  {
374  res = elem (i);
375  absres = std::abs (res);
376  }
377 
378  return res;
379 }
380 
381 Complex
383 {
384  octave_idx_type len = length ();
385  if (len == 0)
386  return Complex (0.0);
387 
388  Complex res = elem (0);
389  double absres = std::abs (res);
390 
391  for (octave_idx_type i = 1; i < len; i++)
392  if (std::abs (elem (i)) > absres)
393  {
394  res = elem (i);
395  absres = std::abs (res);
396  }
397 
398  return res;
399 }
400 
401 // i/o
402 
403 std::ostream&
404 operator << (std::ostream& os, const ComplexRowVector& a)
405 {
406 // int field_width = os.precision () + 7;
407  for (octave_idx_type i = 0; i < a.length (); i++)
408  os << " " /* setw (field_width) */ << a.elem (i);
409  return os;
410 }
411 
412 std::istream&
413 operator >> (std::istream& is, ComplexRowVector& a)
414 {
415  octave_idx_type len = a.length ();
416 
417  if (len > 0)
418  {
419  Complex tmp;
420  for (octave_idx_type i = 0; i < len; i++)
421  {
422  is >> tmp;
423  if (is)
424  a.elem (i) = tmp;
425  else
426  break;
427  }
428  }
429  return is;
430 }
431 
432 // row vector by column vector -> scalar
433 
434 // row vector by column vector -> scalar
435 
436 Complex
438 {
439  ComplexColumnVector tmp (a);
440  return v * tmp;
441 }
442 
443 Complex
445 {
446  Complex retval (0.0, 0.0);
447 
448  octave_idx_type len = v.length ();
449 
450  octave_idx_type a_len = a.length ();
451 
452  if (len != a_len)
453  gripe_nonconformant ("operator *", len, a_len);
454  else if (len != 0)
455  F77_FUNC (xzdotu, XZDOTU) (len, v.data (), 1, a.data (), 1, retval);
456 
457  return retval;
458 }
459 
460 // other operations
461 
463 linspace (const Complex& x1, const Complex& x2, octave_idx_type n)
464 {
465  if (n < 1) n = 1;
466 
467  NoAlias<ComplexRowVector> retval (n);
468 
469  Complex delta = (x2 - x1) / (n - 1.0);
470  retval(0) = x1;
471  for (octave_idx_type i = 1; i < n-1; i++)
472  retval(i) = x1 + static_cast<double> (i)*delta;
473  retval(n-1) = x2;
474 
475  return retval;
476 }