Mercurial > octave
annotate libinterp/corefcn/hess.cc @ 21267:f5b8c3aca5f8
make better use of templates for Hessenberg decomposition
* liboctave/numeric/hess.h, liboctave/numeric/hess.cc:
New files generated from HESS.h, HESS.cc, CmplxHESS.h,
CmplxHESS.cc, dbleHESS.h, dbleHESS.cc, fCmplxHESS.h,
fCmplxHESS.cc, floatHESS.h, and floatHESS.cc and making them
templates.
* liboctave/numeric/module.mk: Update.
* libinterp/corefcn/hess.cc, mx-defs.h, mx-ext.h: Use new template
classes and header file.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 15 Feb 2016 22:33:45 -0500 |
parents | fcac5dbbf9ed |
children | 40de9f8f23a6 |
rev | line source |
---|---|
2928 | 1 /* |
2 | |
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19040
diff
changeset
|
3 Copyright (C) 1996-2015 John W. Eaton |
2928 | 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 | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
2928 | 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 | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
2928 | 20 |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
21200
fcac5dbbf9ed
maint: Indent #ifdef blocks in libinterp.
Rik <rik@octave.org>
parents:
21129
diff
changeset
|
24 # include <config.h> |
2928 | 25 #endif |
26 | |
21267
f5b8c3aca5f8
make better use of templates for Hessenberg decomposition
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
27 #include "hess.h" |
2928 | 28 |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
29 #include "defun.h" |
2928 | 30 #include "error.h" |
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
21024
diff
changeset
|
31 #include "errwarn.h" |
20940
48b2ad5ee801
maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents:
20939
diff
changeset
|
32 #include "ovl.h" |
2928 | 33 #include "utils.h" |
34 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
35 DEFUN (hess, args, nargout, |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
36 "-*- texinfo -*-\n\ |
20853
1142cf6abc0d
2015 Code Sprint: remove class of function from docstring for all C++ files.
Rik <rik@octave.org>
parents:
20838
diff
changeset
|
37 @deftypefn {} {@var{H} =} hess (@var{A})\n\ |
1142cf6abc0d
2015 Code Sprint: remove class of function from docstring for all C++ files.
Rik <rik@octave.org>
parents:
20838
diff
changeset
|
38 @deftypefnx {} {[@var{P}, @var{H}] =} hess (@var{A})\n\ |
3372 | 39 @cindex Hessenberg decomposition\n\ |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
40 Compute the Hessenberg decomposition of the matrix @var{A}.\n\ |
3372 | 41 \n\ |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
42 The Hessenberg decomposition is\n\ |
3372 | 43 @tex\n\ |
44 $$\n\ | |
45 A = PHP^T\n\ | |
46 $$\n\ | |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
47 where $P$ is a square unitary matrix ($P^TP = I$), and $H$\n\ |
20823
c132fbc748da
doc: Fix incorrect description of Hessenberg decomposition (bug #46622).
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
48 is upper Hessenberg ($H_{i,j} = 0, \\forall i > j+1$).\n\ |
3372 | 49 @end tex\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
50 @ifnottex\n\ |
12584
7ef7e20057fa
Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents:
11553
diff
changeset
|
51 @code{@var{P} * @var{H} * @var{P}' = @var{A}} where @var{P} is a square\n\ |
7ef7e20057fa
Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents:
11553
diff
changeset
|
52 unitary matrix (@code{@var{P}' * @var{P} = I}, using complex-conjugate\n\ |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
53 transposition) and @var{H} is upper Hessenberg\n\ |
20823
c132fbc748da
doc: Fix incorrect description of Hessenberg decomposition (bug #46622).
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
54 (@code{@var{H}(i, j) = 0 forall i > j+1)}.\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
55 @end ifnottex\n\ |
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
56 \n\ |
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
57 The Hessenberg decomposition is usually used as the first step in an\n\ |
19040
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
58 eigenvalue computation, but has other applications as well\n\ |
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
59 (see @nospell{Golub, Nash, and Van Loan},\n\ |
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
60 IEEE Transactions on Automatic Control, 1979).\n\ |
16920
53eaa83e4181
doc: Add seealso links between various factorization forms.
Rik <rik@octave.org>
parents:
15195
diff
changeset
|
61 @seealso{eig, chol, lu, qr, qz, schur, svd}\n\ |
3372 | 62 @end deftypefn") |
2928 | 63 { |
20909
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
64 if (args.length () != 1) |
20801
a542a9bf177e
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20555
diff
changeset
|
65 print_usage (); |
2928 | 66 |
67 octave_value arg = args(0); | |
68 | |
5275 | 69 octave_idx_type nr = arg.rows (); |
70 octave_idx_type nc = arg.columns (); | |
2928 | 71 |
72 int arg_is_empty = empty_arg ("hess", nr, nc); | |
73 | |
74 if (arg_is_empty < 0) | |
20941
a4f5da7c5463
maint: Replace "octave_value_list ()" with "ovl ()".
Rik <rik@octave.org>
parents:
20940
diff
changeset
|
75 return ovl (); |
2928 | 76 else if (arg_is_empty > 0) |
77 return octave_value_list (2, Matrix ()); | |
78 | |
79 if (nr != nc) | |
21110
3d0d84305600
Use err_square_matrix_required more widely.
Rik <rik@octave.org>
parents:
21100
diff
changeset
|
80 err_square_matrix_required ("hess", "A"); |
2928 | 81 |
20939
b17fda023ca6
maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents:
20909
diff
changeset
|
82 octave_value_list retval; |
b17fda023ca6
maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents:
20909
diff
changeset
|
83 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
84 if (arg.is_single_type ()) |
2928 | 85 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
86 if (arg.is_real_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
87 { |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
88 FloatMatrix tmp = arg.float_matrix_value (); |
2928 | 89 |
21267
f5b8c3aca5f8
make better use of templates for Hessenberg decomposition
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
90 hess<FloatMatrix> result (tmp); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
91 |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
92 if (nargout <= 1) |
20884
f1b2a2dbc0e1
2015 Code Sprint: use ovl () in C++ files.
José Luis García Pallero <jgpallero@gmail.com>
parents:
20853
diff
changeset
|
93 retval = ovl (result.hess_matrix ()); |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
94 else |
20909
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
95 retval = ovl (result.unitary_hess_matrix (), |
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
96 result.hess_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
97 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
98 else if (arg.is_complex_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
99 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
100 FloatComplexMatrix ctmp = arg.float_complex_matrix_value (); |
2928 | 101 |
21267
f5b8c3aca5f8
make better use of templates for Hessenberg decomposition
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
102 hess<FloatComplexMatrix> result (ctmp); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
103 |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
104 if (nargout <= 1) |
20884
f1b2a2dbc0e1
2015 Code Sprint: use ovl () in C++ files.
José Luis García Pallero <jgpallero@gmail.com>
parents:
20853
diff
changeset
|
105 retval = ovl (result.hess_matrix ()); |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
106 else |
20909
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
107 retval = ovl (result.unitary_hess_matrix (), |
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
108 result.hess_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
109 } |
2928 | 110 } |
111 else | |
112 { | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
113 if (arg.is_real_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
114 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
115 Matrix tmp = arg.matrix_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
116 |
21267
f5b8c3aca5f8
make better use of templates for Hessenberg decomposition
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
117 hess<Matrix> result (tmp); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
118 |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
119 if (nargout <= 1) |
20884
f1b2a2dbc0e1
2015 Code Sprint: use ovl () in C++ files.
José Luis García Pallero <jgpallero@gmail.com>
parents:
20853
diff
changeset
|
120 retval = ovl (result.hess_matrix ()); |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
121 else |
20909
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
122 retval = ovl (result.unitary_hess_matrix (), |
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
123 result.hess_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
124 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
125 else if (arg.is_complex_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
126 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
127 ComplexMatrix ctmp = arg.complex_matrix_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
128 |
21267
f5b8c3aca5f8
make better use of templates for Hessenberg decomposition
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
129 hess<ComplexMatrix> result (ctmp); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
130 |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
131 if (nargout <= 1) |
20884
f1b2a2dbc0e1
2015 Code Sprint: use ovl () in C++ files.
José Luis García Pallero <jgpallero@gmail.com>
parents:
20853
diff
changeset
|
132 retval = ovl (result.hess_matrix ()); |
20555
f90c8372b7ba
eliminate many more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
19697
diff
changeset
|
133 else |
20909
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
134 retval = ovl (result.unitary_hess_matrix (), |
03e4ddd49396
omit unnecessary nargout checks
John W. Eaton <jwe@octave.org>
parents:
20884
diff
changeset
|
135 result.hess_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
136 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
137 else |
21129
228b65504557
maint: Eliminate useless statements after err_XXX.
Rik <rik@octave.org>
parents:
21110
diff
changeset
|
138 err_wrong_type_arg ("hess", arg); |
2928 | 139 } |
140 | |
141 return retval; | |
142 } | |
143 | |
144 /* | |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
145 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
146 %! a = [1, 2, 3; 5, 4, 6; 8, 7, 9]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
147 %! [p, h] = hess (a); |
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
148 %! assert (p * h * p', a, sqrt (eps)); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
149 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
150 %!test |
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
151 %! a = single ([1, 2, 3; 5, 4, 6; 8, 7, 9]); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
152 %! [p, h] = hess (a); |
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
153 %! assert (p * h * p', a, sqrt (eps ("single"))); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
154 |
14501
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
155 %!error hess () |
60e5cf354d80
Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
156 %!error hess ([1, 2; 3, 4], 2) |
21110
3d0d84305600
Use err_square_matrix_required more widely.
Rik <rik@octave.org>
parents:
21100
diff
changeset
|
157 %!error <must be a square matrix> hess ([1, 2; 3, 4; 5, 6]) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
158 */ |