Mercurial > hg > octave-nkf

comparison liboctave/CColVector.cc @ 458:38cb88095913

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-06-06 00:41:10 by jwe] Initial revision
author jwe
date Mon, 06 Jun 1994 00:41:10 +0000
parents
children 2ca256b77602
comparison
equal deleted inserted replaced
457:3d4b4f0fa5ba 458:38cb88095913
1 // ColumnVector manipulations. -*- C++ -*-
2 /*
3
4 Copyright (C) 1992, 1993, 1994 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 2, or (at your option) any
11 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, write to the Free
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27
28 #if defined (__GNUG__)
29 #pragma implementation
30 #endif
31
32 #include <iostream.h>
33
34 #include <Complex.h>
35
36 #include "mx-base.h"
37 #include "mx-inlines.cc"
38 #include "f77-uscore.h"
39 #include "lo-error.h"
40
41 // Fortran functions we call.
42
43 extern "C"
44 {
45 int F77_FCN (zgemm) (const char*, const char*, const int*,
46 const int*, const int*, const Complex*,
47 const Complex*, const int*, const Complex*,
48 const int*, const Complex*, Complex*, const int*,
49 long, long);
50 }
51
52 /*
53 * Complex Column Vector class
54 */
55
56 #define KLUDGE_VECTORS
57 #define TYPE Complex
58 #define KL_VEC_TYPE ComplexColumnVector
59 #include "mx-kludge.cc"
60 #undef KLUDGE_VECTORS
61 #undef TYPE
62 #undef KL_VEC_TYPE
63
64 ComplexColumnVector::ComplexColumnVector (const ColumnVector& a)
65 : Array<Complex> (a.length ())
66 {
67 for (int i = 0; i < length (); i++)
68 elem (i) = a.elem (i);
69 }
70
71 #if 0
72 ComplexColumnVector&
73 ComplexColumnVector::resize (int n)
74 {
75 if (n < 0)
76 {
77 (*current_liboctave_error_handler)
78 ("can't resize to negative dimension");
79 return *this;
80 }
81
82 Complex *new_data = (Complex *) NULL;
83 if (n > 0)
84 {
85 new_data = new Complex [n];
86 int min_len = len < n ? len : n;
87
88 for (int i = 0; i < min_len; i++)
89 new_data[i] = data[i];
90 }
91
92 delete [] data;
93 len = n;
94 data = new_data;
95
96 return *this;
97 }
98
99 ComplexColumnVector&
100 ComplexColumnVector::resize (int n, double val)
101 {
102 int old_len = len;
103 resize (n);
104 for (int i = old_len; i < len; i++)
105 data[i] = val;
106
107 return *this;
108 }
109
110 ComplexColumnVector&
111 ComplexColumnVector::resize (int n, const Complex& val)
112 {
113 int old_len = len;
114 resize (n);
115 for (int i = old_len; i < len; i++)
116 data[i] = val;
117
118 return *this;
119 }
120 #endif
121
122 int
123 ComplexColumnVector::operator == (const ComplexColumnVector& a) const
124 {
125 int len = length ();
126 if (len != a.length ())
127 return 0;
128 return equal (data (), a.data (), len);
129 }
130
131 int
132 ComplexColumnVector::operator != (const ComplexColumnVector& a) const
133 {
134 return !(*this == a);
135 }
136
137 // destructive insert/delete/reorder operations
138
139 ComplexColumnVector&
140 ComplexColumnVector::insert (const ColumnVector& a, int r)
141 {
142 int a_len = a.length ();
143 if (r < 0 || r + a_len - 1 > length ())
144 {
145 (*current_liboctave_error_handler) ("range error for insert");
146 return *this;
147 }
148
149 for (int i = 0; i < a_len; i++)
150 elem (r+i) = a.elem (i);
151
152 return *this;
153 }
154
155 ComplexColumnVector&
156 ComplexColumnVector::insert (const ComplexColumnVector& a, int r)
157 {
158 int a_len = a.length ();
159 if (r < 0 || r + a_len - 1 > length ())
160 {
161 (*current_liboctave_error_handler) ("range error for insert");
162 return *this;
163 }
164
165 for (int i = 0; i < a_len; i++)
166 elem (r+i) = a.elem (i);
167
168 return *this;
169 }
170
171 ComplexColumnVector&
172 ComplexColumnVector::fill (double val)
173 {
174 int len = length ();
175 if (len > 0)
176 for (int i = 0; i < len; i++)
177 elem (i) = val;
178 return *this;
179 }
180
181 ComplexColumnVector&
182 ComplexColumnVector::fill (const Complex& val)
183 {
184 int len = length ();
185 if (len > 0)
186 for (int i = 0; i < len; i++)
187 elem (i) = val;
188 return *this;
189 }
190
191 ComplexColumnVector&
192 ComplexColumnVector::fill (double val, int r1, int r2)
193 {
194 int len = length ();
195 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
196 {
197 (*current_liboctave_error_handler) ("range error for fill");
198 return *this;
199 }
200
201 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
202
203 for (int i = r1; i <= r2; i++)
204 elem (i) = val;
205
206 return *this;
207 }
208
209 ComplexColumnVector&
210 ComplexColumnVector::fill (const Complex& val, int r1, int r2)
211 {
212 int len = length ();
213 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
214 {
215 (*current_liboctave_error_handler) ("range error for fill");
216 return *this;
217 }
218
219 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
220
221 for (int i = r1; i <= r2; i++)
222 elem (i) = val;
223
224 return *this;
225 }
226
227 ComplexColumnVector
228 ComplexColumnVector::stack (const ColumnVector& a) const
229 {
230 int len = length ();
231 int nr_insert = len;
232 ComplexColumnVector retval (len + a.length ());
233 retval.insert (*this, 0);
234 retval.insert (a, nr_insert);
235 return retval;
236 }
237
238 ComplexColumnVector
239 ComplexColumnVector::stack (const ComplexColumnVector& a) const
240 {
241 int len = length ();
242 int nr_insert = len;
243 ComplexColumnVector retval (len + a.length ());
244 retval.insert (*this, 0);
245 retval.insert (a, nr_insert);
246 return retval;
247 }
248
249 ComplexRowVector
250 ComplexColumnVector::hermitian (void) const
251 {
252 int len = length ();
253 return ComplexRowVector (conj_dup (data (), len), len);
254 }
255
256 ComplexRowVector
257 ComplexColumnVector::transpose (void) const
258 {
259 int len = length ();
260 return ComplexRowVector (dup (data (), len), len);
261 }
262
263 ColumnVector
264 real (const ComplexColumnVector& a)
265 {
266 int a_len = a.length ();
267 ColumnVector retval;
268 if (a_len > 0)
269 retval = ColumnVector (real_dup (a.data (), a_len), a_len);
270 return retval;
271 }
272
273 ColumnVector
274 imag (const ComplexColumnVector& a)
275 {
276 int a_len = a.length ();
277 ColumnVector retval;
278 if (a_len > 0)
279 retval = ColumnVector (imag_dup (a.data (), a_len), a_len);
280 return retval;
281 }
282
283 ComplexColumnVector
284 conj (const ComplexColumnVector& a)
285 {
286 int a_len = a.length ();
287 ComplexColumnVector retval;
288 if (a_len > 0)
289 retval = ComplexColumnVector (conj_dup (a.data (), a_len), a_len);
290 return retval;
291 }
292
293 // resize is the destructive equivalent for this one
294
295 ComplexColumnVector
296 ComplexColumnVector::extract (int r1, int r2) const
297 {
298 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
299
300 int new_r = r2 - r1 + 1;
301
302 ComplexColumnVector result (new_r);
303
304 for (int i = 0; i < new_r; i++)
305 result.elem (i) = elem (r1+i);
306
307 return result;
308 }
309
310 // column vector by column vector -> column vector operations
311
312 ComplexColumnVector&
313 ComplexColumnVector::operator += (const ColumnVector& a)
314 {
315 int len = length ();
316 if (len != a.length ())
317 {
318 (*current_liboctave_error_handler)
319 ("nonconformant vector += operation attempted");
320 return *this;
321 }
322
323 if (len == 0)
324 return *this;
325
326 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
327
328 add2 (d, a.data (), len);
329 return *this;
330 }
331
332 ComplexColumnVector&
333 ComplexColumnVector::operator -= (const ColumnVector& a)
334 {
335 int len = length ();
336 if (len != a.length ())
337 {
338 (*current_liboctave_error_handler)
339 ("nonconformant vector -= operation attempted");
340 return *this;
341 }
342
343 if (len == 0)
344 return *this;
345
346 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
347
348 subtract2 (d, a.data (), len);
349 return *this;
350 }
351
352 ComplexColumnVector&
353 ComplexColumnVector::operator += (const ComplexColumnVector& a)
354 {
355 int len = length ();
356
357 if (len != a.length ())
358 {
359 (*current_liboctave_error_handler)
360 ("nonconformant vector += operation attempted");
361 return *this;
362 }
363
364 if (len == 0)
365 return *this;
366
367 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
368
369 add2 (d, a.data (), len);
370 return *this;
371 }
372
373 ComplexColumnVector&
374 ComplexColumnVector::operator -= (const ComplexColumnVector& a)
375 {
376 int len = length ();
377 if (len != a.length ())
378 {
379 (*current_liboctave_error_handler)
380 ("nonconformant vector -= operation attempted");
381 return *this;
382 }
383
384 if (len == 0)
385 return *this;
386
387 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
388
389 subtract2 (d, a.data (), len);
390 return *this;
391 }
392
393 // column vector by scalar -> column vector operations
394
395 ComplexColumnVector
396 operator + (const ComplexColumnVector& v, double s)
397 {
398 int len = v.length ();
399 return ComplexColumnVector (add (v.data (), len, s), len);
400 }
401
402 ComplexColumnVector
403 operator - (const ComplexColumnVector& v, double s)
404 {
405 int len = v.length ();
406 return ComplexColumnVector (subtract (v.data (), len, s), len);
407 }
408
409 ComplexColumnVector
410 operator * (const ComplexColumnVector& v, double s)
411 {
412 int len = v.length ();
413 return ComplexColumnVector (multiply (v.data (), len, s), len);
414 }
415
416 ComplexColumnVector
417 operator / (const ComplexColumnVector& v, double s)
418 {
419 int len = v.length ();
420 return ComplexColumnVector (divide (v.data (), len, s), len);
421 }
422
423 // scalar by column vector -> column vector operations
424
425 ComplexColumnVector
426 operator + (double s, const ComplexColumnVector& a)
427 {
428 int a_len = a.length ();
429 return ComplexColumnVector (add (a.data (), a_len, s), a_len);
430 }
431
432 ComplexColumnVector
433 operator - (double s, const ComplexColumnVector& a)
434 {
435 int a_len = a.length ();
436 return ComplexColumnVector (subtract (s, a.data (), a_len), a_len);
437 }
438
439 ComplexColumnVector
440 operator * (double s, const ComplexColumnVector& a)
441 {
442 int a_len = a.length ();
443 return ComplexColumnVector (multiply (a.data (), a_len, s), a_len);
444 }
445
446 ComplexColumnVector
447 operator / (double s, const ComplexColumnVector& a)
448 {
449 int a_len = a.length ();
450 return ComplexColumnVector (divide (s, a.data (), a_len), a_len);
451 }
452
453 // column vector by row vector -> matrix operations
454
455 ComplexMatrix
456 operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
457 {
458 int len = v.length ();
459 int a_len = a.length ();
460 if (len != a_len)
461 {
462 (*current_liboctave_error_handler)
463 ("nonconformant vector multiplication attempted");
464 return ComplexMatrix ();
465 }
466
467 if (len == 0)
468 return ComplexMatrix (len, len, 0.0);
469
470 char transa = 'N';
471 char transb = 'N';
472 Complex alpha (1.0);
473 Complex beta (0.0);
474 int anr = 1;
475
476 Complex *c = new Complex [len * a_len];
477
478 F77_FCN (zgemm) (&transa, &transb, &len, &a_len, &anr, &alpha,
479 v.data (), &len, a.data (), &anr, &beta, c, &len,
480 1L, 1L);
481
482 return ComplexMatrix (c, len, a_len);
483 }
484
485 // column vector by column vector -> column vector operations
486
487 ComplexColumnVector
488 operator + (const ComplexColumnVector& v, const ColumnVector& a)
489 {
490 int len = v.length ();
491 if (len != a.length ())
492 {
493 (*current_liboctave_error_handler)
494 ("nonconformant vector addition attempted");
495 return ComplexColumnVector ();
496 }
497
498 if (len == 0)
499 return ComplexColumnVector (0);
500
501 return ComplexColumnVector (add (v.data (), a.data (), len), len);
502 }
503
504 ComplexColumnVector
505 operator - (const ComplexColumnVector& v, const ColumnVector& a)
506 {
507 int len = v.length ();
508 if (len != a.length ())
509 {
510 (*current_liboctave_error_handler)
511 ("nonconformant vector subtraction attempted");
512 return ComplexColumnVector ();
513 }
514
515 if (len == 0)
516 return ComplexColumnVector (0);
517
518 return ComplexColumnVector (subtract (v.data (), a.data (), len), len);
519 }
520
521 ComplexColumnVector
522 product (const ComplexColumnVector& v, const ColumnVector& a)
523 {
524 int len = v.length ();
525 if (len != a.length ())
526 {
527 (*current_liboctave_error_handler)
528 ("nonconformant vector product attempted");
529 return ComplexColumnVector ();
530 }
531
532 if (len == 0)
533 return ComplexColumnVector (0);
534
535 return ComplexColumnVector (multiply (v.data (), a.data (), len), len);
536 }
537
538 ComplexColumnVector
539 quotient (const ComplexColumnVector& v, const ColumnVector& a)
540 {
541 int len = v.length ();
542 if (len != a.length ())
543 {
544 (*current_liboctave_error_handler)
545 ("nonconformant vector quotient attempted");
546 return ComplexColumnVector ();
547 }
548
549 if (len == 0)
550 return ComplexColumnVector (0);
551
552 return ComplexColumnVector (divide (v.data (), a.data (), len), len);
553 }
554
555 // other operations
556
557 ComplexColumnVector
558 map (c_c_Mapper f, const ComplexColumnVector& a)
559 {
560 ComplexColumnVector b (a);
561 b.map (f);
562 return b;
563 }
564
565 ColumnVector
566 map (d_c_Mapper f, const ComplexColumnVector& a)
567 {
568 int a_len = a.length ();
569 ColumnVector b (a_len);
570 for (int i = 0; i < a_len; i++)
571 b.elem (i) = f (a.elem (i));
572 return b;
573 }
574
575 void
576 ComplexColumnVector::map (c_c_Mapper f)
577 {
578 for (int i = 0; i < length (); i++)
579 elem (i) = f (elem (i));
580 }
581
582 Complex
583 ComplexColumnVector::min (void) const
584 {
585 int len = length ();
586 if (len == 0)
587 return 0.0;
588
589 Complex res = elem (0);
590 double absres = abs (res);
591
592 for (int i = 1; i < len; i++)
593 if (abs (elem (i)) < absres)
594 {
595 res = elem (i);
596 absres = abs (res);
597 }
598
599 return res;
600 }
601
602 Complex
603 ComplexColumnVector::max (void) const
604 {
605 int len = length ();
606 if (len == 0)
607 return 0.0;
608
609 Complex res = elem (0);
610 double absres = abs (res);
611
612 for (int i = 1; i < len; i++)
613 if (abs (elem (i)) > absres)
614 {
615 res = elem (i);
616 absres = abs (res);
617 }
618
619 return res;
620 }
621
622 // i/o
623
624 ostream&
625 operator << (ostream& os, const ComplexColumnVector& a)
626 {
627 // int field_width = os.precision () + 7;
628 for (int i = 0; i < a.length (); i++)
629 os << /* setw (field_width) << */ a.elem (i) << "\n";
630 return os;
631 }
632
633 istream&
634 operator >> (istream& is, ComplexColumnVector& a)
635 {
636 int len = a.length();
637
638 if (len < 1)
639 is.clear (ios::badbit);
640 else
641 {
642 double tmp;
643 for (int i = 0; i < len; i++)
644 {
645 is >> tmp;
646 if (is)
647 a.elem (i) = tmp;
648 else
649 break;
650 }
651 }
652 }
653
654 /*
655 ;;; Local Variables: ***
656 ;;; mode: C++ ***
657 ;;; page-delimiter: "^/\\*" ***
658 ;;; End: ***
659 */