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&,
74  F77_CONST_CHAR_ARG_DECL,
75  F77_CONST_CHAR_ARG_DECL,
76  const octave_idx_type&, const octave_idx_type&,
77  const octave_idx_type&, const octave_idx_type&,
78  octave_idx_type&
79  F77_CHAR_ARG_LEN_DECL
80  F77_CHAR_ARG_LEN_DECL);
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&
87  F77_CHAR_ARG_LEN_DECL);
88 
89  F77_RET_T
90  F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL,
91  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&
96  F77_CHAR_ARG_LEN_DECL
97  F77_CHAR_ARG_LEN_DECL);
98 
99 
100  F77_RET_T
101  F77_FUNC (sgemm, SGEMM) (F77_CONST_CHAR_ARG_DECL,
102  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&
108  F77_CHAR_ARG_LEN_DECL
109  F77_CHAR_ARG_LEN_DECL);
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&
118  F77_CHAR_ARG_LEN_DECL);
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,
127  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&
131  F77_CHAR_ARG_LEN_DECL
132  F77_CHAR_ARG_LEN_DECL);
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&
146  F77_CHAR_ARG_LEN_DECL);
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&
158  F77_CHAR_ARG_LEN_DECL);
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&
181  F77_CHAR_ARG_LEN_DECL);
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&
189  F77_CHAR_ARG_LEN_DECL);
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&
195  F77_CHAR_ARG_LEN_DECL);
196 
197  F77_RET_T
198  F77_FUNC (strtri, STRTRI) (F77_CONST_CHAR_ARG_DECL,
199  F77_CONST_CHAR_ARG_DECL,
200  const octave_idx_type&, const float*,
201  const octave_idx_type&, octave_idx_type&
202  F77_CHAR_ARG_LEN_DECL
203  F77_CHAR_ARG_LEN_DECL);
204  F77_RET_T
205  F77_FUNC (strcon, STRCON) (F77_CONST_CHAR_ARG_DECL,
206  F77_CONST_CHAR_ARG_DECL,
207  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&
211  F77_CHAR_ARG_LEN_DECL
212  F77_CHAR_ARG_LEN_DECL
213  F77_CHAR_ARG_LEN_DECL);
214  F77_RET_T
215  F77_FUNC (strtrs, STRTRS) (F77_CONST_CHAR_ARG_DECL,
216  F77_CONST_CHAR_ARG_DECL,
217  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&
222  F77_CHAR_ARG_LEN_DECL
223  F77_CHAR_ARG_LEN_DECL
224  F77_CHAR_ARG_LEN_DECL);
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,
232  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&
238  F77_CHAR_ARG_LEN_DECL
239  F77_CHAR_ARG_LEN_DECL);
240 
241  F77_RET_T
242  F77_FUNC (xslange, XSLANGE) (F77_CONST_CHAR_ARG_DECL,
243  const octave_idx_type&,
244  const octave_idx_type&, const float*,
245  const octave_idx_type&, float*, float&
246  F77_CHAR_ARG_LEN_DECL);
247 }
248 
249 // Matrix class.
250 
251 FloatMatrix::FloatMatrix (const FloatRowVector& rv)
252  : FloatNDArray (rv)
253 {
254 }
255 
256 FloatMatrix::FloatMatrix (const FloatColumnVector& cv)
257  : FloatNDArray (cv)
258 {
259 }
260 
261 FloatMatrix::FloatMatrix (const FloatDiagMatrix& a)
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 
268 FloatMatrix::FloatMatrix (const MDiagArray2<float>& a)
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 
275 FloatMatrix::FloatMatrix (const DiagArray2<float>& a)
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 
282 FloatMatrix::FloatMatrix (const PermMatrix& a)
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 
293 FloatMatrix::FloatMatrix (const boolMatrix& a)
294  : FloatNDArray (a)
295 {
296 }
297 
298 FloatMatrix::FloatMatrix (const charMatrix& a)
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
307 FloatMatrix::operator == (const FloatMatrix& a) const
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
316 FloatMatrix::operator != (const FloatMatrix& a) const
317 {
318  return !(*this == a);
319 }
320 
321 bool
322 FloatMatrix::is_symmetric (void) const
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 
337 FloatMatrix&
338 FloatMatrix::insert (const FloatMatrix& a,
339  octave_idx_type r, octave_idx_type c)
340 {
341  FloatNDArray::insert (a, r, c);
342  return *this;
343 }
344 
345 FloatMatrix&
346 FloatMatrix::insert (const FloatRowVector& a,
347  octave_idx_type r, octave_idx_type c)
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 
368 FloatMatrix&
369 FloatMatrix::insert (const FloatColumnVector& a,
370  octave_idx_type r, octave_idx_type c)
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 
391 FloatMatrix&
392 FloatMatrix::insert (const FloatDiagMatrix& a,
393  octave_idx_type r, octave_idx_type c)
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 
419 FloatMatrix&
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 
437 FloatMatrix&
438 FloatMatrix::fill (float val, octave_idx_type r1, octave_idx_type c1,
439  octave_idx_type r2, octave_idx_type c2)
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 
466 FloatMatrix
467 FloatMatrix::append (const FloatMatrix& a) const
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 
484 FloatMatrix
485 FloatMatrix::append (const FloatRowVector& a) const
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 
502 FloatMatrix
503 FloatMatrix::append (const FloatColumnVector& a) const
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 
520 FloatMatrix
521 FloatMatrix::append (const FloatDiagMatrix& a) const
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 
538 FloatMatrix
539 FloatMatrix::stack (const FloatMatrix& a) const
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 
557 FloatMatrix
558 FloatMatrix::stack (const FloatRowVector& a) const
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 
576 FloatMatrix
577 FloatMatrix::stack (const FloatColumnVector& a) const
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 
595 FloatMatrix
596 FloatMatrix::stack (const FloatDiagMatrix& a) const
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 
614 FloatMatrix
615 real (const FloatComplexMatrix& a)
616 {
617  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_real);
618 }
619 
620 FloatMatrix
621 imag (const FloatComplexMatrix& a)
622 {
623  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_imag);
624 }
625 
626 FloatMatrix
627 FloatMatrix::extract (octave_idx_type r1, octave_idx_type c1,
628  octave_idx_type r2, octave_idx_type c2) const
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 
636 FloatMatrix
637 FloatMatrix::extract_n (octave_idx_type r1, octave_idx_type c1,
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 
645 FloatRowVector
646 FloatMatrix::row (octave_idx_type i) const
647 {
648  return index (idx_vector (i), idx_vector::colon);
649 }
650 
651 FloatColumnVector
652 FloatMatrix::column (octave_idx_type i) const
653 {
654  return index (idx_vector::colon, idx_vector (i));
655 }
656 
657 FloatMatrix
658 FloatMatrix::inverse (void) const
659 {
660  octave_idx_type info;
661  float rcon;
662  MatrixType mattype (*this);
663  return inverse (mattype, info, rcon, 0, 0);
664 }
665 
666 FloatMatrix
667 FloatMatrix::inverse (octave_idx_type& info) const
668 {
669  float rcon;
670  MatrixType mattype (*this);
671  return inverse (mattype, info, rcon, 0, 0);
672 }
673 
674 FloatMatrix
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 
682 FloatMatrix
683 FloatMatrix::inverse (MatrixType& mattype) const
684 {
685  octave_idx_type info;
686  float rcon;
687  return inverse (mattype, info, rcon, 0, 0);
688 }
689 
690 FloatMatrix
691 FloatMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const
692 {
693  float rcon;
694  return inverse (mattype, info, rcon, 0, 0);
695 }
696 
697 FloatMatrix
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);
732  OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 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 
755 FloatMatrix
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 
837 FloatMatrix
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 
876 FloatMatrix
877 FloatMatrix::pseudo_inverse (float tol) const
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 
885  FloatColumnVector sigma = S.extract_diag ();
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 
915 FloatComplexMatrix
916 FloatMatrix::fourier (void) const
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 
944 FloatComplexMatrix
945 FloatMatrix::ifourier (void) const
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 
974 FloatComplexMatrix
975 FloatMatrix::fourier2d (void) const
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 
986 FloatComplexMatrix
987 FloatMatrix::ifourier2d (void) const
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 
1020 FloatComplexMatrix
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 
1061 FloatComplexMatrix
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 
1105 FloatComplexMatrix
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 
1171 FloatComplexMatrix
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
1243 FloatMatrix::determinant (void) const
1244 {
1245  octave_idx_type info;
1246  float rcon;
1247  return determinant (info, rcon, 0);
1248 }
1249 
1250 FloatDET
1251 FloatMatrix::determinant (octave_idx_type& info) const
1252 {
1253  float rcon;
1254  return determinant (info, rcon, 0);
1255 }
1256 
1257 FloatDET
1258 FloatMatrix::determinant (octave_idx_type& info, float& rcon,
1259  int calc_cond) const
1260 {
1261  MatrixType mattype (*this);
1262  return determinant (mattype, info, rcon, calc_cond);
1263 }
1264 
1265 FloatDET
1266 FloatMatrix::determinant (MatrixType& mattype,
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
1403 FloatMatrix::rcond (void) const
1404 {
1405  MatrixType mattype (*this);
1406  return rcond (mattype);
1407 }
1408 
1409 float
1410 FloatMatrix::rcond (MatrixType &mattype) const
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 
1572 FloatMatrix
1573 FloatMatrix::utsolve (MatrixType &mattype, const FloatMatrix& b,
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 ())
1584  (*current_liboctave_error_handler)
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 
1668 FloatMatrix
1669 FloatMatrix::ltsolve (MatrixType &mattype, const FloatMatrix& b,
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 ())
1680  (*current_liboctave_error_handler)
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 
1763 FloatMatrix
1764 FloatMatrix::fsolve (MatrixType &mattype, const FloatMatrix& b,
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 ())
1775  (*current_liboctave_error_handler)
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
1888  gripe_singular_matrix ();
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 
1944 FloatMatrix
1945 FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b) const
1946 {
1947  octave_idx_type info;
1948  float rcon;
1949  return solve (typ, b, info, rcon, 0);
1950 }
1951 
1952 FloatMatrix
1953 FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b,
1954  octave_idx_type& info) const
1955 {
1956  float rcon;
1957  return solve (typ, b, info, rcon, 0);
1958 }
1959 
1960 FloatMatrix
1961 FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b,
1962  octave_idx_type& info, float& rcon) const
1963 {
1964  return solve (typ, b, info, rcon, 0);
1965 }
1966 
1967 FloatMatrix
1968 FloatMatrix::solve (MatrixType &mattype, const FloatMatrix& b,
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 
2005 FloatComplexMatrix
2006 FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b) const
2007 {
2008  octave_idx_type info;
2009  float rcon;
2010  return solve (typ, b, info, rcon, 0);
2011 }
2012 
2013 FloatComplexMatrix
2014 FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b,
2015  octave_idx_type& info) const
2016 {
2017  float rcon;
2018  return solve (typ, b, info, rcon, 0);
2019 }
2020 
2021 FloatComplexMatrix
2022 FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b,
2023  octave_idx_type& info,
2024  float& rcon) const
2025 {
2026  return solve (typ, b, info, rcon, 0);
2027 }
2028 
2029 static FloatMatrix
2030 stack_complex_matrix (const FloatComplexMatrix& cm)
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
2047 unstack_complex_matrix (const FloatMatrix& sm)
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 
2060 FloatComplexMatrix
2061 FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b,
2062  octave_idx_type& info,
2063  float& rcon, solve_singularity_handler sing_handler,
2064  bool singular_fallback, blas_trans_type transt) const
2065 {
2066  FloatMatrix tmp = stack_complex_matrix (b);
2067  tmp = solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt);
2068  return unstack_complex_matrix (tmp);
2069 }
2070 
2071 FloatColumnVector
2072 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b) const
2073 {
2074  octave_idx_type info; float rcon;
2075  return solve (typ, b, info, rcon);
2076 }
2077 
2078 FloatColumnVector
2079 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b,
2080  octave_idx_type& info) const
2081 {
2082  float rcon;
2083  return solve (typ, b, info, rcon);
2084 }
2085 
2086 FloatColumnVector
2087 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b,
2088  octave_idx_type& info,
2089  float& rcon) const
2090 {
2091  return solve (typ, b, info, rcon, 0);
2092 }
2093 
2094 FloatColumnVector
2095 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b,
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 
2105 FloatComplexColumnVector
2106 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b) const
2107 {
2108  FloatComplexMatrix tmp (*this);
2109  return tmp.solve (typ, b);
2110 }
2111 
2112 FloatComplexColumnVector
2113 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
2114  octave_idx_type& info) const
2115 {
2116  FloatComplexMatrix tmp (*this);
2117  return tmp.solve (typ, b, info);
2118 }
2119 
2120 FloatComplexColumnVector
2121 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
2122  octave_idx_type& info, float& rcon) const
2123 {
2124  FloatComplexMatrix tmp (*this);
2125  return tmp.solve (typ, b, info, rcon);
2126 }
2127 
2128 FloatComplexColumnVector
2129 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
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 
2138 FloatMatrix
2139 FloatMatrix::solve (const FloatMatrix& b) const
2140 {
2141  octave_idx_type info;
2142  float rcon;
2143  return solve (b, info, rcon, 0);
2144 }
2145 
2146 FloatMatrix
2147 FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info) const
2148 {
2149  float rcon;
2150  return solve (b, info, rcon, 0);
2151 }
2152 
2153 FloatMatrix
2154 FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info,
2155  float& rcon) const
2156 {
2157  return solve (b, info, rcon, 0);
2158 }
2159 
2160 FloatMatrix
2161 FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info,
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 
2169 FloatComplexMatrix
2170 FloatMatrix::solve (const FloatComplexMatrix& b) const
2171 {
2172  FloatComplexMatrix tmp (*this);
2173  return tmp.solve (b);
2174 }
2175 
2176 FloatComplexMatrix
2177 FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info) const
2178 {
2179  FloatComplexMatrix tmp (*this);
2180  return tmp.solve (b, info);
2181 }
2182 
2183 FloatComplexMatrix
2184 FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info,
2185  float& rcon) const
2186 {
2187  FloatComplexMatrix tmp (*this);
2188  return tmp.solve (b, info, rcon);
2189 }
2190 
2191 FloatComplexMatrix
2192 FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info,
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 
2201 FloatColumnVector
2202 FloatMatrix::solve (const FloatColumnVector& b) const
2203 {
2204  octave_idx_type info; float rcon;
2205  return solve (b, info, rcon);
2206 }
2207 
2208 FloatColumnVector
2209 FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info) const
2210 {
2211  float rcon;
2212  return solve (b, info, rcon);
2213 }
2214 
2215 FloatColumnVector
2216 FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info,
2217  float& rcon) const
2218 {
2219  return solve (b, info, rcon, 0);
2220 }
2221 
2222 FloatColumnVector
2223 FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info,
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 
2231 FloatComplexColumnVector
2232 FloatMatrix::solve (const FloatComplexColumnVector& b) const
2233 {
2234  FloatComplexMatrix tmp (*this);
2235  return tmp.solve (b);
2236 }
2237 
2238 FloatComplexColumnVector
2239 FloatMatrix::solve (const FloatComplexColumnVector& b,
2240  octave_idx_type& info) const
2241 {
2242  FloatComplexMatrix tmp (*this);
2243  return tmp.solve (b, info);
2244 }
2245 
2246 FloatComplexColumnVector
2247 FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info,
2248  float& rcon) const
2249 {
2250  FloatComplexMatrix tmp (*this);
2251  return tmp.solve (b, info, rcon);
2252 }
2253 
2254 FloatComplexColumnVector
2255 FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info,
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 
2263 FloatMatrix
2264 FloatMatrix::lssolve (const FloatMatrix& b) const
2265 {
2266  octave_idx_type info;
2267  octave_idx_type rank;
2268  float rcon;
2269  return lssolve (b, info, rank, rcon);
2270 }
2271 
2272 FloatMatrix
2273 FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info) const
2274 {
2275  octave_idx_type rank;
2276  float rcon;
2277  return lssolve (b, info, rank, rcon);
2278 }
2279 
2280 FloatMatrix
2281 FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info,
2282  octave_idx_type& rank) const
2283 {
2284  float rcon;
2285  return lssolve (b, info, rank, rcon);
2286 }
2287 
2288 FloatMatrix
2289 FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info,
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 ())
2300  (*current_liboctave_error_handler)
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 
2422 FloatComplexMatrix
2423 FloatMatrix::lssolve (const FloatComplexMatrix& b) const
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 
2432 FloatComplexMatrix
2433 FloatMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info) const
2434 {
2435  FloatComplexMatrix tmp (*this);
2436  octave_idx_type rank;
2437  float rcon;
2438  return tmp.lssolve (b, info, rank, rcon);
2439 }
2440 
2441 FloatComplexMatrix
2442 FloatMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info,
2443  octave_idx_type& rank) const
2444 {
2445  FloatComplexMatrix tmp (*this);
2446  float rcon;
2447  return tmp.lssolve (b, info, rank, rcon);
2448 }
2449 
2450 FloatComplexMatrix
2451 FloatMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info,
2452  octave_idx_type& rank, float& rcon) const
2453 {
2454  FloatComplexMatrix tmp (*this);
2455  return tmp.lssolve (b, info, rank, rcon);
2456 }
2457 
2458 FloatColumnVector
2459 FloatMatrix::lssolve (const FloatColumnVector& b) const
2460 {
2461  octave_idx_type info;
2462  octave_idx_type rank;
2463  float rcon;
2464  return lssolve (b, info, rank, rcon);
2465 }
2466 
2467 FloatColumnVector
2468 FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info) const
2469 {
2470  octave_idx_type rank;
2471  float rcon;
2472  return lssolve (b, info, rank, rcon);
2473 }
2474 
2475 FloatColumnVector
2476 FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info,
2477  octave_idx_type& rank) const
2478 {
2479  float rcon;
2480  return lssolve (b, info, rank, rcon);
2481 }
2482 
2483 FloatColumnVector
2484 FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info,
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 ())
2495  (*current_liboctave_error_handler)
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 
2576 FloatComplexColumnVector
2577 FloatMatrix::lssolve (const FloatComplexColumnVector& b) const
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 
2586 FloatComplexColumnVector
2587 FloatMatrix::lssolve (const FloatComplexColumnVector& b,
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 
2596 FloatComplexColumnVector
2597 FloatMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info,
2598  octave_idx_type& rank) const
2599 {
2600  FloatComplexMatrix tmp (*this);
2601  float rcon;
2602  return tmp.lssolve (b, info, rank, rcon);
2603 }
2604 
2605 FloatComplexColumnVector
2606 FloatMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info,
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&
2614 FloatMatrix::operator += (const FloatDiagMatrix& a)
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&
2635 FloatMatrix::operator -= (const FloatDiagMatrix& a)
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 
2657 FloatMatrix
2658 operator * (const FloatColumnVector& v, const FloatRowVector& a)
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 
2684 FloatMatrix
2685 FloatMatrix::cumprod (int dim) const
2686 {
2687  return FloatNDArray::cumprod (dim);
2688 }
2689 
2690 FloatMatrix
2691 FloatMatrix::cumsum (int dim) const
2692 {
2693  return FloatNDArray::cumsum (dim);
2694 }
2695 
2696 FloatMatrix
2697 FloatMatrix::prod (int dim) const
2698 {
2699  return FloatNDArray::prod (dim);
2700 }
2701 
2702 FloatMatrix
2703 FloatMatrix::sum (int dim) const
2704 {
2705  return FloatNDArray::sum (dim);
2706 }
2707 
2708 FloatMatrix
2709 FloatMatrix::sumsq (int dim) const
2710 {
2711  return FloatNDArray::sumsq (dim);
2712 }
2713 
2714 FloatMatrix
2715 FloatMatrix::abs (void) const
2716 {
2717  return FloatNDArray::abs ();
2718 }
2719 
2720 FloatMatrix
2721 FloatMatrix::diag (octave_idx_type k) const
2722 {
2723  return FloatNDArray::diag (k);
2724 }
2725 
2726 FloatDiagMatrix
2727 FloatMatrix::diag (octave_idx_type m, octave_idx_type n) const
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
2737  (*current_liboctave_error_handler)
2738  ("diag: expecting vector argument");
2739 
2740  return retval;
2741 }
2742 
2743 FloatColumnVector
2744 FloatMatrix::row_min (void) const
2745 {
2746  Array<octave_idx_type> dummy_idx;
2747  return row_min (dummy_idx);
2748 }
2749 
2750 FloatColumnVector
2751 FloatMatrix::row_min (Array<octave_idx_type>& idx_arg) const
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 
2798 FloatColumnVector
2799 FloatMatrix::row_max (void) const
2800 {
2801  Array<octave_idx_type> dummy_idx;
2802  return row_max (dummy_idx);
2803 }
2804 
2805 FloatColumnVector
2806 FloatMatrix::row_max (Array<octave_idx_type>& idx_arg) const
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 
2853 FloatRowVector
2854 FloatMatrix::column_min (void) const
2855 {
2856  Array<octave_idx_type> dummy_idx;
2857  return column_min (dummy_idx);
2858 }
2859 
2860 FloatRowVector
2861 FloatMatrix::column_min (Array<octave_idx_type>& idx_arg) const
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 
2908 FloatRowVector
2909 FloatMatrix::column_max (void) const
2910 {
2911  Array<octave_idx_type> dummy_idx;
2912  return column_max (dummy_idx);
2913 }
2914 
2915 FloatRowVector
2916 FloatMatrix::column_max (Array<octave_idx_type>& idx_arg) const
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 
3003 FloatMatrix
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 
3020 FloatMatrix
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
3092 get_blas_trans_arg (bool trans)
3093 {
3094  return trans ? 'T' : 'N';
3095 }
3096 
3097 // the general GEMM operation
3098 
3099 FloatMatrix
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 
3187 FloatMatrix
3188 operator * (const FloatMatrix& a, const FloatMatrix& b)
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 
3199 FloatMatrix
3200 min (float d, const FloatMatrix& m)
3201 {
3202  octave_idx_type nr = m.rows ();
3203  octave_idx_type nc = m.columns ();
3204 
3205  EMPTY_RETURN_CHECK (FloatMatrix);
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 
3219 FloatMatrix
3220 min (const FloatMatrix& m, float d)
3221 {
3222  octave_idx_type nr = m.rows ();
3223  octave_idx_type nc = m.columns ();
3224 
3225  EMPTY_RETURN_CHECK (FloatMatrix);
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 
3239 FloatMatrix
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 
3252  EMPTY_RETURN_CHECK (FloatMatrix);
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 
3266 FloatMatrix
3267 max (float d, const FloatMatrix& m)
3268 {
3269  octave_idx_type nr = m.rows ();
3270  octave_idx_type nc = m.columns ();
3271 
3272  EMPTY_RETURN_CHECK (FloatMatrix);
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 
3286 FloatMatrix
3287 max (const FloatMatrix& m, float d)
3288 {
3289  octave_idx_type nr = m.rows ();
3290  octave_idx_type nc = m.columns ();
3291 
3292  EMPTY_RETURN_CHECK (FloatMatrix);
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 
3306 FloatMatrix
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 
3319  EMPTY_RETURN_CHECK (FloatMatrix);
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 
3333 FloatMatrix linspace (const FloatColumnVector& x1,
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)
3343  (*current_liboctave_error_handler)
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 
3367 MS_CMP_OPS (FloatMatrix, float)
3368 MS_BOOL_OPS (FloatMatrix, float)
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)
xgemm
FloatMatrix xgemm(const FloatMatrix &a, const FloatMatrix &b, blas_trans_type transa, blas_trans_type transb)
Definition: fMatrix.cc:3100
FloatMatrix::extract_n
FloatMatrix extract_n(octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const
Definition: fMatrix.cc:637
EMPTY_RETURN_CHECK
#define EMPTY_RETURN_CHECK(T)
Definition: fMatrix.cc:3195
F77_CHAR_ARG_LEN
#define F77_CHAR_ARG_LEN(l)
Definition: f77-fcn.h:253
FloatMatrix::operator-=
FloatMatrix & operator-=(const FloatDiagMatrix &a)
Definition: fMatrix.cc:2635
FloatComplexMatrix::solve
FloatComplexMatrix solve(MatrixType &typ, const FloatMatrix &b) const
Definition: fCMatrix.cc:2298
FloatMatrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
Definition: fMatrix.h:134
FloatMatrix::cumprod
FloatMatrix cumprod(int dim=-1) const
Definition: fMatrix.cc:2685
FloatMatrix::finverse
FloatMatrix finverse(MatrixType &mattype, octave_idx_type &info, float &rcon, int force, int calc_cond) const
Definition: fMatrix.cc:756
operator*
FloatMatrix operator*(const FloatColumnVector &v, const FloatRowVector &a)
Definition: fMatrix.cc:2658
gripe_nonconformant
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Definition: lo-array-gripes.cc:56
MM_BOOL_OPS
#define MM_BOOL_OPS(M1, M2)
Definition: mx-op-defs.h:210
FloatMatrix::abs
FloatMatrix abs(void) const
Definition: fMatrix.cc:2715
octave_fftw::fft
static int fft(const double *in, Complex *out, size_t npts, size_t nsamples=1, octave_idx_type stride=1, octave_idx_type dist=-1)
Definition: oct-fftw.cc:827
DiagArray2::elem
T elem(octave_idx_type r, octave_idx_type c) const
Definition: DiagArray2.h:110
operator<<
std::ostream & operator<<(std::ostream &os, const FloatMatrix &a)
Definition: fMatrix.cc:2964
FloatMatrix::insert
FloatMatrix & insert(const FloatMatrix &a, octave_idx_type r, octave_idx_type c)
Definition: fMatrix.cc:338
min
FloatMatrix min(float d, const FloatMatrix &m)
Definition: fMatrix.cc:3200
Givens
FloatMatrix Givens(float x, float y)
Definition: fMatrix.cc:3004
xisnan
bool xisnan(double x)
Definition: lo-mappers.cc:144
PermMatrix::rows
octave_idx_type rows(void) const
Definition: PermMatrix.h:55
FloatMatrix::sumsq
FloatMatrix sumsq(int dim=-1) const
Definition: fMatrix.cc:2709
octave_Float_Inf
#define octave_Float_Inf
Definition: lo-ieee.h:40
FloatMatrix::operator+=
FloatMatrix & operator+=(const FloatDiagMatrix &a)
Definition: fMatrix.cc:2614
idx_vector::colon
static const idx_vector colon
Definition: idx-vector.h:492
nn
static octave_idx_type nn
Definition: DASPK.cc:71
xslange
subroutine xslange(norm, m, n, a, lda, work, retval)
Definition: xslange.f:1
FloatNDArray::cumsum
FloatNDArray cumsum(int dim=-1) const
Definition: fNDArray.cc:619
FloatMatrix::fsolve
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
FloatMatrix transpose(void) const
Definition: fMatrix.h:118
mx_inline_real
void mx_inline_real(size_t n, T *r, const std::complex< T > *x)
Definition: mx-inlines.cc:251
SM_BOOL_OPS
#define SM_BOOL_OPS(S, M)
Definition: mx-op-defs.h:167
FloatMatrix::determinant
FloatDET determinant(void) const
Definition: fMatrix.cc:1243
FloatColumnVector::resize
void resize(octave_idx_type n, const float &rfv=0)
Definition: fColVector.h:109
FloatMatrix::fill
FloatMatrix & fill(float val)
Definition: fMatrix.cc:420
xmax
Complex xmax(const Complex &x, const Complex &y)
Definition: lo-mappers.cc:269
FloatMatrix::column_min
FloatRowVector column_min(void) const
Definition: fMatrix.cc:2854
operator>>
std::istream & operator>>(std::istream &is, FloatMatrix &a)
Definition: fMatrix.cc:2979
FloatNDArray::diag
FloatNDArray diag(octave_idx_type k=0) const
Definition: fNDArray.cc:831
DiagArray2::rows
octave_idx_type rows(void) const
Definition: DiagArray2.h:86
FloatMatrix::diag
FloatMatrix diag(octave_idx_type k=0) const
Definition: fMatrix.cc:2721
get_blas_trans_arg
static char get_blas_trans_arg(bool trans)
Definition: fMatrix.cc:3092
FloatMatrix::row_min
FloatColumnVector row_min(void) const
Definition: fMatrix.cc:2744
FloatMatrix::utsolve
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::FloatMatrix
FloatMatrix(void)
Definition: fMatrix.h:46
FloatNDArray::sumsq
FloatNDArray sumsq(int dim=-1) const
Definition: fNDArray.cc:649
xmin
Complex xmin(const Complex &x, const Complex &y)
Definition: lo-mappers.cc:263
FloatCHOL::rcond
float rcond(void) const
Definition: floatCHOL.h:67
xilaenv
subroutine xilaenv(ispec, name, opts, n1, n2, n3, n4, retval)
Definition: xilaenv.f:1
Array< float >::elem
float & elem(octave_idx_type n)
Definition: Array.h:380
MatrixType::is_hermitian
bool is_hermitian(void) const
Definition: MatrixType.h:122
FloatMatrix::solve
FloatMatrix solve(MatrixType &typ, const FloatMatrix &b) const
Definition: fMatrix.cc:1945
FloatMatrix::row
FloatRowVector row(octave_idx_type i) const
Definition: fMatrix.cc:646
octave_Float_NaN
#define octave_Float_NaN
Definition: lo-ieee.h:46
FloatMatrix::is_symmetric
bool is_symmetric(void) const
Definition: fMatrix.cc:322
FloatComplexMatrix::lssolve
FloatComplexMatrix lssolve(const FloatMatrix &b) const
Definition: fCMatrix.cc:2597
cfftb
subroutine cfftb(n, c, wsave)
Definition: cfftb.f:1
FloatMatrix::ifourier
FloatComplexMatrix ifourier(void) const
Definition: fMatrix.cc:945
FloatMatrix::tinverse
FloatMatrix tinverse(MatrixType &mattype, octave_idx_type &info, float &rcon, int force, int calc_cond) const
Definition: fMatrix.cc:698
xlog2
double xlog2(double x)
Definition: lo-mappers.cc:93
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
MS_CMP_OPS
#define MS_CMP_OPS(M, S)
Definition: mx-op-defs.h:107
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:313
d
F77_RET_T const double const double double * d
Definition: __pchip_deriv__.cc:38
FloatMatrix::extract
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
FloatNDArray abs(void) const
Definition: fNDArray.cc:792
FloatNDArray::sum
FloatNDArray sum(int dim=-1) const
Definition: fNDArray.cc:637
FloatMatrix::cumsum
FloatMatrix cumsum(int dim=-1) const
Definition: fMatrix.cc:2691
octave_fftw::fftNd
static int fftNd(const double *, Complex *, const int, const dim_vector &)
Definition: oct-fftw.cc:889
F77_CONST_CHAR_ARG2
#define F77_CONST_CHAR_ARG2(x, l)
Definition: f77-fcn.h:251
FloatMatrix::stack
FloatMatrix stack(const FloatMatrix &a) const
Definition: fMatrix.cc:539
mx_inline_imag
void mx_inline_imag(size_t n, T *r, const std::complex< T > *x)
Definition: mx-inlines.cc:254
FloatDET
base_det< float > FloatDET
Definition: DET.h:86
r2
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type r2
Definition: sparse-dmsolve.cc:160
cffti
subroutine cffti(n, wsave)
Definition: cffti.f:1
current_liboctave_error_handler
liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
MatrixType::type
int type(bool quiet=true)
Definition: MatrixType.cc:963
FloatMatrix::fourier
FloatComplexMatrix fourier(void) const
Definition: fMatrix.cc:916
octave_Inf
#define octave_Inf
Definition: lo-ieee.h:31
FloatSVD::left_singular_matrix
FloatMatrix left_singular_matrix(void) const
Definition: floatSVD.cc:58
Array< float >::make_unique
void make_unique(void)
Definition: Array.h:104
imag
FloatMatrix imag(const FloatComplexMatrix &a)
Definition: fMatrix.cc:621
PermMatrix::col_perm_vec
const Array< octave_idx_type > & col_perm_vec(void) const
Definition: PermMatrix.h:72
FloatSVD::right_singular_matrix
FloatMatrix right_singular_matrix(void) const
Definition: floatSVD.cc:71
real
FloatMatrix real(const FloatComplexMatrix &a)
Definition: fMatrix.cc:615
base_det
Definition: DET.h:31
FloatMatrix::ifourier2d
FloatComplexMatrix ifourier2d(void) const
Definition: fMatrix.cc:987
FloatColumnVector::extract
FloatColumnVector extract(octave_idx_type r1, octave_idx_type r2) const
Definition: fColVector.cc:169
octave_fftw::ifft
static int ifft(const Complex *in, Complex *out, size_t npts, size_t nsamples=1, octave_idx_type stride=1, octave_idx_type dist=-1)
Definition: oct-fftw.cc:865
Array< float >::data
const float * data(void) const
Definition: Array.h:479
FloatSVD::singular_values
FloatDiagMatrix singular_values(void) const
Definition: floatSVD.h:75
norm
double norm(const ColumnVector &v)
Definition: graphics.cc:5328
F77_CONST_CHAR_ARG_DECL
F77_RET_T F77_CONST_CHAR_ARG_DECL
Definition: fMatrix.cc:74
Array::resize
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1033
cfftf
subroutine cfftf(n, c, wsave)
Definition: cfftf.f:1
base_det::square
base_det square() const
Definition: DET.h:69
MatrixType::mark_as_unsymmetric
void mark_as_unsymmetric(void)
Definition: MatrixType.cc:1224
Array< float >::is_square
bool is_square(void) const
Definition: Array.h:470
xnorm
OCTAVE_API double xnorm(const ColumnVector &x, double p)
Definition: oct-norm.cc:536
gripe_singular_matrix
void gripe_singular_matrix(double rcond)
Definition: lo-array-gripes.cc:174
xsdot
subroutine xsdot(n, dx, incx, dy, incy, retval)
Definition: xsdot.f:1
FloatNDArray::prod
FloatNDArray prod(int dim=-1) const
Definition: fNDArray.cc:625
FloatDiagMatrix::extract_diag
FloatColumnVector extract_diag(octave_idx_type k=0) const
Definition: fDiagMatrix.h:106
F77_RET_T
#define F77_RET_T
Definition: f77-fcn.h:264
FloatMatrix::column
FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:652
FloatMatrix::row_max
FloatColumnVector row_max(void) const
Definition: fMatrix.cc:2799
MatrixType::mark_as_rectangular
void mark_as_rectangular(void)
Definition: MatrixType.h:155
FloatMatrix::FloatComplexMatrix
friend class FloatComplexMatrix
Definition: fMatrix.h:116
octave_fftw::ifftNd
static int ifftNd(const Complex *, Complex *, const int, const dim_vector &)
Definition: oct-fftw.cc:933
Array< float >::xelem
float & xelem(octave_idx_type n)
Definition: Array.h:353
FloatMatrix::sum
FloatMatrix sum(int dim=-1) const
Definition: fMatrix.cc:2703
DiagArray2::cols
octave_idx_type cols(void) const
Definition: DiagArray2.h:87
FloatSCHUR::schur_matrix
FloatMatrix schur_matrix(void) const
Definition: floatSCHUR.h:71
Array< float >::length
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
FloatMatrix::column_max
FloatRowVector column_max(void) const
Definition: fMatrix.cc:2909
FloatMatrix::inverse
FloatMatrix inverse(void) const
Definition: fMatrix.cc:658
FloatNDArray::cumprod
FloatNDArray cumprod(int dim=-1) const
Definition: fNDArray.cc:613
NoAlias
This is a simple wrapper template that will subclass an Array type or any later type derived from ...
Definition: Array.h:762
octave_write_float
void octave_write_float(std::ostream &os, float d)
Definition: lo-utils.cc:409
V
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
FloatSCHUR::unitary_matrix
FloatMatrix unitary_matrix(void) const
Definition: floatSCHUR.h:73
F77_FUNC
F77_RET_T F77_FUNC(xilaenv, XILAENV)(const octave_idx_type &
octave_NaN
#define octave_NaN
Definition: lo-ieee.h:37
FloatRowVector::resize
void resize(octave_idx_type n, const float &rfv=0)
Definition: fRowVector.h:99
FloatCHOL::inverse
FloatMatrix inverse(void) const
Definition: floatCHOL.cc:187
unstack_complex_matrix
static FloatComplexMatrix unstack_complex_matrix(const FloatMatrix &sm)
Definition: fMatrix.cc:2047
FloatMatrix::lssolve
FloatMatrix lssolve(const FloatMatrix &b) const
Definition: fMatrix.cc:2264
FloatRowVector::max
float max(void) const
Definition: fRowVector.cc:246
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5281
FloatMatrix::prod
FloatMatrix prod(int dim=-1) const
Definition: fMatrix.cc:2697
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
MS_BOOL_OPS
#define MS_BOOL_OPS(M, S)
Definition: mx-op-defs.h:124
MM_CMP_OPS
#define MM_CMP_OPS(M1, M2)
Definition: mx-op-defs.h:193
get_blas_char
char get_blas_char(blas_trans_type transt)
Definition: mx-defs.h:136
octave_idx_type
FloatMatrix::operator!=
bool operator!=(const FloatMatrix &a) const
Definition: fMatrix.cc:316
FloatNDArray::insert
FloatNDArray & insert(const FloatNDArray &a, octave_idx_type r, octave_idx_type c)
Definition: fNDArray.cc:776
FloatMatrix::rcond
float rcond(void) const
Definition: fMatrix.cc:1403
idx_vector
Definition: idx-vector.h:50
FloatMatrix::ltsolve
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
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
blas_trans_type
blas_trans_type
Definition: mx-defs.h:128
c1
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type octave_idx_type c1
Definition: sparse-dmsolve.cc:160
r1
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type r1
Definition: sparse-dmsolve.cc:160
Array::fortran_vec
const T * fortran_vec(void) const
Definition: Array.h:481
Sylvester
FloatMatrix Sylvester(const FloatMatrix &a, const FloatMatrix &b, const FloatMatrix &c)
Definition: fMatrix.cc:3021
mx_inline_equal
bool mx_inline_equal(size_t n, const T1 *x, const T2 *y)
Definition: mx-inlines.cc:441
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:321
FloatComplexSVD::sigma
FloatDiagMatrix sigma
Definition: fCmplxSVD.h:90
FloatMatrix::append
FloatMatrix append(const FloatMatrix &a) const
Definition: fMatrix.cc:467
dim_vector
Definition: dim-vector.h:53
stack_complex_matrix
static FloatMatrix stack_complex_matrix(const FloatComplexMatrix &cm)
Definition: fMatrix.cc:2030
FloatMatrix::fourier2d
FloatComplexMatrix fourier2d(void) const
Definition: fMatrix.cc:975
Array< float >::columns
octave_idx_type columns(void) const
Definition: Array.h:322
DiagArray2::length
octave_idx_type length(void) const
Definition: DiagArray2.h:92
Array< float >::index
Array< float > index(const idx_vector &i) const
Indexing without resizing.
F77_CHAR_ARG_LEN_DECL
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
SM_CMP_OPS
#define SM_CMP_OPS(S, M)
Definition: mx-op-defs.h:150
x
F77_RET_T const double * x
Definition: __pchip_deriv__.cc:38
linspace
FloatMatrix linspace(const FloatColumnVector &x1, const FloatColumnVector &x2, octave_idx_type n)
Definition: fMatrix.cc:3333
max
FloatMatrix max(float d, const FloatMatrix &m)
Definition: fMatrix.cc:3267
FloatMatrix::pseudo_inverse
FloatMatrix pseudo_inverse(float tol=0.0) const
Definition: fMatrix.cc:877
FloatMatrix::operator==
bool operator==(const FloatMatrix &a) const
Definition: fMatrix.cc:307
Generated on Thu Jun 4 2015 23:30:23 for GNU Octave by   doxygen 1.8.8