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