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