GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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