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