GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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