Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/det.cc @ 8337:e02242c54c49

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reuse matrix type detected in det
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 19 Nov 2008 16:55:47 +0100
parents 64cf956a109c
children c3f7e2549abb
comparison
equal deleted inserted replaced
8336:9813c07ca946 8337:e02242c54c49
30 #include "defun-dld.h" 30 #include "defun-dld.h"
31 #include "error.h" 31 #include "error.h"
32 #include "gripes.h" 32 #include "gripes.h"
33 #include "oct-obj.h" 33 #include "oct-obj.h"
34 #include "utils.h" 34 #include "utils.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
41 #define MAYBE_CAST(VAR, CLASS) \
42 const CLASS *VAR = arg.type_id () == CLASS::static_type_id () ? \
43 dynamic_cast<const CLASS *> (&arg.get_rep ()) : 0
35 44
36 DEFUN_DLD (det, args, , 45 DEFUN_DLD (det, args, ,
37 "-*- texinfo -*-\n\ 46 "-*- texinfo -*-\n\
38 @deftypefn {Loadable Function} {[@var{d}, @var{rcond}] =} det (@var{a})\n\ 47 @deftypefn {Loadable Function} {[@var{d}, @var{rcond}] =} det (@var{a})\n\
39 Compute the determinant of @var{a} using @sc{Lapack} for full and UMFPACK\n\ 48 Compute the determinant of @var{a} using @sc{Lapack} for full and UMFPACK\n\
66 if (arg_is_empty < 0) 75 if (arg_is_empty < 0)
67 return retval; 76 return retval;
68 if (arg_is_empty > 0) 77 if (arg_is_empty > 0)
69 return octave_value (Matrix (1, 1, 1.0)); 78 return octave_value (Matrix (1, 1, 1.0));
70 79
80
71 if (nr != nc) 81 if (nr != nc)
72 { 82 {
73 gripe_square_matrix_required ("det"); 83 gripe_square_matrix_required ("det");
74 return retval; 84 return retval;
75 } 85 }
84 // Always compute rcond, so we can detect numerically 94 // Always compute rcond, so we can detect numerically
85 // singular matrices. 95 // singular matrices.
86 FloatMatrix m = arg.float_matrix_value (); 96 FloatMatrix m = arg.float_matrix_value ();
87 if (! error_state) 97 if (! error_state)
88 { 98 {
89 FloatDET det = m.determinant (info, rcond); 99 MAYBE_CAST (rep, octave_float_matrix);
100 MatrixType mtype = rep ? rep -> matrix_type () : MatrixType ();
101 FloatDET det = m.determinant (mtype, info, rcond);
90 retval(1) = rcond; 102 retval(1) = rcond;
91 retval(0) = info == -1 ? static_cast<float>(0.0) : det.value (); 103 retval(0) = info == -1 ? static_cast<float>(0.0) : det.value ();
104 if (rep) rep->matrix_type (mtype);
92 } 105 }
93 } 106 }
94 else if (arg.is_complex_type ()) 107 else if (arg.is_complex_type ())
95 { 108 {
96 octave_idx_type info; 109 octave_idx_type info;
98 // Always compute rcond, so we can detect numerically 111 // Always compute rcond, so we can detect numerically
99 // singular matrices. 112 // singular matrices.
100 FloatComplexMatrix m = arg.float_complex_matrix_value (); 113 FloatComplexMatrix m = arg.float_complex_matrix_value ();
101 if (! error_state) 114 if (! error_state)
102 { 115 {
103 FloatComplexDET det = m.determinant (info, rcond); 116 MAYBE_CAST (rep, octave_float_complex_matrix);
117 MatrixType mtype = rep ? rep -> matrix_type () : MatrixType ();
118 FloatComplexDET det = m.determinant (mtype, info, rcond);
104 retval(1) = rcond; 119 retval(1) = rcond;
105 retval(0) = info == -1 ? FloatComplex (0.0) : det.value (); 120 retval(0) = info == -1 ? FloatComplex (0.0) : det.value ();
121 if (rep) rep->matrix_type (mtype);
106 } 122 }
107 } 123 }
108 } 124 }
109 else 125 else
110 { 126 {
127 else 143 else
128 { 144 {
129 Matrix m = arg.matrix_value (); 145 Matrix m = arg.matrix_value ();
130 if (! error_state) 146 if (! error_state)
131 { 147 {
132 DET det = m.determinant (info, rcond); 148 MAYBE_CAST (rep, octave_matrix);
149 MatrixType mtype = rep ? rep -> matrix_type () : MatrixType ();
150 DET det = m.determinant (mtype, info, rcond);
133 retval(1) = rcond; 151 retval(1) = rcond;
134 retval(0) = info == -1 ? 0.0 : det.value (); 152 retval(0) = info == -1 ? 0.0 : det.value ();
153 if (rep) rep->matrix_type (mtype);
135 } 154 }
136 } 155 }
137 } 156 }
138 else if (arg.is_complex_type ()) 157 else if (arg.is_complex_type ())
139 { 158 {
154 else 173 else
155 { 174 {
156 ComplexMatrix m = arg.complex_matrix_value (); 175 ComplexMatrix m = arg.complex_matrix_value ();
157 if (! error_state) 176 if (! error_state)
158 { 177 {
159 ComplexDET det = m.determinant (info, rcond); 178 MAYBE_CAST (rep, octave_complex_matrix);
179 MatrixType mtype = rep ? rep -> matrix_type () : MatrixType ();
180 ComplexDET det = m.determinant (mtype, info, rcond);
160 retval(1) = rcond; 181 retval(1) = rcond;
161 retval(0) = info == -1 ? Complex (0.0) : det.value (); 182 retval(0) = info == -1 ? Complex (0.0) : det.value ();
183 if (rep) rep->matrix_type (mtype);
162 } 184 }
163 } 185 }
164 } 186 }
165 else 187 else
166 gripe_wrong_type_arg ("det", arg); 188 gripe_wrong_type_arg ("det", arg);