Mercurial > octave-libgccjit

comparison liboctave/fCRowVector.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 // RowVector manipulations.
2 /*
3
4 Copyright (C) 1994, 1995, 1996, 1997, 1999, 2000, 2001, 2002, 2003,
5 2004, 2005, 2006, 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 (cgemv, CGEMV) (F77_CONST_CHAR_ARG_DECL,
45 const octave_idx_type&, const octave_idx_type&, const FloatComplex&,
46 const FloatComplex*, const octave_idx_type&, const FloatComplex*,
47 const octave_idx_type&, const FloatComplex&, FloatComplex*, const octave_idx_type&
48 F77_CHAR_ARG_LEN_DECL);
49
50 F77_RET_T
51 F77_FUNC (xcdotu, XCDOTU) (const octave_idx_type&, const FloatComplex*, const octave_idx_type&,
52 const FloatComplex*, const octave_idx_type&, FloatComplex&);
53 }
54
55 // FloatComplex Row Vector class
56
57 FloatComplexRowVector::FloatComplexRowVector (const FloatRowVector& a)
58 : MArray<FloatComplex> (a.length ())
59 {
60 for (octave_idx_type i = 0; i < length (); i++)
61 elem (i) = a.elem (i);
62 }
63
64 bool
65 FloatComplexRowVector::operator == (const FloatComplexRowVector& a) const
66 {
67 octave_idx_type len = length ();
68 if (len != a.length ())
69 return 0;
70 return mx_inline_equal (data (), a.data (), len);
71 }
72
73 bool
74 FloatComplexRowVector::operator != (const FloatComplexRowVector& a) const
75 {
76 return !(*this == a);
77 }
78
79 // destructive insert/delete/reorder operations
80
81 FloatComplexRowVector&
82 FloatComplexRowVector::insert (const FloatRowVector& a, octave_idx_type c)
83 {
84 octave_idx_type a_len = a.length ();
85
86 if (c < 0 || c + a_len > length ())
87 {
88 (*current_liboctave_error_handler) ("range error for insert");
89 return *this;
90 }
91
92 if (a_len > 0)
93 {
94 make_unique ();
95
96 for (octave_idx_type i = 0; i < a_len; i++)
97 xelem (c+i) = a.elem (i);
98 }
99
100 return *this;
101 }
102
103 FloatComplexRowVector&
104 FloatComplexRowVector::insert (const FloatComplexRowVector& a, octave_idx_type c)
105 {
106 octave_idx_type a_len = a.length ();
107
108 if (c < 0 || c + a_len > length ())
109 {
110 (*current_liboctave_error_handler) ("range error for insert");
111 return *this;
112 }
113
114 if (a_len > 0)
115 {
116 make_unique ();
117
118 for (octave_idx_type i = 0; i < a_len; i++)
119 xelem (c+i) = a.elem (i);
120 }
121
122 return *this;
123 }
124
125 FloatComplexRowVector&
126 FloatComplexRowVector::fill (float val)
127 {
128 octave_idx_type len = length ();
129
130 if (len > 0)
131 {
132 make_unique ();
133
134 for (octave_idx_type i = 0; i < len; i++)
135 xelem (i) = val;
136 }
137
138 return *this;
139 }
140
141 FloatComplexRowVector&
142 FloatComplexRowVector::fill (const FloatComplex& val)
143 {
144 octave_idx_type len = length ();
145
146 if (len > 0)
147 {
148 make_unique ();
149
150 for (octave_idx_type i = 0; i < len; i++)
151 xelem (i) = val;
152 }
153
154 return *this;
155 }
156
157 FloatComplexRowVector&
158 FloatComplexRowVector::fill (float val, octave_idx_type c1, octave_idx_type c2)
159 {
160 octave_idx_type len = length ();
161
162 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
163 {
164 (*current_liboctave_error_handler) ("range error for fill");
165 return *this;
166 }
167
168 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
169
170 if (c2 >= c1)
171 {
172 make_unique ();
173
174 for (octave_idx_type i = c1; i <= c2; i++)
175 xelem (i) = val;
176 }
177
178 return *this;
179 }
180
181 FloatComplexRowVector&
182 FloatComplexRowVector::fill (const FloatComplex& val, octave_idx_type c1, octave_idx_type c2)
183 {
184 octave_idx_type len = length ();
185
186 if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
187 {
188 (*current_liboctave_error_handler) ("range error for fill");
189 return *this;
190 }
191
192 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
193
194 if (c2 >= c1)
195 {
196 make_unique ();
197
198 for (octave_idx_type i = c1; i <= c2; i++)
199 xelem (i) = val;
200 }
201
202 return *this;
203 }
204
205 FloatComplexRowVector
206 FloatComplexRowVector::append (const FloatRowVector& a) const
207 {
208 octave_idx_type len = length ();
209 octave_idx_type nc_insert = len;
210 FloatComplexRowVector retval (len + a.length ());
211 retval.insert (*this, 0);
212 retval.insert (a, nc_insert);
213 return retval;
214 }
215
216 FloatComplexRowVector
217 FloatComplexRowVector::append (const FloatComplexRowVector& a) const
218 {
219 octave_idx_type len = length ();
220 octave_idx_type nc_insert = len;
221 FloatComplexRowVector retval (len + a.length ());
222 retval.insert (*this, 0);
223 retval.insert (a, nc_insert);
224 return retval;
225 }
226
227 FloatComplexColumnVector
228 FloatComplexRowVector::hermitian (void) const
229 {
230 return MArray<FloatComplex>::hermitian (std::conj);
231 }
232
233 FloatComplexColumnVector
234 FloatComplexRowVector::transpose (void) const
235 {
236 return MArray<FloatComplex>::transpose ();
237 }
238
239 FloatComplexRowVector
240 conj (const FloatComplexRowVector& a)
241 {
242 octave_idx_type a_len = a.length ();
243 FloatComplexRowVector retval;
244 if (a_len > 0)
245 retval = FloatComplexRowVector (mx_inline_conj_dup (a.data (), a_len), a_len);
246 return retval;
247 }
248
249 // resize is the destructive equivalent for this one
250
251 FloatComplexRowVector
252 FloatComplexRowVector::extract (octave_idx_type c1, octave_idx_type c2) const
253 {
254 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
255
256 octave_idx_type new_c = c2 - c1 + 1;
257
258 FloatComplexRowVector result (new_c);
259
260 for (octave_idx_type i = 0; i < new_c; i++)
261 result.elem (i) = elem (c1+i);
262
263 return result;
264 }
265
266 FloatComplexRowVector
267 FloatComplexRowVector::extract_n (octave_idx_type r1, octave_idx_type n) const
268 {
269 FloatComplexRowVector 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 // row vector by row vector -> row vector operations
278
279 FloatComplexRowVector&
280 FloatComplexRowVector::operator += (const FloatRowVector& a)
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 FloatComplex *d = fortran_vec (); // Ensures only one reference to my privates!
296
297 mx_inline_add2 (d, a.data (), len);
298 return *this;
299 }
300
301 FloatComplexRowVector&
302 FloatComplexRowVector::operator -= (const FloatRowVector& a)
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 FloatComplex *d = fortran_vec (); // Ensures only one reference to my privates!
318
319 mx_inline_subtract2 (d, a.data (), len);
320 return *this;
321 }
322
323 // row vector by matrix -> row vector
324
325 FloatComplexRowVector
326 operator * (const FloatComplexRowVector& v, const FloatComplexMatrix& a)
327 {
328 FloatComplexRowVector retval;
329
330 octave_idx_type len = v.length ();
331
332 octave_idx_type a_nr = a.rows ();
333 octave_idx_type a_nc = a.cols ();
334
335 if (a_nr != len)
336 gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
337 else
338 {
339 if (len == 0)
340 retval.resize (a_nc, 0.0);
341 else
342 {
343 // Transpose A to form A'*x == (x'*A)'
344
345 octave_idx_type ld = a_nr;
346
347 retval.resize (a_nc);
348 FloatComplex *y = retval.fortran_vec ();
349
350 F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 ("T", 1),
351 a_nr, a_nc, 1.0, a.data (),
352 ld, v.data (), 1, 0.0, y, 1
353 F77_CHAR_ARG_LEN (1)));
354 }
355 }
356
357 return retval;
358 }
359
360 FloatComplexRowVector
361 operator * (const FloatRowVector& v, const FloatComplexMatrix& a)
362 {
363 FloatComplexRowVector tmp (v);
364 return tmp * a;
365 }
366
367 // other operations
368
369 FloatRowVector
370 FloatComplexRowVector::map (dmapper fcn) const
371 {
372 return MArray<FloatComplex>::map<float> (func_ptr (fcn));
373 }
374
375 FloatComplexRowVector
376 FloatComplexRowVector::map (cmapper fcn) const
377 {
378 return MArray<FloatComplex>::map<FloatComplex> (func_ptr (fcn));
379 }
380
381 FloatComplex
382 FloatComplexRowVector::min (void) const
383 {
384 octave_idx_type len = length ();
385 if (len == 0)
386 return FloatComplex (0.0);
387
388 FloatComplex res = elem (0);
389 float 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 FloatComplex
402 FloatComplexRowVector::max (void) const
403 {
404 octave_idx_type len = length ();
405 if (len == 0)
406 return FloatComplex (0.0);
407
408 FloatComplex res = elem (0);
409 float absres = std::abs (res);
410
411 for (octave_idx_type i = 1; i < len; i++)
412 if (std::abs (elem (i)) > absres)
413 {
414 res = elem (i);
415 absres = std::abs (res);
416 }
417
418 return res;
419 }
420
421 // i/o
422
423 std::ostream&
424 operator << (std::ostream& os, const FloatComplexRowVector& a)
425 {
426 // int field_width = os.precision () + 7;
427 for (octave_idx_type i = 0; i < a.length (); i++)
428 os << " " /* setw (field_width) */ << a.elem (i);
429 return os;
430 }
431
432 std::istream&
433 operator >> (std::istream& is, FloatComplexRowVector& a)
434 {
435 octave_idx_type len = a.length();
436
437 if (len < 1)
438 is.clear (std::ios::badbit);
439 else
440 {
441 FloatComplex tmp;
442 for (octave_idx_type i = 0; i < len; i++)
443 {
444 is >> tmp;
445 if (is)
446 a.elem (i) = tmp;
447 else
448 break;
449 }
450 }
451 return is;
452 }
453
454 // row vector by column vector -> scalar
455
456 // row vector by column vector -> scalar
457
458 FloatComplex
459 operator * (const FloatComplexRowVector& v, const FloatColumnVector& a)
460 {
461 FloatComplexColumnVector tmp (a);
462 return v * tmp;
463 }
464
465 FloatComplex
466 operator * (const FloatComplexRowVector& v, const FloatComplexColumnVector& a)
467 {
468 FloatComplex retval (0.0, 0.0);
469
470 octave_idx_type len = v.length ();
471
472 octave_idx_type a_len = a.length ();
473
474 if (len != a_len)
475 gripe_nonconformant ("operator *", len, a_len);
476 else if (len != 0)
477 F77_FUNC (xcdotu, XCDOTU) (len, v.data (), 1, a.data (), 1, retval);
478
479 return retval;
480 }
481
482 // other operations
483
484 FloatComplexRowVector
485 linspace (const FloatComplex& x1, const FloatComplex& x2, octave_idx_type n)
486 {
487 FloatComplexRowVector retval;
488
489 if (n > 0)
490 {
491 retval.resize (n);
492 FloatComplex delta = (x2 - x1) / static_cast<float> (n - 1.0);
493 retval.elem (0) = x1;
494 for (octave_idx_type i = 1; i < n-1; i++)
495 retval.elem (i) = x1 + static_cast<float> (1.0) * i * delta;
496 retval.elem (n-1) = x2;
497 }
498 else
499 {
500 retval.resize (1);
501 retval.elem (0) = x2;
502 }
503
504 return retval;
505 }
506
507 /*
508 ;;; Local Variables: ***
509 ;;; mode: C++ ***
510 ;;; End: ***
511 */