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