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