File: | libinterp/corefcn/inv.cc |
Location: | line 218, column 32 |
Description: | The left operand of '==' is a garbage value |
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_H1 | |||
24 | #include <config.h> | |||
25 | #endif | |||
26 | ||||
27 | #include "defun.h" | |||
28 | #include "error.h" | |||
29 | #include "gripes.h" | |||
30 | #include "oct-obj.h" | |||
31 | #include "ops.h" | |||
32 | #include "ov-re-diag.h" | |||
33 | #include "ov-cx-diag.h" | |||
34 | #include "ov-flt-re-diag.h" | |||
35 | #include "ov-flt-cx-diag.h" | |||
36 | #include "ov-perm.h" | |||
37 | #include "utils.h" | |||
38 | ||||
39 | DEFUN (inv, args, nargout,octave_value_list Finv (const octave_value_list& args, int nargout) | |||
40 | "-*- texinfo -*-\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
41 | @deftypefn {Built-in Function} {@var{x} =} inv (@var{A})\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
42 | @deftypefnx {Built-in Function} {[@var{x}, @var{rcond}] =} inv (@var{A})\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
43 | Compute the inverse of the square matrix @var{A}. Return an estimate\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
44 | of the reciprocal condition number if requested, otherwise warn of an\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
45 | ill-conditioned matrix if the reciprocal condition number is small.\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
46 | \n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
47 | In general it is best to avoid calculating the inverse of a matrix\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
48 | directly. For example, it is both faster and more accurate to solve\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
49 | systems of equations (@var{A}*@math{x} = @math{b}) with\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
50 | @code{@var{y} = @var{A} \\ @math{b}}, rather than\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
51 | @code{@var{y} = inv (@var{A}) * @math{b}}.\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
52 | \n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
53 | If called with a sparse matrix, then in general @var{x} will be a full\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
54 | matrix requiring significantly more storage. Avoid forming the inverse\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
55 | of a sparse matrix if possible.\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
56 | @seealso{ldivide, rdivide}\n\octave_value_list Finv (const octave_value_list& args, int nargout) | |||
57 | @end deftypefn")octave_value_list Finv (const octave_value_list& args, int nargout) | |||
58 | { | |||
59 | octave_value_list retval; | |||
60 | ||||
61 | int nargin = args.length (); | |||
62 | ||||
63 | if (nargin != 1) | |||
| ||||
64 | { | |||
65 | print_usage (); | |||
66 | return retval; | |||
67 | } | |||
68 | ||||
69 | octave_value arg = args(0); | |||
70 | ||||
71 | octave_idx_type nr = arg.rows (); | |||
72 | octave_idx_type nc = arg.columns (); | |||
73 | ||||
74 | int arg_is_empty = empty_arg ("inverse", nr, nc); | |||
75 | ||||
76 | if (arg_is_empty < 0) | |||
77 | return retval; | |||
78 | else if (arg_is_empty > 0) | |||
79 | return octave_value (Matrix ()); | |||
80 | ||||
81 | if (nr != nc) | |||
82 | { | |||
83 | gripe_square_matrix_required ("inverse"); | |||
84 | return retval; | |||
85 | } | |||
86 | ||||
87 | octave_value result; | |||
88 | octave_idx_type info; | |||
89 | double rcond = 0.0; | |||
90 | float frcond = 0.0; | |||
91 | bool isfloat = arg.is_single_type (); | |||
92 | ||||
93 | if (arg.is_diag_matrix ()) | |||
94 | { | |||
95 | rcond = 1.0; | |||
96 | frcond = 1.0f; | |||
97 | if (arg.is_complex_type ()) | |||
98 | { | |||
99 | if (isfloat) | |||
100 | { | |||
101 | result = arg.float_complex_diag_matrix_value ().inverse (info); | |||
102 | if (nargout > 1) | |||
103 | frcond = arg.float_complex_diag_matrix_value ().rcond (); | |||
104 | } | |||
105 | else | |||
106 | { | |||
107 | result = arg.complex_diag_matrix_value ().inverse (info); | |||
108 | if (nargout > 1) | |||
109 | rcond = arg.complex_diag_matrix_value ().rcond (); | |||
110 | } | |||
111 | } | |||
112 | else | |||
113 | { | |||
114 | if (isfloat) | |||
115 | { | |||
116 | result = arg.float_diag_matrix_value ().inverse (info); | |||
117 | if (nargout > 1) | |||
118 | frcond = arg.float_diag_matrix_value ().rcond (); | |||
119 | } | |||
120 | else | |||
121 | { | |||
122 | result = arg.diag_matrix_value ().inverse (info); | |||
123 | if (nargout > 1) | |||
124 | rcond = arg.diag_matrix_value ().rcond (); | |||
125 | } | |||
126 | } | |||
127 | } | |||
128 | else if (arg.is_perm_matrix ()) | |||
129 | { | |||
130 | rcond = 1.0; | |||
131 | info = 0; | |||
132 | result = arg.perm_matrix_value ().inverse (); | |||
133 | } | |||
134 | else if (isfloat) | |||
135 | { | |||
136 | if (arg.is_real_type ()) | |||
137 | { | |||
138 | FloatMatrix m = arg.float_matrix_value (); | |||
139 | if (! error_state) | |||
140 | { | |||
141 | MatrixType mattyp = args(0).matrix_type (); | |||
142 | result = m.inverse (mattyp, info, frcond, 1); | |||
143 | args(0).matrix_type (mattyp); | |||
144 | } | |||
145 | } | |||
146 | else if (arg.is_complex_type ()) | |||
147 | { | |||
148 | FloatComplexMatrix m = arg.float_complex_matrix_value (); | |||
149 | if (! error_state) | |||
150 | { | |||
151 | MatrixType mattyp = args(0).matrix_type (); | |||
152 | result = m.inverse (mattyp, info, frcond, 1); | |||
153 | args(0).matrix_type (mattyp); | |||
154 | } | |||
155 | } | |||
156 | } | |||
157 | else | |||
158 | { | |||
159 | if (arg.is_real_type ()) | |||
160 | { | |||
161 | if (arg.is_sparse_type ()) | |||
162 | { | |||
163 | SparseMatrix m = arg.sparse_matrix_value (); | |||
164 | if (! error_state) | |||
165 | { | |||
166 | MatrixType mattyp = args(0).matrix_type (); | |||
167 | result = m.inverse (mattyp, info, rcond, 1); | |||
168 | args(0).matrix_type (mattyp); | |||
169 | } | |||
170 | } | |||
171 | else | |||
172 | { | |||
173 | Matrix m = arg.matrix_value (); | |||
174 | if (! error_state) | |||
175 | { | |||
176 | MatrixType mattyp = args(0).matrix_type (); | |||
177 | result = m.inverse (mattyp, info, rcond, 1); | |||
178 | args(0).matrix_type (mattyp); | |||
179 | } | |||
180 | } | |||
181 | } | |||
182 | else if (arg.is_complex_type ()) | |||
183 | { | |||
184 | if (arg.is_sparse_type ()) | |||
185 | { | |||
186 | SparseComplexMatrix m = arg.sparse_complex_matrix_value (); | |||
187 | if (! error_state) | |||
188 | { | |||
189 | MatrixType mattyp = args(0).matrix_type (); | |||
190 | result = m.inverse (mattyp, info, rcond, 1); | |||
191 | args(0).matrix_type (mattyp); | |||
192 | } | |||
193 | } | |||
194 | else | |||
195 | { | |||
196 | ComplexMatrix m = arg.complex_matrix_value (); | |||
197 | if (! error_state) | |||
198 | { | |||
199 | MatrixType mattyp = args(0).matrix_type (); | |||
200 | result = m.inverse (mattyp, info, rcond, 1); | |||
201 | args(0).matrix_type (mattyp); | |||
202 | } | |||
203 | } | |||
204 | } | |||
205 | else | |||
206 | gripe_wrong_type_arg ("inv", arg); | |||
207 | } | |||
208 | ||||
209 | if (! error_state) | |||
210 | { | |||
211 | if (nargout > 1) | |||
212 | retval(1) = isfloat ? octave_value (frcond) : octave_value (rcond); | |||
213 | ||||
214 | retval(0) = result; | |||
215 | ||||
216 | volatile double xrcond = rcond; | |||
217 | xrcond += 1.0; | |||
218 | if (nargout < 2 && (info == -1 || xrcond == 1.0)) | |||
| ||||
219 | warning ("inverse: matrix singular to machine precision, rcond = %g", | |||
220 | rcond); | |||
221 | } | |||
222 | ||||
223 | return retval; | |||
224 | } | |||
225 | ||||
226 | /* | |||
227 | %!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], sqrt (eps)) | |||
228 | %!assert (inv (single ([1, 2; 3, 4])), single ([-2, 1; 1.5, -0.5]), sqrt (eps ("single"))) | |||
229 | ||||
230 | %!error inv () | |||
231 | %!error inv ([1, 2; 3, 4], 2) | |||
232 | %!error <argument must be a square matrix> inv ([1, 2; 3, 4; 5, 6]) | |||
233 | */ | |||
234 | ||||
235 | // FIXME: this should really be done with an alias, but | |||
236 | // alias_builtin() won't do the right thing if we are actually using | |||
237 | // dynamic linking. | |||
238 | ||||
239 | DEFUN (inverse, args, nargout,octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
240 | "-*- texinfo -*-\n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
241 | @deftypefn {Built-in Function} {@var{x} =} inverse (@var{A})\n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
242 | @deftypefnx {Built-in Function} {[@var{x}, @var{rcond}] =} inverse (@var{A})\n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
243 | Compute the inverse of the square matrix @var{A}.\n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
244 | \n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
245 | This is an alias for @code{inv}.\n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
246 | @seealso{inv}\n\octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
247 | @end deftypefn")octave_value_list Finverse (const octave_value_list& args , int nargout) | |||
248 | { | |||
249 | return Finv (args, nargout); | |||
250 | } |