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