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