1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include <istream>
31 #include <ostream>
32
33 #include "Array-util.h"
34 #include "lo-blas-proto.h"
35 #include "lo-error.h"
36 #include "mx-base.h"
37 #include "mx-inlines.cc"
38 #include "oct-cmplx.h"
39
40 // FloatComplex Column Vector class
41
FloatComplexColumnVector(const FloatColumnVector & a)42 FloatComplexColumnVector::FloatComplexColumnVector (const FloatColumnVector& a)
43 : MArray<FloatComplex> (a)
44 { }
45
46 bool
operator ==(const FloatComplexColumnVector & a) const47 FloatComplexColumnVector::operator == (const FloatComplexColumnVector& a) const
48 {
49 octave_idx_type len = numel ();
50 if (len != a.numel ())
51 return 0;
52 return mx_inline_equal (len, data (), a.data ());
53 }
54
55 bool
operator !=(const FloatComplexColumnVector & a) const56 FloatComplexColumnVector::operator != (const FloatComplexColumnVector& a) const
57 {
58 return !(*this == a);
59 }
60
61 // destructive insert/delete/reorder operations
62
63 FloatComplexColumnVector&
insert(const FloatColumnVector & a,octave_idx_type r)64 FloatComplexColumnVector::insert (const FloatColumnVector& a, octave_idx_type r)
65 {
66 octave_idx_type a_len = a.numel ();
67
68 if (r < 0 || r + a_len > numel ())
69 (*current_liboctave_error_handler) ("range error for insert");
70
71 if (a_len > 0)
72 {
73 make_unique ();
74
75 for (octave_idx_type i = 0; i < a_len; i++)
76 xelem (r+i) = a.elem (i);
77 }
78
79 return *this;
80 }
81
82 FloatComplexColumnVector&
insert(const FloatComplexColumnVector & a,octave_idx_type r)83 FloatComplexColumnVector::insert (const FloatComplexColumnVector& a,
84 octave_idx_type r)
85 {
86 octave_idx_type a_len = a.numel ();
87
88 if (r < 0 || r + a_len > numel ())
89 (*current_liboctave_error_handler) ("range error for insert");
90
91 if (a_len > 0)
92 {
93 make_unique ();
94
95 for (octave_idx_type i = 0; i < a_len; i++)
96 xelem (r+i) = a.elem (i);
97 }
98
99 return *this;
100 }
101
102 FloatComplexColumnVector&
fill(float val)103 FloatComplexColumnVector::fill (float val)
104 {
105 octave_idx_type len = numel ();
106
107 if (len > 0)
108 {
109 make_unique ();
110
111 for (octave_idx_type i = 0; i < len; i++)
112 xelem (i) = val;
113 }
114
115 return *this;
116 }
117
118 FloatComplexColumnVector&
fill(const FloatComplex & val)119 FloatComplexColumnVector::fill (const FloatComplex& val)
120 {
121 octave_idx_type len = numel ();
122
123 if (len > 0)
124 {
125 make_unique ();
126
127 for (octave_idx_type i = 0; i < len; i++)
128 xelem (i) = val;
129 }
130
131 return *this;
132 }
133
134 FloatComplexColumnVector&
fill(float val,octave_idx_type r1,octave_idx_type r2)135 FloatComplexColumnVector::fill (float val,
136 octave_idx_type r1, octave_idx_type r2)
137 {
138 octave_idx_type len = numel ();
139
140 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
141 (*current_liboctave_error_handler) ("range error for fill");
142
143 if (r1 > r2) { std::swap (r1, r2); }
144
145 if (r2 >= r1)
146 {
147 make_unique ();
148
149 for (octave_idx_type i = r1; i <= r2; i++)
150 xelem (i) = val;
151 }
152
153 return *this;
154 }
155
156 FloatComplexColumnVector&
fill(const FloatComplex & val,octave_idx_type r1,octave_idx_type r2)157 FloatComplexColumnVector::fill (const FloatComplex& val,
158 octave_idx_type r1, octave_idx_type r2)
159 {
160 octave_idx_type len = numel ();
161
162 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
163 (*current_liboctave_error_handler) ("range error for fill");
164
165 if (r1 > r2) { std::swap (r1, r2); }
166
167 if (r2 >= r1)
168 {
169 make_unique ();
170
171 for (octave_idx_type i = r1; i <= r2; i++)
172 xelem (i) = val;
173 }
174
175 return *this;
176 }
177
178 FloatComplexColumnVector
stack(const FloatColumnVector & a) const179 FloatComplexColumnVector::stack (const FloatColumnVector& a) const
180 {
181 octave_idx_type len = numel ();
182 octave_idx_type nr_insert = len;
183 FloatComplexColumnVector retval (len + a.numel ());
184 retval.insert (*this, 0);
185 retval.insert (a, nr_insert);
186 return retval;
187 }
188
189 FloatComplexColumnVector
stack(const FloatComplexColumnVector & a) const190 FloatComplexColumnVector::stack (const FloatComplexColumnVector& a) const
191 {
192 octave_idx_type len = numel ();
193 octave_idx_type nr_insert = len;
194 FloatComplexColumnVector retval (len + a.numel ());
195 retval.insert (*this, 0);
196 retval.insert (a, nr_insert);
197 return retval;
198 }
199
200 FloatComplexRowVector
hermitian(void) const201 FloatComplexColumnVector::hermitian (void) const
202 {
203 return MArray<FloatComplex>::hermitian (std::conj);
204 }
205
206 FloatComplexRowVector
transpose(void) const207 FloatComplexColumnVector::transpose (void) const
208 {
209 return MArray<FloatComplex>::transpose ();
210 }
211
212 FloatColumnVector
abs(void) const213 FloatComplexColumnVector::abs (void) const
214 {
215 return do_mx_unary_map<float, FloatComplex, std::abs> (*this);
216 }
217
218 FloatComplexColumnVector
conj(const FloatComplexColumnVector & a)219 conj (const FloatComplexColumnVector& a)
220 {
221 return do_mx_unary_map<FloatComplex, FloatComplex, std::conj<float>> (a);
222 }
223
224 // resize is the destructive equivalent for this one
225
226 FloatComplexColumnVector
extract(octave_idx_type r1,octave_idx_type r2) const227 FloatComplexColumnVector::extract (octave_idx_type r1, octave_idx_type r2) const
228 {
229 if (r1 > r2) { std::swap (r1, r2); }
230
231 octave_idx_type new_r = r2 - r1 + 1;
232
233 FloatComplexColumnVector result (new_r);
234
235 for (octave_idx_type i = 0; i < new_r; i++)
236 result.elem (i) = elem (r1+i);
237
238 return result;
239 }
240
241 FloatComplexColumnVector
extract_n(octave_idx_type r1,octave_idx_type n) const242 FloatComplexColumnVector::extract_n (octave_idx_type r1,
243 octave_idx_type n) const
244 {
245 FloatComplexColumnVector result (n);
246
247 for (octave_idx_type i = 0; i < n; i++)
248 result.elem (i) = elem (r1+i);
249
250 return result;
251 }
252
253 // column vector by column vector -> column vector operations
254
255 FloatComplexColumnVector&
operator +=(const FloatColumnVector & a)256 FloatComplexColumnVector::operator += (const FloatColumnVector& a)
257 {
258 octave_idx_type len = numel ();
259
260 octave_idx_type a_len = a.numel ();
261
262 if (len != a_len)
263 octave::err_nonconformant ("operator +=", len, a_len);
264
265 if (len == 0)
266 return *this;
267
268 FloatComplex *d = fortran_vec (); // Ensures only 1 reference to my privates!
269
270 mx_inline_add2 (len, d, a.data ());
271 return *this;
272 }
273
274 FloatComplexColumnVector&
operator -=(const FloatColumnVector & a)275 FloatComplexColumnVector::operator -= (const FloatColumnVector& a)
276 {
277 octave_idx_type len = numel ();
278
279 octave_idx_type a_len = a.numel ();
280
281 if (len != a_len)
282 octave::err_nonconformant ("operator -=", len, a_len);
283
284 if (len == 0)
285 return *this;
286
287 FloatComplex *d = fortran_vec (); // Ensures only 1 reference to my privates!
288
289 mx_inline_sub2 (len, d, a.data ());
290 return *this;
291 }
292
293 // matrix by column vector -> column vector operations
294
295 FloatComplexColumnVector
operator *(const FloatComplexMatrix & m,const FloatColumnVector & a)296 operator * (const FloatComplexMatrix& m, const FloatColumnVector& a)
297 {
298 FloatComplexColumnVector tmp (a);
299 return m * tmp;
300 }
301
302 FloatComplexColumnVector
operator *(const FloatComplexMatrix & m,const FloatComplexColumnVector & a)303 operator * (const FloatComplexMatrix& m, const FloatComplexColumnVector& a)
304 {
305 FloatComplexColumnVector retval;
306
307 F77_INT nr = octave::to_f77_int (m.rows ());
308 F77_INT nc = octave::to_f77_int (m.cols ());
309
310 F77_INT a_len = octave::to_f77_int (a.numel ());
311
312 if (nc != a_len)
313 octave::err_nonconformant ("operator *", nr, nc, a_len, 1);
314
315 retval.clear (nr);
316
317 if (nr != 0)
318 {
319 if (nc == 0)
320 retval.fill (0.0);
321 else
322 {
323 FloatComplex *y = retval.fortran_vec ();
324
325 F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
326 nr, nc, 1.0f, F77_CONST_CMPLX_ARG (m.data ()), nr,
327 F77_CONST_CMPLX_ARG (a.data ()), 1, 0.0f, F77_CMPLX_ARG (y), 1
328 F77_CHAR_ARG_LEN (1)));
329 }
330 }
331
332 return retval;
333 }
334
335 // matrix by column vector -> column vector operations
336
337 FloatComplexColumnVector
operator *(const FloatMatrix & m,const FloatComplexColumnVector & a)338 operator * (const FloatMatrix& m, const FloatComplexColumnVector& a)
339 {
340 FloatComplexMatrix tmp (m);
341 return tmp * a;
342 }
343
344 // diagonal matrix by column vector -> column vector operations
345
346 FloatComplexColumnVector
operator *(const FloatDiagMatrix & m,const FloatComplexColumnVector & a)347 operator * (const FloatDiagMatrix& m, const FloatComplexColumnVector& a)
348 {
349 F77_INT nr = octave::to_f77_int (m.rows ());
350 F77_INT nc = octave::to_f77_int (m.cols ());
351
352 F77_INT a_len = octave::to_f77_int (a.numel ());
353
354 if (nc != a_len)
355 octave::err_nonconformant ("operator *", nr, nc, a_len, 1);
356
357 if (nc == 0 || nr == 0)
358 return FloatComplexColumnVector (0);
359
360 FloatComplexColumnVector result (nr);
361
362 for (octave_idx_type i = 0; i < a_len; i++)
363 result.elem (i) = a.elem (i) * m.elem (i, i);
364
365 for (octave_idx_type i = a_len; i < nr; i++)
366 result.elem (i) = 0.0;
367
368 return result;
369 }
370
371 FloatComplexColumnVector
operator *(const FloatComplexDiagMatrix & m,const FloatColumnVector & a)372 operator * (const FloatComplexDiagMatrix& m, const FloatColumnVector& a)
373 {
374 F77_INT nr = octave::to_f77_int (m.rows ());
375 F77_INT nc = octave::to_f77_int (m.cols ());
376
377 F77_INT a_len = octave::to_f77_int (a.numel ());
378
379 if (nc != a_len)
380 octave::err_nonconformant ("operator *", nr, nc, a_len, 1);
381
382 if (nc == 0 || nr == 0)
383 return FloatComplexColumnVector (0);
384
385 FloatComplexColumnVector result (nr);
386
387 for (octave_idx_type i = 0; i < a_len; i++)
388 result.elem (i) = a.elem (i) * m.elem (i, i);
389
390 for (octave_idx_type i = a_len; i < nr; i++)
391 result.elem (i) = 0.0;
392
393 return result;
394 }
395
396 FloatComplexColumnVector
operator *(const FloatComplexDiagMatrix & m,const FloatComplexColumnVector & a)397 operator * (const FloatComplexDiagMatrix& m, const FloatComplexColumnVector& a)
398 {
399 F77_INT nr = octave::to_f77_int (m.rows ());
400 F77_INT nc = octave::to_f77_int (m.cols ());
401
402 F77_INT a_len = octave::to_f77_int (a.numel ());
403
404 if (nc != a_len)
405 octave::err_nonconformant ("operator *", nr, nc, a_len, 1);
406
407 if (nc == 0 || nr == 0)
408 return FloatComplexColumnVector (0);
409
410 FloatComplexColumnVector result (nr);
411
412 for (octave_idx_type i = 0; i < a_len; i++)
413 result.elem (i) = a.elem (i) * m.elem (i, i);
414
415 for (octave_idx_type i = a_len; i < nr; i++)
416 result.elem (i) = 0.0;
417
418 return result;
419 }
420
421 // other operations
422
423 FloatComplex
min(void) const424 FloatComplexColumnVector::min (void) const
425 {
426 octave_idx_type len = numel ();
427 if (len == 0)
428 return 0.0;
429
430 FloatComplex res = elem (0);
431 float absres = std::abs (res);
432
433 for (octave_idx_type i = 1; i < len; i++)
434 if (std::abs (elem (i)) < absres)
435 {
436 res = elem (i);
437 absres = std::abs (res);
438 }
439
440 return res;
441 }
442
443 FloatComplex
max(void) const444 FloatComplexColumnVector::max (void) const
445 {
446 octave_idx_type len = numel ();
447 if (len == 0)
448 return 0.0;
449
450 FloatComplex res = elem (0);
451 float absres = std::abs (res);
452
453 for (octave_idx_type i = 1; i < len; i++)
454 if (std::abs (elem (i)) > absres)
455 {
456 res = elem (i);
457 absres = std::abs (res);
458 }
459
460 return res;
461 }
462
463 // i/o
464
465 std::ostream&
operator <<(std::ostream & os,const FloatComplexColumnVector & a)466 operator << (std::ostream& os, const FloatComplexColumnVector& a)
467 {
468 // int field_width = os.precision () + 7;
469 for (octave_idx_type i = 0; i < a.numel (); i++)
470 os << /* setw (field_width) << */ a.elem (i) << "\n";
471 return os;
472 }
473
474 std::istream&
operator >>(std::istream & is,FloatComplexColumnVector & a)475 operator >> (std::istream& is, FloatComplexColumnVector& a)
476 {
477 octave_idx_type len = a.numel ();
478
479 if (len > 0)
480 {
481 float tmp;
482 for (octave_idx_type i = 0; i < len; i++)
483 {
484 is >> tmp;
485 if (is)
486 a.elem (i) = tmp;
487 else
488 break;
489 }
490 }
491 return is;
492 }
493