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