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