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