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