GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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