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