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