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