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