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