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