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