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