GNU Octave  3.8.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-2013 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 Routines from @sc{lapack} are used for full matrices and code from\n\
59 @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 ? static_cast<float>(0.0) : 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 */