GNU Octave  4.0.0
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
det.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "DET.h"
28 
29 #include "defun.h"
30 #include "error.h"
31 #include "gripes.h"
32 #include "oct-obj.h"
33 #include "utils.h"
34 #include "ops.h"
35 
36 #include "ov-re-mat.h"
37 #include "ov-cx-mat.h"
38 #include "ov-flt-re-mat.h"
39 #include "ov-flt-cx-mat.h"
40 #include "ov-re-diag.h"
41 #include "ov-cx-diag.h"
42 #include "ov-flt-re-diag.h"
43 #include "ov-flt-cx-diag.h"
44 #include "ov-perm.h"
45 
46 #define MAYBE_CAST(VAR, CLASS) \
47  const CLASS *VAR = arg.type_id () == CLASS::static_type_id () ? \
48  dynamic_cast<const CLASS *> (&arg.get_rep ()) : 0
49 
50 DEFUN (det, args, nargout,
51  "-*- texinfo -*-\n\
52 @deftypefn {Built-in Function} {} det (@var{A})\n\
53 @deftypefnx {Built-in Function} {[@var{d}, @var{rcond}] =} det (@var{A})\n\
54 Compute the determinant of @var{A}.\n\
55 \n\
56 Return an estimate of the reciprocal condition number if requested.\n\
57 \n\
58 Programming Notes: Routines from @sc{lapack} are used for full matrices and\n\
59 code from @sc{umfpack} is used for sparse matrices.\n\
60 \n\
61 The determinant should not be used to check a matrix for singularity.\n\
62 For that, use any of the condition number functions: @code{cond},\n\
63 @code{condest}, @code{rcond}.\n\
64 @seealso{cond, condest, rcond}\n\
65 @end deftypefn")
66 {
67  octave_value_list retval;
68 
69  int nargin = args.length ();
70 
71  if (nargin != 1)
72  {
73  print_usage ();
74  return retval;
75  }
76 
77  octave_value arg = args(0);
78 
79  octave_idx_type nr = arg.rows ();
80  octave_idx_type nc = arg.columns ();
81 
82  if (nr == 0 && nc == 0)
83  {
84  retval(0) = 1.0;
85  return retval;
86  }
87 
88  int arg_is_empty = empty_arg ("det", nr, nc);
89  if (arg_is_empty < 0)
90  return retval;
91  if (arg_is_empty > 0)
92  return octave_value (Matrix (1, 1, 1.0));
93 
94 
95  if (nr != nc)
96  {
98  return retval;
99  }
100 
101  bool isfloat = arg.is_single_type ();
102 
103  if (arg.is_diag_matrix ())
104  {
105  if (arg.is_complex_type ())
106  {
107  if (isfloat)
108  {
109  retval(0) = arg.float_complex_diag_matrix_value ()
110  .determinant ().value ();
111  if (nargout > 1)
112  retval(1) = arg.float_complex_diag_matrix_value ().rcond ();
113  }
114  else
115  {
116  retval(0) = arg.complex_diag_matrix_value ()
117  .determinant ().value ();
118  if (nargout > 1)
119  retval(1) = arg.complex_diag_matrix_value ().rcond ();
120  }
121  }
122  else
123  {
124  if (isfloat)
125  {
126  retval(0) = arg.float_diag_matrix_value ()
127  .determinant ().value ();
128  if (nargout > 1)
129  retval(1) = arg.float_diag_matrix_value ().rcond ();
130  }
131  else
132  {
133  retval(0) = arg.diag_matrix_value ().determinant ().value ();
134  if (nargout > 1)
135  retval(1) = arg.diag_matrix_value ().rcond ();
136  }
137  }
138  }
139  else if (arg.is_perm_matrix ())
140  {
141  retval(0) = static_cast<double> (arg.perm_matrix_value ().determinant ());
142  if (nargout > 1)
143  retval(1) = 1.0;
144  }
145  else if (arg.is_single_type ())
146  {
147  if (arg.is_real_type ())
148  {
149  octave_idx_type info;
150  float rcond = 0.0;
151  // Always compute rcond, so we can detect numerically
152  // singular matrices.
153  FloatMatrix m = arg.float_matrix_value ();
154  if (! error_state)
155  {
157  MatrixType mtype = rep ? rep -> matrix_type () : MatrixType ();
158  FloatDET det = m.determinant (mtype, info, rcond);
159  retval(1) = rcond;
160  retval(0) = info == -1 ? 0.0f : det.value ();
161  if (rep) rep->matrix_type (mtype);
162  }
163  }
164  else if (arg.is_complex_type ())
165  {
166  octave_idx_type info;
167  float rcond = 0.0;
168  // Always compute rcond, so we can detect numerically
169  // singular matrices.
171  if (! error_state)
172  {
174  MatrixType mtype = rep ? rep -> matrix_type () : MatrixType ();
175  FloatComplexDET det = m.determinant (mtype, info, rcond);
176  retval(1) = rcond;
177  retval(0) = info == -1 ? FloatComplex (0.0) : det.value ();
178  if (rep) rep->matrix_type (mtype);
179  }
180  }
181  }
182  else
183  {
184  if (arg.is_real_type ())
185  {
186  octave_idx_type info;
187  double rcond = 0.0;
188  // Always compute rcond, so we can detect numerically
189  // singular matrices.
190  if (arg.is_sparse_type ())
191  {
193  if (! error_state)
194  {
195  DET det = m.determinant (info, rcond);
196  retval(1) = rcond;
197  retval(0) = info == -1 ? 0.0 : det.value ();
198  }
199  }
200  else
201  {
202  Matrix m = arg.matrix_value ();
203  if (! error_state)
204  {
205  MAYBE_CAST (rep, octave_matrix);
206  MatrixType mtype = rep ? rep -> matrix_type ()
207  : MatrixType ();
208  DET det = m.determinant (mtype, info, rcond);
209  retval(1) = rcond;
210  retval(0) = info == -1 ? 0.0 : det.value ();
211  if (rep) rep->matrix_type (mtype);
212  }
213  }
214  }
215  else if (arg.is_complex_type ())
216  {
217  octave_idx_type info;
218  double rcond = 0.0;
219  // Always compute rcond, so we can detect numerically
220  // singular matrices.
221  if (arg.is_sparse_type ())
222  {
224  if (! error_state)
225  {
226  ComplexDET det = m.determinant (info, rcond);
227  retval(1) = rcond;
228  retval(0) = info == -1 ? Complex (0.0) : det.value ();
229  }
230  }
231  else
232  {
234  if (! error_state)
235  {
237  MatrixType mtype = rep ? rep -> matrix_type ()
238  : MatrixType ();
239  ComplexDET det = m.determinant (mtype, info, rcond);
240  retval(1) = rcond;
241  retval(0) = info == -1 ? Complex (0.0) : det.value ();
242  if (rep) rep->matrix_type (mtype);
243  }
244  }
245  }
246  else
247  gripe_wrong_type_arg ("det", arg);
248  }
249  return retval;
250 }
251 
252 /*
253 %!assert (det ([1, 2; 3, 4]), -2, 10*eps)
254 %!assert (det (single ([1, 2; 3, 4])), single (-2), 10*eps ("single"))
255 %!error det ()
256 %!error det (1, 2)
257 %!error <argument must be a square matrix> det ([1, 2; 3, 4; 5, 6])
258 */
FloatComplexDiagMatrix float_complex_diag_matrix_value(bool force=false) const
Definition: ov.h:840
bool is_real_type(void) const
Definition: ov.h:651
ComplexDET determinant(void) const
Definition: CMatrix.cc:1590
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
octave_idx_type rows(void) const
Definition: ov.h:473
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
float rcond(void) const
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
FloatDET determinant(void) const
Definition: fMatrix.cc:1243
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
bool is_perm_matrix(void) const
Definition: ov.h:559
DET determinant(void) const
Definition: dDiagMatrix.cc:347
FloatDET determinant(void) const
Definition: fDiagMatrix.cc:347
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:246
void gripe_square_matrix_required(const char *name)
Definition: gripes.cc:81
T value() const
Definition: DET.h:66
ComplexDiagMatrix complex_diag_matrix_value(bool force=false) const
Definition: ov.h:836
octave_idx_type determinant(void) const
Definition: PermMatrix.cc:138
octave_idx_type columns(void) const
Definition: ov.h:475
DET determinant(void) const
Definition: dSparse.cc:1246
bool is_sparse_type(void) const
Definition: ov.h:666
FloatDiagMatrix float_diag_matrix_value(bool force=false) const
Definition: ov.h:833
Definition: DET.h:31
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
ComplexDET determinant(void) const
Definition: CSparse.cc:1155
#define MAYBE_CAST(VAR, CLASS)
Definition: det.cc:46
float rcond(void) const
Definition: fDiagMatrix.cc:366
Definition: dMatrix.h:35
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
FloatComplexDET determinant(void) const
double arg(double x)
Definition: lo-mappers.h:37
double rcond(void) const
Definition: dDiagMatrix.cc:366
DET determinant(void) const
Definition: dMatrix.cc:1234
double rcond(void) const
Definition: CDiagMatrix.cc:551
DiagMatrix diag_matrix_value(bool force=false) const
Definition: ov.h:830
ComplexDET determinant(void) const
Definition: CDiagMatrix.cc:532
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
PermMatrix perm_matrix_value(void) const
Definition: ov.h:843
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:820
FloatComplexDET determinant(void) const
Definition: fCMatrix.cc:1594
std::complex< double > Complex
Definition: oct-cmplx.h:29
bool is_single_type(void) const
Definition: ov.h:611
bool is_diag_matrix(void) const
Definition: ov.h:556
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))