Mercurial > hg > octave-lyh

comparison liboctave/fColVector.cc @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children eb63fbe60fab
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 // ColumnVector manipulations.
2 /*
3
4 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2004,
5 2005, 2007 John W. Eaton
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 (sgemv, SGEMV) (F77_CONST_CHAR_ARG_DECL,
45 const octave_idx_type&, const octave_idx_type&, const float&,
46 const float*, const octave_idx_type&, const float*,
47 const octave_idx_type&, const float&, float*,
48 const octave_idx_type&
49 F77_CHAR_ARG_LEN_DECL);
50 }
51
52 // Column Vector class.
53
54 bool
55 FloatColumnVector::operator == (const FloatColumnVector& a) const
56 {
57 octave_idx_type len = length ();
58 if (len != a.length ())
59 return 0;
60 return mx_inline_equal (data (), a.data (), len);
61 }
62
63 bool
64 FloatColumnVector::operator != (const FloatColumnVector& a) const
65 {
66 return !(*this == a);
67 }
68
69 FloatColumnVector&
70 FloatColumnVector::insert (const FloatColumnVector& a, octave_idx_type r)
71 {
72 octave_idx_type a_len = a.length ();
73
74 if (r < 0 || r + a_len > length ())
75 {
76 (*current_liboctave_error_handler) ("range error for insert");
77 return *this;
78 }
79
80 if (a_len > 0)
81 {
82 make_unique ();
83
84 for (octave_idx_type i = 0; i < a_len; i++)
85 xelem (r+i) = a.elem (i);
86 }
87
88 return *this;
89 }
90
91 FloatColumnVector&
92 FloatColumnVector::fill (float val)
93 {
94 octave_idx_type len = length ();
95
96 if (len > 0)
97 {
98 make_unique ();
99
100 for (octave_idx_type i = 0; i < len; i++)
101 xelem (i) = val;
102 }
103
104 return *this;
105 }
106
107 FloatColumnVector&
108 FloatColumnVector::fill (float val, octave_idx_type r1, octave_idx_type r2)
109 {
110 octave_idx_type len = length ();
111
112 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
113 {
114 (*current_liboctave_error_handler) ("range error for fill");
115 return *this;
116 }
117
118 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
119
120 if (r2 >= r1)
121 {
122 make_unique ();
123
124 for (octave_idx_type i = r1; i <= r2; i++)
125 xelem (i) = val;
126 }
127
128 return *this;
129 }
130
131 FloatColumnVector
132 FloatColumnVector::stack (const FloatColumnVector& a) const
133 {
134 octave_idx_type len = length ();
135 octave_idx_type nr_insert = len;
136 FloatColumnVector retval (len + a.length ());
137 retval.insert (*this, 0);
138 retval.insert (a, nr_insert);
139 return retval;
140 }
141
142 FloatRowVector
143 FloatColumnVector::transpose (void) const
144 {
145 return MArray<float>::transpose();
146 }
147
148 FloatColumnVector
149 real (const FloatComplexColumnVector& a)
150 {
151 octave_idx_type a_len = a.length ();
152 FloatColumnVector retval;
153 if (a_len > 0)
154 retval = FloatColumnVector (mx_inline_real_dup (a.data (), a_len), a_len);
155 return retval;
156 }
157
158 FloatColumnVector
159 imag (const FloatComplexColumnVector& a)
160 {
161 octave_idx_type a_len = a.length ();
162 FloatColumnVector retval;
163 if (a_len > 0)
164 retval = FloatColumnVector (mx_inline_imag_dup (a.data (), a_len), a_len);
165 return retval;
166 }
167
168 // resize is the destructive equivalent for this one
169
170 FloatColumnVector
171 FloatColumnVector::extract (octave_idx_type r1, octave_idx_type r2) const
172 {
173 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
174
175 octave_idx_type new_r = r2 - r1 + 1;
176
177 FloatColumnVector result (new_r);
178
179 for (octave_idx_type i = 0; i < new_r; i++)
180 result.xelem (i) = elem (r1+i);
181
182 return result;
183 }
184
185 FloatColumnVector
186 FloatColumnVector::extract_n (octave_idx_type r1, octave_idx_type n) const
187 {
188 FloatColumnVector result (n);
189
190 for (octave_idx_type i = 0; i < n; i++)
191 result.xelem (i) = elem (r1+i);
192
193 return result;
194 }
195
196 // matrix by column vector -> column vector operations
197
198 FloatColumnVector
199 operator * (const FloatMatrix& m, const FloatColumnVector& a)
200 {
201 FloatColumnVector retval;
202
203 octave_idx_type nr = m.rows ();
204 octave_idx_type nc = m.cols ();
205
206 octave_idx_type a_len = a.length ();
207
208 if (nc != a_len)
209 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
210 else
211 {
212 if (nr == 0 || nc == 0)
213 retval.resize (nr, 0.0);
214 else
215 {
216 octave_idx_type ld = nr;
217
218 retval.resize (nr);
219 float *y = retval.fortran_vec ();
220
221 F77_XFCN (sgemv, SGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
222 nr, nc, 1.0, m.data (), ld,
223 a.data (), 1, 0.0, y, 1
224 F77_CHAR_ARG_LEN (1)));
225 }
226 }
227
228 return retval;
229 }
230
231 // diagonal matrix by column vector -> column vector operations
232
233 FloatColumnVector
234 operator * (const FloatDiagMatrix& m, const FloatColumnVector& a)
235 {
236 FloatColumnVector retval;
237
238 octave_idx_type nr = m.rows ();
239 octave_idx_type nc = m.cols ();
240
241 octave_idx_type a_len = a.length ();
242
243 if (nc != a_len)
244 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
245 else
246 {
247 if (nr == 0 || nc == 0)
248 retval.resize (nr, 0.0);
249 else
250 {
251 retval.resize (nr);
252
253 for (octave_idx_type i = 0; i < a_len; i++)
254 retval.elem (i) = a.elem (i) * m.elem (i, i);
255
256 for (octave_idx_type i = a_len; i < nr; i++)
257 retval.elem (i) = 0.0;
258 }
259 }
260
261 return retval;
262 }
263
264 // other operations
265
266 FloatColumnVector
267 FloatColumnVector::map (dmapper fcn) const
268 {
269 return MArray<float>::map<float> (func_ptr (fcn));
270 }
271
272 FloatComplexColumnVector
273 FloatColumnVector::map (cmapper fcn) const
274 {
275 return MArray<float>::map<FloatComplex> (func_ptr (fcn));
276 }
277
278 float
279 FloatColumnVector::min (void) const
280 {
281 octave_idx_type len = length ();
282 if (len == 0)
283 return 0.0;
284
285 float res = elem (0);
286
287 for (octave_idx_type i = 1; i < len; i++)
288 if (elem (i) < res)
289 res = elem (i);
290
291 return res;
292 }
293
294 float
295 FloatColumnVector::max (void) const
296 {
297 octave_idx_type len = length ();
298 if (len == 0)
299 return 0.0;
300
301 float res = elem (0);
302
303 for (octave_idx_type i = 1; i < len; i++)
304 if (elem (i) > res)
305 res = elem (i);
306
307 return res;
308 }
309
310 std::ostream&
311 operator << (std::ostream& os, const FloatColumnVector& a)
312 {
313 // int field_width = os.precision () + 7;
314 for (octave_idx_type i = 0; i < a.length (); i++)
315 os << /* setw (field_width) << */ a.elem (i) << "\n";
316 return os;
317 }
318
319 std::istream&
320 operator >> (std::istream& is, FloatColumnVector& a)
321 {
322 octave_idx_type len = a.length();
323
324 if (len < 1)
325 is.clear (std::ios::badbit);
326 else
327 {
328 float tmp;
329 for (octave_idx_type i = 0; i < len; i++)
330 {
331 is >> tmp;
332 if (is)
333 a.elem (i) = tmp;
334 else
335 break;
336 }
337 }
338 return is;
339 }
340
341 /*
342 ;;; Local Variables: ***
343 ;;; mode: C++ ***
344 ;;; End: ***
345 */