GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
inv.cc
Go to the documentation of this file.
1 /*
2
3 Copyright (C) 1996-2017 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
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
20
21 */
22
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26
27 #include "defun.h"
28 #include "error.h"
29 #include "errwarn.h"
30 #include "ovl.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
38 DEFUN (inv, args, nargout,
39  doc: /* -*- texinfo -*-
40 @deftypefn {} {@var{x} =} inv (@var{A})
41 @deftypefnx {} {[@var{x}, @var{rcond}] =} inv (@var{A})
42 Compute the inverse of the square matrix @var{A}.
43
44 Return an estimate of the reciprocal condition number if requested,
45 otherwise warn of an ill-conditioned matrix if the reciprocal condition
46 number is small.
47
48 In general it is best to avoid calculating the inverse of a matrix directly.
49 For example, it is both faster and more accurate to solve systems of
50 equations (@var{A}*@math{x} = @math{b}) with
51 @code{@var{y} = @var{A} \ @math{b}}, rather than
52 @code{@var{y} = inv (@var{A}) * @math{b}}.
53
54 If called with a sparse matrix, then in general @var{x} will be a full
55 matrix requiring significantly more storage. Avoid forming the inverse of a
56 sparse matrix if possible.
57 @seealso{ldivide, rdivide}
58 @end deftypefn */)
59 {
60  if (args.length () != 1)
61  print_usage ();
62
63  octave_value arg = args(0);
64
65  if (arg.is_empty ())
66  return ovl (Matrix ());
67
68  if (arg.rows () != arg.columns ())
69  err_square_matrix_required ("inverse", "A");
70
72  octave_idx_type info;
73  double rcond = 0.0;
74  float frcond = 0.0;
75  bool isfloat = arg.is_single_type ();
76
77  if (arg.is_diag_matrix ())
78  {
79  rcond = 1.0;
80  frcond = 1.0f;
81  if (arg.is_complex_type ())
82  {
83  if (isfloat)
84  {
85  result = arg.float_complex_diag_matrix_value ().inverse (info);
86  if (nargout > 1)
87  frcond = arg.float_complex_diag_matrix_value ().rcond ();
88  }
89  else
90  {
91  result = arg.complex_diag_matrix_value ().inverse (info);
92  if (nargout > 1)
93  rcond = arg.complex_diag_matrix_value ().rcond ();
94  }
95  }
96  else
97  {
98  if (isfloat)
99  {
100  result = arg.float_diag_matrix_value ().inverse (info);
101  if (nargout > 1)
102  frcond = arg.float_diag_matrix_value ().rcond ();
103  }
104  else
105  {
106  result = arg.diag_matrix_value ().inverse (info);
107  if (nargout > 1)
108  rcond = arg.diag_matrix_value ().rcond ();
109  }
110  }
111  }
112  else if (arg.is_perm_matrix ())
113  {
114  rcond = 1.0;
115  info = 0;
116  result = arg.perm_matrix_value ().inverse ();
117  }
118  else if (isfloat)
119  {
120  if (arg.is_real_type ())
121  {
123
124  MatrixType mattyp = args(0).matrix_type ();
125  result = m.inverse (mattyp, info, frcond, 1);
126  args(0).matrix_type (mattyp);
127  }
128  else if (arg.is_complex_type ())
129  {
131
132  MatrixType mattyp = args(0).matrix_type ();
133  result = m.inverse (mattyp, info, frcond, 1);
134  args(0).matrix_type (mattyp);
135  }
136  }
137  else
138  {
139  if (arg.is_real_type ())
140  {
141  if (arg.is_sparse_type ())
142  {
144
145  MatrixType mattyp = args(0).matrix_type ();
146  result = m.inverse (mattyp, info, rcond, 1);
147  args(0).matrix_type (mattyp);
148  }
149  else
150  {
151  Matrix m = arg.matrix_value ();
152
153  MatrixType mattyp = args(0).matrix_type ();
154  result = m.inverse (mattyp, info, rcond, 1);
155  args(0).matrix_type (mattyp);
156  }
157  }
158  else if (arg.is_complex_type ())
159  {
160  if (arg.is_sparse_type ())
161  {
163
164  MatrixType mattyp = args(0).matrix_type ();
165  result = m.inverse (mattyp, info, rcond, 1);
166  args(0).matrix_type (mattyp);
167  }
168  else
169  {
171
172  MatrixType mattyp = args(0).matrix_type ();
173  result = m.inverse (mattyp, info, rcond, 1);
174  args(0).matrix_type (mattyp);
175  }
176  }
177  else
178  err_wrong_type_arg ("inv", arg);
179  }
180
181  octave_value_list retval (nargout > 1 ? 2 : 1);
182
183  retval(0) = result;
184  if (nargout > 1)
185  retval(1) = isfloat ? octave_value (frcond) : octave_value (rcond);
186
187  bool rcond_plus_one_eq_one = false;
188
189  if (isfloat)
190  {
191  volatile float xrcond = frcond;
192  rcond_plus_one_eq_one = xrcond + 1.0F == 1.0F;
193  }
194  else
195  {
196  volatile double xrcond = rcond;
197  rcond_plus_one_eq_one = xrcond + 1.0 == 1.0;
198  }
199
200  if (nargout < 2 && (info == -1 || rcond_plus_one_eq_one))
201  octave::warn_singular_matrix (isfloat ? frcond : rcond);
202
203  return retval;
204 }
205
206 /*
207 %!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], sqrt (eps))
208 %!assert (inv (single ([1, 2; 3, 4])), single ([-2, 1; 1.5, -0.5]), sqrt (eps ("single")))
209
210 %!error inv ()
211 %!error inv ([1, 2; 3, 4], 2)
212 %!error <must be a square matrix> inv ([1, 2; 3, 4; 5, 6])
213
214 %!test
215 %! [xinv, rcond] = inv (single ([1,2;3,4]));
216 %! assert (isa (xinv, 'single'));
217 %! assert (isa (rcond, 'single'));
218
219 %!test
220 %! [xinv, rcond] = inv ([1,2;3,4]);
221 %! assert (isa (xinv, 'double'));
222 %! assert (isa (rcond, 'double'));
223 */
224
225 // FIXME: this should really be done with an alias, but
226 // alias_builtin() won't do the right thing if we are actually using
228
229 DEFUN (inverse, args, nargout,
230  doc: /* -*- texinfo -*-
231 @deftypefn {} {@var{x} =} inverse (@var{A})
232 @deftypefnx {} {[@var{x}, @var{rcond}] =} inverse (@var{A})
233 Compute the inverse of the square matrix @var{A}.
234
235 This is an alias for @code{inv}.
236 @seealso{inv}
237 @end deftypefn */)
238 {
239  return Finv (args, nargout);
240 }
FloatComplexDiagMatrix float_complex_diag_matrix_value(bool force=false) const
Definition: ov.h:854
bool is_real_type(void) const
Definition: ov.h:667
octave_idx_type rows(void) const
Definition: ov.h:489
ComplexMatrix inverse(void) const
Definition: CMatrix.cc:717
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:809
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
float rcond(void) const
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
SparseComplexMatrix inverse(void) const
Definition: CSparse.cc:661
SparseMatrix inverse(void) const
Definition: dSparse.cc:741
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
bool is_perm_matrix(void) const
Definition: ov.h:575
DiagMatrix inverse(void) const
Definition: dDiagMatrix.cc:226
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:112
FloatComplexMatrix inverse(void) const
Definition: fCMatrix.cc:720
ComplexDiagMatrix complex_diag_matrix_value(bool force=false) const
Definition: ov.h:850
octave_value arg
Definition: pr-output.cc:3440
JNIEnv void * args
Definition: ov-java.cc:67
octave_idx_type columns(void) const
Definition: ov.h:491
PermMatrix inverse(void) const
Definition: PermMatrix.cc:123
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:935
bool is_sparse_type(void) const
Definition: ov.h:682
OCTAVE_EXPORT octave_value_list isfloat
Definition: data.cc:3327
FloatDiagMatrix float_diag_matrix_value(bool force=false) const
Definition: ov.h:847
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
bool is_complex_type(void) const
Definition: ov.h:670
OCTAVE_EXPORT octave_value_list Finv(const octave_value_list &args, int nargout) ar
Definition: inv.cc:58
octave_value retval
Definition: data.cc:6294
float rcond(void) const
Definition: fDiagMatrix.cc:322
Definition: dMatrix.h:37
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:838
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:787
void err_wrong_type_arg(const char *name, const char *s)
Definition: errwarn.cc:156
With real return the complex result
Definition: data.cc:3375
double rcond(void) const
Definition: dDiagMatrix.cc:322
double rcond(void) const
Definition: CDiagMatrix.cc:486
DiagMatrix diag_matrix_value(bool force=false) const
Definition: ov.h:844
FloatDiagMatrix inverse(void) const
Definition: fDiagMatrix.cc:226
bool is_empty(void) const
Definition: ov.h:542
FloatMatrix inverse(void) const
Definition: fMatrix.cc:435
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:805
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:790
PermMatrix perm_matrix_value(void) const
Definition: ov.h:857
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:834
ComplexDiagMatrix inverse(octave_idx_type &info) const
Definition: CDiagMatrix.cc:310
bool is_single_type(void) const
Definition: ov.h:627
FloatComplexDiagMatrix inverse(octave_idx_type &info) const
void warn_singular_matrix(double rcond)
bool is_diag_matrix(void) const
Definition: ov.h:572
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:429