GNU Octave 10.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-2025 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 "boolMatrix.h"
40#include "byte-swap.h"
41#include "chMatrix.h"
42#include "chol.h"
43#include "fCColVector.h"
44#include "fCMatrix.h"
45#include "fColVector.h"
46#include "fDiagMatrix.h"
47#include "fMatrix.h"
48#include "fMatrix.h"
49#include "fNDArray.h"
50#include "fRowVector.h"
51#include "lo-blas-proto.h"
52#include "lo-error.h"
53#include "lo-ieee.h"
54#include "lo-lapack-proto.h"
55#include "lo-mappers.h"
56#include "lo-utils.h"
57#include "mx-fdm-fm.h"
58#include "mx-fm-fdm.h"
59#include "mx-inlines.cc"
60#include "mx-op-defs.h"
61#include "oct-cmplx.h"
62#include "oct-fftw.h"
63#include "oct-locbuf.h"
64#include "oct-norm.h"
65#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// FIXME: Do these really belong here? Maybe they should be in a base class?
2380
2381FloatMatrix
2382FloatMatrix::cumprod (int dim) const
2383{
2384 return FloatNDArray::cumprod (dim);
2385}
2386
2387FloatMatrix
2388FloatMatrix::cumsum (int dim) const
2389{
2390 return FloatNDArray::cumsum (dim);
2391}
2392
2393FloatMatrix
2394FloatMatrix::prod (int dim) const
2395{
2396 return FloatNDArray::prod (dim);
2397}
2398
2399FloatMatrix
2400FloatMatrix::sum (int dim) const
2401{
2402 return FloatNDArray::sum (dim);
2403}
2404
2405FloatMatrix
2406FloatMatrix::sumsq (int dim) const
2407{
2408 return FloatNDArray::sumsq (dim);
2409}
2410
2411FloatMatrix
2412FloatMatrix::abs () const
2413{
2414 return FloatNDArray::abs ();
2415}
2416
2417FloatMatrix
2418FloatMatrix::diag (octave_idx_type k) const
2419{
2420 return FloatNDArray::diag (k);
2421}
2422
2423FloatDiagMatrix
2424FloatMatrix::diag (octave_idx_type m, octave_idx_type n) const
2425{
2426 FloatDiagMatrix retval;
2427
2428 octave_idx_type nr = rows ();
2429 octave_idx_type nc = cols ();
2430
2431 if (nr == 1 || nc == 1)
2432 retval = FloatDiagMatrix (*this, m, n);
2433 else
2434 (*current_liboctave_error_handler) ("diag: expecting vector argument");
2435
2436 return retval;
2437}
2438
2439FloatColumnVector
2440FloatMatrix::row_min () const
2441{
2442 Array<octave_idx_type> dummy_idx;
2443 return row_min (dummy_idx);
2444}
2445
2446FloatColumnVector
2447FloatMatrix::row_min (Array<octave_idx_type>& idx_arg) const
2448{
2449 FloatColumnVector result;
2450
2451 octave_idx_type nr = rows ();
2452 octave_idx_type nc = cols ();
2453
2454 if (nr > 0 && nc > 0)
2455 {
2456 result.resize (nr);
2457 idx_arg.resize (dim_vector (nr, 1));
2458
2459 for (octave_idx_type i = 0; i < nr; i++)
2460 {
2461 octave_idx_type idx_j;
2462
2463 float tmp_min = octave::numeric_limits<float>::NaN ();
2464
2465 for (idx_j = 0; idx_j < nc; idx_j++)
2466 {
2467 tmp_min = elem (i, idx_j);
2468
2469 if (! octave::math::isnan (tmp_min))
2470 break;
2471 }
2472
2473 for (octave_idx_type j = idx_j+1; j < nc; j++)
2474 {
2475 float tmp = elem (i, j);
2476
2477 if (octave::math::isnan (tmp))
2478 continue;
2479 else if (tmp < tmp_min)
2480 {
2481 idx_j = j;
2482 tmp_min = tmp;
2483 }
2484 }
2485
2486 result.elem (i) = tmp_min;
2487 idx_arg.elem (i) = (octave::math::isnan (tmp_min) ? 0 : idx_j);
2488 }
2489 }
2490
2491 return result;
2492}
2493
2494FloatColumnVector
2495FloatMatrix::row_max () const
2496{
2497 Array<octave_idx_type> dummy_idx;
2498 return row_max (dummy_idx);
2499}
2500
2501FloatColumnVector
2502FloatMatrix::row_max (Array<octave_idx_type>& idx_arg) const
2503{
2504 FloatColumnVector result;
2505
2506 octave_idx_type nr = rows ();
2507 octave_idx_type nc = cols ();
2508
2509 if (nr > 0 && nc > 0)
2510 {
2511 result.resize (nr);
2512 idx_arg.resize (dim_vector (nr, 1));
2513
2514 for (octave_idx_type i = 0; i < nr; i++)
2515 {
2516 octave_idx_type idx_j;
2517
2518 float tmp_max = octave::numeric_limits<float>::NaN ();
2519
2520 for (idx_j = 0; idx_j < nc; idx_j++)
2521 {
2522 tmp_max = elem (i, idx_j);
2523
2524 if (! octave::math::isnan (tmp_max))
2525 break;
2526 }
2527
2528 for (octave_idx_type j = idx_j+1; j < nc; j++)
2529 {
2530 float tmp = elem (i, j);
2531
2532 if (octave::math::isnan (tmp))
2533 continue;
2534 else if (tmp > tmp_max)
2535 {
2536 idx_j = j;
2537 tmp_max = tmp;
2538 }
2539 }
2540
2541 result.elem (i) = tmp_max;
2542 idx_arg.elem (i) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
2543 }
2544 }
2545
2546 return result;
2547}
2548
2549FloatRowVector
2550FloatMatrix::column_min () const
2551{
2552 Array<octave_idx_type> dummy_idx;
2553 return column_min (dummy_idx);
2554}
2555
2556FloatRowVector
2557FloatMatrix::column_min (Array<octave_idx_type>& idx_arg) const
2558{
2559 FloatRowVector result;
2560
2561 octave_idx_type nr = rows ();
2562 octave_idx_type nc = cols ();
2563
2564 if (nr > 0 && nc > 0)
2565 {
2566 result.resize (nc);
2567 idx_arg.resize (dim_vector (1, nc));
2568
2569 for (octave_idx_type j = 0; j < nc; j++)
2570 {
2571 octave_idx_type idx_i;
2572
2573 float tmp_min = octave::numeric_limits<float>::NaN ();
2574
2575 for (idx_i = 0; idx_i < nr; idx_i++)
2576 {
2577 tmp_min = elem (idx_i, j);
2578
2579 if (! octave::math::isnan (tmp_min))
2580 break;
2581 }
2582
2583 for (octave_idx_type i = idx_i+1; i < nr; i++)
2584 {
2585 float tmp = elem (i, j);
2586
2587 if (octave::math::isnan (tmp))
2588 continue;
2589 else if (tmp < tmp_min)
2590 {
2591 idx_i = i;
2592 tmp_min = tmp;
2593 }
2594 }
2595
2596 result.elem (j) = tmp_min;
2597 idx_arg.elem (j) = (octave::math::isnan (tmp_min) ? 0 : idx_i);
2598 }
2599 }
2600
2601 return result;
2602}
2603
2604FloatRowVector
2605FloatMatrix::column_max () const
2606{
2607 Array<octave_idx_type> dummy_idx;
2608 return column_max (dummy_idx);
2609}
2610
2611FloatRowVector
2612FloatMatrix::column_max (Array<octave_idx_type>& idx_arg) const
2613{
2614 FloatRowVector result;
2615
2616 octave_idx_type nr = rows ();
2617 octave_idx_type nc = cols ();
2618
2619 if (nr > 0 && nc > 0)
2620 {
2621 result.resize (nc);
2622 idx_arg.resize (dim_vector (1, nc));
2623
2624 for (octave_idx_type j = 0; j < nc; j++)
2625 {
2626 octave_idx_type idx_i;
2627
2628 float tmp_max = octave::numeric_limits<float>::NaN ();
2629
2630 for (idx_i = 0; idx_i < nr; idx_i++)
2631 {
2632 tmp_max = elem (idx_i, j);
2633
2634 if (! octave::math::isnan (tmp_max))
2635 break;
2636 }
2637
2638 for (octave_idx_type i = idx_i+1; i < nr; i++)
2639 {
2640 float tmp = elem (i, j);
2641
2642 if (octave::math::isnan (tmp))
2643 continue;
2644 else if (tmp > tmp_max)
2645 {
2646 idx_i = i;
2647 tmp_max = tmp;
2648 }
2649 }
2650
2651 result.elem (j) = tmp_max;
2652 idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_i);
2653 }
2654 }
2655
2656 return result;
2657}
2658
2659std::ostream&
2660operator << (std::ostream& os, const FloatMatrix& a)
2661{
2662 for (octave_idx_type i = 0; i < a.rows (); i++)
2663 {
2664 for (octave_idx_type j = 0; j < a.cols (); j++)
2665 {
2666 os << ' ';
2667 octave::write_value<float> (os, a.elem (i, j));
2668 }
2669 os << "\n";
2670 }
2671 return os;
2672}
2673
2674std::istream&
2675operator >> (std::istream& is, FloatMatrix& a)
2676{
2677 octave_idx_type nr = a.rows ();
2678 octave_idx_type nc = a.cols ();
2679
2680 if (nr > 0 && nc > 0)
2681 {
2682 float tmp;
2683 for (octave_idx_type i = 0; i < nr; i++)
2684 for (octave_idx_type j = 0; j < nc; j++)
2685 {
2686 tmp = octave::read_value<float> (is);
2687 if (is)
2688 a.elem (i, j) = tmp;
2689 else
2690 return is;
2691 }
2692 }
2693
2694 return is;
2695}
2696
2697FloatMatrix
2698Givens (float x, float y)
2699{
2700 float cc, s, temp_r;
2701
2702 F77_FUNC (slartg, SLARTG) (x, y, cc, s, temp_r);
2703
2704 FloatMatrix g (2, 2);
2705
2706 g.elem (0, 0) = cc;
2707 g.elem (1, 1) = cc;
2708 g.elem (0, 1) = s;
2709 g.elem (1, 0) = -s;
2710
2711 return g;
2712}
2713
2714FloatMatrix
2715Sylvester (const FloatMatrix& a, const FloatMatrix& b, const FloatMatrix& c)
2716{
2717 FloatMatrix retval;
2718
2719 // FIXME: need to check that a, b, and c are all the same size.
2720
2721 // Compute Schur decompositions.
2722
2723 octave::math::schur<FloatMatrix> as (a, "U");
2724 octave::math::schur<FloatMatrix> bs (b, "U");
2725
2726 // Transform c to new coordinates.
2727
2728 FloatMatrix ua = as.unitary_schur_matrix ();
2729 FloatMatrix sch_a = as.schur_matrix ();
2730
2731 FloatMatrix ub = bs.unitary_schur_matrix ();
2732 FloatMatrix sch_b = bs.schur_matrix ();
2733
2734 FloatMatrix cx = ua.transpose () * c * ub;
2735
2736 // Solve the sylvester equation, back-transform, and return the
2737 // solution.
2738
2739 F77_INT a_nr = octave::to_f77_int (a.rows ());
2740 F77_INT b_nr = octave::to_f77_int (b.rows ());
2741
2742 float scale;
2743 F77_INT info;
2744
2745 float *pa = sch_a.rwdata ();
2746 float *pb = sch_b.rwdata ();
2747 float *px = cx.rwdata ();
2748
2749 F77_XFCN (strsyl, STRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
2750 F77_CONST_CHAR_ARG2 ("N", 1),
2751 1, a_nr, b_nr, pa, a_nr, pb,
2752 b_nr, px, a_nr, scale, info
2753 F77_CHAR_ARG_LEN (1)
2754 F77_CHAR_ARG_LEN (1)));
2755
2756 // FIXME: check info?
2757
2758 retval = ua*cx*ub.transpose ();
2759
2760 return retval;
2761}
2762
2763// matrix by matrix -> matrix operations
2764
2765/*
2766## Simple Dot Product, Matrix-Vector, and Matrix-Matrix Unit tests
2767%!assert (single ([1 2 3]) * single ([ 4 ; 5 ; 6]), single (32), 5e-7)
2768%!assert (single ([1 2 ; 3 4]) * single ([5 ; 6]), single ([17 ; 39]), 5e-7)
2769%!assert (single ([1 2 ; 3 4]) * single ([5 6 ; 7 8]), single ([19 22; 43 50]), 5e-7)
2770
2771## Test some simple identities
2772%!shared M, cv, rv
2773%! M = single (randn (10,10));
2774%! cv = single (randn (10,1));
2775%! rv = single (randn (1,10));
2776%!assert ([M*cv,M*cv], M*[cv,cv], 5e-6)
2777%!assert ([M'*cv,M'*cv], M'*[cv,cv], 5e-6)
2778%!assert ([rv*M;rv*M], [rv;rv]*M, 5e-6)
2779%!assert ([rv*M';rv*M'], [rv;rv]*M', 5e-6)
2780%!assert (2*rv*cv, [rv,rv]*[cv;cv], 5e-6)
2781
2782*/
2783
2784static char
2785get_blas_trans_arg (bool trans)
2786{
2787 return trans ? 'T' : 'N';
2788}
2789
2790// the general GEMM operation
2791
2792FloatMatrix
2793xgemm (const FloatMatrix& a, const FloatMatrix& b,
2794 blas_trans_type transa, blas_trans_type transb)
2795{
2796 FloatMatrix retval;
2797
2798 bool tra = transa != blas_no_trans;
2799 bool trb = transb != blas_no_trans;
2800
2801 F77_INT a_nr = octave::to_f77_int (tra ? a.cols () : a.rows ());
2802 F77_INT a_nc = octave::to_f77_int (tra ? a.rows () : a.cols ());
2803
2804 F77_INT b_nr = octave::to_f77_int (trb ? b.cols () : b.rows ());
2805 F77_INT b_nc = octave::to_f77_int (trb ? b.rows () : b.cols ());
2806
2807 if (a_nc != b_nr)
2808 octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
2809
2810 if (a_nr == 0 || a_nc == 0 || b_nc == 0)
2811 retval = FloatMatrix (a_nr, b_nc, 0.0);
2812 else if (a.data () == b.data () && a_nr == b_nc && tra != trb)
2813 {
2814 F77_INT lda = octave::to_f77_int (a.rows ());
2815
2816 retval = FloatMatrix (a_nr, b_nc);
2817 float *c = retval.rwdata ();
2818
2819 const char ctra = get_blas_trans_arg (tra);
2820 F77_XFCN (ssyrk, SSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),
2821 F77_CONST_CHAR_ARG2 (&ctra, 1),
2822 a_nr, a_nc, 1.0,
2823 a.data (), lda, 0.0, c, a_nr
2824 F77_CHAR_ARG_LEN (1)
2825 F77_CHAR_ARG_LEN (1)));
2826 for (int j = 0; j < a_nr; j++)
2827 for (int i = 0; i < j; i++)
2828 retval.xelem (j, i) = retval.xelem (i, j);
2829
2830 }
2831 else
2832 {
2833 F77_INT lda = octave::to_f77_int (a.rows ());
2834 F77_INT tda = octave::to_f77_int (a.cols ());
2835 F77_INT ldb = octave::to_f77_int (b.rows ());
2836 F77_INT tdb = octave::to_f77_int (b.cols ());
2837
2838 retval = FloatMatrix (a_nr, b_nc);
2839 float *c = retval.rwdata ();
2840
2841 if (b_nc == 1)
2842 {
2843 if (a_nr == 1)
2844 F77_FUNC (xsdot, XSDOT) (a_nc, a.data (), 1, b.data (), 1, *c);
2845 else
2846 {
2847 const char ctra = get_blas_trans_arg (tra);
2848 F77_XFCN (sgemv, SGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),
2849 lda, tda, 1.0, a.data (), lda,
2850 b.data (), 1, 0.0, c, 1
2851 F77_CHAR_ARG_LEN (1)));
2852 }
2853 }
2854 else if (a_nr == 1)
2855 {
2856 const char crevtrb = get_blas_trans_arg (! trb);
2857 F77_XFCN (sgemv, SGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),
2858 ldb, tdb, 1.0, b.data (), ldb,
2859 a.data (), 1, 0.0, c, 1
2860 F77_CHAR_ARG_LEN (1)));
2861 }
2862 else
2863 {
2864 const char ctra = get_blas_trans_arg (tra);
2865 const char ctrb = get_blas_trans_arg (trb);
2866 F77_XFCN (sgemm, SGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),
2867 F77_CONST_CHAR_ARG2 (&ctrb, 1),
2868 a_nr, b_nc, a_nc, 1.0, a.data (),
2869 lda, b.data (), ldb, 0.0, c, a_nr
2870 F77_CHAR_ARG_LEN (1)
2871 F77_CHAR_ARG_LEN (1)));
2872 }
2873 }
2874
2875 return retval;
2876}
2877
2878FloatMatrix
2879operator * (const FloatMatrix& a, const FloatMatrix& b)
2880{
2881 return xgemm (a, b);
2882}
2883
2884// FIXME: it would be nice to share code among the min/max functions below.
2885
2886#define EMPTY_RETURN_CHECK(T) \
2887 if (nr == 0 || nc == 0) \
2888 return T (nr, nc);
2889
2890FloatMatrix
2891min (float d, const FloatMatrix& m)
2892{
2893 octave_idx_type nr = m.rows ();
2894 octave_idx_type nc = m.columns ();
2895
2896 EMPTY_RETURN_CHECK (FloatMatrix);
2897
2898 FloatMatrix result (nr, nc);
2899
2900 for (octave_idx_type j = 0; j < nc; j++)
2901 for (octave_idx_type i = 0; i < nr; i++)
2902 {
2903 octave_quit ();
2904 result(i, j) = octave::math::min (d, m(i, j));
2905 }
2906
2907 return result;
2908}
2909
2910FloatMatrix
2911min (const FloatMatrix& m, float d)
2912{
2913 octave_idx_type nr = m.rows ();
2914 octave_idx_type nc = m.columns ();
2915
2916 EMPTY_RETURN_CHECK (FloatMatrix);
2917
2918 FloatMatrix result (nr, nc);
2919
2920 for (octave_idx_type j = 0; j < nc; j++)
2921 for (octave_idx_type i = 0; i < nr; i++)
2922 {
2923 octave_quit ();
2924 result(i, j) = octave::math::min (m(i, j), d);
2925 }
2926
2927 return result;
2928}
2929
2930FloatMatrix
2931min (const FloatMatrix& a, const FloatMatrix& b)
2932{
2933 octave_idx_type nr = a.rows ();
2934 octave_idx_type nc = a.columns ();
2935
2936 if (nr != b.rows () || nc != b.columns ())
2937 (*current_liboctave_error_handler)
2938 ("two-arg min requires same size arguments");
2939
2940 EMPTY_RETURN_CHECK (FloatMatrix);
2941
2942 FloatMatrix result (nr, nc);
2943
2944 for (octave_idx_type j = 0; j < nc; j++)
2945 for (octave_idx_type i = 0; i < nr; i++)
2946 {
2947 octave_quit ();
2948 result(i, j) = octave::math::min (a(i, j), b(i, j));
2949 }
2950
2951 return result;
2952}
2953
2954FloatMatrix
2955max (float d, const FloatMatrix& m)
2956{
2957 octave_idx_type nr = m.rows ();
2958 octave_idx_type nc = m.columns ();
2959
2960 EMPTY_RETURN_CHECK (FloatMatrix);
2961
2962 FloatMatrix result (nr, nc);
2963
2964 for (octave_idx_type j = 0; j < nc; j++)
2965 for (octave_idx_type i = 0; i < nr; i++)
2966 {
2967 octave_quit ();
2968 result(i, j) = octave::math::max (d, m(i, j));
2969 }
2970
2971 return result;
2972}
2973
2974FloatMatrix
2975max (const FloatMatrix& m, float d)
2976{
2977 octave_idx_type nr = m.rows ();
2978 octave_idx_type nc = m.columns ();
2979
2980 EMPTY_RETURN_CHECK (FloatMatrix);
2981
2982 FloatMatrix result (nr, nc);
2983
2984 for (octave_idx_type j = 0; j < nc; j++)
2985 for (octave_idx_type i = 0; i < nr; i++)
2986 {
2987 octave_quit ();
2988 result(i, j) = octave::math::max (m(i, j), d);
2989 }
2990
2991 return result;
2992}
2993
2994FloatMatrix
2995max (const FloatMatrix& a, const FloatMatrix& b)
2996{
2997 octave_idx_type nr = a.rows ();
2998 octave_idx_type nc = a.columns ();
2999
3000 if (nr != b.rows () || nc != b.columns ())
3001 (*current_liboctave_error_handler)
3002 ("two-arg max requires same size arguments");
3003
3004 EMPTY_RETURN_CHECK (FloatMatrix);
3005
3006 FloatMatrix result (nr, nc);
3007
3008 for (octave_idx_type j = 0; j < nc; j++)
3009 for (octave_idx_type i = 0; i < nr; i++)
3010 {
3011 octave_quit ();
3012 result(i, j) = octave::math::max (a(i, j), b(i, j));
3013 }
3014
3015 return result;
3016}
3017
3018FloatMatrix
3019linspace (const FloatColumnVector& x1,
3020 const FloatColumnVector& x2,
3021 octave_idx_type n)
3022
3023{
3024 octave_idx_type m = x1.numel ();
3025
3026 if (x2.numel () != m)
3027 (*current_liboctave_error_handler)
3028 ("linspace: vectors must be of equal length");
3029
3030 FloatMatrix retval;
3031
3032 if (n < 1)
3033 {
3034 retval.clear (m, 0);
3035 return retval;
3036 }
3037
3038 retval.clear (m, n);
3039 for (octave_idx_type i = 0; i < m; i++)
3040 retval.insert (linspace (x1(i), x2(i), n), i, 0);
3041
3042 return retval;
3043}
3044
3045MS_CMP_OPS (FloatMatrix, float)
3046MS_BOOL_OPS (FloatMatrix, float)
3047
3048SM_CMP_OPS (float, FloatMatrix)
3049SM_BOOL_OPS (float, FloatMatrix)
3050
3051MM_CMP_OPS (FloatMatrix, FloatMatrix)
3052MM_BOOL_OPS (FloatMatrix, FloatMatrix)
FloatDET
base_det< float > FloatDET
Definition DET.h:86
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:525
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.h:649
Array::clear
void clear()
Definition Array-base.cc:110
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:563
Array::rows
octave_idx_type rows() const
Definition Array.h:463
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.h:475
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::make_unique
void make_unique()
Definition Array.h:220
Array::data
const T * data() const
Size of the specified dimension.
Definition Array.h:665
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
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::sumsq
FloatMatrix sumsq(int dim=-1) const
Definition fMatrix.cc:2406
FloatMatrix::sum
FloatMatrix sum(int dim=-1) const
Definition fMatrix.cc:2400
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:2495
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:2412
FloatMatrix::fill
FloatMatrix & fill(float val)
Definition fMatrix.cc:228
FloatMatrix::inverse
FloatMatrix inverse() const
Definition fMatrix.cc:466
FloatMatrix::prod
FloatMatrix prod(int dim=-1) const
Definition fMatrix.cc:2394
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::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::cumsum
FloatMatrix cumsum(int dim=-1) const
Definition fMatrix.cc:2388
FloatMatrix::fourier
FloatComplexMatrix fourier() const
Definition fMatrix.cc:749
FloatMatrix::column_min
FloatRowVector column_min() const
Definition fMatrix.cc:2550
FloatMatrix::append
FloatMatrix append(const FloatMatrix &a) const
Definition fMatrix.cc:272
FloatMatrix::column_max
FloatRowVector column_max() const
Definition fMatrix.cc:2605
FloatMatrix::diag
FloatMatrix diag(octave_idx_type k=0) const
Definition fMatrix.cc:2418
FloatMatrix::operator!=
bool operator!=(const FloatMatrix &a) const
Definition fMatrix.cc:133
FloatMatrix::row_min
FloatColumnVector row_min() const
Definition fMatrix.cc:2440
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::cumprod
FloatMatrix cumprod(int dim=-1) const
Definition fMatrix.cc:2382
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
FloatNDArray::sum
FloatNDArray sum(int dim=-1) const
Definition fNDArray.cc:388
FloatNDArray::prod
FloatNDArray prod(int dim=-1) const
Definition fNDArray.cc:376
FloatNDArray::abs
FloatNDArray abs() const
Definition fNDArray.cc:538
FloatNDArray::cumprod
FloatNDArray cumprod(int dim=-1) const
Definition fNDArray.cc:364
FloatNDArray::insert
FloatNDArray & insert(const FloatNDArray &a, octave_idx_type r, octave_idx_type c)
Definition fNDArray.cc:522
FloatNDArray::cumsum
FloatNDArray cumsum(int dim=-1) const
Definition fNDArray.cc:370
FloatNDArray::sumsq
FloatNDArray sumsq(int dim=-1) const
Definition fNDArray.cc:400
FloatNDArray::diag
FloatNDArray diag(octave_idx_type k=0) const
Definition fNDArray.cc:577
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:37
chol::rcond
COND_T rcond() const
Definition chol.h:62
chol::inverse
T inverse() const
Definition chol.cc:250
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
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
linspace
FloatMatrix linspace(const FloatColumnVector &x1, const FloatColumnVector &x2, octave_idx_type n)
Definition fMatrix.cc:3019
min
FloatMatrix min(float d, const FloatMatrix &m)
Definition fMatrix.cc:2891
EMPTY_RETURN_CHECK
#define EMPTY_RETURN_CHECK(T)
Definition fMatrix.cc:2886
real
FloatMatrix real(const FloatComplexMatrix &a)
Definition fMatrix.cc:392
imag
FloatMatrix imag(const FloatComplexMatrix &a)
Definition fMatrix.cc:398
max
FloatMatrix max(float d, const FloatMatrix &m)
Definition fMatrix.cc:2955
operator<<
std::ostream & operator<<(std::ostream &os, const FloatMatrix &a)
Definition fMatrix.cc:2660
Givens
FloatMatrix Givens(float x, float y)
Definition fMatrix.cc:2698
Sylvester
FloatMatrix Sylvester(const FloatMatrix &a, const FloatMatrix &b, const FloatMatrix &c)
Definition fMatrix.cc:2715
xgemm
FloatMatrix xgemm(const FloatMatrix &a, const FloatMatrix &b, blas_trans_type transa, blas_trans_type transb)
Definition fMatrix.cc:2793
operator>>
std::istream & operator>>(std::istream &is, FloatMatrix &a)
Definition fMatrix.cc:2675
norm
double norm(const ColumnVector &v)
Definition graphics.cc:5778
scale
void scale(Matrix &m, double x, double y, double z)
Definition graphics.cc:5731
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
mx_inline_equal
bool mx_inline_equal(std::size_t n, const T1 *x, const T2 *y)
Definition mx-inlines.cc:577
mx_inline_real
void mx_inline_real(std::size_t n, T *r, const std::complex< T > *x)
Definition mx-inlines.cc:325
mx_inline_imag
void mx_inline_imag(std::size_t n, T *r, const std::complex< T > *x)
Definition mx-inlines.cc:333
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
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