GNU Octave  3.8.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
eig.cc
Go to the documentation of this file.
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_H
24 #include <config.h>
25 #endif
26 
27 #include "EIG.h"
28 #include "fEIG.h"
29 
30 #include "defun.h"
31 #include "error.h"
32 #include "gripes.h"
33 #include "oct-obj.h"
34 #include "utils.h"
35 
36 DEFUN (eig, args, nargout,
37  "-*- texinfo -*-\n\
38 @deftypefn {Built-in Function} {@var{lambda} =} eig (@var{A})\n\
39 @deftypefnx {Built-in Function} {@var{lambda} =} eig (@var{A}, @var{B})\n\
40 @deftypefnx {Built-in Function} {[@var{V}, @var{lambda}] =} eig (@var{A})\n\
41 @deftypefnx {Built-in Function} {[@var{V}, @var{lambda}] =} eig (@var{A}, @var{B})\n\
42 Compute the eigenvalues (and optionally the eigenvectors) of a matrix\n\
43 or a pair of matrices\n\
44 \n\
45 The algorithm used depends on whether there are one or two input\n\
46 matrices, if they are real or complex and if they are symmetric\n\
47 (Hermitian if complex) or non-symmetric.\n\
48 \n\
49 The eigenvalues returned by @code{eig} are not ordered.\n\
50 @seealso{eigs, svd}\n\
51 @end deftypefn")
52 {
53  octave_value_list retval;
54 
55  int nargin = args.length ();
56 
57  if (nargin > 2 || nargin == 0 || nargout > 2)
58  {
59  print_usage ();
60  return retval;
61  }
62 
63  octave_value arg_a, arg_b;
64 
65  octave_idx_type nr_a = 0, nr_b = 0;
66  octave_idx_type nc_a = 0, nc_b = 0;
67 
68  arg_a = args(0);
69  nr_a = arg_a.rows ();
70  nc_a = arg_a.columns ();
71 
72  int arg_is_empty = empty_arg ("eig", nr_a, nc_a);
73  if (arg_is_empty < 0)
74  return retval;
75  else if (arg_is_empty > 0)
76  return octave_value_list (2, Matrix ());
77 
78  if (!(arg_a.is_single_type () || arg_a.is_double_type ()))
79  {
80  gripe_wrong_type_arg ("eig", arg_a);
81  return retval;
82  }
83 
84  if (nargin == 2)
85  {
86  arg_b = args(1);
87  nr_b = arg_b.rows ();
88  nc_b = arg_b.columns ();
89 
90  arg_is_empty = empty_arg ("eig", nr_b, nc_b);
91  if (arg_is_empty < 0)
92  return retval;
93  else if (arg_is_empty > 0)
94  return octave_value_list (2, Matrix ());
95 
96  if (!(arg_b.is_single_type () || arg_b.is_double_type ()))
97  {
98  gripe_wrong_type_arg ("eig", arg_b);
99  return retval;
100  }
101  }
102 
103  if (nr_a != nc_a)
104  {
106  return retval;
107  }
108 
109  if (nargin == 2 && nr_b != nc_b)
110  {
112  return retval;
113  }
114 
115  Matrix tmp_a, tmp_b;
116  ComplexMatrix ctmp_a, ctmp_b;
117  FloatMatrix ftmp_a, ftmp_b;
118  FloatComplexMatrix fctmp_a, fctmp_b;
119 
120  if (arg_a.is_single_type ())
121  {
122  FloatEIG result;
123 
124  if (nargin == 1)
125  {
126  if (arg_a.is_real_type ())
127  {
128  ftmp_a = arg_a.float_matrix_value ();
129 
130  if (error_state)
131  return retval;
132  else
133  result = FloatEIG (ftmp_a, nargout > 1);
134  }
135  else
136  {
137  fctmp_a = arg_a.float_complex_matrix_value ();
138 
139  if (error_state)
140  return retval;
141  else
142  result = FloatEIG (fctmp_a, nargout > 1);
143  }
144  }
145  else if (nargin == 2)
146  {
147  if (arg_a.is_real_type () && arg_b.is_real_type ())
148  {
149  ftmp_a = arg_a.float_matrix_value ();
150  ftmp_b = arg_b.float_matrix_value ();
151 
152  if (error_state)
153  return retval;
154  else
155  result = FloatEIG (ftmp_a, ftmp_b, nargout > 1);
156  }
157  else
158  {
159  fctmp_a = arg_a.float_complex_matrix_value ();
160  fctmp_b = arg_b.float_complex_matrix_value ();
161 
162  if (error_state)
163  return retval;
164  else
165  result = FloatEIG (fctmp_a, fctmp_b, nargout > 1);
166  }
167  }
168 
169  if (! error_state)
170  {
171  if (nargout == 0 || nargout == 1)
172  {
173  retval(0) = result.eigenvalues ();
174  }
175  else
176  {
177  // Blame it on Matlab.
178 
180 
181  retval(1) = d;
182  retval(0) = result.eigenvectors ();
183  }
184  }
185  }
186  else
187  {
188  EIG result;
189 
190  if (nargin == 1)
191  {
192  if (arg_a.is_real_type ())
193  {
194  tmp_a = arg_a.matrix_value ();
195 
196  if (error_state)
197  return retval;
198  else
199  result = EIG (tmp_a, nargout > 1);
200  }
201  else
202  {
203  ctmp_a = arg_a.complex_matrix_value ();
204 
205  if (error_state)
206  return retval;
207  else
208  result = EIG (ctmp_a, nargout > 1);
209  }
210  }
211  else if (nargin == 2)
212  {
213  if (arg_a.is_real_type () && arg_b.is_real_type ())
214  {
215  tmp_a = arg_a.matrix_value ();
216  tmp_b = arg_b.matrix_value ();
217 
218  if (error_state)
219  return retval;
220  else
221  result = EIG (tmp_a, tmp_b, nargout > 1);
222  }
223  else
224  {
225  ctmp_a = arg_a.complex_matrix_value ();
226  ctmp_b = arg_b.complex_matrix_value ();
227 
228  if (error_state)
229  return retval;
230  else
231  result = EIG (ctmp_a, ctmp_b, nargout > 1);
232  }
233  }
234 
235  if (! error_state)
236  {
237  if (nargout == 0 || nargout == 1)
238  {
239  retval(0) = result.eigenvalues ();
240  }
241  else
242  {
243  // Blame it on Matlab.
244 
245  ComplexDiagMatrix d (result.eigenvalues ());
246 
247  retval(1) = d;
248  retval(0) = result.eigenvectors ();
249  }
250  }
251  }
252 
253  return retval;
254 }
255 
256 /*
257 %!assert (eig ([1, 2; 2, 1]), [-1; 3], sqrt (eps))
258 
259 %!test
260 %! [v, d] = eig ([1, 2; 2, 1]);
261 %! x = 1 / sqrt (2);
262 %! assert (d, [-1, 0; 0, 3], sqrt (eps));
263 %! assert (v, [-x, x; x, x], sqrt (eps));
264 
265 %!assert (eig (single ([1, 2; 2, 1])), single ([-1; 3]), sqrt (eps ("single")))
266 
267 %!test
268 %! [v, d] = eig (single ([1, 2; 2, 1]));
269 %! x = single (1 / sqrt (2));
270 %! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single")));
271 %! assert (v, [-x, x; x, x], sqrt (eps ("single")));
272 
273 %!test
274 %! A = [1, 2; -1, 1]; B = [3, 3; 1, 2];
275 %! [v, d] = eig (A, B);
276 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
277 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
278 
279 %!test
280 %! A = single ([1, 2; -1, 1]); B = single ([3, 3; 1, 2]);
281 %! [v, d] = eig (A, B);
282 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
283 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
284 
285 %!test
286 %! A = [1, 2; 2, 1]; B = [3, -2; -2, 3];
287 %! [v, d] = eig (A, B);
288 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
289 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
290 
291 %!test
292 %! A = single ([1, 2; 2, 1]); B = single ([3, -2; -2, 3]);
293 %! [v, d] = eig (A, B);
294 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
295 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
296 
297 %!test
298 %! A = [1+3i, 2+i; 2-i, 1+3i]; B = [5+9i, 2+i; 2-i, 5+9i];
299 %! [v, d] = eig (A, B);
300 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
301 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
302 
303 %!test
304 %! A = single ([1+3i, 2+i; 2-i, 1+3i]); B = single ([5+9i, 2+i; 2-i, 5+9i]);
305 %! [v, d] = eig (A, B);
306 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
307 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
308 
309 %!test
310 %! A = [1+3i, 2+3i; 3-8i, 8+3i]; B = [8+i, 3+i; 4-9i, 3+i];
311 %! [v, d] = eig (A, B);
312 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
313 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
314 
315 %!test
316 %! A = single ([1+3i, 2+3i; 3-8i, 8+3i]); B = single ([8+i, 3+i; 4-9i, 3+i]);
317 %! [v, d] = eig (A, B);
318 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
319 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
320 
321 %!test
322 %! A = [1, 2; 3, 8]; B = [8, 3; 4, 3];
323 %! [v, d] = eig (A, B);
324 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
325 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
326 
327 %!error eig ()
328 %!error eig ([1, 2; 3, 4], [4, 3; 2, 1], 1)
329 %!error <EIG requires same size matrices> eig ([1, 2; 3, 4], 2)
330 %!error <argument must be a square matrix> eig ([1, 2; 3, 4; 5, 6])
331 %!error <wrong type argument> eig ("abcd")
332 %!error <wrong type argument> eig ([1 2 ; 2 3], "abcd")
333 %!error <wrong type argument> eig (false, [1 2 ; 2 3])
334 */