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