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