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