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