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