1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1996-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 #if defined (HAVE_CONFIG_H) 27 # include "config.h" 28 #endif 29 30 #include "defun.h" 31 #include "error.h" 32 #include "errwarn.h" 33 #include "ovl.h" 34 #include "ops.h" 35 #include "ov-re-diag.h" 36 #include "ov-cx-diag.h" 37 #include "ov-flt-re-diag.h" 38 #include "ov-flt-cx-diag.h" 39 #include "ov-perm.h" 40 41 DEFUN (inv, args, nargout, 42 doc: /* -*- texinfo -*- 43 @deftypefn {} {@var{x} =} inv (@var{A}) 44 @deftypefnx {} {[@var{x}, @var{rcond}] =} inv (@var{A}) 45 @deftypefnx {} {[@dots{}] =} inverse (@dots{}) 46 Compute the inverse of the square matrix @var{A}. 47 48 Return an estimate of the reciprocal condition number if requested, 49 otherwise warn of an ill-conditioned matrix if the reciprocal condition 50 number is small. 51 52 In general it is best to avoid calculating the inverse of a matrix directly. 53 For example, it is both faster and more accurate to solve systems of 54 equations (@var{A}*@math{x} = @math{b}) with 55 @code{@var{y} = @var{A} \ @math{b}}, rather than 56 @code{@var{y} = inv (@var{A}) * @math{b}}. 57 58 If called with a sparse matrix, then in general @var{x} will be a full 59 matrix requiring significantly more storage. Avoid forming the inverse of a 60 sparse matrix if possible. 61 62 @code{inverse} is an alias and may be used identically in place of @code{inv}. 63 @seealso{ldivide, rdivide, pinv} 64 @end deftypefn */) 65 { 66 if (args.length () != 1) 67 print_usage (); 68 69 octave_value arg = args(0); 70 71 if (arg.isempty ()) 72 return ovl (Matrix ()); 73 74 if (arg.rows () != arg.columns ()) 75 err_square_matrix_required ("inverse", "A"); 76 77 octave_value result; 78 octave_idx_type info; 79 double rcond = 0.0; 80 float frcond = 0.0; 81 bool isfloat = arg.is_single_type (); 82 83 if (arg.is_diag_matrix ()) 84 { 85 rcond = 1.0; 86 frcond = 1.0f; 87 if (arg.iscomplex ()) 88 { 89 if (isfloat) 90 { 91 result = arg.float_complex_diag_matrix_value ().inverse (info); 92 if (info == -1) 93 frcond = 0.0f; 94 else if (nargout > 1) 95 frcond = arg.float_complex_diag_matrix_value ().rcond (); 96 } 97 else 98 { 99 result = arg.complex_diag_matrix_value ().inverse (info); 100 if (info == -1) 101 rcond = 0.0; 102 else if (nargout > 1) 103 rcond = arg.complex_diag_matrix_value ().rcond (); 104 } 105 } 106 else 107 { 108 if (isfloat) 109 { 110 result = arg.float_diag_matrix_value ().inverse (info); 111 if (info == -1) 112 frcond = 0.0f; 113 else if (nargout > 1) 114 frcond = arg.float_diag_matrix_value ().rcond (); 115 } 116 else 117 { 118 result = arg.diag_matrix_value ().inverse (info); 119 if (info == -1) 120 rcond = 0.0; 121 else if (nargout > 1) 122 rcond = arg.diag_matrix_value ().rcond (); 123 } 124 } 125 } 126 else if (arg.is_perm_matrix ()) 127 { 128 rcond = 1.0; 129 info = 0; 130 result = arg.perm_matrix_value ().inverse (); 131 } 132 else if (isfloat) 133 { 134 if (arg.isreal ()) 135 { 136 FloatMatrix m = arg.float_matrix_value (); 137 138 MatrixType mattyp = args(0).matrix_type (); 139 result = m.inverse (mattyp, info, frcond, 1); 140 args(0).matrix_type (mattyp); 141 } 142 else if (arg.iscomplex ()) 143 { 144 FloatComplexMatrix m = arg.float_complex_matrix_value (); 145 146 MatrixType mattyp = args(0).matrix_type (); 147 result = m.inverse (mattyp, info, frcond, 1); 148 args(0).matrix_type (mattyp); 149 } 150 } 151 else 152 { 153 if (arg.isreal ()) 154 { 155 if (arg.issparse ()) 156 { 157 SparseMatrix m = arg.sparse_matrix_value (); 158 159 MatrixType mattyp = args(0).matrix_type (); 160 result = m.inverse (mattyp, info, rcond, 1); 161 args(0).matrix_type (mattyp); 162 } 163 else 164 { 165 Matrix m = arg.matrix_value (); 166 167 MatrixType mattyp = args(0).matrix_type (); 168 result = m.inverse (mattyp, info, rcond, 1); 169 args(0).matrix_type (mattyp); 170 } 171 } 172 else if (arg.iscomplex ()) 173 { 174 if (arg.issparse ()) 175 { 176 SparseComplexMatrix m = arg.sparse_complex_matrix_value (); 177 178 MatrixType mattyp = args(0).matrix_type (); 179 result = m.inverse (mattyp, info, rcond, 1); 180 args(0).matrix_type (mattyp); 181 } 182 else 183 { 184 ComplexMatrix m = arg.complex_matrix_value (); 185 186 MatrixType mattyp = args(0).matrix_type (); 187 result = m.inverse (mattyp, info, rcond, 1); 188 args(0).matrix_type (mattyp); 189 } 190 } 191 else 192 err_wrong_type_arg ("inv", arg); 193 } 194 195 octave_value_list retval (nargout > 1 ? 2 : 1); 196 197 retval(0) = result; 198 if (nargout > 1) 199 retval(1) = (isfloat ? octave_value (frcond) : octave_value (rcond)); 200 201 bool rcond_plus_one_eq_one = false; 202 203 if (isfloat) 204 { 205 volatile float xrcond = frcond; 206 rcond_plus_one_eq_one = xrcond + 1.0f == 1.0f; 207 } 208 else 209 { 210 volatile double xrcond = rcond; 211 rcond_plus_one_eq_one = xrcond + 1.0 == 1.0; 212 } 213 214 if (nargout < 2 && (info == -1 || rcond_plus_one_eq_one)) 215 octave::warn_singular_matrix (isfloat ? frcond : rcond); 216 217 return retval; 218 } 219 220 /* 221 %!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], sqrt (eps)) 222 %!assert (inv (single ([1, 2; 3, 4])), single ([-2, 1; 1.5, -0.5]), sqrt (eps ("single"))) 223 224 ## Test special inputs 225 %!assert (inv (zeros (2,0)), []) 226 %!warning <matrix singular> assert (inv (Inf), 0) 227 %!warning <matrix singular> assert (inv (-Inf), -0) 228 %!warning <matrix singular> assert (inv (single (Inf)), single (0)) 229 %!warning <matrix singular> assert (inv (complex (1, Inf)), 0) 230 %!warning <matrix singular> assert (inv (single (complex (1,Inf))), single (0)) 231 232 %!test 233 %! [xinv, rcond] = inv (single ([1,2;3,4])); 234 %! assert (isa (xinv, "single")); 235 %! assert (isa (rcond, "single")); 236 237 %!test 238 %! [xinv, rcond] = inv ([1,2;3,4]); 239 %! assert (isa (xinv, "double")); 240 %! assert (isa (rcond, "double")); 241 242 %!testif HAVE_UMFPACK <*56232> 243 %! fail ("A = inv (sparse ([1, 2;0 ,0]))", "warning", "matrix singular"); 244 %! assert (A, sparse ([Inf, Inf; 0, 0])); 245 246 %!testif HAVE_UMFPACK <*56232> 247 %! fail ("A = inv (sparse ([1i, 2;0 ,0]))", "warning", "matrix singular"); 248 %! assert (A, sparse ([Inf, Inf; 0, 0])); 249 250 %!test 251 %! fail ("A = inv (diag ([1, 0, 1]))", "warning", "matrix singular"); 252 %! assert (A, diag ([Inf, Inf, Inf])); 253 254 %!error <inverse of the null matrix not defined> inv (diag ([0, 0])) 255 %!error <inverse of the null matrix not defined> inv (diag (complex ([0, 0]))) 256 257 %!testif HAVE_UMFPACK <*56232> 258 %! fail ("A = inv (sparse ([1, 0, 0; 0, 0, 0; 0, 0, 1]))", "warning", "matrix singular"); 259 %! assert (A, sparse ([Inf, 0, 0; 0, 0, 0; 0, 0, Inf])); 260 261 %!error inv () 262 %!error inv ([1, 2; 3, 4], 2) 263 %!error <must be a square matrix> inv ([1, 2; 3, 4; 5, 6]) 264 %!error <inverse of the null matrix not defined> inv (sparse (2, 2, 0)) 265 */ 266 267 DEFALIAS (inverse, inv); 268