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