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
hess.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 "CmplxHESS.h"
28 #include "dbleHESS.h"
29 #include "fCmplxHESS.h"
30 #include "floatHESS.h"
31 
32 #include "defun.h"
33 #include "error.h"
34 #include "gripes.h"
35 #include "oct-obj.h"
36 #include "utils.h"
37 
38 DEFUN (hess, args, nargout,
39  "-*- texinfo -*-\n\
40 @deftypefn {Built-in Function} {@var{H} =} hess (@var{A})\n\
41 @deftypefnx {Built-in Function} {[@var{P}, @var{H}] =} hess (@var{A})\n\
42 @cindex Hessenberg decomposition\n\
43 Compute the Hessenberg decomposition of the matrix @var{A}.\n\
44 \n\
45 The Hessenberg decomposition is\n\
46 @tex\n\
47 $$\n\
48 A = PHP^T\n\
49 $$\n\
50 where $P$ is a square unitary matrix ($P^TP = I$), and $H$\n\
51 is upper Hessenberg ($H_{i,j} = 0, \\forall i \\ge j+1$).\n\
52 @end tex\n\
53 @ifnottex\n\
54 @code{@var{P} * @var{H} * @var{P}' = @var{A}} where @var{P} is a square\n\
55 unitary matrix (@code{@var{P}' * @var{P} = I}, using complex-conjugate\n\
56 transposition) and @var{H} is upper Hessenberg\n\
57 (@code{@var{H}(i, j) = 0 forall i >= j+1)}.\n\
58 @end ifnottex\n\
59 \n\
60 The Hessenberg decomposition is usually used as the first step in an\n\
61 eigenvalue computation, but has other applications as well\n\
62 (see @nospell{Golub, Nash, and Van Loan},\n\
63 IEEE Transactions on Automatic Control, 1979).\n\
64 @seealso{eig, chol, lu, qr, qz, schur, svd}\n\
65 @end deftypefn")
66 {
67  octave_value_list retval;
68 
69  int nargin = args.length ();
70 
71  if (nargin != 1 || nargout > 2)
72  {
73  print_usage ();
74  return retval;
75  }
76 
77  octave_value arg = args(0);
78 
79  octave_idx_type nr = arg.rows ();
80  octave_idx_type nc = arg.columns ();
81 
82  int arg_is_empty = empty_arg ("hess", nr, nc);
83 
84  if (arg_is_empty < 0)
85  return retval;
86  else if (arg_is_empty > 0)
87  return octave_value_list (2, Matrix ());
88 
89  if (nr != nc)
90  {
92  return retval;
93  }
94 
95  if (arg.is_single_type ())
96  {
97  if (arg.is_real_type ())
98  {
99  FloatMatrix tmp = arg.float_matrix_value ();
100 
101  if (! error_state)
102  {
103  FloatHESS result (tmp);
104 
105  if (nargout <= 1)
106  retval(0) = result.hess_matrix ();
107  else
108  {
109  retval(1) = result.hess_matrix ();
110  retval(0) = result.unitary_hess_matrix ();
111  }
112  }
113  }
114  else if (arg.is_complex_type ())
115  {
117 
118  if (! error_state)
119  {
120  FloatComplexHESS result (ctmp);
121 
122  if (nargout <= 1)
123  retval(0) = result.hess_matrix ();
124  else
125  {
126  retval(1) = result.hess_matrix ();
127  retval(0) = result.unitary_hess_matrix ();
128  }
129  }
130  }
131  }
132  else
133  {
134  if (arg.is_real_type ())
135  {
136  Matrix tmp = arg.matrix_value ();
137 
138  if (! error_state)
139  {
140  HESS result (tmp);
141 
142  if (nargout <= 1)
143  retval(0) = result.hess_matrix ();
144  else
145  {
146  retval(1) = result.hess_matrix ();
147  retval(0) = result.unitary_hess_matrix ();
148  }
149  }
150  }
151  else if (arg.is_complex_type ())
152  {
153  ComplexMatrix ctmp = arg.complex_matrix_value ();
154 
155  if (! error_state)
156  {
157  ComplexHESS result (ctmp);
158 
159  if (nargout <= 1)
160  retval(0) = result.hess_matrix ();
161  else
162  {
163  retval(1) = result.hess_matrix ();
164  retval(0) = result.unitary_hess_matrix ();
165  }
166  }
167  }
168  else
169  {
170  gripe_wrong_type_arg ("hess", arg);
171  }
172  }
173 
174  return retval;
175 }
176 
177 /*
178 %!test
179 %! a = [1, 2, 3; 5, 4, 6; 8, 7, 9];
180 %! [p, h] = hess (a);
181 %! assert (p * h * p', a, sqrt (eps));
182 
183 %!test
184 %! a = single ([1, 2, 3; 5, 4, 6; 8, 7, 9]);
185 %! [p, h] = hess (a);
186 %! assert (p * h * p', a, sqrt (eps ("single")));
187 
188 %!error hess ()
189 %!error hess ([1, 2; 3, 4], 2)
190 %!error <argument must be a square matrix> hess ([1, 2; 3, 4; 5, 6])
191 */
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_idx_type rows(void) const
Definition: ov.h:473
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:246
void gripe_square_matrix_required(const char *name)
Definition: gripes.cc:81
FloatMatrix unitary_hess_matrix(void) const
Definition: floatHESS.h:67
octave_idx_type columns(void) const
Definition: ov.h:475
FloatMatrix hess_matrix(void) const
Definition: floatHESS.h:65
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
Definition: dMatrix.h:35
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
Matrix unitary_hess_matrix(void) const
Definition: dbleHESS.h:63
double arg(double x)
Definition: lo-mappers.h:37
FloatComplexMatrix hess_matrix(void) const
Definition: fCmplxHESS.h:65
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
ComplexMatrix hess_matrix(void) const
Definition: CmplxHESS.h:65
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
Definition: dbleHESS.h:30
FloatComplexMatrix unitary_hess_matrix(void) const
Definition: fCmplxHESS.h:67
ComplexMatrix unitary_hess_matrix(void) const
Definition: CmplxHESS.h:67
bool is_single_type(void) const
Definition: ov.h:611
Matrix hess_matrix(void) const
Definition: dbleHESS.h:61