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