1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1996-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 #if defined (HAVE_CONFIG_H) 27 # include "config.h" 28 #endif 29 30 #include "DET.h" 31 32 #include "defun.h" 33 #include "error.h" 34 #include "errwarn.h" 35 #include "ovl.h" 36 #include "ops.h" 37 38 #include "ov-re-mat.h" 39 #include "ov-cx-mat.h" 40 #include "ov-flt-re-mat.h" 41 #include "ov-flt-cx-mat.h" 42 #include "ov-re-diag.h" 43 #include "ov-cx-diag.h" 44 #include "ov-flt-re-diag.h" 45 #include "ov-flt-cx-diag.h" 46 #include "ov-perm.h" 47 48 #define MAYBE_CAST(VAR, CLASS) \ 49 const CLASS *VAR = (arg.type_id () == CLASS::static_type_id () \ 50 ? dynamic_cast<const CLASS *> (&arg.get_rep ()) \ 51 : nullptr) 52 53 DEFUN (det, args, nargout, 54 doc: /* -*- texinfo -*- 55 @deftypefn {} {} det (@var{A}) 56 @deftypefnx {} {[@var{d}, @var{rcond}] =} det (@var{A}) 57 Compute the determinant of @var{A}. 58 59 Return an estimate of the reciprocal condition number if requested. 60 61 Programming Notes: Routines from @sc{lapack} are used for full matrices and 62 code from @sc{umfpack} is used for sparse matrices. 63 64 The determinant should not be used to check a matrix for singularity. 65 For that, use any of the condition number functions: @code{cond}, 66 @code{condest}, @code{rcond}. 67 @seealso{cond, condest, rcond} 68 @end deftypefn */) 69 { 70 if (args.length () != 1) 71 print_usage (); 72 73 octave_value arg = args(0); 74 75 if (arg.isempty ()) 76 return ovl (1.0); 77 78 if (arg.rows () != arg.columns ()) 79 err_square_matrix_required ("det", "A"); 80 81 octave_value_list retval (2); 82 83 bool isfloat = arg.is_single_type (); 84 85 if (arg.is_diag_matrix ()) 86 { 87 if (nargout <= 1) 88 retval.resize (1); 89 90 if (arg.iscomplex ()) 91 { 92 if (isfloat) 93 { 94 retval(0) = arg.float_complex_diag_matrix_value () 95 .determinant ().value (); 96 if (nargout > 1) 97 retval(1) = arg.float_complex_diag_matrix_value ().rcond (); 98 } 99 else 100 { 101 retval(0) = arg.complex_diag_matrix_value () 102 .determinant ().value (); 103 if (nargout > 1) 104 retval(1) = arg.complex_diag_matrix_value ().rcond (); 105 } 106 } 107 else 108 { 109 if (isfloat) 110 { 111 retval(0) = arg.float_diag_matrix_value () 112 .determinant ().value (); 113 if (nargout > 1) 114 retval(1) = arg.float_diag_matrix_value ().rcond (); 115 } 116 else 117 { 118 retval(0) = arg.diag_matrix_value ().determinant ().value (); 119 if (nargout > 1) 120 retval(1) = arg.diag_matrix_value ().rcond (); 121 } 122 } 123 } 124 else if (arg.is_perm_matrix ()) 125 { 126 if (nargout <= 1) 127 retval.resize (1); 128 129 retval(0) = static_cast<double> (arg.perm_matrix_value ().determinant ()); 130 if (nargout > 1) 131 retval(1) = 1.0; 132 } 133 else if (arg.is_single_type ()) 134 { 135 if (arg.isreal ()) 136 { 137 octave_idx_type info; 138 float rcond = 0.0; 139 // Always compute rcond, so we can detect singular matrices. 140 FloatMatrix m = arg.float_matrix_value (); 141 142 MAYBE_CAST (rep, octave_float_matrix); 143 MatrixType mtype = (rep ? rep -> matrix_type () : MatrixType ()); 144 FloatDET det = m.determinant (mtype, info, rcond); 145 retval(0) = (info == -1 ? 0.0f : det.value ()); 146 retval(1) = rcond; 147 if (rep) 148 rep->matrix_type (mtype); 149 } 150 else if (arg.iscomplex ()) 151 { 152 octave_idx_type info; 153 float rcond = 0.0; 154 // Always compute rcond, so we can detect singular matrices. 155 FloatComplexMatrix m = arg.float_complex_matrix_value (); 156 157 MAYBE_CAST (rep, octave_float_complex_matrix); 158 MatrixType mtype = (rep ? rep -> matrix_type () : MatrixType ()); 159 FloatComplexDET det = m.determinant (mtype, info, rcond); 160 retval(0) = (info == -1 ? FloatComplex (0.0) : det.value ()); 161 retval(1) = rcond; 162 if (rep) 163 rep->matrix_type (mtype); 164 } 165 } 166 else 167 { 168 if (arg.isreal ()) 169 { 170 octave_idx_type info; 171 double rcond = 0.0; 172 // Always compute rcond, so we can detect singular matrices. 173 if (arg.issparse ()) 174 { 175 SparseMatrix m = arg.sparse_matrix_value (); 176 177 DET det = m.determinant (info, rcond); 178 retval(0) = (info == -1 ? 0.0 : det.value ()); 179 retval(1) = rcond; 180 } 181 else 182 { 183 Matrix m = arg.matrix_value (); 184 185 MAYBE_CAST (rep, octave_matrix); 186 MatrixType mtype = (rep ? rep -> matrix_type () 187 : MatrixType ()); 188 DET det = m.determinant (mtype, info, rcond); 189 retval(0) = (info == -1 ? 0.0 : det.value ()); 190 retval(1) = rcond; 191 if (rep) 192 rep->matrix_type (mtype); 193 } 194 } 195 else if (arg.iscomplex ()) 196 { 197 octave_idx_type info; 198 double rcond = 0.0; 199 // Always compute rcond, so we can detect singular matrices. 200 if (arg.issparse ()) 201 { 202 SparseComplexMatrix m = arg.sparse_complex_matrix_value (); 203 204 ComplexDET det = m.determinant (info, rcond); 205 retval(0) = (info == -1 ? Complex (0.0) : det.value ()); 206 retval(1) = rcond; 207 } 208 else 209 { 210 ComplexMatrix m = arg.complex_matrix_value (); 211 212 MAYBE_CAST (rep, octave_complex_matrix); 213 MatrixType mtype = (rep ? rep -> matrix_type () 214 : MatrixType ()); 215 ComplexDET det = m.determinant (mtype, info, rcond); 216 retval(0) = (info == -1 ? Complex (0.0) : det.value ()); 217 retval(1) = rcond; 218 if (rep) 219 rep->matrix_type (mtype); 220 } 221 } 222 else 223 err_wrong_type_arg ("det", arg); 224 } 225 226 return retval; 227 } 228 229 /* 230 %!assert (det ([1, 2; 3, 4]), -2, 10*eps) 231 %!assert (det (single ([1, 2; 3, 4])), single (-2), 10*eps ("single")) 232 %!assert (det (eye (2000)), 1) 233 %!error det () 234 %!error det (1, 2) 235 %!error <must be a square matrix> det ([1, 2; 3, 4; 5, 6]) 236 */ 237