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
dColVector.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 (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
46  const double&, const double*,
47  const octave_idx_type&, const double*,
48  const octave_idx_type&, const double&, double*,
49  const octave_idx_type&
51 }
52 
53 // Column Vector class.
54 
55 bool
57 {
58  octave_idx_type len = length ();
59  if (len != a.length ())
60  return 0;
61  return mx_inline_equal (len, data (), a.data ());
62 }
63 
64 bool
66 {
67  return !(*this == a);
68 }
69 
72 {
73  octave_idx_type a_len = a.length ();
74 
75  if (r < 0 || r + a_len > length ())
76  {
77  (*current_liboctave_error_handler) ("range error for insert");
78  return *this;
79  }
80 
81  if (a_len > 0)
82  {
83  make_unique ();
84 
85  for (octave_idx_type i = 0; i < a_len; i++)
86  xelem (r+i) = a.elem (i);
87  }
88 
89  return *this;
90 }
91 
93 ColumnVector::fill (double val)
94 {
95  octave_idx_type len = length ();
96 
97  if (len > 0)
98  {
99  make_unique ();
100 
101  for (octave_idx_type i = 0; i < len; i++)
102  xelem (i) = val;
103  }
104 
105  return *this;
106 }
107 
110 {
111  octave_idx_type len = length ();
112 
113  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
114  {
115  (*current_liboctave_error_handler) ("range error for fill");
116  return *this;
117  }
118 
119  if (r1 > r2) { std::swap (r1, r2); }
120 
121  if (r2 >= r1)
122  {
123  make_unique ();
124 
125  for (octave_idx_type i = r1; i <= r2; i++)
126  xelem (i) = val;
127  }
128 
129  return *this;
130 }
131 
134 {
135  octave_idx_type len = length ();
136  octave_idx_type nr_insert = len;
137  ColumnVector retval (len + a.length ());
138  retval.insert (*this, 0);
139  retval.insert (a, nr_insert);
140  return retval;
141 }
142 
143 RowVector
145 {
146  return MArray<double>::transpose ();
147 }
148 
150 ColumnVector::abs (void) const
151 {
152  return do_mx_unary_map<double, double, std::abs> (*this);
153 }
154 
157 {
158  return do_mx_unary_op<double, Complex> (a, mx_inline_real);
159 }
160 
163 {
164  return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
165 }
166 
167 // resize is the destructive equivalent for this one
168 
171 {
172  if (r1 > r2) { std::swap (r1, r2); }
173 
174  octave_idx_type new_r = r2 - r1 + 1;
175 
176  ColumnVector result (new_r);
177 
178  for (octave_idx_type i = 0; i < new_r; i++)
179  result.xelem (i) = elem (r1+i);
180 
181  return result;
182 }
183 
186 {
187  ColumnVector result (n);
188 
189  for (octave_idx_type i = 0; i < n; i++)
190  result.xelem (i) = elem (r1+i);
191 
192  return result;
193 }
194 
195 // matrix by column vector -> column vector operations
196 
198 operator * (const Matrix& m, const ColumnVector& a)
199 {
200  ColumnVector retval;
201 
202  octave_idx_type nr = m.rows ();
203  octave_idx_type nc = m.cols ();
204 
205  octave_idx_type a_len = a.length ();
206 
207  if (nc != a_len)
208  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
209  else
210  {
211  retval.clear (nr);
212 
213  if (nr != 0)
214  {
215  if (nc == 0)
216  retval.fill (0.0);
217  else
218  {
219  double *y = retval.fortran_vec ();
220 
221  F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
222  nr, nc, 1.0, m.data (), nr,
223  a.data (), 1, 0.0, y, 1
224  F77_CHAR_ARG_LEN (1)));
225  }
226  }
227  }
228 
229  return retval;
230 }
231 
232 // diagonal matrix by column vector -> column vector operations
233 
236 {
237  ColumnVector retval;
238 
239  octave_idx_type nr = m.rows ();
240  octave_idx_type nc = m.cols ();
241 
242  octave_idx_type a_len = a.length ();
243 
244  if (nc != a_len)
245  gripe_nonconformant ("operator *", nr, nc, a_len, 1);
246  else
247  {
248  if (nr == 0 || nc == 0)
249  retval.resize (nr, 0.0);
250  else
251  {
252  retval.resize (nr);
253 
254  for (octave_idx_type i = 0; i < a_len; i++)
255  retval.elem (i) = a.elem (i) * m.elem (i, i);
256 
257  for (octave_idx_type i = a_len; i < nr; i++)
258  retval.elem (i) = 0.0;
259  }
260  }
261 
262  return retval;
263 }
264 
265 // other operations
266 
267 double
268 ColumnVector::min (void) const
269 {
270  octave_idx_type len = length ();
271  if (len == 0)
272  return 0.0;
273 
274  double res = elem (0);
275 
276  for (octave_idx_type i = 1; i < len; i++)
277  if (elem (i) < res)
278  res = elem (i);
279 
280  return res;
281 }
282 
283 double
284 ColumnVector::max (void) const
285 {
286  octave_idx_type len = length ();
287  if (len == 0)
288  return 0.0;
289 
290  double res = elem (0);
291 
292  for (octave_idx_type i = 1; i < len; i++)
293  if (elem (i) > res)
294  res = elem (i);
295 
296  return res;
297 }
298 
299 std::ostream&
300 operator << (std::ostream& os, const ColumnVector& a)
301 {
302 // int field_width = os.precision () + 7;
303  for (octave_idx_type i = 0; i < a.length (); i++)
304  os << /* setw (field_width) << */ a.elem (i) << "\n";
305  return os;
306 }
307 
308 std::istream&
309 operator >> (std::istream& is, ColumnVector& a)
310 {
311  octave_idx_type len = a.length ();
312 
313  if (len > 0)
314  {
315  double tmp;
316  for (octave_idx_type i = 0; i < len; i++)
317  {
318  is >> tmp;
319  if (is)
320  a.elem (i) = tmp;
321  else
322  break;
323  }
324  }
325  return is;
326 }