GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
fCMatrix.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1994-2025 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 <algorithm>
31#include <cmath>
32#include <complex>
33#include <istream>
34#include <limits>
35#include <ostream>
36
37#include "Array-util.h"
38#include "DET.h"
39#include "boolMatrix.h"
40#include "chMatrix.h"
41#include "chol.h"
42#include "fCColVector.h"
43#include "fCDiagMatrix.h"
44#include "fCMatrix.h"
45#include "fCNDArray.h"
46#include "fCRowVector.h"
47#include "lo-blas-proto.h"
48#include "lo-error.h"
49#include "lo-ieee.h"
50#include "lo-lapack-proto.h"
51#include "lo-mappers.h"
52#include "lo-utils.h"
53#include "mx-fcm-fdm.h"
54#include "mx-fcm-fs.h"
55#include "mx-fdm-fcm.h"
56#include "mx-inlines.cc"
57#include "mx-op-defs.h"
58#include "oct-cmplx.h"
59#include "oct-fftw.h"
60#include "oct-locbuf.h"
61#include "oct-norm.h"
62#include "schur.h"
63#include "svd.h"
64
65static const FloatComplex FloatComplex_NaN_result (octave::numeric_limits<float>::NaN (),
66 octave::numeric_limits<float>::NaN ());
67
68// FloatComplex Matrix class
69
73
77
81
82FloatComplexMatrix::FloatComplexMatrix (const FloatDiagMatrix& a)
83 : FloatComplexNDArray (a.dims (), 0.0)
84{
85 for (octave_idx_type i = 0; i < a.length (); i++)
86 elem (i, i) = a.elem (i, i);
87}
88
89FloatComplexMatrix::FloatComplexMatrix (const MDiagArray2<float>& a)
90 : FloatComplexNDArray (a.dims (), 0.0)
91{
92 for (octave_idx_type i = 0; i < a.length (); i++)
93 elem (i, i) = a.elem (i, i);
94}
95
96FloatComplexMatrix::FloatComplexMatrix (const DiagArray2<float>& a)
97 : FloatComplexNDArray (a.dims (), 0.0)
98{
99 for (octave_idx_type i = 0; i < a.length (); i++)
100 elem (i, i) = a.elem (i, i);
101}
102
106
110
111FloatComplexMatrix::FloatComplexMatrix (const FloatComplexDiagMatrix& a)
112 : FloatComplexNDArray (a.dims (), 0.0)
113{
114 for (octave_idx_type i = 0; i < a.length (); i++)
115 elem (i, i) = a.elem (i, i);
116}
117
118FloatComplexMatrix::FloatComplexMatrix (const MDiagArray2<FloatComplex>& a)
119 : FloatComplexNDArray (a.dims (), 0.0)
120{
121 for (octave_idx_type i = 0; i < a.length (); i++)
122 elem (i, i) = a.elem (i, i);
123}
124
125FloatComplexMatrix::FloatComplexMatrix (const DiagArray2<FloatComplex>& a)
126 : FloatComplexNDArray (a.dims (), 0.0)
127{
128 for (octave_idx_type i = 0; i < a.length (); i++)
129 elem (i, i) = a.elem (i, i);
130}
131
132// FIXME: could we use a templated mixed-type copy function
133// here?
134
138
139FloatComplexMatrix::FloatComplexMatrix (const charMatrix& a)
140 : FloatComplexNDArray (a.dims (), 0.0)
141{
142 for (octave_idx_type i = 0; i < a.rows (); i++)
143 for (octave_idx_type j = 0; j < a.cols (); j++)
144 elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
145}
146
147FloatComplexMatrix::FloatComplexMatrix (const FloatMatrix& re,
148 const FloatMatrix& im)
149 : FloatComplexNDArray (re.dims ())
150{
151 if (im.rows () != rows () || im.cols () != cols ())
152 (*current_liboctave_error_handler) ("complex: internal error");
153
154 octave_idx_type nel = numel ();
155 for (octave_idx_type i = 0; i < nel; i++)
156 xelem (i) = FloatComplex (re(i), im(i));
157}
158
159bool
160FloatComplexMatrix::operator == (const FloatComplexMatrix& a) const
161{
162 if (rows () != a.rows () || cols () != a.cols ())
163 return false;
164
165 return mx_inline_equal (numel (), data (), a.data ());
166}
167
168bool
169FloatComplexMatrix::operator != (const FloatComplexMatrix& a) const
170{
171 return !(*this == a);
172}
173
174bool
175FloatComplexMatrix::ishermitian () const
176{
177 octave_idx_type nr = rows ();
178 octave_idx_type nc = cols ();
179
180 if (issquare () && nr > 0)
181 {
182 for (octave_idx_type i = 0; i < nr; i++)
183 for (octave_idx_type j = i; j < nc; j++)
184 if (elem (i, j) != conj (elem (j, i)))
185 return false;
186
187 return true;
188 }
189
190 return false;
191}
192
193// destructive insert/delete/reorder operations
194
195FloatComplexMatrix&
196FloatComplexMatrix::insert (const FloatMatrix& a,
197 octave_idx_type r, octave_idx_type c)
198{
199 octave_idx_type a_nr = a.rows ();
200 octave_idx_type a_nc = a.cols ();
201
202 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
203 (*current_liboctave_error_handler) ("range error for insert");
204
205 if (a_nr >0 && a_nc > 0)
206 {
207 make_unique ();
208
209 for (octave_idx_type j = 0; j < a_nc; j++)
210 for (octave_idx_type i = 0; i < a_nr; i++)
211 xelem (r+i, c+j) = a.elem (i, j);
212 }
213
214 return *this;
215}
216
217FloatComplexMatrix&
218FloatComplexMatrix::insert (const FloatRowVector& a,
219 octave_idx_type r, octave_idx_type c)
220{
221 octave_idx_type a_len = a.numel ();
222
223 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
224 (*current_liboctave_error_handler) ("range error for insert");
225
226 if (a_len > 0)
227 {
228 make_unique ();
229
230 for (octave_idx_type i = 0; i < a_len; i++)
231 xelem (r, c+i) = a.elem (i);
232 }
233
234 return *this;
235}
236
237FloatComplexMatrix&
238FloatComplexMatrix::insert (const FloatColumnVector& a,
239 octave_idx_type r, octave_idx_type c)
240{
241 octave_idx_type a_len = a.numel ();
242
243 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
244 (*current_liboctave_error_handler) ("range error for insert");
245
246 if (a_len > 0)
247 {
248 make_unique ();
249
250 for (octave_idx_type i = 0; i < a_len; i++)
251 xelem (r+i, c) = a.elem (i);
252 }
253
254 return *this;
255}
256
257FloatComplexMatrix&
258FloatComplexMatrix::insert (const FloatDiagMatrix& a,
259 octave_idx_type r, octave_idx_type c)
260{
261 octave_idx_type a_nr = a.rows ();
262 octave_idx_type a_nc = a.cols ();
263
264 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
265 (*current_liboctave_error_handler) ("range error for insert");
266
267 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
268
269 octave_idx_type a_len = a.length ();
270
271 if (a_len > 0)
272 {
273 make_unique ();
274
275 for (octave_idx_type i = 0; i < a_len; i++)
276 xelem (r+i, c+i) = a.elem (i, i);
277 }
278
279 return *this;
280}
281
282FloatComplexMatrix&
289
290FloatComplexMatrix&
291FloatComplexMatrix::insert (const FloatComplexRowVector& a,
292 octave_idx_type r, octave_idx_type c)
293{
294 octave_idx_type a_len = a.numel ();
295 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
296 (*current_liboctave_error_handler) ("range error for insert");
297
298 for (octave_idx_type i = 0; i < a_len; i++)
299 elem (r, c+i) = a.elem (i);
300
301 return *this;
302}
303
304FloatComplexMatrix&
305FloatComplexMatrix::insert (const FloatComplexColumnVector& a,
306 octave_idx_type r, octave_idx_type c)
307{
308 octave_idx_type a_len = a.numel ();
309
310 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
311 (*current_liboctave_error_handler) ("range error for insert");
312
313 if (a_len > 0)
314 {
315 make_unique ();
316
317 for (octave_idx_type i = 0; i < a_len; i++)
318 xelem (r+i, c) = a.elem (i);
319 }
320
321 return *this;
322}
323
324FloatComplexMatrix&
325FloatComplexMatrix::insert (const FloatComplexDiagMatrix& a,
326 octave_idx_type r, octave_idx_type c)
327{
328 octave_idx_type a_nr = a.rows ();
329 octave_idx_type a_nc = a.cols ();
330
331 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
332 (*current_liboctave_error_handler) ("range error for insert");
333
334 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
335
336 octave_idx_type a_len = a.length ();
337
338 if (a_len > 0)
339 {
340 make_unique ();
341
342 for (octave_idx_type i = 0; i < a_len; i++)
343 xelem (r+i, c+i) = a.elem (i, i);
344 }
345
346 return *this;
347}
348
349FloatComplexMatrix&
350FloatComplexMatrix::fill (float val)
351{
352 octave_idx_type nr = rows ();
353 octave_idx_type nc = cols ();
354
355 if (nr > 0 && nc > 0)
356 {
357 make_unique ();
358
359 for (octave_idx_type j = 0; j < nc; j++)
360 for (octave_idx_type i = 0; i < nr; i++)
361 xelem (i, j) = val;
362 }
363
364 return *this;
365}
366
367FloatComplexMatrix&
368FloatComplexMatrix::fill (const FloatComplex& val)
369{
370 octave_idx_type nr = rows ();
371 octave_idx_type nc = cols ();
372
373 if (nr > 0 && nc > 0)
374 {
375 make_unique ();
376
377 for (octave_idx_type j = 0; j < nc; j++)
378 for (octave_idx_type i = 0; i < nr; i++)
379 xelem (i, j) = val;
380 }
381
382 return *this;
383}
384
385FloatComplexMatrix&
386FloatComplexMatrix::fill (float val, octave_idx_type r1, octave_idx_type c1,
387 octave_idx_type r2, octave_idx_type c2)
388{
389 octave_idx_type nr = rows ();
390 octave_idx_type nc = cols ();
391
392 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
393 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
394 (*current_liboctave_error_handler) ("range error for fill");
395
396 if (r1 > r2) { std::swap (r1, r2); }
397 if (c1 > c2) { std::swap (c1, c2); }
398
399 if (r2 >= r1 && c2 >= c1)
400 {
401 make_unique ();
402
403 for (octave_idx_type j = c1; j <= c2; j++)
404 for (octave_idx_type i = r1; i <= r2; i++)
405 xelem (i, j) = val;
406 }
407
408 return *this;
409}
410
411FloatComplexMatrix&
412FloatComplexMatrix::fill (const FloatComplex& val,
413 octave_idx_type r1, octave_idx_type c1,
414 octave_idx_type r2, octave_idx_type c2)
415{
416 octave_idx_type nr = rows ();
417 octave_idx_type nc = cols ();
418
419 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
420 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
421 (*current_liboctave_error_handler) ("range error for fill");
422
423 if (r1 > r2) { std::swap (r1, r2); }
424 if (c1 > c2) { std::swap (c1, c2); }
425
426 if (r2 >= r1 && c2 >=c1)
427 {
428 make_unique ();
429
430 for (octave_idx_type j = c1; j <= c2; j++)
431 for (octave_idx_type i = r1; i <= r2; i++)
432 xelem (i, j) = val;
433 }
434
435 return *this;
436}
437
438FloatComplexMatrix
439FloatComplexMatrix::append (const FloatMatrix& a) const
440{
441 octave_idx_type nr = rows ();
442 octave_idx_type nc = cols ();
443 if (nr != a.rows ())
444 (*current_liboctave_error_handler) ("row dimension mismatch for append");
445
446 octave_idx_type nc_insert = nc;
447 FloatComplexMatrix retval (nr, nc + a.cols ());
448 retval.insert (*this, 0, 0);
449 retval.insert (a, 0, nc_insert);
450 return retval;
451}
452
453FloatComplexMatrix
454FloatComplexMatrix::append (const FloatRowVector& a) const
455{
456 octave_idx_type nr = rows ();
457 octave_idx_type nc = cols ();
458 if (nr != 1)
459 (*current_liboctave_error_handler) ("row dimension mismatch for append");
460
461 octave_idx_type nc_insert = nc;
462 FloatComplexMatrix retval (nr, nc + a.numel ());
463 retval.insert (*this, 0, 0);
464 retval.insert (a, 0, nc_insert);
465 return retval;
466}
467
468FloatComplexMatrix
469FloatComplexMatrix::append (const FloatColumnVector& a) const
470{
471 octave_idx_type nr = rows ();
472 octave_idx_type nc = cols ();
473 if (nr != a.numel ())
474 (*current_liboctave_error_handler) ("row dimension mismatch for append");
475
476 octave_idx_type nc_insert = nc;
477 FloatComplexMatrix retval (nr, nc + 1);
478 retval.insert (*this, 0, 0);
479 retval.insert (a, 0, nc_insert);
480 return retval;
481}
482
483FloatComplexMatrix
484FloatComplexMatrix::append (const FloatDiagMatrix& a) const
485{
486 octave_idx_type nr = rows ();
487 octave_idx_type nc = cols ();
488 if (nr != a.rows ())
489 (*current_liboctave_error_handler) ("row dimension mismatch for append");
490
491 octave_idx_type nc_insert = nc;
492 FloatComplexMatrix retval (nr, nc + a.cols ());
493 retval.insert (*this, 0, 0);
494 retval.insert (a, 0, nc_insert);
495 return retval;
496}
497
498FloatComplexMatrix
499FloatComplexMatrix::append (const FloatComplexMatrix& a) const
500{
501 octave_idx_type nr = rows ();
502 octave_idx_type nc = cols ();
503 if (nr != a.rows ())
504 (*current_liboctave_error_handler) ("row dimension mismatch for append");
505
506 octave_idx_type nc_insert = nc;
507 FloatComplexMatrix retval (nr, nc + a.cols ());
508 retval.insert (*this, 0, 0);
509 retval.insert (a, 0, nc_insert);
510 return retval;
511}
512
513FloatComplexMatrix
514FloatComplexMatrix::append (const FloatComplexRowVector& a) const
515{
516 octave_idx_type nr = rows ();
517 octave_idx_type nc = cols ();
518 if (nr != 1)
519 (*current_liboctave_error_handler) ("row dimension mismatch for append");
520
521 octave_idx_type nc_insert = nc;
522 FloatComplexMatrix retval (nr, nc + a.numel ());
523 retval.insert (*this, 0, 0);
524 retval.insert (a, 0, nc_insert);
525 return retval;
526}
527
528FloatComplexMatrix
529FloatComplexMatrix::append (const FloatComplexColumnVector& a) const
530{
531 octave_idx_type nr = rows ();
532 octave_idx_type nc = cols ();
533 if (nr != a.numel ())
534 (*current_liboctave_error_handler) ("row dimension mismatch for append");
535
536 octave_idx_type nc_insert = nc;
537 FloatComplexMatrix retval (nr, nc + 1);
538 retval.insert (*this, 0, 0);
539 retval.insert (a, 0, nc_insert);
540 return retval;
541}
542
543FloatComplexMatrix
544FloatComplexMatrix::append (const FloatComplexDiagMatrix& a) const
545{
546 octave_idx_type nr = rows ();
547 octave_idx_type nc = cols ();
548 if (nr != a.rows ())
549 (*current_liboctave_error_handler) ("row dimension mismatch for append");
550
551 octave_idx_type nc_insert = nc;
552 FloatComplexMatrix retval (nr, nc + a.cols ());
553 retval.insert (*this, 0, 0);
554 retval.insert (a, 0, nc_insert);
555 return retval;
556}
557
558FloatComplexMatrix
559FloatComplexMatrix::stack (const FloatMatrix& a) const
560{
561 octave_idx_type nr = rows ();
562 octave_idx_type nc = cols ();
563 if (nc != a.cols ())
564 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
565
566 octave_idx_type nr_insert = nr;
567 FloatComplexMatrix retval (nr + a.rows (), nc);
568 retval.insert (*this, 0, 0);
569 retval.insert (a, nr_insert, 0);
570 return retval;
571}
572
573FloatComplexMatrix
574FloatComplexMatrix::stack (const FloatRowVector& a) const
575{
576 octave_idx_type nr = rows ();
577 octave_idx_type nc = cols ();
578 if (nc != a.numel ())
579 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
580
581 octave_idx_type nr_insert = nr;
582 FloatComplexMatrix retval (nr + 1, nc);
583 retval.insert (*this, 0, 0);
584 retval.insert (a, nr_insert, 0);
585 return retval;
586}
587
588FloatComplexMatrix
589FloatComplexMatrix::stack (const FloatColumnVector& a) const
590{
591 octave_idx_type nr = rows ();
592 octave_idx_type nc = cols ();
593 if (nc != 1)
594 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
595
596 octave_idx_type nr_insert = nr;
597 FloatComplexMatrix retval (nr + a.numel (), nc);
598 retval.insert (*this, 0, 0);
599 retval.insert (a, nr_insert, 0);
600 return retval;
601}
602
603FloatComplexMatrix
604FloatComplexMatrix::stack (const FloatDiagMatrix& a) const
605{
606 octave_idx_type nr = rows ();
607 octave_idx_type nc = cols ();
608 if (nc != a.cols ())
609 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
610
611 octave_idx_type nr_insert = nr;
612 FloatComplexMatrix retval (nr + a.rows (), nc);
613 retval.insert (*this, 0, 0);
614 retval.insert (a, nr_insert, 0);
615 return retval;
616}
617
618FloatComplexMatrix
619FloatComplexMatrix::stack (const FloatComplexMatrix& a) const
620{
621 octave_idx_type nr = rows ();
622 octave_idx_type nc = cols ();
623 if (nc != a.cols ())
624 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
625
626 octave_idx_type nr_insert = nr;
627 FloatComplexMatrix retval (nr + a.rows (), nc);
628 retval.insert (*this, 0, 0);
629 retval.insert (a, nr_insert, 0);
630 return retval;
631}
632
633FloatComplexMatrix
634FloatComplexMatrix::stack (const FloatComplexRowVector& a) const
635{
636 octave_idx_type nr = rows ();
637 octave_idx_type nc = cols ();
638 if (nc != a.numel ())
639 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
640
641 octave_idx_type nr_insert = nr;
642 FloatComplexMatrix retval (nr + 1, nc);
643 retval.insert (*this, 0, 0);
644 retval.insert (a, nr_insert, 0);
645 return retval;
646}
647
648FloatComplexMatrix
649FloatComplexMatrix::stack (const FloatComplexColumnVector& a) const
650{
651 octave_idx_type nr = rows ();
652 octave_idx_type nc = cols ();
653 if (nc != 1)
654 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
655
656 octave_idx_type nr_insert = nr;
657 FloatComplexMatrix retval (nr + a.numel (), nc);
658 retval.insert (*this, 0, 0);
659 retval.insert (a, nr_insert, 0);
660 return retval;
661}
662
663FloatComplexMatrix
664FloatComplexMatrix::stack (const FloatComplexDiagMatrix& a) const
665{
666 octave_idx_type nr = rows ();
667 octave_idx_type nc = cols ();
668 if (nc != a.cols ())
669 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
670
671 octave_idx_type nr_insert = nr;
672 FloatComplexMatrix retval (nr + a.rows (), nc);
673 retval.insert (*this, 0, 0);
674 retval.insert (a, nr_insert, 0);
675 return retval;
676}
677
678FloatComplexMatrix
679conj (const FloatComplexMatrix& a)
680{
681 return do_mx_unary_map<FloatComplex, FloatComplex, std::conj<float>> (a);
682}
683
684// resize is the destructive equivalent for this one
685
686FloatComplexMatrix
687FloatComplexMatrix::extract (octave_idx_type r1, octave_idx_type c1,
688 octave_idx_type r2, octave_idx_type c2) const
689{
690 if (r1 > r2) { std::swap (r1, r2); }
691 if (c1 > c2) { std::swap (c1, c2); }
692
693 return index (octave::idx_vector (r1, r2+1), octave::idx_vector (c1, c2+1));
694}
695
696FloatComplexMatrix
697FloatComplexMatrix::extract_n (octave_idx_type r1, octave_idx_type c1,
698 octave_idx_type nr, octave_idx_type nc) const
699{
700 return index (octave::idx_vector (r1, r1 + nr), octave::idx_vector (c1, c1 + nc));
701}
702
703// extract row or column i.
704
705FloatComplexRowVector
706FloatComplexMatrix::row (octave_idx_type i) const
707{
708 return index (octave::idx_vector (i), octave::idx_vector::colon);
709}
710
711FloatComplexColumnVector
712FloatComplexMatrix::column (octave_idx_type i) const
713{
714 return index (octave::idx_vector::colon, octave::idx_vector (i));
715}
716
717// Local function to calculate the 1-norm.
718static
719float
720norm1 (const FloatComplexMatrix& a)
721{
722 float anorm = 0.0;
723 FloatRowVector colsum = a.abs ().sum ().row (0);
724
725 for (octave_idx_type i = 0; i < colsum.numel (); i++)
726 {
727 float sum = colsum.xelem (i);
728 if (octave::math::isinf (sum) || octave::math::isnan (sum))
729 {
730 anorm = sum; // Pass Inf or NaN to output
731 break;
732 }
733 else
734 anorm = std::max (anorm, sum);
735 }
736
737 return anorm;
738}
739
740// Local function to check if matrix is singular based on rcond.
741static inline
742bool
743is_singular (const float rcond)
744{
745 return (std::abs (rcond) <= std::numeric_limits<float>::epsilon ());
746}
747
748FloatComplexMatrix
749FloatComplexMatrix::inverse () const
750{
751 octave_idx_type info;
752 float rcon;
753 MatrixType mattype (*this);
754 return inverse (mattype, info, rcon, false, false);
755}
756
757FloatComplexMatrix
758FloatComplexMatrix::inverse (octave_idx_type& info) const
759{
760 float rcon;
761 MatrixType mattype (*this);
762 return inverse (mattype, info, rcon, false, false);
763}
764
765FloatComplexMatrix
766FloatComplexMatrix::inverse (octave_idx_type& info, float& rcon, bool force,
767 bool calc_cond) const
768{
769 MatrixType mattype (*this);
770 return inverse (mattype, info, rcon, force, calc_cond);
771}
772
773FloatComplexMatrix
774FloatComplexMatrix::inverse (MatrixType& mattype) const
775{
776 octave_idx_type info;
777 float rcon;
778 return inverse (mattype, info, rcon, false, false);
779}
780
781FloatComplexMatrix
782FloatComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const
783{
784 float rcon;
785 return inverse (mattype, info, rcon, false, false);
786}
787
788FloatComplexMatrix
789FloatComplexMatrix::tinverse (MatrixType& mattype, octave_idx_type& info,
790 float& rcon, bool force, bool calc_cond) const
791{
792 FloatComplexMatrix retval;
793
794 F77_INT nr = octave::to_f77_int (rows ());
795 F77_INT nc = octave::to_f77_int (cols ());
796
797 if (nr != nc || nr == 0 || nc == 0)
798 (*current_liboctave_error_handler) ("inverse requires square matrix");
799
800 int typ = mattype.type ();
801 char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
802 char udiag = 'N';
803 retval = *this;
804 FloatComplex *tmp_data = retval.rwdata ();
805
806 F77_INT tmp_info = 0;
807
808 F77_XFCN (ctrtri, CTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
809 F77_CONST_CHAR_ARG2 (&udiag, 1),
810 nr, F77_CMPLX_ARG (tmp_data), nr, tmp_info
811 F77_CHAR_ARG_LEN (1)
812 F77_CHAR_ARG_LEN (1)));
813
814 info = tmp_info;
815
816 // Throw away extra info LAPACK gives so as to not change output.
817 rcon = 0.0;
818 if (info != 0)
819 info = -1;
820 else if (calc_cond)
821 {
822 F77_INT ztrcon_info = 0;
823 char job = '1';
824
825 OCTAVE_LOCAL_BUFFER (FloatComplex, cwork, 2*nr);
826 OCTAVE_LOCAL_BUFFER (float, rwork, nr);
827
828 F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
829 F77_CONST_CHAR_ARG2 (&uplo, 1),
830 F77_CONST_CHAR_ARG2 (&udiag, 1),
831 nr, F77_CMPLX_ARG (tmp_data), nr, rcon,
832 F77_CMPLX_ARG (cwork), rwork, ztrcon_info
833 F77_CHAR_ARG_LEN (1)
834 F77_CHAR_ARG_LEN (1)
835 F77_CHAR_ARG_LEN (1)));
836
837 if (ztrcon_info != 0)
838 info = -1;
839 }
840
841 if (info == -1 && ! force)
842 retval = *this; // Restore matrix contents.
843
844 return retval;
845}
846
847FloatComplexMatrix
848FloatComplexMatrix::finverse (MatrixType& mattype, octave_idx_type& info,
849 float& rcon, bool force, bool calc_cond) const
850{
851 FloatComplexMatrix retval;
852
853 F77_INT nr = octave::to_f77_int (rows ());
854 F77_INT nc = octave::to_f77_int (cols ());
855
856 if (nr != nc)
857 (*current_liboctave_error_handler) ("inverse requires square matrix");
858
859 Array<F77_INT> ipvt (dim_vector (nr, 1));
860 F77_INT *pipvt = ipvt.rwdata ();
861
862 retval = *this;
863 FloatComplex *tmp_data = retval.rwdata ();
864
865 Array<FloatComplex> z (dim_vector (1, 1));
866 F77_INT lwork = -1;
867
868 // Query the optimum work array size.
869
870 F77_INT tmp_info = 0;
871
872 F77_XFCN (cgetri, CGETRI, (nc, F77_CMPLX_ARG (tmp_data), nr, pipvt,
873 F77_CMPLX_ARG (z.rwdata ()), lwork,
874 tmp_info));
875
876 lwork = static_cast<F77_INT> (std::real (z(0)));
877 lwork = (lwork < 2 * nc ? 2 * nc : lwork);
878 z.resize (dim_vector (lwork, 1));
879 FloatComplex *pz = z.rwdata ();
880
881 info = 0;
882 tmp_info = 0;
883
884 // Calculate norm of the matrix (always, see bug #45577) for later use.
885 float anorm = norm1 (retval);
886
887 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
888 // and bug #46330, segfault with matrices containing Inf & NaN
889 if (octave::math::isnan (anorm) || octave::math::isinf (anorm))
890 info = -1;
891 else
892 {
893 F77_XFCN (cgetrf, CGETRF, (nc, nc, F77_CMPLX_ARG (tmp_data), nr, pipvt,
894 tmp_info));
895
896 info = tmp_info;
897 }
898
899 // Throw away extra info LAPACK gives so as to not change output.
900 rcon = 0.0;
901 if (info != 0)
902 info = -1;
903 else if (calc_cond)
904 {
905 if (octave::math::isnan (anorm))
906 rcon = octave::numeric_limits<float>::NaN ();
907 else
908 {
909 F77_INT cgecon_info = 0;
910
911 // Now calculate the condition number for non-singular matrix.
912 char job = '1';
913 Array<float> rz (dim_vector (2 * nc, 1));
914 float *prz = rz.rwdata ();
915 F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
916 nc, F77_CMPLX_ARG (tmp_data), nr, anorm,
917 rcon, F77_CMPLX_ARG (pz), prz, cgecon_info
918 F77_CHAR_ARG_LEN (1)));
919
920 if (cgecon_info != 0)
921 info = -1;
922 }
923 }
924
925 if ((info == -1 && ! force)
926 || octave::math::isnan (anorm) || octave::math::isinf (anorm))
927 retval = *this; // Restore contents.
928 else
929 {
930 F77_INT zgetri_info = 0;
931
932 F77_XFCN (cgetri, CGETRI, (nc, F77_CMPLX_ARG (tmp_data), nr, pipvt,
933 F77_CMPLX_ARG (pz), lwork, zgetri_info));
934
935 if (zgetri_info != 0)
936 info = -1;
937 }
938
939 if (info != 0)
940 mattype.mark_as_rectangular ();
941
942 return retval;
943}
944
945FloatComplexMatrix
946FloatComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info,
947 float& rcon, bool force, bool calc_cond) const
948{
949 int typ = mattype.type (false);
950 FloatComplexMatrix ret;
951
952 if (typ == MatrixType::Unknown)
953 typ = mattype.type (*this);
954
955 if (typ == MatrixType::Diagonal) // a scalar is classified as Diagonal.
956 {
957 FloatComplex scalar = this->elem (0);
958 float real = std::real (scalar);
959 float imag = std::imag (scalar);
960
961 if (real == 0 && imag == 0)
962 ret = FloatComplexMatrix (1, 1,
963 FloatComplex (octave::numeric_limits<float>::Inf (), 0.0));
964 else
965 ret = FloatComplex (1, 0) / (*this);
966
967 if (calc_cond)
968 {
969 if (octave::math::isfinite (real) && octave::math::isfinite (imag)
970 && (real != 0 || imag != 0))
971 rcon = 1.0f;
972 else if (octave::math::isinf (real) || octave::math::isinf (imag)
973 || (real == 0 && imag == 0))
974 rcon = 0.0f;
975 else
976 rcon = octave::numeric_limits<float>::NaN ();
977 }
978 }
979 else if (typ == MatrixType::Upper || typ == MatrixType::Lower)
980 ret = tinverse (mattype, info, rcon, force, calc_cond);
981 else
982 {
983 if (mattype.ishermitian ())
984 {
985 octave::math::chol<FloatComplexMatrix> chol (*this, info, true, calc_cond);
986 if (info == 0)
987 {
988 if (calc_cond)
989 rcon = chol.rcond ();
990 else
991 rcon = 1.0;
992 ret = chol.inverse ();
993 }
994 else
995 mattype.mark_as_unsymmetric ();
996 }
997
998 if (! mattype.ishermitian ())
999 ret = finverse (mattype, info, rcon, force, calc_cond);
1000
1001 if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0)
1002 {
1003 ret = FloatComplexMatrix (rows (), columns (),
1004 FloatComplex (octave::numeric_limits<float>::Inf (), 0.0));
1005 }
1006 }
1007
1008 return ret;
1009}
1010
1011FloatComplexMatrix
1012FloatComplexMatrix::pseudo_inverse (float tol) const
1013{
1014 FloatComplexMatrix retval;
1015
1016 octave::math::svd<FloatComplexMatrix> result
1017 (*this, octave::math::svd<FloatComplexMatrix>::Type::economy);
1018
1019 FloatDiagMatrix S = result.singular_values ();
1020 FloatComplexMatrix U = result.left_singular_matrix ();
1021 FloatComplexMatrix V = result.right_singular_matrix ();
1022
1023 FloatColumnVector sigma = S.extract_diag ();
1024
1025 octave_idx_type r = sigma.numel () - 1;
1026 octave_idx_type nr = rows ();
1027 octave_idx_type nc = cols ();
1028
1029 if (tol <= 0.0)
1030 {
1031 tol = std::max (nr, nc) * sigma.elem (0)
1032 * std::numeric_limits<float>::epsilon ();
1033
1034 if (tol == 0)
1035 tol = std::numeric_limits<float>::min ();
1036 }
1037
1038 while (r >= 0 && sigma.elem (r) < tol)
1039 r--;
1040
1041 if (r < 0)
1042 retval = FloatComplexMatrix (nc, nr, 0.0);
1043 else
1044 {
1045 FloatComplexMatrix Ur = U.extract (0, 0, nr-1, r);
1046 FloatDiagMatrix D = FloatDiagMatrix (sigma.extract (0, r)).inverse ();
1047 FloatComplexMatrix Vr = V.extract (0, 0, nc-1, r);
1048 retval = Vr * D * Ur.hermitian ();
1049 }
1050
1051 return retval;
1052}
1053
1054#if defined (HAVE_FFTW)
1055
1056FloatComplexMatrix
1057FloatComplexMatrix::fourier () const
1058{
1059 std::size_t nr = rows ();
1060 std::size_t nc = cols ();
1061
1062 FloatComplexMatrix retval (nr, nc);
1063
1064 std::size_t npts, nsamples;
1065
1066 if (nr == 1 || nc == 1)
1067 {
1068 npts = (nr > nc ? nr : nc);
1069 nsamples = 1;
1070 }
1071 else
1072 {
1073 npts = nr;
1074 nsamples = nc;
1075 }
1076
1077 const FloatComplex *in (data ());
1078 FloatComplex *out (retval.rwdata ());
1079
1080 octave::fftw::fft (in, out, npts, nsamples);
1081
1082 return retval;
1083}
1084
1085FloatComplexMatrix
1086FloatComplexMatrix::ifourier () const
1087{
1088 std::size_t nr = rows ();
1089 std::size_t nc = cols ();
1090
1091 FloatComplexMatrix retval (nr, nc);
1092
1093 std::size_t npts, nsamples;
1094
1095 if (nr == 1 || nc == 1)
1096 {
1097 npts = (nr > nc ? nr : nc);
1098 nsamples = 1;
1099 }
1100 else
1101 {
1102 npts = nr;
1103 nsamples = nc;
1104 }
1105
1106 const FloatComplex *in (data ());
1107 FloatComplex *out (retval.rwdata ());
1108
1109 octave::fftw::ifft (in, out, npts, nsamples);
1110
1111 return retval;
1112}
1113
1114FloatComplexMatrix
1115FloatComplexMatrix::fourier2d () const
1116{
1117 dim_vector dv (rows (), cols ());
1118
1119 FloatComplexMatrix retval (rows (), cols ());
1120 const FloatComplex *in (data ());
1121 FloatComplex *out (retval.rwdata ());
1122
1123 octave::fftw::fftNd (in, out, 2, dv);
1124
1125 return retval;
1126}
1127
1128FloatComplexMatrix
1129FloatComplexMatrix::ifourier2d () const
1130{
1131 dim_vector dv (rows (), cols ());
1132
1133 FloatComplexMatrix retval (rows (), cols ());
1134 const FloatComplex *in (data ());
1135 FloatComplex *out (retval.rwdata ());
1136
1137 octave::fftw::ifftNd (in, out, 2, dv);
1138
1139 return retval;
1140}
1141
1142#else
1143
1144FloatComplexMatrix
1145FloatComplexMatrix::fourier () const
1146{
1147 (*current_liboctave_error_handler)
1148 ("support for FFTW was unavailable or disabled when liboctave was built");
1149
1150 return FloatComplexMatrix ();
1151}
1152
1153FloatComplexMatrix
1154FloatComplexMatrix::ifourier () const
1155{
1156 (*current_liboctave_error_handler)
1157 ("support for FFTW was unavailable or disabled when liboctave was built");
1158
1159 return FloatComplexMatrix ();
1160}
1161
1162FloatComplexMatrix
1163FloatComplexMatrix::fourier2d () const
1164{
1165 (*current_liboctave_error_handler)
1166 ("support for FFTW was unavailable or disabled when liboctave was built");
1167
1168 return FloatComplexMatrix ();
1169}
1170
1171FloatComplexMatrix
1172FloatComplexMatrix::ifourier2d () const
1173{
1174 (*current_liboctave_error_handler)
1175 ("support for FFTW was unavailable or disabled when liboctave was built");
1176
1177 return FloatComplexMatrix ();
1178}
1179
1180#endif
1181
1182FloatComplexDET
1183FloatComplexMatrix::determinant () const
1184{
1185 octave_idx_type info;
1186 float rcon;
1187 return determinant (info, rcon, 0);
1188}
1189
1190FloatComplexDET
1191FloatComplexMatrix::determinant (octave_idx_type& info) const
1192{
1193 float rcon;
1194 return determinant (info, rcon, 0);
1195}
1196
1197FloatComplexDET
1198FloatComplexMatrix::determinant (octave_idx_type& info, float& rcon,
1199 bool calc_cond) const
1200{
1201 MatrixType mattype (*this);
1202 return determinant (mattype, info, rcon, calc_cond);
1203}
1204
1205FloatComplexDET
1206FloatComplexMatrix::determinant (MatrixType& mattype,
1207 octave_idx_type& info, float& rcon,
1208 bool calc_cond) const
1209{
1210 FloatComplexDET retval (1.0);
1211
1212 info = 0;
1213 rcon = 0.0;
1214
1215 F77_INT nr = octave::to_f77_int (rows ());
1216 F77_INT nc = octave::to_f77_int (cols ());
1217
1218 if (nr != nc)
1219 (*current_liboctave_error_handler) ("matrix must be square");
1220
1221 int typ = mattype.type ();
1222
1223 // Even though the matrix is marked as singular (Rectangular), we may
1224 // still get a useful number from the LU factorization, because it always
1225 // completes.
1226
1227 if (typ == MatrixType::Unknown)
1228 typ = mattype.type (*this);
1229 else if (typ == MatrixType::Rectangular)
1230 typ = MatrixType::Full;
1231
1232 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1233 {
1234 for (F77_INT i = 0; i < nc; i++)
1235 retval *= elem (i, i);
1236 }
1237 else if (typ == MatrixType::Hermitian)
1238 {
1239 FloatComplexMatrix atmp = *this;
1240 FloatComplex *tmp_data = atmp.rwdata ();
1241
1242 float anorm;
1243 if (calc_cond)
1244 anorm = norm1 (*this);
1245
1246 F77_INT tmp_info = 0;
1247
1248 char job = 'L';
1249 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1250 F77_CMPLX_ARG (tmp_data), nr, tmp_info
1251 F77_CHAR_ARG_LEN (1)));
1252
1253 info = tmp_info;
1254
1255 if (info != 0)
1256 {
1257 rcon = 0.0;
1258 mattype.mark_as_unsymmetric ();
1259 typ = MatrixType::Full;
1260 }
1261 else
1262 {
1263 if (calc_cond)
1264 {
1265 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1266 FloatComplex *pz = z.rwdata ();
1267 Array<float> rz (dim_vector (nc, 1));
1268 float *prz = rz.rwdata ();
1269
1270 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1271 nr, F77_CMPLX_ARG (tmp_data), nr, anorm,
1272 rcon, F77_CMPLX_ARG (pz), prz, tmp_info
1273 F77_CHAR_ARG_LEN (1)));
1274
1275 info = tmp_info;
1276
1277 if (info != 0)
1278 rcon = 0.0;
1279 }
1280
1281 for (F77_INT i = 0; i < nc; i++)
1282 retval *= atmp(i, i);
1283
1284 retval = retval.square ();
1285 }
1286 }
1287 else if (typ != MatrixType::Full)
1288 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1289
1290 if (typ == MatrixType::Full)
1291 {
1292 Array<F77_INT> ipvt (dim_vector (nr, 1));
1293 F77_INT *pipvt = ipvt.rwdata ();
1294
1295 FloatComplexMatrix atmp = *this;
1296 FloatComplex *tmp_data = atmp.rwdata ();
1297
1298 info = 0;
1299
1300 // Calculate norm of the matrix (always, see bug #45577) for later use.
1301 float anorm = norm1 (*this);
1302
1303 F77_INT tmp_info = 0;
1304
1305 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1306 if (octave::math::isnan (anorm))
1307 info = -1;
1308 else
1309 {
1310 F77_XFCN (cgetrf, CGETRF, (nr, nr, F77_CMPLX_ARG (tmp_data), nr, pipvt,
1311 tmp_info));
1312
1313 info = tmp_info;
1314 }
1315
1316 // Throw away extra info LAPACK gives so as to not change output.
1317 rcon = 0.0;
1318 if (info != 0)
1319 {
1320 info = -1;
1321 retval = FloatComplexDET ();
1322 }
1323 else
1324 {
1325 if (calc_cond)
1326 {
1327 // Now calc the condition number for non-singular matrix.
1328 char job = '1';
1329 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1330 FloatComplex *pz = z.rwdata ();
1331 Array<float> rz (dim_vector (2 * nc, 1));
1332 float *prz = rz.rwdata ();
1333
1334 F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1335 nc, F77_CMPLX_ARG (tmp_data), nr, anorm,
1336 rcon, F77_CMPLX_ARG (pz), prz, tmp_info
1337 F77_CHAR_ARG_LEN (1)));
1338
1339 info = tmp_info;
1340 }
1341
1342 if (info != 0)
1343 {
1344 info = -1;
1345 retval = FloatComplexDET ();
1346 }
1347 else
1348 {
1349 for (F77_INT i = 0; i < nc; i++)
1350 {
1351 FloatComplex c = atmp(i, i);
1352 retval *= (ipvt(i) != (i+1)) ? -c : c;
1353 }
1354 }
1355 }
1356 }
1357
1358 return retval;
1359}
1360
1361float
1362FloatComplexMatrix::rcond () const
1363{
1364 MatrixType mattype (*this);
1365 return rcond (mattype);
1366}
1367
1368float
1369FloatComplexMatrix::rcond (MatrixType& mattype) const
1370{
1371 float rcon = octave::numeric_limits<float>::NaN ();
1372 F77_INT nr = octave::to_f77_int (rows ());
1373 F77_INT nc = octave::to_f77_int (cols ());
1374
1375 if (nr != nc)
1376 (*current_liboctave_error_handler) ("matrix must be square");
1377
1378 if (nr == 0 || nc == 0)
1379 rcon = octave::numeric_limits<float>::Inf ();
1380 else
1381 {
1382 int typ = mattype.type ();
1383
1384 if (typ == MatrixType::Unknown)
1385 typ = mattype.type (*this);
1386
1387 // Only calculate the condition number for LU/Cholesky
1388 if (typ == MatrixType::Upper)
1389 {
1390 const FloatComplex *tmp_data = data ();
1391 F77_INT info = 0;
1392 char norm = '1';
1393 char uplo = 'U';
1394 char dia = 'N';
1395
1396 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1397 FloatComplex *pz = z.rwdata ();
1398 Array<float> rz (dim_vector (nc, 1));
1399 float *prz = rz.rwdata ();
1400
1401 F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1402 F77_CONST_CHAR_ARG2 (&uplo, 1),
1403 F77_CONST_CHAR_ARG2 (&dia, 1),
1404 nr, F77_CONST_CMPLX_ARG (tmp_data), nr, rcon,
1405 F77_CMPLX_ARG (pz), prz, info
1406 F77_CHAR_ARG_LEN (1)
1407 F77_CHAR_ARG_LEN (1)
1408 F77_CHAR_ARG_LEN (1)));
1409
1410 if (info != 0)
1411 rcon = 0;
1412 }
1413 else if (typ == MatrixType::Permuted_Upper)
1414 (*current_liboctave_error_handler)
1415 ("permuted triangular matrix not implemented");
1416 else if (typ == MatrixType::Lower)
1417 {
1418 const FloatComplex *tmp_data = data ();
1419 F77_INT info = 0;
1420 char norm = '1';
1421 char uplo = 'L';
1422 char dia = 'N';
1423
1424 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1425 FloatComplex *pz = z.rwdata ();
1426 Array<float> rz (dim_vector (nc, 1));
1427 float *prz = rz.rwdata ();
1428
1429 F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1430 F77_CONST_CHAR_ARG2 (&uplo, 1),
1431 F77_CONST_CHAR_ARG2 (&dia, 1),
1432 nr, F77_CONST_CMPLX_ARG (tmp_data), nr, rcon,
1433 F77_CMPLX_ARG (pz), prz, info
1434 F77_CHAR_ARG_LEN (1)
1435 F77_CHAR_ARG_LEN (1)
1436 F77_CHAR_ARG_LEN (1)));
1437
1438 if (info != 0)
1439 rcon = 0.0;
1440 }
1441 else if (typ == MatrixType::Permuted_Lower)
1442 (*current_liboctave_error_handler)
1443 ("permuted triangular matrix not implemented");
1444 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1445 {
1446 float anorm = -1.0;
1447
1448 if (typ == MatrixType::Hermitian)
1449 {
1450 F77_INT info = 0;
1451 char job = 'L';
1452
1453 FloatComplexMatrix atmp = *this;
1454 FloatComplex *tmp_data = atmp.rwdata ();
1455
1456 anorm = norm1 (atmp);
1457
1458 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1459 F77_CMPLX_ARG (tmp_data), nr, info
1460 F77_CHAR_ARG_LEN (1)));
1461
1462 if (info != 0)
1463 {
1464 rcon = 0.0;
1465
1466 mattype.mark_as_unsymmetric ();
1467 typ = MatrixType::Full;
1468 }
1469 else
1470 {
1471 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1472 FloatComplex *pz = z.rwdata ();
1473 Array<float> rz (dim_vector (nc, 1));
1474 float *prz = rz.rwdata ();
1475
1476 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1477 nr, F77_CMPLX_ARG (tmp_data), nr, anorm,
1478 rcon, F77_CMPLX_ARG (pz), prz, info
1479 F77_CHAR_ARG_LEN (1)));
1480
1481 if (info != 0)
1482 rcon = 0.0;
1483 }
1484 }
1485
1486 if (typ == MatrixType::Full)
1487 {
1488 F77_INT info = 0;
1489
1490 FloatComplexMatrix atmp = *this;
1491 FloatComplex *tmp_data = atmp.rwdata ();
1492
1493 Array<F77_INT> ipvt (dim_vector (nr, 1));
1494 F77_INT *pipvt = ipvt.rwdata ();
1495
1496 if (anorm < 0.0)
1497 anorm = norm1 (atmp);
1498
1499 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1500 FloatComplex *pz = z.rwdata ();
1501 Array<float> rz (dim_vector (2 * nc, 1));
1502 float *prz = rz.rwdata ();
1503
1504 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1505 if (octave::math::isnan (anorm))
1506 info = -1;
1507 else
1508 F77_XFCN (cgetrf, CGETRF, (nr, nr, F77_CMPLX_ARG (tmp_data),
1509 nr, pipvt, info));
1510
1511 if (info != 0)
1512 {
1513 rcon = 0.0;
1514 mattype.mark_as_rectangular ();
1515 }
1516 else
1517 {
1518 char job = '1';
1519 F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1520 nc, F77_CMPLX_ARG (tmp_data), nr, anorm,
1521 rcon, F77_CMPLX_ARG (pz), prz, info
1522 F77_CHAR_ARG_LEN (1)));
1523
1524 if (info != 0)
1525 rcon = 0.0;
1526 }
1527 }
1528 }
1529 else
1530 rcon = 0.0;
1531 }
1532
1533 return rcon;
1534}
1535
1536FloatComplexMatrix
1537FloatComplexMatrix::utsolve (MatrixType& mattype, const FloatComplexMatrix& b,
1538 octave_idx_type& info, float& rcon,
1539 solve_singularity_handler sing_handler,
1540 bool calc_cond, blas_trans_type transt) const
1541{
1542 FloatComplexMatrix retval;
1543
1544 F77_INT nr = octave::to_f77_int (rows ());
1545 F77_INT nc = octave::to_f77_int (cols ());
1546
1547 F77_INT b_nr = octave::to_f77_int (b.rows ());
1548 F77_INT b_nc = octave::to_f77_int (b.cols ());
1549
1550 if (nr != b_nr)
1551 (*current_liboctave_error_handler)
1552 ("matrix dimension mismatch solution of linear equations");
1553
1554 if (nr == 0 || nc == 0 || b_nc == 0)
1555 retval = FloatComplexMatrix (nc, b_nc, FloatComplex (0.0, 0.0));
1556 else
1557 {
1558 int typ = mattype.type ();
1559
1560 if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
1561 {
1562 rcon = 1.0;
1563 info = 0;
1564
1565 if (typ == MatrixType::Permuted_Upper)
1566 (*current_liboctave_error_handler)
1567 ("permuted triangular matrix not implemented");
1568 else
1569 {
1570 const FloatComplex *tmp_data = data ();
1571
1572 retval = b;
1573 FloatComplex *result = retval.rwdata ();
1574
1575 char uplo = 'U';
1576 char trans = get_blas_char (transt);
1577 char dia = 'N';
1578
1579 F77_INT tmp_info = 0;
1580
1581 F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1582 F77_CONST_CHAR_ARG2 (&trans, 1),
1583 F77_CONST_CHAR_ARG2 (&dia, 1),
1584 nr, b_nc, F77_CONST_CMPLX_ARG (tmp_data), nr,
1585 F77_CMPLX_ARG (result), nr, tmp_info
1586 F77_CHAR_ARG_LEN (1)
1587 F77_CHAR_ARG_LEN (1)
1588 F77_CHAR_ARG_LEN (1)));
1589
1590 info = tmp_info;
1591
1592 if (calc_cond)
1593 {
1594 char norm = '1';
1595 uplo = 'U';
1596 dia = 'N';
1597
1598 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1599 FloatComplex *pz = z.rwdata ();
1600 Array<float> rz (dim_vector (nc, 1));
1601 float *prz = rz.rwdata ();
1602
1603 F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1604 F77_CONST_CHAR_ARG2 (&uplo, 1),
1605 F77_CONST_CHAR_ARG2 (&dia, 1),
1606 nr, F77_CONST_CMPLX_ARG (tmp_data), nr, rcon,
1607 F77_CMPLX_ARG (pz), prz, tmp_info
1608 F77_CHAR_ARG_LEN (1)
1609 F77_CHAR_ARG_LEN (1)
1610 F77_CHAR_ARG_LEN (1)));
1611
1612 info = tmp_info;
1613
1614 if (info != 0)
1615 info = -2;
1616
1617 if (is_singular (rcon) || octave::math::isnan (rcon))
1618 {
1619 info = -2;
1620
1621 if (sing_handler)
1622 sing_handler (rcon);
1623 else
1624 octave::warn_singular_matrix (rcon);
1625 }
1626 }
1627 }
1628 }
1629 else
1630 (*current_liboctave_error_handler) ("incorrect matrix type");
1631 }
1632
1633 return retval;
1634}
1635
1636FloatComplexMatrix
1637FloatComplexMatrix::ltsolve (MatrixType& mattype, const FloatComplexMatrix& b,
1638 octave_idx_type& info, float& rcon,
1639 solve_singularity_handler sing_handler,
1640 bool calc_cond, blas_trans_type transt) const
1641{
1642 FloatComplexMatrix retval;
1643
1644 F77_INT nr = octave::to_f77_int (rows ());
1645 F77_INT nc = octave::to_f77_int (cols ());
1646
1647 F77_INT b_nr = octave::to_f77_int (b.rows ());
1648 F77_INT b_nc = octave::to_f77_int (b.cols ());
1649
1650 if (nr != b_nr)
1651 (*current_liboctave_error_handler)
1652 ("matrix dimension mismatch solution of linear equations");
1653
1654 if (nr == 0 || nc == 0 || b_nc == 0)
1655 retval = FloatComplexMatrix (nc, b_nc, FloatComplex (0.0, 0.0));
1656 else
1657 {
1658 int typ = mattype.type ();
1659
1660 if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
1661 {
1662 rcon = 1.0;
1663 info = 0;
1664
1665 if (typ == MatrixType::Permuted_Lower)
1666 (*current_liboctave_error_handler)
1667 ("permuted triangular matrix not implemented");
1668 else
1669 {
1670 const FloatComplex *tmp_data = data ();
1671
1672 retval = b;
1673 FloatComplex *result = retval.rwdata ();
1674
1675 char uplo = 'L';
1676 char trans = get_blas_char (transt);
1677 char dia = 'N';
1678
1679 F77_INT tmp_info = 0;
1680
1681 F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1682 F77_CONST_CHAR_ARG2 (&trans, 1),
1683 F77_CONST_CHAR_ARG2 (&dia, 1),
1684 nr, b_nc, F77_CONST_CMPLX_ARG (tmp_data), nr,
1685 F77_CMPLX_ARG (result), nr, tmp_info
1686 F77_CHAR_ARG_LEN (1)
1687 F77_CHAR_ARG_LEN (1)
1688 F77_CHAR_ARG_LEN (1)));
1689
1690 info = tmp_info;
1691
1692 if (calc_cond)
1693 {
1694 char norm = '1';
1695 uplo = 'L';
1696 dia = 'N';
1697
1698 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1699 FloatComplex *pz = z.rwdata ();
1700 Array<float> rz (dim_vector (nc, 1));
1701 float *prz = rz.rwdata ();
1702
1703 F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1704 F77_CONST_CHAR_ARG2 (&uplo, 1),
1705 F77_CONST_CHAR_ARG2 (&dia, 1),
1706 nr, F77_CONST_CMPLX_ARG (tmp_data), nr, rcon,
1707 F77_CMPLX_ARG (pz), prz, tmp_info
1708 F77_CHAR_ARG_LEN (1)
1709 F77_CHAR_ARG_LEN (1)
1710 F77_CHAR_ARG_LEN (1)));
1711
1712 info = tmp_info;
1713
1714 if (info != 0)
1715 info = -2;
1716
1717 if (is_singular (rcon) || octave::math::isnan (rcon))
1718 {
1719 info = -2;
1720
1721 if (sing_handler)
1722 sing_handler (rcon);
1723 else
1724 octave::warn_singular_matrix (rcon);
1725 }
1726 }
1727 }
1728 }
1729 else
1730 (*current_liboctave_error_handler) ("incorrect matrix type");
1731 }
1732
1733 return retval;
1734}
1735
1736FloatComplexMatrix
1737FloatComplexMatrix::fsolve (MatrixType& mattype, const FloatComplexMatrix& b,
1738 octave_idx_type& info, float& rcon,
1739 solve_singularity_handler sing_handler,
1740 bool calc_cond) const
1741{
1742 FloatComplexMatrix retval;
1743
1744 F77_INT nr = octave::to_f77_int (rows ());
1745 F77_INT nc = octave::to_f77_int (cols ());
1746
1747 F77_INT b_nr = octave::to_f77_int (b.rows ());
1748 F77_INT b_nc = octave::to_f77_int (b.cols ());
1749
1750 if (nr != nc || nr != b_nr)
1751 (*current_liboctave_error_handler)
1752 ("matrix dimension mismatch solution of linear equations");
1753
1754 if (nr == 0 || b_nc == 0)
1755 retval = FloatComplexMatrix (nc, b_nc, FloatComplex (0.0, 0.0));
1756 else
1757 {
1758 int typ = mattype.type ();
1759
1760 // Calculate the norm of the matrix for later use when determining rcon.
1761 float anorm = -1.0;
1762
1763 if (typ == MatrixType::Hermitian)
1764 {
1765 info = 0;
1766 char job = 'L';
1767
1768 FloatComplexMatrix atmp = *this;
1769 FloatComplex *tmp_data = atmp.rwdata ();
1770
1771 // The norm of the matrix for later use when determining rcon.
1772 if (calc_cond)
1773 anorm = norm1 (atmp);
1774
1775 F77_INT tmp_info = 0;
1776
1777 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1778 F77_CMPLX_ARG (tmp_data), nr, tmp_info
1779 F77_CHAR_ARG_LEN (1)));
1780
1781 info = tmp_info;
1782
1783 // Throw away extra info LAPACK gives so as to not change output.
1784 rcon = 0.0;
1785 if (info != 0)
1786 {
1787 info = -2;
1788
1789 mattype.mark_as_unsymmetric ();
1790 typ = MatrixType::Full;
1791 }
1792 else
1793 {
1794 if (calc_cond)
1795 {
1796 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1797 FloatComplex *pz = z.rwdata ();
1798 Array<float> rz (dim_vector (nc, 1));
1799 float *prz = rz.rwdata ();
1800
1801 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1802 nr, F77_CMPLX_ARG (tmp_data), nr, anorm,
1803 rcon, F77_CMPLX_ARG (pz), prz, tmp_info
1804 F77_CHAR_ARG_LEN (1)));
1805
1806 info = tmp_info;
1807
1808 if (info != 0)
1809 info = -2;
1810
1811 if (is_singular (rcon) || octave::math::isnan (rcon))
1812 {
1813 info = -2;
1814
1815 if (sing_handler)
1816 sing_handler (rcon);
1817 else
1818 octave::warn_singular_matrix (rcon);
1819 }
1820 }
1821
1822 if (info == 0)
1823 {
1824 retval = b;
1825 FloatComplex *result = retval.rwdata ();
1826
1827 F77_XFCN (cpotrs, CPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1828 nr, b_nc, F77_CMPLX_ARG (tmp_data), nr,
1829 F77_CMPLX_ARG (result), b_nr, tmp_info
1830 F77_CHAR_ARG_LEN (1)));
1831
1832 info = tmp_info;
1833 }
1834 else
1835 {
1836 mattype.mark_as_unsymmetric ();
1837 typ = MatrixType::Full;
1838 }
1839 }
1840 }
1841
1842 if (typ == MatrixType::Full)
1843 {
1844 info = 0;
1845
1846 Array<F77_INT> ipvt (dim_vector (nr, 1));
1847 F77_INT *pipvt = ipvt.rwdata ();
1848
1849 FloatComplexMatrix atmp = *this;
1850 FloatComplex *tmp_data = atmp.rwdata ();
1851
1852 Array<FloatComplex> z (dim_vector (2 * nc, 1));
1853 FloatComplex *pz = z.rwdata ();
1854 Array<float> rz (dim_vector (2 * nc, 1));
1855 float *prz = rz.rwdata ();
1856
1857 // Calculate the norm of the matrix, for later use.
1858 if (calc_cond && anorm < 0.0)
1859 anorm = norm1 (atmp);
1860
1861 F77_INT tmp_info = 0;
1862
1863 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1864 // and bug #46330, segfault with matrices containing Inf & NaN
1865 if (octave::math::isnan (anorm) || octave::math::isinf (anorm))
1866 info = -2;
1867 else
1868 {
1869 F77_XFCN (cgetrf, CGETRF, (nr, nr, F77_CMPLX_ARG (tmp_data),
1870 nr, pipvt, tmp_info));
1871
1872 info = tmp_info;
1873 }
1874
1875 // Throw away extra info LAPACK gives so as to not change output.
1876 rcon = 0.0;
1877 if (info != 0)
1878 {
1879 info = -2;
1880
1881 if (sing_handler)
1882 sing_handler (rcon);
1883 else
1884 octave::warn_singular_matrix ();
1885
1886 mattype.mark_as_rectangular ();
1887 }
1888 else
1889 {
1890 if (calc_cond)
1891 {
1892 // Calculate the condition number for non-singular matrix.
1893 char job = '1';
1894 F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1895 nc, F77_CMPLX_ARG (tmp_data), nr, anorm,
1896 rcon, F77_CMPLX_ARG (pz), prz, tmp_info
1897 F77_CHAR_ARG_LEN (1)));
1898
1899 info = tmp_info;
1900
1901 if (info != 0)
1902 info = -2;
1903
1904 if (is_singular (rcon) || octave::math::isnan (rcon))
1905 {
1906 if (sing_handler)
1907 sing_handler (rcon);
1908 else
1909 octave::warn_singular_matrix (rcon);
1910 }
1911 }
1912
1913 if (info == 0)
1914 {
1915 retval = b;
1916 FloatComplex *result = retval.rwdata ();
1917
1918 char job = 'N';
1919 F77_XFCN (cgetrs, CGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1920 nr, b_nc, F77_CMPLX_ARG (tmp_data), nr,
1921 pipvt, F77_CMPLX_ARG (result), b_nr, tmp_info
1922 F77_CHAR_ARG_LEN (1)));
1923
1924 info = tmp_info;
1925 }
1926 else
1927 mattype.mark_as_rectangular ();
1928 }
1929 }
1930
1931 if (octave::math::isinf (anorm))
1932 {
1933 retval = FloatComplexMatrix (b_nr, b_nc,
1934 FloatComplex (0, 0));
1935 mattype.mark_as_full ();
1936 }
1937 }
1938
1939 return retval;
1940}
1941
1942FloatComplexMatrix
1943FloatComplexMatrix::solve (MatrixType& mattype, const FloatMatrix& b) const
1944{
1945 octave_idx_type info;
1946 float rcon;
1947 return solve (mattype, b, info, rcon, nullptr);
1948}
1949
1950FloatComplexMatrix
1951FloatComplexMatrix::solve (MatrixType& mattype, const FloatMatrix& b,
1952 octave_idx_type& info) const
1953{
1954 float rcon;
1955 return solve (mattype, b, info, rcon, nullptr);
1956}
1957
1958FloatComplexMatrix
1959FloatComplexMatrix::solve (MatrixType& mattype, const FloatMatrix& b,
1960 octave_idx_type& info, float& rcon) const
1961{
1962 return solve (mattype, b, info, rcon, nullptr);
1963}
1964
1965FloatComplexMatrix
1966FloatComplexMatrix::solve (MatrixType& mattype, const FloatMatrix& b,
1967 octave_idx_type& info, float& rcon,
1968 solve_singularity_handler sing_handler,
1969 bool singular_fallback, blas_trans_type transt) const
1970{
1971 FloatComplexMatrix tmp (b);
1972 return solve (mattype, tmp, info, rcon, sing_handler, singular_fallback,
1973 transt);
1974}
1975
1976FloatComplexMatrix
1977FloatComplexMatrix::solve (MatrixType& mattype,
1978 const FloatComplexMatrix& b) const
1979{
1980 octave_idx_type info;
1981 float rcon;
1982 return solve (mattype, b, info, rcon, nullptr);
1983}
1984
1985FloatComplexMatrix
1986FloatComplexMatrix::solve (MatrixType& mattype, const FloatComplexMatrix& b,
1987 octave_idx_type& info) const
1988{
1989 float rcon;
1990 return solve (mattype, b, info, rcon, nullptr);
1991}
1992
1993FloatComplexMatrix
1994FloatComplexMatrix::solve (MatrixType& mattype, const FloatComplexMatrix& b,
1995 octave_idx_type& info, float& rcon) const
1996{
1997 return solve (mattype, b, info, rcon, nullptr);
1998}
1999
2000FloatComplexMatrix
2001FloatComplexMatrix::solve (MatrixType& mattype, const FloatComplexMatrix& b,
2002 octave_idx_type& info, float& rcon,
2003 solve_singularity_handler sing_handler,
2004 bool singular_fallback,
2005 blas_trans_type transt) const
2006{
2007 FloatComplexMatrix retval;
2008 int typ = mattype.type ();
2009
2010 if (typ == MatrixType::Unknown)
2011 typ = mattype.type (*this);
2012
2013 // Only calculate the condition number for LU/Cholesky
2014 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
2015 retval = utsolve (mattype, b, info, rcon, sing_handler, true, transt);
2016 else if (typ == MatrixType::Lower
2017 || typ == MatrixType::Permuted_Lower)
2018 retval = ltsolve (mattype, b, info, rcon, sing_handler, true, transt);
2019 else if (transt == blas_trans)
2020 return transpose ().solve (mattype, b, info, rcon, sing_handler,
2021 singular_fallback);
2022 else if (transt == blas_conj_trans)
2023 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler,
2024 singular_fallback);
2025 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
2026 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
2027 else if (typ != MatrixType::Rectangular)
2028 (*current_liboctave_error_handler) ("unknown matrix type");
2029
2030 // Rectangular or one of the above solvers flags a singular matrix
2031 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
2032 {
2033 octave_idx_type rank;
2034 retval = lssolve (b, info, rank, rcon);
2035 }
2036
2037 return retval;
2038}
2039
2040FloatComplexColumnVector
2041FloatComplexMatrix::solve (MatrixType& mattype,
2042 const FloatColumnVector& b) const
2043{
2044 octave_idx_type info;
2045 float rcon;
2046 return solve (mattype, FloatComplexColumnVector (b), info, rcon, nullptr);
2047}
2048
2049FloatComplexColumnVector
2050FloatComplexMatrix::solve (MatrixType& mattype, const FloatColumnVector& b,
2051 octave_idx_type& info) const
2052{
2053 float rcon;
2054 return solve (mattype, FloatComplexColumnVector (b), info, rcon, nullptr);
2055}
2056
2057FloatComplexColumnVector
2058FloatComplexMatrix::solve (MatrixType& mattype, const FloatColumnVector& b,
2059 octave_idx_type& info, float& rcon) const
2060{
2061 return solve (mattype, FloatComplexColumnVector (b), info, rcon, nullptr);
2062}
2063
2064FloatComplexColumnVector
2065FloatComplexMatrix::solve (MatrixType& mattype, const FloatColumnVector& b,
2066 octave_idx_type& info, float& rcon,
2067 solve_singularity_handler sing_handler,
2068 blas_trans_type transt) const
2069{
2070 return solve (mattype, FloatComplexColumnVector (b), info, rcon,
2071 sing_handler, transt);
2072}
2073
2074FloatComplexColumnVector
2075FloatComplexMatrix::solve (MatrixType& mattype,
2076 const FloatComplexColumnVector& b) const
2077{
2078 octave_idx_type info;
2079 float rcon;
2080 return solve (mattype, b, info, rcon, nullptr);
2081}
2082
2083FloatComplexColumnVector
2084FloatComplexMatrix::solve (MatrixType& mattype,
2085 const FloatComplexColumnVector& b,
2086 octave_idx_type& info) const
2087{
2088 float rcon;
2089 return solve (mattype, b, info, rcon, nullptr);
2090}
2091
2092FloatComplexColumnVector
2093FloatComplexMatrix::solve (MatrixType& mattype,
2094 const FloatComplexColumnVector& b,
2095 octave_idx_type& info, float& rcon) const
2096{
2097 return solve (mattype, b, info, rcon, nullptr);
2098}
2099
2100FloatComplexColumnVector
2101FloatComplexMatrix::solve (MatrixType& mattype,
2102 const FloatComplexColumnVector& b,
2103 octave_idx_type& info, float& rcon,
2104 solve_singularity_handler sing_handler,
2105 blas_trans_type transt) const
2106{
2107
2108 FloatComplexMatrix tmp (b);
2109 tmp = solve (mattype, tmp, info, rcon, sing_handler, true, transt);
2110 return tmp.column (static_cast<octave_idx_type> (0));
2111}
2112
2113FloatComplexMatrix
2114FloatComplexMatrix::solve (const FloatMatrix& b) const
2115{
2116 octave_idx_type info;
2117 float rcon;
2118 return solve (b, info, rcon, nullptr);
2119}
2120
2121FloatComplexMatrix
2122FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info) const
2123{
2124 float rcon;
2125 return solve (b, info, rcon, nullptr);
2126}
2127
2128FloatComplexMatrix
2129FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info,
2130 float& rcon) const
2131{
2132 return solve (b, info, rcon, nullptr);
2133}
2134
2135FloatComplexMatrix
2136FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info,
2137 float& rcon,
2138 solve_singularity_handler sing_handler,
2139 blas_trans_type transt) const
2140{
2141 FloatComplexMatrix tmp (b);
2142 return solve (tmp, info, rcon, sing_handler, transt);
2143}
2144
2145FloatComplexMatrix
2146FloatComplexMatrix::solve (const FloatComplexMatrix& b) const
2147{
2148 octave_idx_type info;
2149 float rcon;
2150 return solve (b, info, rcon, nullptr);
2151}
2152
2153FloatComplexMatrix
2154FloatComplexMatrix::solve (const FloatComplexMatrix& b,
2155 octave_idx_type& info) const
2156{
2157 float rcon;
2158 return solve (b, info, rcon, nullptr);
2159}
2160
2161FloatComplexMatrix
2162FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info,
2163 float& rcon) const
2164{
2165 return solve (b, info, rcon, nullptr);
2166}
2167
2168FloatComplexMatrix
2169FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info,
2170 float& rcon,
2171 solve_singularity_handler sing_handler,
2172 blas_trans_type transt) const
2173{
2174 MatrixType mattype (*this);
2175 return solve (mattype, b, info, rcon, sing_handler, true, transt);
2176}
2177
2178FloatComplexColumnVector
2179FloatComplexMatrix::solve (const FloatColumnVector& b) const
2180{
2181 octave_idx_type info;
2182 float rcon;
2183 return solve (FloatComplexColumnVector (b), info, rcon, nullptr);
2184}
2185
2186FloatComplexColumnVector
2187FloatComplexMatrix::solve (const FloatColumnVector& b,
2188 octave_idx_type& info) const
2189{
2190 float rcon;
2191 return solve (FloatComplexColumnVector (b), info, rcon, nullptr);
2192}
2193
2194FloatComplexColumnVector
2195FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info,
2196 float& rcon) const
2197{
2198 return solve (FloatComplexColumnVector (b), info, rcon, nullptr);
2199}
2200
2201FloatComplexColumnVector
2202FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info,
2203 float& rcon,
2204 solve_singularity_handler sing_handler,
2205 blas_trans_type transt) const
2206{
2207 return solve (FloatComplexColumnVector (b), info, rcon, sing_handler, transt);
2208}
2209
2210FloatComplexColumnVector
2211FloatComplexMatrix::solve (const FloatComplexColumnVector& b) const
2212{
2213 octave_idx_type info;
2214 float rcon;
2215 return solve (b, info, rcon, nullptr);
2216}
2217
2218FloatComplexColumnVector
2219FloatComplexMatrix::solve (const FloatComplexColumnVector& b,
2220 octave_idx_type& info) const
2221{
2222 float rcon;
2223 return solve (b, info, rcon, nullptr);
2224}
2225
2226FloatComplexColumnVector
2227FloatComplexMatrix::solve (const FloatComplexColumnVector& b,
2228 octave_idx_type& info,
2229 float& rcon) const
2230{
2231 return solve (b, info, rcon, nullptr);
2232}
2233
2234FloatComplexColumnVector
2235FloatComplexMatrix::solve (const FloatComplexColumnVector& b,
2236 octave_idx_type& info,
2237 float& rcon,
2238 solve_singularity_handler sing_handler,
2239 blas_trans_type transt) const
2240{
2241 MatrixType mattype (*this);
2242 return solve (mattype, b, info, rcon, sing_handler, transt);
2243}
2244
2245FloatComplexMatrix
2246FloatComplexMatrix::lssolve (const FloatMatrix& b) const
2247{
2248 octave_idx_type info;
2249 octave_idx_type rank;
2250 float rcon;
2251 return lssolve (FloatComplexMatrix (b), info, rank, rcon);
2252}
2253
2254FloatComplexMatrix
2255FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info) const
2256{
2257 octave_idx_type rank;
2258 float rcon;
2259 return lssolve (FloatComplexMatrix (b), info, rank, rcon);
2260}
2261
2262FloatComplexMatrix
2263FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info,
2264 octave_idx_type& rank) const
2265{
2266 float rcon;
2267 return lssolve (FloatComplexMatrix (b), info, rank, rcon);
2268}
2269
2270FloatComplexMatrix
2271FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info,
2272 octave_idx_type& rank, float& rcon) const
2273{
2274 return lssolve (FloatComplexMatrix (b), info, rank, rcon);
2275}
2276
2277FloatComplexMatrix
2278FloatComplexMatrix::lssolve (const FloatComplexMatrix& b) const
2279{
2280 octave_idx_type info;
2281 octave_idx_type rank;
2282 float rcon;
2283 return lssolve (b, info, rank, rcon);
2284}
2285
2286FloatComplexMatrix
2287FloatComplexMatrix::lssolve (const FloatComplexMatrix& b,
2288 octave_idx_type& info) const
2289{
2290 octave_idx_type rank;
2291 float rcon;
2292 return lssolve (b, info, rank, rcon);
2293}
2294
2295FloatComplexMatrix
2296FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info,
2297 octave_idx_type& rank) const
2298{
2299 float rcon;
2300 return lssolve (b, info, rank, rcon);
2301}
2302
2303FloatComplexMatrix
2304FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info,
2305 octave_idx_type& rank, float& rcon) const
2306{
2307 FloatComplexMatrix retval;
2308
2309 F77_INT m = octave::to_f77_int (rows ());
2310 F77_INT n = octave::to_f77_int (cols ());
2311
2312 F77_INT b_nr = octave::to_f77_int (b.rows ());
2313 F77_INT b_nc = octave::to_f77_int (b.cols ());
2314 F77_INT nrhs = b_nc; // alias for code readability
2315
2316 if (m != b_nr)
2317 (*current_liboctave_error_handler)
2318 ("matrix dimension mismatch solution of linear equations");
2319
2320 if (m == 0 || n == 0 || b_nc == 0)
2321 retval = FloatComplexMatrix (n, b_nc, FloatComplex (0.0, 0.0));
2322 else
2323 {
2324 F77_INT minmn = (m < n ? m : n);
2325 F77_INT maxmn = (m > n ? m : n);
2326 rcon = -1.0;
2327
2328 if (m != n)
2329 {
2330 retval = FloatComplexMatrix (maxmn, nrhs);
2331
2332 for (F77_INT j = 0; j < nrhs; j++)
2333 for (F77_INT i = 0; i < m; i++)
2334 retval.elem (i, j) = b.elem (i, j);
2335 }
2336 else
2337 retval = b;
2338
2339 FloatComplexMatrix atmp = *this;
2340 FloatComplex *tmp_data = atmp.rwdata ();
2341
2342 FloatComplex *pretval = retval.rwdata ();
2343 Array<float> s (dim_vector (minmn, 1));
2344 float *ps = s.rwdata ();
2345
2346 // Ask ZGELSD what the dimension of WORK should be.
2347 F77_INT lwork = -1;
2348
2349 Array<FloatComplex> work (dim_vector (1, 1));
2350
2351 F77_INT smlsiz;
2352 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("CGELSD", 6),
2353 F77_CONST_CHAR_ARG2 (" ", 1),
2354 0, 0, 0, 0, smlsiz
2355 F77_CHAR_ARG_LEN (6)
2356 F77_CHAR_ARG_LEN (1));
2357
2358 F77_INT mnthr;
2359 F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("CGELSD", 6),
2360 F77_CONST_CHAR_ARG2 (" ", 1),
2361 m, n, nrhs, -1, mnthr
2362 F77_CHAR_ARG_LEN (6)
2363 F77_CHAR_ARG_LEN (1));
2364
2365 // We compute the size of rwork and iwork because ZGELSD in
2366 // older versions of LAPACK does not return them on a query
2367 // call.
2368 float dminmn = static_cast<float> (minmn);
2369 float dsmlsizp1 = static_cast<float> (smlsiz+1);
2370 float tmp = octave::math::log2 (dminmn / dsmlsizp1);
2371
2372 F77_INT nlvl = static_cast<F77_INT> (tmp) + 1;
2373 if (nlvl < 0)
2374 nlvl = 0;
2375
2376 F77_INT lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
2377 + 3*smlsiz*nrhs
2378 + std::max ((smlsiz+1)*(smlsiz+1),
2379 n*(1+nrhs) + 2*nrhs);
2380 if (lrwork < 1)
2381 lrwork = 1;
2382 Array<float> rwork (dim_vector (lrwork, 1));
2383 float *prwork = rwork.rwdata ();
2384
2385 F77_INT liwork = 3 * minmn * nlvl + 11 * minmn;
2386 if (liwork < 1)
2387 liwork = 1;
2388 Array<F77_INT> iwork (dim_vector (liwork, 1));
2389 F77_INT *piwork = iwork.rwdata ();
2390
2391 F77_INT tmp_info = 0;
2392 F77_INT tmp_rank = 0;
2393
2394 F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, F77_CMPLX_ARG (tmp_data), m,
2395 F77_CMPLX_ARG (pretval), maxmn,
2396 ps, rcon, tmp_rank, F77_CMPLX_ARG (work.rwdata ()),
2397 lwork, prwork, piwork, tmp_info));
2398
2399 info = tmp_info;
2400 rank = tmp_rank;
2401
2402 // The workspace query is broken in at least LAPACK 3.0.0
2403 // through 3.1.1 when n >= mnthr. The obtuse formula below
2404 // should provide sufficient workspace for ZGELSD to operate
2405 // efficiently.
2406 if (n > m && n >= mnthr)
2407 {
2408 F77_INT addend = m;
2409
2410 if (2*m-4 > addend)
2411 addend = 2*m-4;
2412
2413 if (nrhs > addend)
2414 addend = nrhs;
2415
2416 if (n-3*m > addend)
2417 addend = n-3*m;
2418
2419 const F77_INT lworkaround = 4*m + m*m + addend;
2420
2421 if (std::real (work(0)) < lworkaround)
2422 work(0) = lworkaround;
2423 }
2424 else if (m >= n)
2425 {
2426 F77_INT lworkaround = 2*m + m*nrhs;
2427
2428 if (std::real (work(0)) < lworkaround)
2429 work(0) = lworkaround;
2430 }
2431
2432 lwork = static_cast<F77_INT> (std::real (work(0)));
2433 work.resize (dim_vector (lwork, 1));
2434
2435 float anorm = norm1 (*this);
2436
2437 if (octave::math::isinf (anorm))
2438 {
2439 rcon = 0.0;
2440 retval = FloatComplexMatrix (n, b_nc, 0.0);
2441 }
2442 else if (octave::math::isnan (anorm))
2443 {
2444 rcon = octave::numeric_limits<float>::NaN ();
2445 retval = FloatComplexMatrix (n, b_nc,
2446 octave::numeric_limits<float>::NaN ());
2447 }
2448 else
2449 {
2450 F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, F77_CMPLX_ARG (tmp_data),
2451 m, F77_CMPLX_ARG (pretval),
2452 maxmn, ps, rcon, tmp_rank,
2453 F77_CMPLX_ARG (work.rwdata ()),
2454 lwork, prwork, piwork, tmp_info));
2455
2456 info = tmp_info;
2457 rank = tmp_rank;
2458
2459 if (s.elem (0) == 0.0)
2460 rcon = 0.0;
2461 else
2462 rcon = s.elem (minmn - 1) / s.elem (0);
2463
2464 retval.resize (n, nrhs);
2465 }
2466 }
2467
2468 return retval;
2469}
2470
2471FloatComplexColumnVector
2472FloatComplexMatrix::lssolve (const FloatColumnVector& b) const
2473{
2474 octave_idx_type info;
2475 octave_idx_type rank;
2476 float rcon;
2477 return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
2478}
2479
2480FloatComplexColumnVector
2481FloatComplexMatrix::lssolve (const FloatColumnVector& b,
2482 octave_idx_type& info) const
2483{
2484 octave_idx_type rank;
2485 float rcon;
2486 return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
2487}
2488
2489FloatComplexColumnVector
2490FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info,
2491 octave_idx_type& rank) const
2492{
2493 float rcon;
2494 return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
2495}
2496
2497FloatComplexColumnVector
2498FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info,
2499 octave_idx_type& rank, float& rcon) const
2500{
2501 return lssolve (FloatComplexColumnVector (b), info, rank, rcon);
2502}
2503
2504FloatComplexColumnVector
2505FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b) const
2506{
2507 octave_idx_type info;
2508 octave_idx_type rank;
2509 float rcon;
2510 return lssolve (b, info, rank, rcon);
2511}
2512
2513FloatComplexColumnVector
2514FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b,
2515 octave_idx_type& info) const
2516{
2517 octave_idx_type rank;
2518 float rcon;
2519 return lssolve (b, info, rank, rcon);
2520}
2521
2522FloatComplexColumnVector
2523FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b,
2524 octave_idx_type& info,
2525 octave_idx_type& rank) const
2526{
2527 float rcon;
2528 return lssolve (b, info, rank, rcon);
2529
2530}
2531
2532FloatComplexColumnVector
2533FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b,
2534 octave_idx_type& info,
2535 octave_idx_type& rank, float& rcon) const
2536{
2537 FloatComplexColumnVector retval;
2538
2539 F77_INT nrhs = 1;
2540
2541 F77_INT m = octave::to_f77_int (rows ());
2542 F77_INT n = octave::to_f77_int (cols ());
2543
2544 F77_INT b_nel = octave::to_f77_int (b.numel ());
2545
2546 if (m != b_nel)
2547 (*current_liboctave_error_handler)
2548 ("matrix dimension mismatch solution of linear equations");
2549
2550 if (m == 0 || n == 0)
2551 retval = FloatComplexColumnVector (n, FloatComplex (0.0, 0.0));
2552 else
2553 {
2554 F77_INT minmn = (m < n ? m : n);
2555 F77_INT maxmn = (m > n ? m : n);
2556 rcon = -1.0;
2557
2558 if (m != n)
2559 {
2560 retval = FloatComplexColumnVector (maxmn);
2561
2562 for (F77_INT i = 0; i < m; i++)
2563 retval.elem (i) = b.elem (i);
2564 }
2565 else
2566 retval = b;
2567
2568 FloatComplexMatrix atmp = *this;
2569 FloatComplex *tmp_data = atmp.rwdata ();
2570
2571 FloatComplex *pretval = retval.rwdata ();
2572 Array<float> s (dim_vector (minmn, 1));
2573 float *ps = s.rwdata ();
2574
2575 // Ask ZGELSD what the dimension of WORK should be.
2576 F77_INT lwork = -1;
2577
2578 Array<FloatComplex> work (dim_vector (1, 1));
2579
2580 F77_INT smlsiz;
2581 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("CGELSD", 6),
2582 F77_CONST_CHAR_ARG2 (" ", 1),
2583 0, 0, 0, 0, smlsiz
2584 F77_CHAR_ARG_LEN (6)
2585 F77_CHAR_ARG_LEN (1));
2586
2587 // We compute the size of rwork and iwork because ZGELSD in
2588 // older versions of LAPACK does not return them on a query
2589 // call.
2590 float dminmn = static_cast<float> (minmn);
2591 float dsmlsizp1 = static_cast<float> (smlsiz+1);
2592 float tmp = octave::math::log2 (dminmn / dsmlsizp1);
2593
2594 F77_INT nlvl = static_cast<F77_INT> (tmp) + 1;
2595 if (nlvl < 0)
2596 nlvl = 0;
2597
2598 F77_INT lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
2599 + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
2600 if (lrwork < 1)
2601 lrwork = 1;
2602 Array<float> rwork (dim_vector (lrwork, 1));
2603 float *prwork = rwork.rwdata ();
2604
2605 F77_INT liwork = 3 * minmn * nlvl + 11 * minmn;
2606 if (liwork < 1)
2607 liwork = 1;
2608 Array<F77_INT> iwork (dim_vector (liwork, 1));
2609 F77_INT *piwork = iwork.rwdata ();
2610
2611 F77_INT tmp_info = 0;
2612 F77_INT tmp_rank = 0;
2613
2614 F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, F77_CMPLX_ARG (tmp_data), m,
2615 F77_CMPLX_ARG (pretval), maxmn,
2616 ps, rcon, tmp_rank, F77_CMPLX_ARG (work.rwdata ()),
2617 lwork, prwork, piwork, tmp_info));
2618
2619 info = tmp_info;
2620 rank = tmp_rank;
2621
2622 lwork = static_cast<F77_INT> (std::real (work(0)));
2623 work.resize (dim_vector (lwork, 1));
2624 rwork.resize (dim_vector (static_cast<F77_INT> (rwork(0)), 1));
2625 iwork.resize (dim_vector (iwork(0), 1));
2626
2627 F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, F77_CMPLX_ARG (tmp_data), m,
2628 F77_CMPLX_ARG (pretval),
2629 maxmn, ps, rcon, tmp_rank,
2630 F77_CMPLX_ARG (work.rwdata ()), lwork,
2631 prwork, piwork, tmp_info));
2632
2633 info = tmp_info;
2634 rank = tmp_rank;
2635
2636 if (rank < minmn)
2637 {
2638 if (s.elem (0) == 0.0)
2639 rcon = 0.0;
2640 else
2641 rcon = s.elem (minmn - 1) / s.elem (0);
2642
2643 retval.resize (n);
2644 }
2645 }
2646
2647 return retval;
2648}
2649
2650// column vector by row vector -> matrix operations
2651
2652FloatComplexMatrix
2653operator * (const FloatColumnVector& v, const FloatComplexRowVector& a)
2654{
2655 FloatComplexColumnVector tmp (v);
2656 return tmp * a;
2657}
2658
2659FloatComplexMatrix
2660operator * (const FloatComplexColumnVector& a, const FloatRowVector& b)
2661{
2662 FloatComplexRowVector tmp (b);
2663 return a * tmp;
2664}
2665
2666FloatComplexMatrix
2667operator * (const FloatComplexColumnVector& v, const FloatComplexRowVector& a)
2668{
2669 FloatComplexMatrix retval;
2670
2671 F77_INT len = octave::to_f77_int (v.numel ());
2672
2673 if (len != 0)
2674 {
2675 F77_INT a_len = octave::to_f77_int (a.numel ());
2676
2677 retval = FloatComplexMatrix (len, a_len);
2678 FloatComplex *c = retval.rwdata ();
2679
2680 F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
2681 F77_CONST_CHAR_ARG2 ("N", 1),
2682 len, a_len, 1, 1.0, F77_CONST_CMPLX_ARG (v.data ()), len,
2683 F77_CONST_CMPLX_ARG (a.data ()), 1, 0.0, F77_CMPLX_ARG (c), len
2684 F77_CHAR_ARG_LEN (1)
2685 F77_CHAR_ARG_LEN (1)));
2686 }
2687
2688 return retval;
2689}
2690
2691// matrix by diagonal matrix -> matrix operations
2692
2693FloatComplexMatrix&
2694FloatComplexMatrix::operator += (const FloatDiagMatrix& a)
2695{
2696 octave_idx_type nr = rows ();
2697 octave_idx_type nc = cols ();
2698
2699 octave_idx_type a_nr = a.rows ();
2700 octave_idx_type a_nc = a.cols ();
2701
2702 if (nr != a_nr || nc != a_nc)
2703 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2704
2705 for (octave_idx_type i = 0; i < a.length (); i++)
2706 elem (i, i) += a.elem (i, i);
2707
2708 return *this;
2709}
2710
2711FloatComplexMatrix&
2712FloatComplexMatrix::operator -= (const FloatDiagMatrix& a)
2713{
2714 octave_idx_type nr = rows ();
2715 octave_idx_type nc = cols ();
2716
2717 octave_idx_type a_nr = a.rows ();
2718 octave_idx_type a_nc = a.cols ();
2719
2720 if (nr != a_nr || nc != a_nc)
2721 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2722
2723 for (octave_idx_type i = 0; i < a.length (); i++)
2724 elem (i, i) -= a.elem (i, i);
2725
2726 return *this;
2727}
2728
2729FloatComplexMatrix&
2730FloatComplexMatrix::operator += (const FloatComplexDiagMatrix& a)
2731{
2732 octave_idx_type nr = rows ();
2733 octave_idx_type nc = cols ();
2734
2735 octave_idx_type a_nr = a.rows ();
2736 octave_idx_type a_nc = a.cols ();
2737
2738 if (nr != a_nr || nc != a_nc)
2739 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2740
2741 for (octave_idx_type i = 0; i < a.length (); i++)
2742 elem (i, i) += a.elem (i, i);
2743
2744 return *this;
2745}
2746
2747FloatComplexMatrix&
2748FloatComplexMatrix::operator -= (const FloatComplexDiagMatrix& a)
2749{
2750 octave_idx_type nr = rows ();
2751 octave_idx_type nc = cols ();
2752
2753 octave_idx_type a_nr = a.rows ();
2754 octave_idx_type a_nc = a.cols ();
2755
2756 if (nr != a_nr || nc != a_nc)
2757 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2758
2759 for (octave_idx_type i = 0; i < a.length (); i++)
2760 elem (i, i) -= a.elem (i, i);
2761
2762 return *this;
2763}
2764
2765// matrix by matrix -> matrix operations
2766
2767FloatComplexMatrix&
2768FloatComplexMatrix::operator += (const FloatMatrix& a)
2769{
2770 octave_idx_type nr = rows ();
2771 octave_idx_type nc = cols ();
2772
2773 octave_idx_type a_nr = a.rows ();
2774 octave_idx_type a_nc = a.cols ();
2775
2776 if (nr != a_nr || nc != a_nc)
2777 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2778
2779 if (nr == 0 || nc == 0)
2780 return *this;
2781
2782 FloatComplex *d = rwdata (); // Ensures only 1 reference to my privates!
2783
2784 mx_inline_add2 (numel (), d, a.data ());
2785 return *this;
2786}
2787
2788FloatComplexMatrix&
2789FloatComplexMatrix::operator -= (const FloatMatrix& a)
2790{
2791 octave_idx_type nr = rows ();
2792 octave_idx_type nc = cols ();
2793
2794 octave_idx_type a_nr = a.rows ();
2795 octave_idx_type a_nc = a.cols ();
2796
2797 if (nr != a_nr || nc != a_nc)
2798 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2799
2800 if (nr == 0 || nc == 0)
2801 return *this;
2802
2803 FloatComplex *d = rwdata (); // Ensures only 1 reference to my privates!
2804
2805 mx_inline_sub2 (numel (), d, a.data ());
2806 return *this;
2807}
2808
2809// unary operations
2810
2811boolMatrix
2812FloatComplexMatrix::all (int dim) const
2813{
2814 return FloatComplexNDArray::all (dim);
2815}
2816
2817boolMatrix
2818FloatComplexMatrix::any (int dim) const
2819{
2820 return FloatComplexNDArray::any (dim);
2821}
2822
2823FloatComplexMatrix
2824FloatComplexMatrix::cumprod (int dim) const
2825{
2826 return FloatComplexNDArray::cumprod (dim);
2827}
2828
2829FloatComplexMatrix
2830FloatComplexMatrix::cumsum (int dim) const
2831{
2832 return FloatComplexNDArray::cumsum (dim);
2833}
2834
2835FloatComplexMatrix
2836FloatComplexMatrix::prod (int dim) const
2837{
2838 return FloatComplexNDArray::prod (dim);
2839}
2840
2841FloatComplexMatrix
2842FloatComplexMatrix::sum (int dim) const
2843{
2844 return FloatComplexNDArray::sum (dim);
2845}
2846
2847FloatComplexMatrix
2848FloatComplexMatrix::sumsq (int dim) const
2849{
2850 return FloatComplexNDArray::sumsq (dim);
2851}
2852
2853FloatMatrix
2854FloatComplexMatrix::abs () const
2855{
2856 return FloatComplexNDArray::abs ();
2857}
2858
2859FloatComplexMatrix
2864
2865FloatComplexDiagMatrix
2866FloatComplexMatrix::diag (octave_idx_type m, octave_idx_type n) const
2867{
2868 FloatComplexDiagMatrix retval;
2869
2870 octave_idx_type nr = rows ();
2871 octave_idx_type nc = cols ();
2872
2873 if (nr == 1 || nc == 1)
2874 retval = FloatComplexDiagMatrix (*this, m, n);
2875 else
2876 (*current_liboctave_error_handler) ("diag: expecting vector argument");
2877
2878 return retval;
2879}
2880
2881bool
2882FloatComplexMatrix::row_is_real_only (octave_idx_type i) const
2883{
2884 bool retval = true;
2885
2886 octave_idx_type nc = columns ();
2887
2888 for (octave_idx_type j = 0; j < nc; j++)
2889 {
2890 if (std::imag (elem (i, j)) != 0.0)
2891 {
2892 retval = false;
2893 break;
2894 }
2895 }
2896
2897 return retval;
2898}
2899
2900bool
2901FloatComplexMatrix::column_is_real_only (octave_idx_type j) const
2902{
2903 bool retval = true;
2904
2905 octave_idx_type nr = rows ();
2906
2907 for (octave_idx_type i = 0; i < nr; i++)
2908 {
2909 if (std::imag (elem (i, j)) != 0.0)
2910 {
2911 retval = false;
2912 break;
2913 }
2914 }
2915
2916 return retval;
2917}
2918
2919FloatComplexColumnVector
2920FloatComplexMatrix::row_min () const
2921{
2922 Array<octave_idx_type> dummy_idx;
2923 return row_min (dummy_idx);
2924}
2925
2926FloatComplexColumnVector
2927FloatComplexMatrix::row_min (Array<octave_idx_type>& idx_arg) const
2928{
2929 FloatComplexColumnVector result;
2930
2931 octave_idx_type nr = rows ();
2932 octave_idx_type nc = cols ();
2933
2934 if (nr > 0 && nc > 0)
2935 {
2936 result.resize (nr);
2937 idx_arg.resize (dim_vector (nr, 1));
2938
2939 for (octave_idx_type i = 0; i < nr; i++)
2940 {
2941 bool real_only = row_is_real_only (i);
2942
2943 octave_idx_type idx_j;
2944
2945 FloatComplex tmp_min;
2946
2947 float abs_min = octave::numeric_limits<float>::NaN ();
2948
2949 for (idx_j = 0; idx_j < nc; idx_j++)
2950 {
2951 tmp_min = elem (i, idx_j);
2952
2953 if (! octave::math::isnan (tmp_min))
2954 {
2955 abs_min = (real_only ? tmp_min.real ()
2956 : std::abs (tmp_min));
2957 break;
2958 }
2959 }
2960
2961 for (octave_idx_type j = idx_j+1; j < nc; j++)
2962 {
2963 FloatComplex tmp = elem (i, j);
2964
2965 if (octave::math::isnan (tmp))
2966 continue;
2967
2968 float abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
2969
2970 if (abs_tmp < abs_min)
2971 {
2972 idx_j = j;
2973 tmp_min = tmp;
2974 abs_min = abs_tmp;
2975 }
2976 }
2977
2978 if (octave::math::isnan (tmp_min))
2979 {
2980 result.elem (i) = FloatComplex_NaN_result;
2981 idx_arg.elem (i) = 0;
2982 }
2983 else
2984 {
2985 result.elem (i) = tmp_min;
2986 idx_arg.elem (i) = idx_j;
2987 }
2988 }
2989 }
2990
2991 return result;
2992}
2993
2994FloatComplexColumnVector
2995FloatComplexMatrix::row_max () const
2996{
2997 Array<octave_idx_type> dummy_idx;
2998 return row_max (dummy_idx);
2999}
3000
3001FloatComplexColumnVector
3002FloatComplexMatrix::row_max (Array<octave_idx_type>& idx_arg) const
3003{
3004 FloatComplexColumnVector result;
3005
3006 octave_idx_type nr = rows ();
3007 octave_idx_type nc = cols ();
3008
3009 if (nr > 0 && nc > 0)
3010 {
3011 result.resize (nr);
3012 idx_arg.resize (dim_vector (nr, 1));
3013
3014 for (octave_idx_type i = 0; i < nr; i++)
3015 {
3016 bool real_only = row_is_real_only (i);
3017
3018 octave_idx_type idx_j;
3019
3020 FloatComplex tmp_max;
3021
3022 float abs_max = octave::numeric_limits<float>::NaN ();
3023
3024 for (idx_j = 0; idx_j < nc; idx_j++)
3025 {
3026 tmp_max = elem (i, idx_j);
3027
3028 if (! octave::math::isnan (tmp_max))
3029 {
3030 abs_max = (real_only ? tmp_max.real ()
3031 : std::abs (tmp_max));
3032 break;
3033 }
3034 }
3035
3036 for (octave_idx_type j = idx_j+1; j < nc; j++)
3037 {
3038 FloatComplex tmp = elem (i, j);
3039
3040 if (octave::math::isnan (tmp))
3041 continue;
3042
3043 float abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
3044
3045 if (abs_tmp > abs_max)
3046 {
3047 idx_j = j;
3048 tmp_max = tmp;
3049 abs_max = abs_tmp;
3050 }
3051 }
3052
3053 if (octave::math::isnan (tmp_max))
3054 {
3055 result.elem (i) = FloatComplex_NaN_result;
3056 idx_arg.elem (i) = 0;
3057 }
3058 else
3059 {
3060 result.elem (i) = tmp_max;
3061 idx_arg.elem (i) = idx_j;
3062 }
3063 }
3064 }
3065
3066 return result;
3067}
3068
3069FloatComplexRowVector
3070FloatComplexMatrix::column_min () const
3071{
3072 Array<octave_idx_type> dummy_idx;
3073 return column_min (dummy_idx);
3074}
3075
3076FloatComplexRowVector
3077FloatComplexMatrix::column_min (Array<octave_idx_type>& idx_arg) const
3078{
3079 FloatComplexRowVector result;
3080
3081 octave_idx_type nr = rows ();
3082 octave_idx_type nc = cols ();
3083
3084 if (nr > 0 && nc > 0)
3085 {
3086 result.resize (nc);
3087 idx_arg.resize (dim_vector (1, nc));
3088
3089 for (octave_idx_type j = 0; j < nc; j++)
3090 {
3091 bool real_only = column_is_real_only (j);
3092
3093 octave_idx_type idx_i;
3094
3095 FloatComplex tmp_min;
3096
3097 float abs_min = octave::numeric_limits<float>::NaN ();
3098
3099 for (idx_i = 0; idx_i < nr; idx_i++)
3100 {
3101 tmp_min = elem (idx_i, j);
3102
3103 if (! octave::math::isnan (tmp_min))
3104 {
3105 abs_min = (real_only ? tmp_min.real ()
3106 : std::abs (tmp_min));
3107 break;
3108 }
3109 }
3110
3111 for (octave_idx_type i = idx_i+1; i < nr; i++)
3112 {
3113 FloatComplex tmp = elem (i, j);
3114
3115 if (octave::math::isnan (tmp))
3116 continue;
3117
3118 float abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
3119
3120 if (abs_tmp < abs_min)
3121 {
3122 idx_i = i;
3123 tmp_min = tmp;
3124 abs_min = abs_tmp;
3125 }
3126 }
3127
3128 if (octave::math::isnan (tmp_min))
3129 {
3130 result.elem (j) = FloatComplex_NaN_result;
3131 idx_arg.elem (j) = 0;
3132 }
3133 else
3134 {
3135 result.elem (j) = tmp_min;
3136 idx_arg.elem (j) = idx_i;
3137 }
3138 }
3139 }
3140
3141 return result;
3142}
3143
3144FloatComplexRowVector
3145FloatComplexMatrix::column_max () const
3146{
3147 Array<octave_idx_type> dummy_idx;
3148 return column_max (dummy_idx);
3149}
3150
3151FloatComplexRowVector
3152FloatComplexMatrix::column_max (Array<octave_idx_type>& idx_arg) const
3153{
3154 FloatComplexRowVector result;
3155
3156 octave_idx_type nr = rows ();
3157 octave_idx_type nc = cols ();
3158
3159 if (nr > 0 && nc > 0)
3160 {
3161 result.resize (nc);
3162 idx_arg.resize (dim_vector (1, nc));
3163
3164 for (octave_idx_type j = 0; j < nc; j++)
3165 {
3166 bool real_only = column_is_real_only (j);
3167
3168 octave_idx_type idx_i;
3169
3170 FloatComplex tmp_max;
3171
3172 float abs_max = octave::numeric_limits<float>::NaN ();
3173
3174 for (idx_i = 0; idx_i < nr; idx_i++)
3175 {
3176 tmp_max = elem (idx_i, j);
3177
3178 if (! octave::math::isnan (tmp_max))
3179 {
3180 abs_max = (real_only ? tmp_max.real ()
3181 : std::abs (tmp_max));
3182 break;
3183 }
3184 }
3185
3186 for (octave_idx_type i = idx_i+1; i < nr; i++)
3187 {
3188 FloatComplex tmp = elem (i, j);
3189
3190 if (octave::math::isnan (tmp))
3191 continue;
3192
3193 float abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
3194
3195 if (abs_tmp > abs_max)
3196 {
3197 idx_i = i;
3198 tmp_max = tmp;
3199 abs_max = abs_tmp;
3200 }
3201 }
3202
3203 if (octave::math::isnan (tmp_max))
3204 {
3205 result.elem (j) = FloatComplex_NaN_result;
3206 idx_arg.elem (j) = 0;
3207 }
3208 else
3209 {
3210 result.elem (j) = tmp_max;
3211 idx_arg.elem (j) = idx_i;
3212 }
3213 }
3214 }
3215
3216 return result;
3217}
3218
3219// i/o
3220
3221std::ostream&
3222operator << (std::ostream& os, const FloatComplexMatrix& a)
3223{
3224 for (octave_idx_type i = 0; i < a.rows (); i++)
3225 {
3226 for (octave_idx_type j = 0; j < a.cols (); j++)
3227 {
3228 os << ' ';
3229 octave::write_value<Complex> (os, a.elem (i, j));
3230 }
3231 os << "\n";
3232 }
3233 return os;
3234}
3235
3236std::istream&
3237operator >> (std::istream& is, FloatComplexMatrix& a)
3238{
3239 octave_idx_type nr = a.rows ();
3240 octave_idx_type nc = a.cols ();
3241
3242 if (nr > 0 && nc > 0)
3243 {
3244 FloatComplex tmp;
3245 for (octave_idx_type i = 0; i < nr; i++)
3246 for (octave_idx_type j = 0; j < nc; j++)
3247 {
3248 tmp = octave::read_value<FloatComplex> (is);
3249 if (is)
3250 a.elem (i, j) = tmp;
3251 else
3252 return is;
3253 }
3254 }
3255
3256 return is;
3257}
3258
3259FloatComplexMatrix
3260Givens (const FloatComplex& x, const FloatComplex& y)
3261{
3262 float cc;
3263 FloatComplex cs, temp_r;
3264
3265 F77_FUNC (clartg, CLARTG) (F77_CONST_CMPLX_ARG (&x), F77_CONST_CMPLX_ARG (&y),
3266 cc, F77_CMPLX_ARG (&cs), F77_CMPLX_ARG (&temp_r));
3267
3268 FloatComplexMatrix g (2, 2);
3269
3270 g.elem (0, 0) = cc;
3271 g.elem (1, 1) = cc;
3272 g.elem (0, 1) = cs;
3273 g.elem (1, 0) = -conj (cs);
3274
3275 return g;
3276}
3277
3278FloatComplexMatrix
3279Sylvester (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
3280 const FloatComplexMatrix& c)
3281{
3282 FloatComplexMatrix retval;
3283
3284 // FIXME: need to check that a, b, and c are all the same
3285 // size.
3286
3287 // Compute Schur decompositions
3288
3289 octave::math::schur<FloatComplexMatrix> as (a, "U");
3290 octave::math::schur<FloatComplexMatrix> bs (b, "U");
3291
3292 // Transform c to new coordinates.
3293
3294 FloatComplexMatrix ua = as.unitary_schur_matrix ();
3295 FloatComplexMatrix sch_a = as.schur_matrix ();
3296
3297 FloatComplexMatrix ub = bs.unitary_schur_matrix ();
3298 FloatComplexMatrix sch_b = bs.schur_matrix ();
3299
3300 FloatComplexMatrix cx = ua.hermitian () * c * ub;
3301
3302 // Solve the sylvester equation, back-transform, and return the
3303 // solution.
3304
3305 F77_INT a_nr = octave::to_f77_int (a.rows ());
3306 F77_INT b_nr = octave::to_f77_int (b.rows ());
3307
3308 float scale;
3309 F77_INT info;
3310
3311 FloatComplex *pa = sch_a.rwdata ();
3312 FloatComplex *pb = sch_b.rwdata ();
3313 FloatComplex *px = cx.rwdata ();
3314
3315 F77_XFCN (ctrsyl, CTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
3316 F77_CONST_CHAR_ARG2 ("N", 1),
3317 1, a_nr, b_nr, F77_CMPLX_ARG (pa), a_nr, F77_CMPLX_ARG (pb),
3318 b_nr, F77_CMPLX_ARG (px), a_nr, scale, info
3319 F77_CHAR_ARG_LEN (1)
3320 F77_CHAR_ARG_LEN (1)));
3321
3322 // FIXME: check info?
3323
3324 retval = ua * cx * ub.hermitian ();
3325
3326 return retval;
3327}
3328
3329FloatComplexMatrix
3330operator * (const FloatComplexMatrix& m, const FloatMatrix& a)
3331{
3332 if (m.columns () > std::min (m.rows (), a.columns ()) / 10)
3333 return FloatComplexMatrix (real (m) * a, imag (m) * a);
3334 else
3335 return m * FloatComplexMatrix (a);
3336}
3337
3338FloatComplexMatrix
3339operator * (const FloatMatrix& m, const FloatComplexMatrix& a)
3340{
3341 if (a.rows () > std::min (m.rows (), a.columns ()) / 10)
3342 return FloatComplexMatrix (m * real (a), m * imag (a));
3343 else
3344 return FloatComplexMatrix (m) * a;
3345}
3346
3347/*
3348
3349## Simple Dot Product, Matrix-Vector, and Matrix-Matrix Unit tests
3350%!assert (single ([1+i 2+i 3+i]) * single ([ 4+i ; 5+i ; 6+i]), single (29+21i), 5e-7)
3351%!assert (single ([1+i 2+i ; 3+i 4+i]) * single ([5+i ; 6+i]), single ([15 + 14i ; 37 + 18i]), 5e-7)
3352%!assert (single ([1+i 2+i ; 3+i 4+i ]) * single ([5+i 6+i ; 7+i 8+i]), single ([17 + 15i 20 + 17i; 41 + 19i 48 + 21i]), 5e-7)
3353%!assert (single ([1 i])*single ([i 0])', single (-i))
3354
3355## Test some simple identities
3356%!shared M, cv, rv
3357%! M = single (randn (10,10))+ i*single (rand (10,10));
3358%! cv = single (randn (10,1))+ i*single (rand (10,1));
3359%! rv = single (randn (1,10))+ i*single (rand (1,10));
3360%!assert ([M*cv,M*cv], M*[cv,cv], 5e-6)
3361%!assert ([M.'*cv,M.'*cv], M.'*[cv,cv], 5e-6)
3362%!assert ([M'*cv,M'*cv], M'*[cv,cv], 5e-6)
3363%!assert ([rv*M;rv*M], [rv;rv]*M, 5e-6)
3364%!assert ([rv*M.';rv*M.'], [rv;rv]*M.', 5e-6)
3365%!assert ([rv*M';rv*M'], [rv;rv]*M', 5e-6)
3366%!assert (2*rv*cv, [rv,rv]*[cv;cv], 5e-6)
3367
3368*/
3369
3370static char
3371get_blas_trans_arg (bool trans, bool conj)
3372{
3373 return trans ? (conj ? 'C' : 'T') : 'N';
3374}
3375
3376// the general GEMM operation
3377
3378FloatComplexMatrix
3379xgemm (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
3380 blas_trans_type transa, blas_trans_type transb)
3381{
3382 FloatComplexMatrix retval;
3383
3384 bool tra = transa != blas_no_trans;
3385 bool trb = transb != blas_no_trans;
3386 bool cja = transa == blas_conj_trans;
3387 bool cjb = transb == blas_conj_trans;
3388
3389 F77_INT a_nr = octave::to_f77_int (tra ? a.cols () : a.rows ());
3390 F77_INT a_nc = octave::to_f77_int (tra ? a.rows () : a.cols ());
3391
3392 F77_INT b_nr = octave::to_f77_int (trb ? b.cols () : b.rows ());
3393 F77_INT b_nc = octave::to_f77_int (trb ? b.rows () : b.cols ());
3394
3395 if (a_nc != b_nr)
3396 octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
3397
3398 if (a_nr == 0 || a_nc == 0 || b_nc == 0)
3399 retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
3400 else if (a.data () == b.data () && a_nr == b_nc && tra != trb)
3401 {
3402 F77_INT lda = octave::to_f77_int (a.rows ());
3403
3404 // FIXME: looking at the reference BLAS, it appears that it
3405 // should not be necessary to initialize the output matrix if
3406 // BETA is 0 in the call to CHERK, but ATLAS appears to
3407 // use the result matrix before zeroing the elements.
3408
3409 retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
3410 FloatComplex *c = retval.rwdata ();
3411
3412 const char ctra = get_blas_trans_arg (tra, cja);
3413 if (cja || cjb)
3414 {
3415 F77_XFCN (cherk, CHERK, (F77_CONST_CHAR_ARG2 ("U", 1),
3416 F77_CONST_CHAR_ARG2 (&ctra, 1),
3417 a_nr, a_nc, 1.0,
3418 F77_CONST_CMPLX_ARG (a.data ()), lda, 0.0, F77_CMPLX_ARG (c), a_nr
3419 F77_CHAR_ARG_LEN (1)
3420 F77_CHAR_ARG_LEN (1)));
3421 for (F77_INT j = 0; j < a_nr; j++)
3422 for (F77_INT i = 0; i < j; i++)
3423 retval.xelem (j, i) = octave::math::conj (retval.xelem (i, j));
3424 }
3425 else
3426 {
3427 F77_XFCN (csyrk, CSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),
3428 F77_CONST_CHAR_ARG2 (&ctra, 1),
3429 a_nr, a_nc, 1.0,
3430 F77_CONST_CMPLX_ARG (a.data ()), lda, 0.0, F77_CMPLX_ARG (c), a_nr
3431 F77_CHAR_ARG_LEN (1)
3432 F77_CHAR_ARG_LEN (1)));
3433 for (F77_INT j = 0; j < a_nr; j++)
3434 for (F77_INT i = 0; i < j; i++)
3435 retval.xelem (j, i) = retval.xelem (i, j);
3436
3437 }
3438
3439 }
3440 else
3441 {
3442 F77_INT lda = octave::to_f77_int (a.rows ());
3443 F77_INT tda = octave::to_f77_int (a.cols ());
3444 F77_INT ldb = octave::to_f77_int (b.rows ());
3445 F77_INT tdb = octave::to_f77_int (b.cols ());
3446
3447 retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
3448 FloatComplex *c = retval.rwdata ();
3449
3450 if (b_nc == 1 && a_nr == 1)
3451 {
3452 if (cja == cjb)
3453 {
3454 F77_FUNC (xcdotu, XCDOTU) (a_nc, F77_CONST_CMPLX_ARG (a.data ()), 1,
3455 F77_CONST_CMPLX_ARG (b.data ()), 1,
3456 F77_CMPLX_ARG (c));
3457 if (cja) *c = octave::math::conj (*c);
3458 }
3459 else if (cja)
3460 F77_FUNC (xcdotc, XCDOTC) (a_nc, F77_CONST_CMPLX_ARG (a.data ()), 1,
3461 F77_CONST_CMPLX_ARG (b.data ()), 1,
3462 F77_CMPLX_ARG (c));
3463 else
3464 F77_FUNC (xcdotc, XCDOTC) (a_nc, F77_CONST_CMPLX_ARG (b.data ()), 1,
3465 F77_CONST_CMPLX_ARG (a.data ()), 1,
3466 F77_CMPLX_ARG (c));
3467 }
3468 else if (b_nc == 1 && ! cjb)
3469 {
3470 const char ctra = get_blas_trans_arg (tra, cja);
3471 F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),
3472 lda, tda, 1.0, F77_CONST_CMPLX_ARG (a.data ()), lda,
3473 F77_CONST_CMPLX_ARG (b.data ()), 1, 0.0, F77_CMPLX_ARG (c), 1
3474 F77_CHAR_ARG_LEN (1)));
3475 }
3476 else if (a_nr == 1 && ! cja && ! cjb)
3477 {
3478 const char crevtrb = get_blas_trans_arg (! trb, cjb);
3479 F77_XFCN (cgemv, CGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),
3480 ldb, tdb, 1.0, F77_CONST_CMPLX_ARG (b.data ()), ldb,
3481 F77_CONST_CMPLX_ARG (a.data ()), 1, 0.0, F77_CMPLX_ARG (c), 1
3482 F77_CHAR_ARG_LEN (1)));
3483 }
3484 else
3485 {
3486 const char ctra = get_blas_trans_arg (tra, cja);
3487 const char ctrb = get_blas_trans_arg (trb, cjb);
3488 F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),
3489 F77_CONST_CHAR_ARG2 (&ctrb, 1),
3490 a_nr, b_nc, a_nc, 1.0, F77_CONST_CMPLX_ARG (a.data ()),
3491 lda, F77_CONST_CMPLX_ARG (b.data ()), ldb, 0.0, F77_CMPLX_ARG (c), a_nr
3492 F77_CHAR_ARG_LEN (1)
3493 F77_CHAR_ARG_LEN (1)));
3494 }
3495 }
3496
3497 return retval;
3498}
3499
3500FloatComplexMatrix
3501operator * (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
3502{
3503 return xgemm (a, b);
3504}
3505
3506// FIXME: it would be nice to share code among the min/max
3507// functions below.
3508
3509#define EMPTY_RETURN_CHECK(T) \
3510 if (nr == 0 || nc == 0) \
3511 return T (nr, nc);
3512
3513FloatComplexMatrix
3514min (const FloatComplex& c, const FloatComplexMatrix& m)
3515{
3516 octave_idx_type nr = m.rows ();
3517 octave_idx_type nc = m.columns ();
3518
3519 EMPTY_RETURN_CHECK (FloatComplexMatrix);
3520
3521 FloatComplexMatrix result (nr, nc);
3522
3523 for (octave_idx_type j = 0; j < nc; j++)
3524 for (octave_idx_type i = 0; i < nr; i++)
3525 {
3526 octave_quit ();
3527 result(i, j) = octave::math::min (c, m(i, j));
3528 }
3529
3530 return result;
3531}
3532
3533FloatComplexMatrix
3534min (const FloatComplexMatrix& m, const FloatComplex& c)
3535{
3536 return min (c, m);
3537}
3538
3539FloatComplexMatrix
3540min (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
3541{
3542 octave_idx_type nr = a.rows ();
3543 octave_idx_type nc = a.columns ();
3544
3545 if (nr != b.rows () || nc != b.columns ())
3546 (*current_liboctave_error_handler)
3547 ("two-arg min requires same size arguments");
3548
3549 EMPTY_RETURN_CHECK (FloatComplexMatrix);
3550
3551 FloatComplexMatrix result (nr, nc);
3552
3553 for (octave_idx_type j = 0; j < nc; j++)
3554 {
3555 bool columns_are_real_only = true;
3556 for (octave_idx_type i = 0; i < nr; i++)
3557 {
3558 octave_quit ();
3559 if (std::imag (a(i, j)) != 0.0 || std::imag (b(i, j)) != 0.0)
3560 {
3561 columns_are_real_only = false;
3562 break;
3563 }
3564 }
3565
3566 if (columns_are_real_only)
3567 {
3568 for (octave_idx_type i = 0; i < nr; i++)
3569 result(i, j) = octave::math::min (std::real (a(i, j)),
3570 std::real (b(i, j)));
3571 }
3572 else
3573 {
3574 for (octave_idx_type i = 0; i < nr; i++)
3575 {
3576 octave_quit ();
3577 result(i, j) = octave::math::min (a(i, j), b(i, j));
3578 }
3579 }
3580 }
3581
3582 return result;
3583}
3584
3585FloatComplexMatrix
3586max (const FloatComplex& c, const FloatComplexMatrix& m)
3587{
3588 octave_idx_type nr = m.rows ();
3589 octave_idx_type nc = m.columns ();
3590
3591 EMPTY_RETURN_CHECK (FloatComplexMatrix);
3592
3593 FloatComplexMatrix result (nr, nc);
3594
3595 for (octave_idx_type j = 0; j < nc; j++)
3596 for (octave_idx_type i = 0; i < nr; i++)
3597 {
3598 octave_quit ();
3599 result(i, j) = octave::math::max (c, m(i, j));
3600 }
3601
3602 return result;
3603}
3604
3605FloatComplexMatrix
3606max (const FloatComplexMatrix& m, const FloatComplex& c)
3607{
3608 return max (c, m);
3609}
3610
3611FloatComplexMatrix
3612max (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
3613{
3614 octave_idx_type nr = a.rows ();
3615 octave_idx_type nc = a.columns ();
3616
3617 if (nr != b.rows () || nc != b.columns ())
3618 (*current_liboctave_error_handler)
3619 ("two-arg max requires same size arguments");
3620
3621 EMPTY_RETURN_CHECK (FloatComplexMatrix);
3622
3623 FloatComplexMatrix result (nr, nc);
3624
3625 for (octave_idx_type j = 0; j < nc; j++)
3626 {
3627 bool columns_are_real_only = true;
3628 for (octave_idx_type i = 0; i < nr; i++)
3629 {
3630 octave_quit ();
3631 if (std::imag (a(i, j)) != 0.0 || std::imag (b(i, j)) != 0.0)
3632 {
3633 columns_are_real_only = false;
3634 break;
3635 }
3636 }
3637
3638 if (columns_are_real_only)
3639 {
3640 for (octave_idx_type i = 0; i < nr; i++)
3641 {
3642 octave_quit ();
3643 result(i, j) = octave::math::max (std::real (a(i, j)),
3644 std::real (b(i, j)));
3645 }
3646 }
3647 else
3648 {
3649 for (octave_idx_type i = 0; i < nr; i++)
3650 {
3651 octave_quit ();
3652 result(i, j) = octave::math::max (a(i, j), b(i, j));
3653 }
3654 }
3655 }
3656
3657 return result;
3658}
3659
3660FloatComplexMatrix
3661linspace (const FloatComplexColumnVector& x1,
3662 const FloatComplexColumnVector& x2,
3663 octave_idx_type n)
3664
3665{
3666 octave_idx_type m = x1.numel ();
3667
3668 if (x2.numel () != m)
3669 (*current_liboctave_error_handler)
3670 ("linspace: vectors must be of equal length");
3671
3672 FloatComplexMatrix retval;
3673
3674 if (n < 1)
3675 {
3676 retval.clear (m, 0);
3677 return retval;
3678 }
3679
3680 retval.clear (m, n);
3681 for (octave_idx_type i = 0; i < m; i++)
3682 retval.insert (linspace (x1(i), x2(i), n), i, 0);
3683
3684 return retval;
3685}
3686
3687MS_CMP_OPS (FloatComplexMatrix, FloatComplex)
3688MS_BOOL_OPS (FloatComplexMatrix, FloatComplex)
3689
3690SM_CMP_OPS (FloatComplex, FloatComplexMatrix)
3691SM_BOOL_OPS (FloatComplex, FloatComplexMatrix)
3692
3693MM_CMP_OPS (FloatComplexMatrix, FloatComplexMatrix)
3694MM_BOOL_OPS (FloatComplexMatrix, FloatComplexMatrix)
FloatComplexDET
base_det< FloatComplex > FloatComplexDET
Definition DET.h:88
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:525
Array::index
Array< T, Alloc > index(const octave::idx_vector &i) const
Indexing without resizing.
Definition Array-base.cc:718
Array::issquare
bool issquare() const
Size of the specified dimension.
Definition Array.h:649
Array::clear
void clear()
Definition Array-base.cc:110
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:563
Array::rows
octave_idx_type rows() const
Definition Array.h:463
Array::resize
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition Array-base.cc:1031
Array::insert
Array< T, Alloc > & insert(const Array< T, Alloc > &a, const Array< octave_idx_type > &idx)
Insert an array into another at a specified position.
Definition Array-base.cc:1631
Array::columns
octave_idx_type columns() const
Definition Array.h:475
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::make_unique
void make_unique()
Definition Array.h:220
Array::data
const T * data() const
Size of the specified dimension.
Definition Array.h:665
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
DiagArray2::rows
octave_idx_type rows() const
Definition DiagArray2.h:86
DiagArray2::length
octave_idx_type length() const
Definition DiagArray2.h:92
DiagArray2::elem
T elem(octave_idx_type r, octave_idx_type c) const
Definition DiagArray2.h:114
DiagArray2::cols
octave_idx_type cols() const
Definition DiagArray2.h:87
FloatColumnVector::extract
FloatColumnVector extract(octave_idx_type r1, octave_idx_type r2) const
Definition fColVector.cc:151
FloatComplexColumnVector::resize
void resize(octave_idx_type n, const FloatComplex &rfv=FloatComplex(0))
Definition fCColVector.h:152
FloatComplexMatrix::row_min
FloatComplexColumnVector row_min() const
Definition fCMatrix.cc:2920
FloatComplexMatrix::hermitian
FloatComplexMatrix hermitian() const
Definition fCMatrix.h:176
FloatComplexMatrix::ifourier
FloatComplexMatrix ifourier() const
Definition fCMatrix.cc:1086
FloatComplexMatrix::row_is_real_only
bool row_is_real_only(octave_idx_type) const
Definition fCMatrix.cc:2882
FloatComplexMatrix::all
boolMatrix all(int dim=-1) const
Definition fCMatrix.cc:2812
FloatComplexMatrix::sumsq
FloatComplexMatrix sumsq(int dim=-1) const
Definition fCMatrix.cc:2848
FloatComplexMatrix::pseudo_inverse
FloatComplexMatrix pseudo_inverse(float tol=0.0) const
Definition fCMatrix.cc:1012
FloatComplexMatrix::extract
FloatComplexMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition fCMatrix.cc:687
FloatComplexMatrix::lssolve
FloatComplexMatrix lssolve(const FloatMatrix &b) const
Definition fCMatrix.cc:2246
FloatComplexMatrix::operator+=
FloatComplexMatrix & operator+=(const FloatDiagMatrix &a)
Definition fCMatrix.cc:2694
FloatComplexMatrix::row
FloatComplexRowVector row(octave_idx_type i) const
Definition fCMatrix.cc:706
FloatComplexMatrix::operator!=
bool operator!=(const FloatComplexMatrix &a) const
Definition fCMatrix.cc:169
FloatComplexMatrix::rcond
float rcond() const
Definition fCMatrix.cc:1362
FloatComplexMatrix::column
FloatComplexColumnVector column(octave_idx_type i) const
Definition fCMatrix.cc:712
FloatComplexMatrix::fourier2d
FloatComplexMatrix fourier2d() const
Definition fCMatrix.cc:1115
FloatComplexMatrix::abs
FloatMatrix abs() const
Definition fCMatrix.cc:2854
FloatComplexMatrix::stack
FloatComplexMatrix stack(const FloatMatrix &a) const
Definition fCMatrix.cc:559
FloatComplexMatrix::cumsum
FloatComplexMatrix cumsum(int dim=-1) const
Definition fCMatrix.cc:2830
FloatComplexMatrix::insert
FloatComplexMatrix & insert(const FloatMatrix &a, octave_idx_type r, octave_idx_type c)
Definition fCMatrix.cc:196
FloatComplexMatrix::operator==
bool operator==(const FloatComplexMatrix &a) const
Definition fCMatrix.cc:160
FloatComplexMatrix::diag
FloatComplexMatrix diag(octave_idx_type k=0) const
Definition fCMatrix.cc:2860
FloatComplexMatrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition fCMatrix.h:199
FloatComplexMatrix::ifourier2d
FloatComplexMatrix ifourier2d() const
Definition fCMatrix.cc:1129
FloatComplexMatrix::transpose
FloatComplexMatrix transpose() const
Definition fCMatrix.h:178
FloatComplexMatrix::fourier
FloatComplexMatrix fourier() const
Definition fCMatrix.cc:1057
FloatComplexMatrix::sum
FloatComplexMatrix sum(int dim=-1) const
Definition fCMatrix.cc:2842
FloatComplexMatrix::FloatComplexMatrix
FloatComplexMatrix()=default
FloatComplexMatrix::conj
friend FloatComplexMatrix conj(const FloatComplexMatrix &a)
Definition fCMatrix.cc:679
FloatComplexMatrix::column_max
FloatComplexRowVector column_max() const
Definition fCMatrix.cc:3145
FloatComplexMatrix::solve
FloatComplexMatrix solve(MatrixType &mattype, const FloatMatrix &b) const
Definition fCMatrix.cc:1943
FloatComplexMatrix::operator-=
FloatComplexMatrix & operator-=(const FloatDiagMatrix &a)
Definition fCMatrix.cc:2712
FloatComplexMatrix::any
boolMatrix any(int dim=-1) const
Definition fCMatrix.cc:2818
FloatComplexMatrix::column_min
FloatComplexRowVector column_min() const
Definition fCMatrix.cc:3070
FloatComplexMatrix::fill
FloatComplexMatrix & fill(float val)
Definition fCMatrix.cc:350
FloatComplexMatrix::determinant
FloatComplexDET determinant() const
Definition fCMatrix.cc:1183
FloatComplexMatrix::column_is_real_only
bool column_is_real_only(octave_idx_type) const
Definition fCMatrix.cc:2901
FloatComplexMatrix::inverse
FloatComplexMatrix inverse() const
Definition fCMatrix.cc:749
FloatComplexMatrix::ishermitian
bool ishermitian() const
Definition fCMatrix.cc:175
FloatComplexMatrix::append
FloatComplexMatrix append(const FloatMatrix &a) const
Definition fCMatrix.cc:439
FloatComplexMatrix::cumprod
FloatComplexMatrix cumprod(int dim=-1) const
Definition fCMatrix.cc:2824
FloatComplexMatrix::row_max
FloatComplexColumnVector row_max() const
Definition fCMatrix.cc:2995
FloatComplexMatrix::prod
FloatComplexMatrix prod(int dim=-1) const
Definition fCMatrix.cc:2836
FloatComplexMatrix::extract_n
FloatComplexMatrix extract_n(octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const
Definition fCMatrix.cc:697
FloatComplexNDArray::sum
FloatComplexNDArray sum(int dim=-1) const
Definition fCNDArray.cc:384
FloatComplexNDArray::cumprod
FloatComplexNDArray cumprod(int dim=-1) const
Definition fCNDArray.cc:358
FloatComplexNDArray::abs
FloatNDArray abs() const
Definition fCNDArray.cc:490
FloatComplexNDArray::sumsq
FloatComplexNDArray sumsq(int dim=-1) const
Definition fCNDArray.cc:396
FloatComplexNDArray::prod
FloatComplexNDArray prod(int dim=-1) const
Definition fCNDArray.cc:372
FloatComplexNDArray::cumsum
FloatComplexNDArray cumsum(int dim=-1) const
Definition fCNDArray.cc:365
FloatComplexNDArray::all
boolNDArray all(int dim=-1) const
Definition fCNDArray.cc:346
FloatComplexNDArray::any
boolNDArray any(int dim=-1) const
Definition fCNDArray.cc:352
FloatComplexNDArray::diag
FloatComplexNDArray diag(octave_idx_type k=0) const
Definition fCNDArray.cc:600
FloatComplexRowVector::resize
void resize(octave_idx_type n, const FloatComplex &rfv=FloatComplex(0))
Definition fCRowVector.h:133
FloatDiagMatrix::inverse
FloatDiagMatrix inverse() const
Definition fDiagMatrix.cc:227
FloatDiagMatrix::extract_diag
FloatColumnVector extract_diag(octave_idx_type k=0) const
Definition fDiagMatrix.h:108
FloatMatrix::sum
FloatMatrix sum(int dim=-1) const
Definition fMatrix.cc:2400
FloatMatrix::row
FloatRowVector row(octave_idx_type i) const
Definition fMatrix.cc:423
MDiagArray2
Template for two dimensional diagonal array with math operators.
Definition MDiagArray2.h:54
MatrixType::ishermitian
bool ishermitian() const
Definition MatrixType.h:121
MatrixType::mark_as_unsymmetric
void mark_as_unsymmetric()
Definition MatrixType.cc:920
MatrixType::mark_as_full
void mark_as_full()
Definition MatrixType.h:152
MatrixType::type
int type(bool quiet=true)
Definition MatrixType.cc:656
MatrixType::mark_as_rectangular
void mark_as_rectangular()
Definition MatrixType.h:154
base_det
Definition DET.h:38
base_det::square
base_det square() const
Definition DET.h:67
chol
Definition chol.h:37
chol::rcond
COND_T rcond() const
Definition chol.h:62
chol::inverse
T inverse() const
Definition chol.cc:250
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
octave_idx_type
real
ColumnVector real(const ComplexColumnVector &a)
Definition dColVector.cc:137
imag
ColumnVector imag(const ComplexColumnVector &a)
Definition dColVector.cc:143
F77_CONST_CMPLX_ARG
#define F77_CONST_CMPLX_ARG(x)
Definition f77-fcn.h:313
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
Sylvester
FloatComplexMatrix Sylvester(const FloatComplexMatrix &a, const FloatComplexMatrix &b, const FloatComplexMatrix &c)
Definition fCMatrix.cc:3279
xgemm
FloatComplexMatrix xgemm(const FloatComplexMatrix &a, const FloatComplexMatrix &b, blas_trans_type transa, blas_trans_type transb)
Definition fCMatrix.cc:3379
operator*
FloatComplexMatrix operator*(const FloatColumnVector &v, const FloatComplexRowVector &a)
Definition fCMatrix.cc:2653
max
FloatComplexMatrix max(const FloatComplex &c, const FloatComplexMatrix &m)
Definition fCMatrix.cc:3586
operator<<
std::ostream & operator<<(std::ostream &os, const FloatComplexMatrix &a)
Definition fCMatrix.cc:3222
linspace
FloatComplexMatrix linspace(const FloatComplexColumnVector &x1, const FloatComplexColumnVector &x2, octave_idx_type n)
Definition fCMatrix.cc:3661
EMPTY_RETURN_CHECK
#define EMPTY_RETURN_CHECK(T)
Definition fCMatrix.cc:3509
min
FloatComplexMatrix min(const FloatComplex &c, const FloatComplexMatrix &m)
Definition fCMatrix.cc:3514
conj
FloatComplexMatrix conj(const FloatComplexMatrix &a)
Definition fCMatrix.cc:679
Givens
FloatComplexMatrix Givens(const FloatComplex &x, const FloatComplex &y)
Definition fCMatrix.cc:3260
operator>>
std::istream & operator>>(std::istream &is, FloatComplexMatrix &a)
Definition fCMatrix.cc:3237
norm
double norm(const ColumnVector &v)
Definition graphics.cc:5778
scale
void scale(Matrix &m, double x, double y, double z)
Definition graphics.cc:5731
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition lo-error.c:41
V
F77_RET_T const F77_INT const F77_INT const F77_INT const F77_DBLE const F77_DBLE F77_INT F77_DBLE * V
Definition lo-lapack-proto.h:904
d
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
Definition lo-slatec-proto.h:39
x
F77_RET_T const F77_DBLE * x
Definition lo-slatec-proto.h:38
blas_trans_type
blas_trans_type
Definition mx-defs.h:80
blas_no_trans
@ blas_no_trans
Definition mx-defs.h:81
blas_conj_trans
@ blas_conj_trans
Definition mx-defs.h:83
blas_trans
@ blas_trans
Definition mx-defs.h:82
get_blas_char
char get_blas_char(blas_trans_type transt)
Definition mx-defs.h:87
mx_inline_sub2
void mx_inline_sub2(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:128
mx_inline_add2
void mx_inline_add2(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:127
mx_inline_equal
bool mx_inline_equal(std::size_t n, const T1 *x, const T2 *y)
Definition mx-inlines.cc:577
MM_BOOL_OPS
#define MM_BOOL_OPS(M1, M2)
Definition mx-op-defs.h:213
MM_CMP_OPS
#define MM_CMP_OPS(M1, M2)
Definition mx-op-defs.h:196
MS_CMP_OPS
#define MS_CMP_OPS(M, S)
Definition mx-op-defs.h:110
SM_CMP_OPS
#define SM_CMP_OPS(S, M)
Definition mx-op-defs.h:153
SM_BOOL_OPS
#define SM_BOOL_OPS(S, M)
Definition mx-op-defs.h:170
MS_BOOL_OPS
#define MS_BOOL_OPS(M, S)
Definition mx-op-defs.h:127
FloatComplex
std::complex< float > FloatComplex
Definition oct-cmplx.h:34
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
xcdotc
subroutine xcdotc(n, zx, incx, zy, incy, retval)
Definition xcdotc.f:2
xcdotu
subroutine xcdotu(n, zx, incx, zy, incy, retval)
Definition xcdotu.f:2
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg
len
F77_RET_T len
Definition xerbla.cc:61
xilaenv
subroutine xilaenv(ispec, name, opts, n1, n2, n3, n4, retval)
Definition xilaenv.f:2