GNU Octave  4.0.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
inv.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 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 "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,
40  "-*- texinfo -*-\n\
41 @deftypefn {Built-in Function} {@var{x} =} inv (@var{A})\n\
42 @deftypefnx {Built-in Function} {[@var{x}, @var{rcond}] =} inv (@var{A})\n\
43 Compute the inverse of the square matrix @var{A}.\n\
44 \n\
45 Return an estimate of the reciprocal condition number if requested,\n\
46 otherwise warn of an ill-conditioned matrix if the reciprocal condition\n\
47 number is small.\n\
48 \n\
49 In general it is best to avoid calculating the inverse of a matrix directly.\n\
50 For example, it is both faster and more accurate to solve systems of\n\
51 equations (@var{A}*@math{x} = @math{b}) with\n\
52 @code{@var{y} = @var{A} \\ @math{b}}, rather than\n\
53 @code{@var{y} = inv (@var{A}) * @math{b}}.\n\
54 \n\
55 If called with a sparse matrix, then in general @var{x} will be a full\n\
56 matrix requiring significantly more storage. Avoid forming the inverse of a\n\
57 sparse matrix if possible.\n\
58 @seealso{ldivide, rdivide}\n\
59 @end deftypefn")
60 {
61  octave_value_list retval;
62 
63  int nargin = args.length ();
64 
65  if (nargin != 1)
66  {
67  print_usage ();
68  return retval;
69  }
70 
71  octave_value arg = args(0);
72 
73  octave_idx_type nr = arg.rows ();
74  octave_idx_type nc = arg.columns ();
75 
76  int arg_is_empty = empty_arg ("inverse", nr, nc);
77 
78  if (arg_is_empty < 0)
79  return retval;
80  else if (arg_is_empty > 0)
81  return octave_value (Matrix ());
82 
83  if (nr != nc)
84  {
85  gripe_square_matrix_required ("inverse");
86  return retval;
87  }
88 
89  octave_value result;
90  octave_idx_type info;
91  double rcond = 0.0;
92  float frcond = 0.0;
93  bool isfloat = arg.is_single_type ();
94 
95  if (arg.is_diag_matrix ())
96  {
97  rcond = 1.0;
98  frcond = 1.0f;
99  if (arg.is_complex_type ())
100  {
101  if (isfloat)
102  {
103  result = arg.float_complex_diag_matrix_value ().inverse (info);
104  if (nargout > 1)
105  frcond = arg.float_complex_diag_matrix_value ().rcond ();
106  }
107  else
108  {
109  result = arg.complex_diag_matrix_value ().inverse (info);
110  if (nargout > 1)
111  rcond = arg.complex_diag_matrix_value ().rcond ();
112  }
113  }
114  else
115  {
116  if (isfloat)
117  {
118  result = arg.float_diag_matrix_value ().inverse (info);
119  if (nargout > 1)
120  frcond = arg.float_diag_matrix_value ().rcond ();
121  }
122  else
123  {
124  result = arg.diag_matrix_value ().inverse (info);
125  if (nargout > 1)
126  rcond = arg.diag_matrix_value ().rcond ();
127  }
128  }
129  }
130  else if (arg.is_perm_matrix ())
131  {
132  rcond = 1.0;
133  info = 0;
134  result = arg.perm_matrix_value ().inverse ();
135  }
136  else if (isfloat)
137  {
138  if (arg.is_real_type ())
139  {
140  FloatMatrix m = arg.float_matrix_value ();
141  if (! error_state)
142  {
143  MatrixType mattyp = args(0).matrix_type ();
144  result = m.inverse (mattyp, info, frcond, 1);
145  args(0).matrix_type (mattyp);
146  }
147  }
148  else if (arg.is_complex_type ())
149  {
151  if (! error_state)
152  {
153  MatrixType mattyp = args(0).matrix_type ();
154  result = m.inverse (mattyp, info, frcond, 1);
155  args(0).matrix_type (mattyp);
156  }
157  }
158  }
159  else
160  {
161  if (arg.is_real_type ())
162  {
163  if (arg.is_sparse_type ())
164  {
166  if (! error_state)
167  {
168  MatrixType mattyp = args(0).matrix_type ();
169  result = m.inverse (mattyp, info, rcond, 1);
170  args(0).matrix_type (mattyp);
171  }
172  }
173  else
174  {
175  Matrix m = arg.matrix_value ();
176  if (! error_state)
177  {
178  MatrixType mattyp = args(0).matrix_type ();
179  result = m.inverse (mattyp, info, rcond, 1);
180  args(0).matrix_type (mattyp);
181  }
182  }
183  }
184  else if (arg.is_complex_type ())
185  {
186  if (arg.is_sparse_type ())
187  {
189  if (! error_state)
190  {
191  MatrixType mattyp = args(0).matrix_type ();
192  result = m.inverse (mattyp, info, rcond, 1);
193  args(0).matrix_type (mattyp);
194  }
195  }
196  else
197  {
199  if (! error_state)
200  {
201  MatrixType mattyp = args(0).matrix_type ();
202  result = m.inverse (mattyp, info, rcond, 1);
203  args(0).matrix_type (mattyp);
204  }
205  }
206  }
207  else
208  gripe_wrong_type_arg ("inv", arg);
209  }
210 
211  if (! error_state)
212  {
213  if (nargout > 1)
214  retval(1) = isfloat ? octave_value (frcond) : octave_value (rcond);
215 
216  retval(0) = result;
217 
218  bool rcond_plus_one_eq_one = false;
219 
220  if (isfloat)
221  {
222  volatile float xrcond = frcond;
223  rcond_plus_one_eq_one = xrcond + 1.0F == 1.0F;
224  }
225  else
226  {
227  volatile double xrcond = rcond;
228  rcond_plus_one_eq_one = xrcond + 1.0 == 1.0;
229  }
230 
231  if (nargout < 2 && (info == -1 || rcond_plus_one_eq_one))
232  gripe_singular_matrix (isfloat ? frcond : rcond);
233  }
234 
235  return retval;
236 }
237 
238 /*
239 %!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], sqrt (eps))
240 %!assert (inv (single ([1, 2; 3, 4])), single ([-2, 1; 1.5, -0.5]), sqrt (eps ("single")))
241 
242 %!error inv ()
243 %!error inv ([1, 2; 3, 4], 2)
244 %!error <argument must be a square matrix> inv ([1, 2; 3, 4; 5, 6])
245 
246 %!test
247 %! [xinv, rcond] = inv (single ([1,2;3,4]));
248 %! assert (isa (xinv, 'single'));
249 %! assert (isa (rcond, 'single'));
250 
251 %!test
252 %! [xinv, rcond] = inv ([1,2;3,4]);
253 %! assert (isa (xinv, 'double'));
254 %! assert (isa (rcond, 'double'));
255 */
256 
257 // FIXME: this should really be done with an alias, but
258 // alias_builtin() won't do the right thing if we are actually using
259 // dynamic linking.
260 
261 DEFUN (inverse, args, nargout,
262  "-*- texinfo -*-\n\
263 @deftypefn {Built-in Function} {@var{x} =} inverse (@var{A})\n\
264 @deftypefnx {Built-in Function} {[@var{x}, @var{rcond}] =} inverse (@var{A})\n\
265 Compute the inverse of the square matrix @var{A}.\n\
266 \n\
267 This is an alias for @code{inv}.\n\
268 @seealso{inv}\n\
269 @end deftypefn")
270 {
271  return Finv (args, nargout);
272 }
FloatComplexDiagMatrix float_complex_diag_matrix_value(bool force=false) const
Definition: ov.h:840
bool is_real_type(void) const
Definition: ov.h:651
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
OCTAVE_EXPORT octave_value_list Finv(const octave_value_list &args, int nargout)
Definition: inv.cc:59
octave_idx_type rows(void) const
Definition: ov.h:473
ComplexMatrix inverse(void) const
Definition: CMatrix.cc:1002
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
float rcond(void) const
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
SparseComplexMatrix inverse(void) const
Definition: CSparse.cc:733
octave_idx_type length(void) const
Definition: oct-obj.h:89
SparseMatrix inverse(void) const
Definition: dSparse.cc:828
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
bool is_perm_matrix(void) const
Definition: ov.h:559
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:246
DiagMatrix inverse(void) const
Definition: dDiagMatrix.cc:259
FloatComplexMatrix inverse(void) const
Definition: fCMatrix.cc:1008
void gripe_square_matrix_required(const char *name)
Definition: gripes.cc:81
ComplexDiagMatrix complex_diag_matrix_value(bool force=false) const
Definition: ov.h:836
octave_idx_type columns(void) const
Definition: ov.h:475
PermMatrix inverse(void) const
Definition: PermMatrix.cc:132
bool is_sparse_type(void) const
Definition: ov.h:666
FloatDiagMatrix float_diag_matrix_value(bool force=false) const
Definition: ov.h:833
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
float rcond(void) const
Definition: fDiagMatrix.cc:366
void gripe_singular_matrix(double rcond)
Definition: dMatrix.h:35
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
double arg(double x)
Definition: lo-mappers.h:37
double rcond(void) const
Definition: dDiagMatrix.cc:366
double rcond(void) const
Definition: CDiagMatrix.cc:551
DiagMatrix diag_matrix_value(bool force=false) const
Definition: ov.h:830
FloatDiagMatrix inverse(void) const
Definition: fDiagMatrix.cc:259
FloatMatrix inverse(void) const
Definition: fMatrix.cc:658
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
PermMatrix perm_matrix_value(void) const
Definition: ov.h:843
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:820
ComplexDiagMatrix inverse(octave_idx_type &info) const
Definition: CDiagMatrix.cc:358
bool is_single_type(void) const
Definition: ov.h:611
FloatComplexDiagMatrix inverse(octave_idx_type &info) const
bool is_diag_matrix(void) const
Definition: ov.h:556
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
Matrix inverse(void) const
Definition: dMatrix.cc:651