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