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