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