3
|
1 // ColumnVector manipulations. -*- C++ -*- |
|
2 /* |
|
3 |
|
4 Copyright (C) 1992, 1993 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 // I\'m not sure how this is supposed to work if the .h file declares |
|
25 // several classes, each of which is defined in a separate file... |
|
26 // |
|
27 // #ifdef __GNUG__ |
|
28 // #pragma implementation "Matrix.h" |
|
29 // #endif |
|
30 |
|
31 #include "Matrix.h" |
|
32 #include "mx-inlines.cc" |
227
|
33 #include "lo-error.h" |
3
|
34 |
|
35 /* |
|
36 * Column Vector class. |
|
37 */ |
|
38 |
|
39 ColumnVector::ColumnVector (int n) |
|
40 { |
|
41 if (n < 0) |
227
|
42 { |
|
43 (*current_liboctave_error_handler) |
|
44 ("can't create column vector with negative dimension"); |
|
45 len = 0; |
|
46 data = (double *) NULL; |
|
47 return; |
|
48 } |
3
|
49 |
|
50 len = n; |
|
51 if (n > 0) |
|
52 data = new double [len]; |
|
53 else |
|
54 data = (double *) NULL; |
|
55 } |
|
56 |
|
57 ColumnVector::ColumnVector (int n, double val) |
|
58 { |
|
59 if (n < 0) |
227
|
60 { |
|
61 (*current_liboctave_error_handler) |
|
62 ("can't create column vector with negative dimension"); |
|
63 len = 0; |
|
64 data = (double *) NULL; |
|
65 return; |
|
66 } |
3
|
67 |
|
68 len = n; |
|
69 if (n > 0) |
|
70 { |
|
71 data = new double [len]; |
|
72 copy (data, len, val); |
|
73 } |
|
74 else |
|
75 data = (double *) NULL; |
|
76 } |
|
77 |
|
78 ColumnVector::ColumnVector (const ColumnVector& a) |
|
79 { |
|
80 len = a.len; |
|
81 if (len > 0) |
|
82 { |
|
83 data = new double [len]; |
|
84 copy (data, a.data, len); |
|
85 } |
|
86 else |
|
87 data = (double *) NULL; |
|
88 } |
|
89 |
|
90 ColumnVector::ColumnVector (double a) |
|
91 { |
|
92 len = 1; |
|
93 data = new double [1]; |
|
94 data[0] = a; |
|
95 } |
|
96 |
|
97 ColumnVector& |
|
98 ColumnVector::operator = (const ColumnVector& a) |
|
99 { |
|
100 if (this != &a) |
|
101 { |
|
102 delete [] data; |
|
103 len = a.len; |
|
104 if (len > 0) |
|
105 { |
|
106 data = new double [len]; |
|
107 copy (data, a.data, len); |
|
108 } |
|
109 else |
|
110 data = (double *) NULL; |
|
111 } |
|
112 return *this; |
|
113 } |
|
114 |
227
|
115 double& |
|
116 ColumnVector::checkelem (int n) |
|
117 { |
|
118 #ifndef NO_RANGE_CHECK |
|
119 if (n < 0 || n >= len) |
|
120 { |
|
121 (*current_liboctave_error_handler) ("range error"); |
|
122 static double foo = 0.0; |
|
123 return foo; |
|
124 } |
|
125 #endif |
|
126 |
|
127 return elem (n); |
|
128 } |
|
129 |
|
130 double |
|
131 ColumnVector::checkelem (int n) const |
|
132 { |
|
133 #ifndef NO_RANGE_CHECK |
|
134 if (n < 0 || n >= len) |
|
135 { |
|
136 (*current_liboctave_error_handler) ("range error"); |
|
137 return 0.0; |
|
138 } |
|
139 #endif |
|
140 |
|
141 return elem (n); |
|
142 } |
|
143 |
3
|
144 ColumnVector& |
|
145 ColumnVector::resize (int n) |
|
146 { |
|
147 if (n < 0) |
227
|
148 { |
|
149 (*current_liboctave_error_handler) |
|
150 ("can't resize to negative dimension"); |
|
151 return *this; |
|
152 } |
3
|
153 |
|
154 double *new_data = (double *) NULL; |
|
155 if (n > 0) |
|
156 { |
|
157 new_data = new double [n]; |
|
158 int min_len = len < n ? len : n; |
|
159 |
|
160 for (int i = 0; i < min_len; i++) |
|
161 new_data[i] = data[i]; |
|
162 } |
|
163 |
|
164 delete [] data; |
|
165 len = n; |
|
166 data = new_data; |
|
167 |
|
168 return *this; |
|
169 } |
|
170 |
|
171 ColumnVector& |
|
172 ColumnVector::resize (int n, double val) |
|
173 { |
|
174 int old_len = len; |
|
175 resize (n); |
|
176 for (int i = old_len; i < len; i++) |
|
177 data[i] = val; |
|
178 |
|
179 return *this; |
|
180 } |
|
181 |
|
182 int |
|
183 ColumnVector::operator == (const ColumnVector& a) const |
|
184 { |
|
185 if (len != a.len) |
|
186 return 0; |
|
187 return equal (data, a.data, len); |
|
188 } |
|
189 |
|
190 int |
|
191 ColumnVector::operator != (const ColumnVector& a) const |
|
192 { |
|
193 if (len != a.len) |
|
194 return 1; |
|
195 return !equal (data, a.data, len); |
|
196 } |
|
197 |
|
198 ColumnVector& |
|
199 ColumnVector::insert (const ColumnVector& a, int r) |
|
200 { |
|
201 if (r < 0 || r + a.len - 1 > len) |
227
|
202 { |
|
203 (*current_liboctave_error_handler) ("range error for insert"); |
|
204 return *this; |
|
205 } |
3
|
206 |
|
207 for (int i = 0; i < a.len; i++) |
|
208 data[r+i] = a.data[i]; |
|
209 |
|
210 return *this; |
|
211 } |
|
212 |
|
213 ColumnVector& |
|
214 ColumnVector::fill (double val) |
|
215 { |
|
216 if (len > 0) |
|
217 copy (data, len, val); |
|
218 return *this; |
|
219 } |
|
220 |
|
221 ColumnVector& |
|
222 ColumnVector::fill (double val, int r1, int r2) |
|
223 { |
|
224 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) |
227
|
225 { |
|
226 (*current_liboctave_error_handler) ("range error for fill"); |
|
227 return *this; |
|
228 } |
3
|
229 |
|
230 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
231 |
|
232 for (int i = r1; i <= r2; i++) |
|
233 data[i] = val; |
|
234 |
|
235 return *this; |
|
236 } |
|
237 |
|
238 ColumnVector |
|
239 ColumnVector::stack (const ColumnVector& a) const |
|
240 { |
|
241 int nr_insert = len; |
|
242 ColumnVector retval (len + a.len); |
|
243 retval.insert (*this, 0); |
|
244 retval.insert (a, nr_insert); |
|
245 return retval; |
|
246 } |
|
247 |
|
248 RowVector |
|
249 ColumnVector::transpose (void) const |
|
250 { |
|
251 return RowVector (dup (data, len), len); |
|
252 } |
|
253 |
|
254 // resize is the destructive equivalent for this one |
|
255 |
|
256 ColumnVector |
|
257 ColumnVector::extract (int r1, int r2) const |
|
258 { |
|
259 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
260 |
|
261 int new_r = r2 - r1 + 1; |
|
262 |
|
263 ColumnVector result (new_r); |
|
264 |
|
265 for (int i = 0; i < new_r; i++) |
|
266 result.data[i] = elem (r1+i); |
|
267 |
|
268 return result; |
|
269 } |
|
270 |
|
271 // column vector by scalar -> column vector operations |
|
272 |
|
273 ColumnVector |
|
274 ColumnVector::operator + (double s) const |
|
275 { |
|
276 return ColumnVector (add (data, len, s), len); |
|
277 } |
|
278 |
|
279 ColumnVector |
|
280 ColumnVector::operator - (double s) const |
|
281 { |
|
282 return ColumnVector (subtract (data, len, s), len); |
|
283 } |
|
284 |
|
285 ColumnVector |
|
286 ColumnVector::operator * (double s) const |
|
287 { |
|
288 return ColumnVector (multiply (data, len, s), len); |
|
289 } |
|
290 |
|
291 ColumnVector |
|
292 ColumnVector::operator / (double s) const |
|
293 { |
|
294 return ColumnVector (divide (data, len, s), len); |
|
295 } |
|
296 |
|
297 // scalar by column vector -> column vector operations |
|
298 |
|
299 ColumnVector |
|
300 operator + (double s, const ColumnVector& a) |
|
301 { |
|
302 return ColumnVector (add (a.data, a.len, s), a.len); |
|
303 } |
|
304 |
|
305 ColumnVector |
|
306 operator - (double s, const ColumnVector& a) |
|
307 { |
|
308 return ColumnVector (subtract (s, a.data, a.len), a.len); |
|
309 } |
|
310 |
|
311 ColumnVector |
|
312 operator * (double s, const ColumnVector& a) |
|
313 { |
|
314 return ColumnVector (multiply (a.data, a.len, s), a.len); |
|
315 } |
|
316 |
|
317 ColumnVector |
|
318 operator / (double s, const ColumnVector& a) |
|
319 { |
|
320 return ColumnVector (divide (s, a.data, a.len), a.len); |
|
321 } |
|
322 |
|
323 ComplexColumnVector |
161
|
324 ColumnVector::operator + (const Complex& s) const |
3
|
325 { |
|
326 return ComplexColumnVector (add (data, len, s), len); |
|
327 } |
|
328 |
|
329 ComplexColumnVector |
161
|
330 ColumnVector::operator - (const Complex& s) const |
3
|
331 { |
|
332 return ComplexColumnVector (subtract (data, len, s), len); |
|
333 } |
|
334 |
|
335 ComplexColumnVector |
161
|
336 ColumnVector::operator * (const Complex& s) const |
3
|
337 { |
|
338 return ComplexColumnVector (multiply (data, len, s), len); |
|
339 } |
|
340 |
|
341 ComplexColumnVector |
161
|
342 ColumnVector::operator / (const Complex& s) const |
3
|
343 { |
|
344 return ComplexColumnVector (divide (data, len, s), len); |
|
345 } |
|
346 |
|
347 // column vector by row vector -> matrix operations |
|
348 |
|
349 Matrix |
|
350 ColumnVector::operator * (const RowVector& a) const |
|
351 { |
|
352 if (len != a.len) |
227
|
353 { |
|
354 (*current_liboctave_error_handler) |
|
355 ("nonconformant vector multiplication attempted"); |
|
356 return Matrix (); |
|
357 } |
3
|
358 |
|
359 if (len == 0) |
|
360 return Matrix (len, len, 0.0); |
|
361 |
|
362 char transa = 'N'; |
|
363 char transb = 'N'; |
|
364 double alpha = 1.0; |
|
365 double beta = 0.0; |
|
366 int anr = 1; |
|
367 int anc = a.len; |
|
368 |
|
369 double *c = new double [len * a.len]; |
|
370 |
|
371 F77_FCN (dgemm) (&transa, &transb, &len, &anc, &anr, &alpha, data, |
|
372 &len, a.data, &anr, &beta, c, &len, 1L, 1L); |
|
373 |
|
374 return Matrix (c, len, a.len); |
|
375 } |
|
376 |
|
377 ComplexMatrix |
|
378 ColumnVector::operator * (const ComplexRowVector& a) const |
|
379 { |
|
380 ComplexColumnVector tmp (*this); |
|
381 return tmp * a; |
|
382 } |
|
383 |
|
384 // column vector by column vector -> column vector operations |
|
385 |
|
386 ColumnVector |
|
387 ColumnVector::operator + (const ColumnVector& a) const |
|
388 { |
|
389 if (len != a.len) |
227
|
390 { |
|
391 (*current_liboctave_error_handler) |
|
392 ("nonconformant vector addition attempted"); |
|
393 return ColumnVector (); |
|
394 } |
3
|
395 |
|
396 if (len == 0) |
|
397 return ColumnVector (0); |
|
398 |
|
399 return ColumnVector (add (data, a.data, len), len); |
|
400 } |
|
401 |
|
402 ColumnVector |
|
403 ColumnVector::operator - (const ColumnVector& a) const |
|
404 { |
|
405 if (len != a.len) |
227
|
406 { |
|
407 (*current_liboctave_error_handler) |
|
408 ("nonconformant vector subtraction attempted"); |
|
409 return ColumnVector (); |
|
410 } |
3
|
411 |
|
412 if (len == 0) |
|
413 return ColumnVector (0); |
|
414 |
|
415 return ColumnVector (subtract (data, a.data, len), len); |
|
416 } |
|
417 |
|
418 ComplexColumnVector |
|
419 ColumnVector::operator + (const ComplexColumnVector& a) const |
|
420 { |
|
421 if (len != a.len) |
227
|
422 { |
|
423 (*current_liboctave_error_handler) |
|
424 ("nonconformant vector addition attempted"); |
|
425 return ComplexColumnVector (); |
|
426 } |
3
|
427 |
|
428 if (len == 0) |
|
429 return ComplexColumnVector (0); |
|
430 |
|
431 return ComplexColumnVector (add (data, a.data, len), len); |
|
432 } |
|
433 |
|
434 ComplexColumnVector |
|
435 ColumnVector::operator - (const ComplexColumnVector& a) const |
|
436 { |
|
437 if (len != a.len) |
227
|
438 { |
|
439 (*current_liboctave_error_handler) |
|
440 ("nonconformant vector subtraction attempted"); |
|
441 return ComplexColumnVector (); |
|
442 } |
3
|
443 |
|
444 if (len == 0) |
|
445 return ComplexColumnVector (0); |
|
446 |
|
447 return ComplexColumnVector (subtract (data, a.data, len), len); |
|
448 } |
|
449 |
|
450 ColumnVector |
|
451 ColumnVector::product (const ColumnVector& a) const |
|
452 { |
|
453 if (len != a.len) |
227
|
454 { |
|
455 (*current_liboctave_error_handler) |
|
456 ("nonconformant vector product attempted"); |
|
457 return ColumnVector (); |
|
458 } |
3
|
459 |
|
460 if (len == 0) |
|
461 return ColumnVector (0); |
|
462 |
|
463 return ColumnVector (multiply (data, a.data, len), len); |
|
464 } |
|
465 |
|
466 ColumnVector |
|
467 ColumnVector::quotient (const ColumnVector& a) const |
|
468 { |
|
469 if (len != a.len) |
227
|
470 { |
|
471 (*current_liboctave_error_handler) |
|
472 ("nonconformant vector quotient attempted"); |
|
473 return ColumnVector (); |
|
474 } |
3
|
475 |
|
476 if (len == 0) |
|
477 return ColumnVector (0); |
|
478 |
|
479 return ColumnVector (divide (data, a.data, len), len); |
|
480 } |
|
481 |
|
482 ComplexColumnVector |
|
483 ColumnVector::product (const ComplexColumnVector& a) const |
|
484 { |
|
485 if (len != a.len) |
227
|
486 { |
|
487 (*current_liboctave_error_handler) |
|
488 ("nonconformant vector product attempted"); |
|
489 return ColumnVector (); |
|
490 } |
3
|
491 |
|
492 if (len == 0) |
|
493 return ComplexColumnVector (0); |
|
494 |
|
495 return ComplexColumnVector (multiply (data, a.data, len), len); |
|
496 } |
|
497 |
|
498 ComplexColumnVector |
|
499 ColumnVector::quotient (const ComplexColumnVector& a) const |
|
500 { |
|
501 if (len != a.len) |
227
|
502 { |
|
503 (*current_liboctave_error_handler) |
|
504 ("nonconformant vector quotient attempted"); |
|
505 return ColumnVector (); |
|
506 } |
3
|
507 |
|
508 if (len == 0) |
|
509 return ComplexColumnVector (0); |
|
510 |
|
511 return ComplexColumnVector (divide (data, a.data, len), len); |
|
512 } |
|
513 |
|
514 ColumnVector& |
|
515 ColumnVector::operator += (const ColumnVector& a) |
|
516 { |
|
517 if (len != a.len) |
227
|
518 { |
|
519 (*current_liboctave_error_handler) |
|
520 ("nonconformant vector += operation attempted"); |
|
521 return ColumnVector (); |
|
522 } |
3
|
523 |
|
524 if (len == 0) |
|
525 return *this; |
|
526 |
|
527 add2 (data, a.data, len); |
|
528 return *this; |
|
529 } |
|
530 |
|
531 ColumnVector& |
|
532 ColumnVector::operator -= (const ColumnVector& a) |
|
533 { |
|
534 if (len != a.len) |
227
|
535 { |
|
536 (*current_liboctave_error_handler) |
|
537 ("nonconformant vector -= operation attempted"); |
|
538 return ColumnVector (); |
|
539 } |
3
|
540 |
|
541 if (len == 0) |
|
542 return *this; |
|
543 |
|
544 subtract2 (data, a.data, len); |
|
545 return *this; |
|
546 } |
|
547 |
|
548 // unary operations |
|
549 |
|
550 ColumnVector |
|
551 ColumnVector::operator - (void) const |
|
552 { |
|
553 if (len == 0) |
|
554 return ColumnVector (0); |
|
555 |
|
556 return ColumnVector (negate (data, len), len); |
|
557 } |
|
558 |
|
559 ColumnVector |
|
560 map (d_d_Mapper f, const ColumnVector& a) |
|
561 { |
|
562 ColumnVector b (a); |
|
563 b.map (f); |
|
564 return b; |
|
565 } |
|
566 |
|
567 void |
|
568 ColumnVector::map (d_d_Mapper f) |
|
569 { |
|
570 for (int i = 0; i < len; i++) |
|
571 data[i] = f (data[i]); |
|
572 } |
|
573 |
|
574 double |
|
575 ColumnVector::min (void) const |
|
576 { |
|
577 if (len == 0) |
|
578 return 0.0; |
|
579 |
|
580 double res = data[0]; |
|
581 |
|
582 for (int i = 1; i < len; i++) |
|
583 if (data[i] < res) |
|
584 res = data[i]; |
|
585 |
|
586 return res; |
|
587 } |
|
588 |
|
589 double |
|
590 ColumnVector::max (void) const |
|
591 { |
|
592 if (len == 0) |
|
593 return 0.0; |
|
594 |
|
595 double res = data[0]; |
|
596 |
|
597 for (int i = 1; i < len; i++) |
|
598 if (data[i] > res) |
|
599 res = data[i]; |
|
600 |
|
601 return res; |
|
602 } |
|
603 |
|
604 ostream& |
|
605 operator << (ostream& os, const ColumnVector& a) |
|
606 { |
|
607 // int field_width = os.precision () + 7; |
|
608 for (int i = 0; i < a.len; i++) |
|
609 os << /* setw (field_width) << */ a.data[i] << "\n"; |
|
610 return os; |
|
611 } |
|
612 |
|
613 /* |
|
614 * Complex Column Vector class |
|
615 */ |
|
616 |
|
617 ComplexColumnVector::ComplexColumnVector (int n) |
|
618 { |
|
619 if (n < 0) |
227
|
620 { |
|
621 (*current_liboctave_error_handler) |
|
622 ("can't create column vector with negative dimension"); |
|
623 len = 0; |
|
624 data = (Complex *) NULL; |
|
625 return; |
|
626 } |
3
|
627 |
|
628 len = n; |
|
629 if (n > 0) |
|
630 data = new Complex [len]; |
|
631 else |
|
632 data = (Complex *) NULL; |
|
633 } |
|
634 |
|
635 ComplexColumnVector::ComplexColumnVector (int n, double val) |
|
636 { |
|
637 if (n < 0) |
227
|
638 { |
|
639 (*current_liboctave_error_handler) |
|
640 ("can't create column vector with negative dimension"); |
|
641 len = 0; |
|
642 data = (Complex *) NULL; |
|
643 return; |
|
644 } |
3
|
645 |
|
646 len = n; |
|
647 if (n > 0) |
|
648 { |
|
649 data = new Complex [len]; |
|
650 copy (data, len, val); |
|
651 } |
|
652 else |
|
653 data = (Complex *) NULL; |
|
654 } |
|
655 |
161
|
656 ComplexColumnVector::ComplexColumnVector (int n, const Complex& val) |
3
|
657 { |
|
658 if (n < 0) |
227
|
659 { |
|
660 (*current_liboctave_error_handler) |
|
661 ("can't create column vector with negative dimension"); |
|
662 len = 0; |
|
663 data = (Complex *) NULL; |
|
664 return; |
|
665 } |
3
|
666 |
|
667 len = n; |
|
668 if (n > 0) |
|
669 { |
|
670 data = new Complex [len]; |
|
671 copy (data, len, val); |
|
672 } |
|
673 else |
|
674 data = (Complex *) NULL; |
|
675 } |
|
676 |
|
677 ComplexColumnVector::ComplexColumnVector (const ColumnVector& a) |
|
678 { |
|
679 len = a.len; |
|
680 if (len > 0) |
|
681 { |
|
682 data = new Complex [len]; |
|
683 copy (data, a.data, len); |
|
684 } |
|
685 else |
|
686 data = (Complex *) NULL; |
|
687 } |
|
688 |
|
689 ComplexColumnVector::ComplexColumnVector (const ComplexColumnVector& a) |
|
690 { |
|
691 len = a.len; |
|
692 if (len > 0) |
|
693 { |
|
694 data = new Complex [len]; |
|
695 copy (data, a.data, len); |
|
696 } |
|
697 else |
|
698 data = (Complex *) NULL; |
|
699 } |
|
700 |
|
701 ComplexColumnVector::ComplexColumnVector (double a) |
|
702 { |
|
703 len = 1; |
|
704 data = new Complex [1]; |
|
705 data[0] = a; |
|
706 } |
|
707 |
161
|
708 ComplexColumnVector::ComplexColumnVector (const Complex& a) |
3
|
709 { |
|
710 len = 1; |
|
711 data = new Complex [1]; |
|
712 data[0] = Complex (a); |
|
713 } |
|
714 |
|
715 ComplexColumnVector& |
|
716 ComplexColumnVector::operator = (const ColumnVector& a) |
|
717 { |
|
718 delete [] data; |
|
719 len = a.len; |
|
720 if (len > 0) |
|
721 { |
|
722 data = new Complex [len]; |
|
723 copy (data, a.data, len); |
|
724 } |
|
725 else |
|
726 data = (Complex *) NULL; |
|
727 |
|
728 return *this; |
|
729 } |
|
730 |
|
731 ComplexColumnVector& |
|
732 ComplexColumnVector::operator = (const ComplexColumnVector& a) |
|
733 { |
|
734 if (this != &a) |
|
735 { |
|
736 delete [] data; |
|
737 len = a.len; |
|
738 if (len > 0) |
|
739 { |
|
740 data = new Complex [len]; |
|
741 copy (data, a.data, len); |
|
742 } |
|
743 else |
|
744 data = (Complex *) NULL; |
|
745 } |
|
746 return *this; |
|
747 } |
|
748 |
227
|
749 Complex& |
|
750 ComplexColumnVector::checkelem (int n) |
|
751 { |
|
752 #ifndef NO_RANGE_CHECK |
|
753 if (n < 0 || n >= len) |
|
754 { |
|
755 (*current_liboctave_error_handler) ("range error"); |
|
756 static Complex foo (0.0); |
|
757 return foo; |
|
758 } |
|
759 #endif |
|
760 |
|
761 return elem (n); |
|
762 } |
|
763 |
|
764 Complex |
|
765 ComplexColumnVector::checkelem (int n) const |
|
766 { |
|
767 #ifndef NO_RANGE_CHECK |
|
768 if (n < 0 || n >= len) |
|
769 { |
|
770 (*current_liboctave_error_handler) ("range error"); |
|
771 return Complex (0.0); |
|
772 } |
|
773 #endif |
|
774 |
|
775 return elem (n); |
|
776 } |
|
777 |
3
|
778 ComplexColumnVector& |
|
779 ComplexColumnVector::resize (int n) |
|
780 { |
|
781 if (n < 0) |
227
|
782 { |
|
783 (*current_liboctave_error_handler) |
|
784 ("can't resize to negative dimension"); |
|
785 return *this; |
|
786 } |
3
|
787 |
|
788 Complex *new_data = (Complex *) NULL; |
|
789 if (n > 0) |
|
790 { |
|
791 new_data = new Complex [n]; |
|
792 int min_len = len < n ? len : n; |
|
793 |
|
794 for (int i = 0; i < min_len; i++) |
|
795 new_data[i] = data[i]; |
|
796 } |
|
797 |
|
798 delete [] data; |
|
799 len = n; |
|
800 data = new_data; |
|
801 |
|
802 return *this; |
|
803 } |
|
804 |
|
805 ComplexColumnVector& |
|
806 ComplexColumnVector::resize (int n, double val) |
|
807 { |
|
808 int old_len = len; |
|
809 resize (n); |
|
810 for (int i = old_len; i < len; i++) |
|
811 data[i] = val; |
|
812 |
|
813 return *this; |
|
814 } |
|
815 |
|
816 ComplexColumnVector& |
161
|
817 ComplexColumnVector::resize (int n, const Complex& val) |
3
|
818 { |
|
819 int old_len = len; |
|
820 resize (n); |
|
821 for (int i = old_len; i < len; i++) |
|
822 data[i] = val; |
|
823 |
|
824 return *this; |
|
825 } |
|
826 |
|
827 int |
|
828 ComplexColumnVector::operator == (const ComplexColumnVector& a) const |
|
829 { |
|
830 if (len != a.len) |
|
831 return 0; |
|
832 return equal (data, a.data, len); |
|
833 } |
|
834 |
|
835 int |
|
836 ComplexColumnVector::operator != (const ComplexColumnVector& a) const |
|
837 { |
|
838 if (len != a.len) |
|
839 return 0; |
|
840 return !equal (data, a.data, len); |
|
841 } |
|
842 |
|
843 // destructive insert/delete/reorder operations |
|
844 |
|
845 ComplexColumnVector& |
|
846 ComplexColumnVector::insert (const ColumnVector& a, int r) |
|
847 { |
|
848 if (r < 0 || r + a.len - 1 > len) |
227
|
849 { |
|
850 (*current_liboctave_error_handler) ("range error for insert"); |
|
851 return *this; |
|
852 } |
3
|
853 |
|
854 for (int i = 0; i < a.len; i++) |
|
855 data[r+i] = a.data[i]; |
|
856 |
|
857 return *this; |
|
858 } |
|
859 |
|
860 ComplexColumnVector& |
|
861 ComplexColumnVector::insert (const ComplexColumnVector& a, int r) |
|
862 { |
|
863 if (r < 0 || r + a.len - 1 > len) |
227
|
864 { |
|
865 (*current_liboctave_error_handler) ("range error for insert"); |
|
866 return *this; |
|
867 } |
3
|
868 |
|
869 for (int i = 0; i < a.len; i++) |
|
870 data[r+i] = a.data[i]; |
|
871 |
|
872 return *this; |
|
873 } |
|
874 |
|
875 ComplexColumnVector& |
|
876 ComplexColumnVector::fill (double val) |
|
877 { |
|
878 if (len > 0) |
|
879 copy (data, len, val); |
|
880 return *this; |
|
881 } |
|
882 |
|
883 ComplexColumnVector& |
161
|
884 ComplexColumnVector::fill (const Complex& val) |
3
|
885 { |
|
886 if (len > 0) |
|
887 copy (data, len, val); |
|
888 return *this; |
|
889 } |
|
890 |
|
891 ComplexColumnVector& |
|
892 ComplexColumnVector::fill (double val, int r1, int r2) |
|
893 { |
|
894 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) |
227
|
895 { |
|
896 (*current_liboctave_error_handler) ("range error for fill"); |
|
897 return *this; |
|
898 } |
3
|
899 |
|
900 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
901 |
|
902 for (int i = r1; i <= r2; i++) |
|
903 data[i] = val; |
|
904 |
|
905 return *this; |
|
906 } |
|
907 |
|
908 ComplexColumnVector& |
161
|
909 ComplexColumnVector::fill (const Complex& val, int r1, int r2) |
3
|
910 { |
|
911 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) |
227
|
912 { |
|
913 (*current_liboctave_error_handler) ("range error for fill"); |
|
914 return *this; |
|
915 } |
3
|
916 |
|
917 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
918 |
|
919 for (int i = r1; i <= r2; i++) |
|
920 data[i] = val; |
|
921 |
|
922 return *this; |
|
923 } |
|
924 |
|
925 ComplexColumnVector |
|
926 ComplexColumnVector::stack (const ColumnVector& a) const |
|
927 { |
|
928 int nr_insert = len; |
|
929 ComplexColumnVector retval (len + a.len); |
|
930 retval.insert (*this, 0); |
|
931 retval.insert (a, nr_insert); |
|
932 return retval; |
|
933 } |
|
934 |
|
935 ComplexColumnVector |
|
936 ComplexColumnVector::stack (const ComplexColumnVector& a) const |
|
937 { |
|
938 int nr_insert = len; |
|
939 ComplexColumnVector retval (len + a.len); |
|
940 retval.insert (*this, 0); |
|
941 retval.insert (a, nr_insert); |
|
942 return retval; |
|
943 } |
|
944 |
|
945 ComplexRowVector |
|
946 ComplexColumnVector::hermitian (void) const |
|
947 { |
|
948 return ComplexRowVector (conj_dup (data, len), len); |
|
949 } |
|
950 |
|
951 ComplexRowVector |
|
952 ComplexColumnVector::transpose (void) const |
|
953 { |
|
954 return ComplexRowVector (dup (data, len), len); |
|
955 } |
|
956 |
|
957 ColumnVector |
|
958 real (const ComplexColumnVector& a) |
|
959 { |
|
960 ColumnVector retval; |
|
961 if (a.len > 0) |
|
962 retval = ColumnVector (real_dup (a.data, a.len), a.len); |
|
963 return retval; |
|
964 } |
|
965 |
|
966 ColumnVector |
|
967 imag (const ComplexColumnVector& a) |
|
968 { |
|
969 ColumnVector retval; |
|
970 if (a.len > 0) |
|
971 retval = ColumnVector (imag_dup (a.data, a.len), a.len); |
|
972 return retval; |
|
973 } |
|
974 |
|
975 ComplexColumnVector |
|
976 conj (const ComplexColumnVector& a) |
|
977 { |
|
978 ComplexColumnVector retval; |
|
979 if (a.len > 0) |
|
980 retval = ComplexColumnVector (conj_dup (a.data, a.len), a.len); |
|
981 return retval; |
|
982 } |
|
983 |
|
984 // resize is the destructive equivalent for this one |
|
985 |
|
986 ComplexColumnVector |
|
987 ComplexColumnVector::extract (int r1, int r2) const |
|
988 { |
|
989 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } |
|
990 |
|
991 int new_r = r2 - r1 + 1; |
|
992 |
|
993 ComplexColumnVector result (new_r); |
|
994 |
|
995 for (int i = 0; i < new_r; i++) |
|
996 result.data[i] = elem (r1+i); |
|
997 |
|
998 return result; |
|
999 } |
|
1000 |
|
1001 // column vector by scalar -> column vector operations |
|
1002 |
|
1003 ComplexColumnVector |
|
1004 ComplexColumnVector::operator + (double s) const |
|
1005 { |
|
1006 return ComplexColumnVector (add (data, len, s), len); |
|
1007 } |
|
1008 |
|
1009 ComplexColumnVector |
|
1010 ComplexColumnVector::operator - (double s) const |
|
1011 { |
|
1012 return ComplexColumnVector (subtract (data, len, s), len); |
|
1013 } |
|
1014 |
|
1015 ComplexColumnVector |
|
1016 ComplexColumnVector::operator * (double s) const |
|
1017 { |
|
1018 return ComplexColumnVector (multiply (data, len, s), len); |
|
1019 } |
|
1020 |
|
1021 ComplexColumnVector |
|
1022 ComplexColumnVector::operator / (double s) const |
|
1023 { |
|
1024 return ComplexColumnVector (divide (data, len, s), len); |
|
1025 } |
|
1026 |
|
1027 ComplexColumnVector |
161
|
1028 ComplexColumnVector::operator + (const Complex& s) const |
3
|
1029 { |
|
1030 return ComplexColumnVector (add (data, len, s), len); |
|
1031 } |
|
1032 |
|
1033 ComplexColumnVector |
161
|
1034 ComplexColumnVector::operator - (const Complex& s) const |
3
|
1035 { |
|
1036 return ComplexColumnVector (subtract (data, len, s), len); |
|
1037 } |
|
1038 |
|
1039 ComplexColumnVector |
161
|
1040 ComplexColumnVector::operator * (const Complex& s) const |
3
|
1041 { |
|
1042 return ComplexColumnVector (multiply (data, len, s), len); |
|
1043 } |
|
1044 |
|
1045 ComplexColumnVector |
161
|
1046 ComplexColumnVector::operator / (const Complex& s) const |
3
|
1047 { |
|
1048 return ComplexColumnVector (divide (data, len, s), len); |
|
1049 } |
|
1050 |
|
1051 // scalar by column vector -> column vector operations |
|
1052 |
|
1053 ComplexColumnVector |
|
1054 operator + (double s, const ComplexColumnVector& a) |
|
1055 { |
|
1056 return ComplexColumnVector (add (a.data, a.len, s), a.len); |
|
1057 } |
|
1058 |
|
1059 ComplexColumnVector |
|
1060 operator - (double s, const ComplexColumnVector& a) |
|
1061 { |
|
1062 return ComplexColumnVector (subtract (s, a.data, a.len), a.len); |
|
1063 } |
|
1064 |
|
1065 ComplexColumnVector |
|
1066 operator * (double s, const ComplexColumnVector& a) |
|
1067 { |
|
1068 return ComplexColumnVector (multiply (a.data, a.len, s), a.len); |
|
1069 } |
|
1070 |
|
1071 ComplexColumnVector |
|
1072 operator / (double s, const ComplexColumnVector& a) |
|
1073 { |
|
1074 return ComplexColumnVector (divide (s, a.data, a.len), a.len); |
|
1075 } |
|
1076 |
|
1077 ComplexColumnVector |
161
|
1078 operator + (const Complex& s, const ComplexColumnVector& a) |
3
|
1079 { |
|
1080 return ComplexColumnVector (add (a.data, a.len, s), a.len); |
|
1081 } |
|
1082 |
|
1083 ComplexColumnVector |
161
|
1084 operator - (const Complex& s, const ComplexColumnVector& a) |
3
|
1085 { |
|
1086 return ComplexColumnVector (subtract (s, a.data, a.len), a.len); |
|
1087 } |
|
1088 |
|
1089 ComplexColumnVector |
161
|
1090 operator * (const Complex& s, const ComplexColumnVector& a) |
3
|
1091 { |
|
1092 return ComplexColumnVector (multiply (a.data, a.len, s), a.len); |
|
1093 } |
|
1094 |
|
1095 ComplexColumnVector |
161
|
1096 operator / (const Complex& s, const ComplexColumnVector& a) |
3
|
1097 { |
|
1098 return ComplexColumnVector (divide (s, a.data, a.len), a.len); |
|
1099 } |
|
1100 |
|
1101 // column vector by row vector -> matrix operations |
|
1102 |
|
1103 ComplexMatrix |
|
1104 ComplexColumnVector::operator * (const RowVector& a) const |
|
1105 { |
|
1106 ComplexRowVector tmp (a); |
|
1107 return *this * tmp; |
|
1108 } |
|
1109 |
|
1110 ComplexMatrix |
|
1111 ComplexColumnVector::operator * (const ComplexRowVector& a) const |
|
1112 { |
|
1113 if (len != a.len) |
227
|
1114 { |
|
1115 (*current_liboctave_error_handler) |
|
1116 ("nonconformant vector multiplication attempted"); |
|
1117 return ComplexMatrix (); |
|
1118 } |
3
|
1119 |
|
1120 if (len == 0) |
|
1121 return ComplexMatrix (len, len, 0.0); |
|
1122 |
|
1123 char transa = 'N'; |
|
1124 char transb = 'N'; |
|
1125 Complex alpha (1.0); |
|
1126 Complex beta (0.0); |
|
1127 int anr = 1; |
|
1128 int anc = a.len; |
|
1129 |
|
1130 Complex *c = new Complex [len * a.len]; |
|
1131 |
|
1132 F77_FCN (zgemm) (&transa, &transb, &len, &anc, &anr, &alpha, data, |
|
1133 &len, a.data, &anr, &beta, c, &len, 1L, 1L); |
|
1134 |
|
1135 return ComplexMatrix (c, len, a.len); |
|
1136 } |
|
1137 |
|
1138 // column vector by column vector -> column vector operations |
|
1139 |
|
1140 ComplexColumnVector |
|
1141 ComplexColumnVector::operator + (const ColumnVector& a) const |
|
1142 { |
|
1143 if (len != a.len) |
227
|
1144 { |
|
1145 (*current_liboctave_error_handler) |
|
1146 ("nonconformant vector addition attempted"); |
|
1147 return ComplexColumnVector (); |
|
1148 } |
3
|
1149 |
|
1150 if (len == 0) |
|
1151 return ComplexColumnVector (0); |
|
1152 |
|
1153 return ComplexColumnVector (add (data, a.data, len), len); |
|
1154 } |
|
1155 |
|
1156 ComplexColumnVector |
|
1157 ComplexColumnVector::operator - (const ColumnVector& a) const |
|
1158 { |
|
1159 if (len != a.len) |
227
|
1160 { |
|
1161 (*current_liboctave_error_handler) |
|
1162 ("nonconformant vector subtraction attempted"); |
|
1163 return ComplexColumnVector (); |
|
1164 } |
3
|
1165 |
|
1166 if (len == 0) |
|
1167 return ComplexColumnVector (0); |
|
1168 |
|
1169 return ComplexColumnVector (subtract (data, a.data, len), len); |
|
1170 } |
|
1171 |
|
1172 ComplexColumnVector |
|
1173 ComplexColumnVector::operator + (const ComplexColumnVector& a) const |
|
1174 { |
|
1175 if (len != a.len) |
227
|
1176 { |
|
1177 (*current_liboctave_error_handler) |
|
1178 ("nonconformant vector addition attempted"); |
|
1179 return ComplexColumnVector (); |
|
1180 } |
3
|
1181 |
|
1182 if (len == 0) |
|
1183 return ComplexColumnVector (0); |
|
1184 |
|
1185 return ComplexColumnVector (add (data, a.data, len), len); |
|
1186 } |
|
1187 |
|
1188 ComplexColumnVector |
|
1189 ComplexColumnVector::operator - (const ComplexColumnVector& a) const |
|
1190 { |
|
1191 if (len != a.len) |
227
|
1192 { |
|
1193 (*current_liboctave_error_handler) |
|
1194 ("nonconformant vector subtraction attempted"); |
|
1195 return ComplexColumnVector (); |
|
1196 } |
3
|
1197 |
|
1198 if (len == 0) |
|
1199 return ComplexColumnVector (0); |
|
1200 |
|
1201 return ComplexColumnVector (subtract (data, a.data, len), len); |
|
1202 } |
|
1203 |
|
1204 ComplexColumnVector |
|
1205 ComplexColumnVector::product (const ColumnVector& a) const |
|
1206 { |
|
1207 if (len != a.len) |
227
|
1208 { |
|
1209 (*current_liboctave_error_handler) |
|
1210 ("nonconformant vector product attempted"); |
|
1211 return ComplexColumnVector (); |
|
1212 } |
3
|
1213 |
|
1214 if (len == 0) |
|
1215 return ComplexColumnVector (0); |
|
1216 |
|
1217 return ComplexColumnVector (multiply (data, a.data, len), len); |
|
1218 } |
|
1219 |
|
1220 ComplexColumnVector |
|
1221 ComplexColumnVector::quotient (const ColumnVector& a) const |
|
1222 { |
|
1223 if (len != a.len) |
227
|
1224 { |
|
1225 (*current_liboctave_error_handler) |
|
1226 ("nonconformant vector quotient attempted"); |
|
1227 return ComplexColumnVector (); |
|
1228 } |
3
|
1229 |
|
1230 if (len == 0) |
|
1231 return ComplexColumnVector (0); |
|
1232 |
|
1233 return ComplexColumnVector (divide (data, a.data, len), len); |
|
1234 } |
|
1235 |
|
1236 ComplexColumnVector |
|
1237 ComplexColumnVector::product (const ComplexColumnVector& a) const |
|
1238 { |
|
1239 if (len != a.len) |
227
|
1240 { |
|
1241 (*current_liboctave_error_handler) |
|
1242 ("nonconformant vector product attempted"); |
|
1243 return ComplexColumnVector (); |
|
1244 } |
3
|
1245 |
|
1246 if (len == 0) |
|
1247 return ComplexColumnVector (0); |
|
1248 |
|
1249 return ComplexColumnVector (multiply (data, a.data, len), len); |
|
1250 } |
|
1251 |
|
1252 ComplexColumnVector |
|
1253 ComplexColumnVector::quotient (const ComplexColumnVector& a) const |
|
1254 { |
|
1255 if (len != a.len) |
227
|
1256 { |
|
1257 (*current_liboctave_error_handler) |
|
1258 ("nonconformant vector quotient attempted"); |
|
1259 return ComplexColumnVector (); |
|
1260 } |
3
|
1261 |
|
1262 if (len == 0) |
|
1263 return ComplexColumnVector (0); |
|
1264 |
|
1265 return ComplexColumnVector (divide (data, a.data, len), len); |
|
1266 } |
|
1267 |
|
1268 ComplexColumnVector& |
|
1269 ComplexColumnVector::operator += (const ColumnVector& a) |
|
1270 { |
|
1271 if (len != a.len) |
227
|
1272 { |
|
1273 (*current_liboctave_error_handler) |
|
1274 ("nonconformant vector += operation attempted"); |
|
1275 return *this; |
|
1276 } |
3
|
1277 |
|
1278 if (len == 0) |
|
1279 return *this; |
|
1280 |
|
1281 add2 (data, a.data, len); |
|
1282 return *this; |
|
1283 } |
|
1284 |
|
1285 ComplexColumnVector& |
|
1286 ComplexColumnVector::operator -= (const ColumnVector& a) |
|
1287 { |
|
1288 if (len != a.len) |
227
|
1289 { |
|
1290 (*current_liboctave_error_handler) |
|
1291 ("nonconformant vector -= operation attempted"); |
|
1292 return *this; |
|
1293 } |
3
|
1294 |
|
1295 if (len == 0) |
|
1296 return *this; |
|
1297 |
|
1298 subtract2 (data, a.data, len); |
|
1299 return *this; |
|
1300 } |
|
1301 |
|
1302 ComplexColumnVector& |
|
1303 ComplexColumnVector::operator += (const ComplexColumnVector& a) |
|
1304 { |
|
1305 if (len != a.len) |
227
|
1306 { |
|
1307 (*current_liboctave_error_handler) |
|
1308 ("nonconformant vector += operation attempted"); |
|
1309 return *this; |
|
1310 } |
3
|
1311 |
|
1312 if (len == 0) |
|
1313 return *this; |
|
1314 |
|
1315 add2 (data, a.data, len); |
|
1316 return *this; |
|
1317 } |
|
1318 |
|
1319 ComplexColumnVector& |
|
1320 ComplexColumnVector::operator -= (const ComplexColumnVector& a) |
|
1321 { |
|
1322 if (len != a.len) |
227
|
1323 { |
|
1324 (*current_liboctave_error_handler) |
|
1325 ("nonconformant vector -= operation attempted"); |
|
1326 return *this; |
|
1327 } |
3
|
1328 |
|
1329 if (len == 0) |
|
1330 return *this; |
|
1331 |
|
1332 subtract2 (data, a.data, len); |
|
1333 return *this; |
|
1334 } |
|
1335 |
|
1336 // unary operations |
|
1337 |
|
1338 ComplexColumnVector |
|
1339 ComplexColumnVector::operator - (void) const |
|
1340 { |
|
1341 if (len == 0) |
|
1342 return ComplexColumnVector (0); |
|
1343 |
|
1344 return ComplexColumnVector (negate (data, len), len); |
|
1345 } |
|
1346 |
|
1347 ComplexColumnVector |
|
1348 map (c_c_Mapper f, const ComplexColumnVector& a) |
|
1349 { |
|
1350 ComplexColumnVector b (a); |
|
1351 b.map (f); |
|
1352 return b; |
|
1353 } |
|
1354 |
|
1355 ColumnVector |
|
1356 map (d_c_Mapper f, const ComplexColumnVector& a) |
|
1357 { |
|
1358 ColumnVector b (a.len); |
|
1359 for (int i = 0; i < a.len; i++) |
|
1360 b.elem (i) = f (a.elem (i)); |
|
1361 return b; |
|
1362 } |
|
1363 |
|
1364 void |
|
1365 ComplexColumnVector::map (c_c_Mapper f) |
|
1366 { |
|
1367 for (int i = 0; i < len; i++) |
|
1368 data[i] = f (data[i]); |
|
1369 } |
|
1370 |
|
1371 Complex |
|
1372 ComplexColumnVector::min (void) const |
|
1373 { |
|
1374 if (len == 0) |
|
1375 return 0.0; |
|
1376 |
|
1377 Complex res = data[0]; |
|
1378 double absres = abs (res); |
|
1379 |
|
1380 for (int i = 1; i < len; i++) |
|
1381 if (abs (data[i]) < absres) |
|
1382 { |
|
1383 res = data[i]; |
|
1384 absres = abs (res); |
|
1385 } |
|
1386 |
|
1387 return res; |
|
1388 } |
|
1389 |
|
1390 Complex |
|
1391 ComplexColumnVector::max (void) const |
|
1392 { |
|
1393 if (len == 0) |
|
1394 return 0.0; |
|
1395 |
|
1396 Complex res = data[0]; |
|
1397 double absres = abs (res); |
|
1398 |
|
1399 for (int i = 1; i < len; i++) |
|
1400 if (abs (data[i]) > absres) |
|
1401 { |
|
1402 res = data[i]; |
|
1403 absres = abs (res); |
|
1404 } |
|
1405 |
|
1406 return res; |
|
1407 } |
|
1408 |
|
1409 // i/o |
|
1410 |
|
1411 ostream& |
|
1412 operator << (ostream& os, const ComplexColumnVector& a) |
|
1413 { |
|
1414 // int field_width = os.precision () + 7; |
|
1415 for (int i = 0; i < a.len; i++) |
|
1416 os << /* setw (field_width) << */ a.data[i] << "\n"; |
|
1417 return os; |
|
1418 } |
|
1419 |
|
1420 /* |
|
1421 ;;; Local Variables: *** |
|
1422 ;;; mode: C++ *** |
|
1423 ;;; page-delimiter: "^/\\*" *** |
|
1424 ;;; End: *** |
|
1425 */ |