Mercurial > hg > octave-lyh

comparison src/DLD-FUNCTIONS/det.cc @ 7505:f5005d9510f4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove dispatched sparse functions and treat in the generic versions of the functions
author David Bateman <dbateman@free.fr>
date Wed, 20 Feb 2008 15:52:11 -0500
parents a1dbe9d80eee
children eb7bdde776f2
comparison
equal deleted inserted replaced
7504:ddcf233d765b 7505:f5005d9510f4
35 #include "utils.h" 35 #include "utils.h"
36 36
37 DEFUN_DLD (det, args, , 37 DEFUN_DLD (det, args, ,
38 "-*- texinfo -*-\n\ 38 "-*- texinfo -*-\n\
39 @deftypefn {Loadable Function} {[@var{d}, @var{rcond}] = } det (@var{a})\n\ 39 @deftypefn {Loadable Function} {[@var{d}, @var{rcond}] = } det (@var{a})\n\
40 Compute the determinant of @var{a} using @sc{Lapack}. Return an estimate\n\ 40 Compute the determinant of @var{a} using @sc{Lapack} for full and UMFPACK\n\
41 of the reciprocal condition number if requested.\n\ 41 for sparse matrices. Return an estimate of the reciprocal condition number\n\
42 if requested.\n\
42 @end deftypefn") 43 @end deftypefn")
43 { 44 {
44 octave_value_list retval; 45 octave_value_list retval;
45 46
46 int nargin = args.length (); 47 int nargin = args.length ();
51 return retval; 52 return retval;
52 } 53 }
53 54
54 octave_value arg = args(0); 55 octave_value arg = args(0);
55 56
56 int nr = arg.rows (); 57 octave_idx_type nr = arg.rows ();
57 int nc = arg.columns (); 58 octave_idx_type nc = arg.columns ();
58 59
59 if (nr == 0 && nc == 0) 60 if (nr == 0 && nc == 0)
60 { 61 {
61 retval(0) = 1.0; 62 retval(0) = 1.0;
62 return retval; 63 return retval;
74 return retval; 75 return retval;
75 } 76 }
76 77
77 if (arg.is_real_type ()) 78 if (arg.is_real_type ())
78 { 79 {
79 Matrix m = arg.matrix_value (); 80 octave_idx_type info;
80 81 double rcond = 0.0;
81 if (! error_state) 82 // Always compute rcond, so we can detect numerically
83 // singular matrices.
84 if (arg.is_sparse_type ())
82 { 85 {
83 // Always compute rcond, so we can detect numerically 86 SparseMatrix m = arg.sparse_matrix_value ();
84 // singular matrices. 87 if (! error_state)
85 88 {
86 octave_idx_type info; 89 DET det = m.determinant (info, rcond);
87 double rcond = 0.0; 90 retval(1) = rcond;
88 DET det = m.determinant (info, rcond); 91 volatile double xrcond = rcond;
89 retval(1) = rcond; 92 xrcond += 1.0;
90 volatile double xrcond = rcond; 93 retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ());
91 xrcond += 1.0; 94 }
92 retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ()); 95 }
96 else
97 {
98 Matrix m = arg.matrix_value ();
99 if (! error_state)
100 {
101 DET det = m.determinant (info, rcond);
102 retval(1) = rcond;
103 volatile double xrcond = rcond;
104 xrcond += 1.0;
105 retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ());
106 }
93 } 107 }
94 } 108 }
95 else if (arg.is_complex_type ()) 109 else if (arg.is_complex_type ())
96 { 110 {
97 ComplexMatrix m = arg.complex_matrix_value (); 111 octave_idx_type info;
112 double rcond = 0.0;
113 // Always compute rcond, so we can detect numerically
114 // singular matrices.
115 if (arg.is_sparse_type ())
116 {
117 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
118 if (! error_state)
119 {
120 ComplexDET det = m.determinant (info, rcond);
121 retval(1) = rcond;
122 volatile double xrcond = rcond;
123 xrcond += 1.0;
124 retval(0) = ((info == -1 || xrcond == 1.0)
125 ? Complex (0.0) : det.value ());
126 }
127 }
128 else
129 {
130 ComplexMatrix m = arg.complex_matrix_value ();
131 if (! error_state)
132 {
133 ComplexDET det = m.determinant (info, rcond);
134 retval(1) = rcond;
135 volatile double xrcond = rcond;
136 xrcond += 1.0;
137 retval(0) = ((info == -1 || xrcond == 1.0)
138 ? Complex (0.0) : det.value ());
98 139
99 if (! error_state) 140 }
100 {
101 // Always compute rcond, so we can detect numerically
102 // singular matrices.
103
104 octave_idx_type info;
105 double rcond = 0.0;
106 ComplexDET det = m.determinant (info, rcond);
107 retval(1) = rcond;
108 volatile double xrcond = rcond;
109 xrcond += 1.0;
110 retval(0) = ((info == -1 || xrcond == 1.0)
111 ? Complex (0.0) : det.value ());
112 } 141 }
113 } 142 }
114 else 143 else
115 gripe_wrong_type_arg ("det", arg); 144 gripe_wrong_type_arg ("det", arg);
116 145