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