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