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