GNU Octave 11.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
hess.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1994-2026 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 "Array-oct.h"
31#include "CMatrix.h"
32#include "dMatrix.h"
33#include "fCMatrix.h"
34#include "fMatrix.h"
35#include "hess.h"
36#include "lapack-proto.h"
37#include "oct-error.h"
38
39OCTAVE_BEGIN_NAMESPACE(octave)
40OCTAVE_BEGIN_NAMESPACE(math)
41
42template <>
43OCTAVE_API octave_idx_type
44hess<Matrix>::init (const Matrix& a)
45{
46 F77_INT a_nr = to_f77_int (a.rows ());
47 F77_INT a_nc = to_f77_int (a.cols ());
48
49 if (a_nr != a_nc)
50 (*current_liboctave_error_handler) ("hess: requires square matrix");
51
52 char job = 'N';
53 char side = 'R';
54
55 F77_INT n = a_nc;
56 F77_INT lwork = 32 * n;
57 F77_INT info;
58 F77_INT ilo;
59 F77_INT ihi;
60
61 m_hess_mat = a;
62 double *h = m_hess_mat.rwdata ();
63
64 Array<double> scale (dim_vector (n, 1));
65 double *pscale = scale.rwdata ();
66
67 F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
68 n, h, n, ilo, ihi, pscale, info
69 F77_CHAR_ARG_LEN (1)));
70
71 Array<double> tau (dim_vector (n-1, 1));
72 double *ptau = tau.rwdata ();
73
74 Array<double> work (dim_vector (lwork, 1));
75 double *pwork = work.rwdata ();
76
77 F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
78 lwork, info));
79
80 m_unitary_hess_mat = m_hess_mat;
81 double *z = m_unitary_hess_mat.rwdata ();
82
83 F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork,
84 lwork, info));
85
86 F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
87 F77_CONST_CHAR_ARG2 (&side, 1),
88 n, ilo, ihi, pscale, n, z,
89 n, info
90 F77_CHAR_ARG_LEN (1)
91 F77_CHAR_ARG_LEN (1)));
92
93 // If someone thinks of a more graceful way of doing
94 // this (or faster for that matter :-)), please let
95 // me know!
96
97 if (n > 2)
98 for (F77_INT j = 0; j < a_nc; j++)
99 for (F77_INT i = j+2; i < a_nr; i++)
100 m_hess_mat.elem (i, j) = 0;
101
102 return info;
103}
104
105template <>
106OCTAVE_API octave_idx_type
107hess<FloatMatrix>::init (const FloatMatrix& a)
108{
109 F77_INT a_nr = to_f77_int (a.rows ());
110 F77_INT a_nc = to_f77_int (a.cols ());
111
112 if (a_nr != a_nc)
113 (*current_liboctave_error_handler) ("hess: requires square matrix");
114
115 char job = 'N';
116 char side = 'R';
117
118 F77_INT n = a_nc;
119 F77_INT lwork = 32 * n;
120 F77_INT info;
121 F77_INT ilo;
122 F77_INT ihi;
123
124 m_hess_mat = a;
125 float *h = m_hess_mat.rwdata ();
126
127 Array<float> scale (dim_vector (n, 1));
128 float *pscale = scale.rwdata ();
129
130 F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
131 n, h, n, ilo, ihi, pscale, info
132 F77_CHAR_ARG_LEN (1)));
133
134 Array<float> tau (dim_vector (n-1, 1));
135 float *ptau = tau.rwdata ();
136
137 Array<float> work (dim_vector (lwork, 1));
138 float *pwork = work.rwdata ();
139
140 F77_XFCN (sgehrd, SGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
141 lwork, info));
142
143 m_unitary_hess_mat = m_hess_mat;
144 float *z = m_unitary_hess_mat.rwdata ();
145
146 F77_XFCN (sorghr, SORGHR, (n, ilo, ihi, z, n, ptau, pwork,
147 lwork, info));
148
149 F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
150 F77_CONST_CHAR_ARG2 (&side, 1),
151 n, ilo, ihi, pscale, n, z,
152 n, info
153 F77_CHAR_ARG_LEN (1)
154 F77_CHAR_ARG_LEN (1)));
155
156 // If someone thinks of a more graceful way of doing
157 // this (or faster for that matter :-)), please let
158 // me know!
159
160 if (n > 2)
161 for (F77_INT j = 0; j < a_nc; j++)
162 for (F77_INT i = j+2; i < a_nr; i++)
163 m_hess_mat.elem (i, j) = 0;
164
165 return info;
166}
167
168template <>
169OCTAVE_API octave_idx_type
170hess<ComplexMatrix>::init (const ComplexMatrix& a)
171{
172 F77_INT a_nr = to_f77_int (a.rows ());
173 F77_INT a_nc = to_f77_int (a.cols ());
174
175 if (a_nr != a_nc)
176 (*current_liboctave_error_handler) ("hess: requires square matrix");
177
178 char job = 'N';
179 char side = 'R';
180
181 F77_INT n = a_nc;
182 F77_INT lwork = 32 * n;
183 F77_INT info;
184 F77_INT ilo;
185 F77_INT ihi;
186
187 m_hess_mat = a;
188 Complex *h = m_hess_mat.rwdata ();
189
190 Array<double> scale (dim_vector (n, 1));
191 double *pscale = scale.rwdata ();
192
193 F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
194 n, F77_DBLE_CMPLX_ARG (h), n, ilo, ihi, pscale, info
195 F77_CHAR_ARG_LEN (1)));
196
197 Array<Complex> tau (dim_vector (n-1, 1));
198 Complex *ptau = tau.rwdata ();
199
200 Array<Complex> work (dim_vector (lwork, 1));
201 Complex *pwork = work.rwdata ();
202
203 F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, F77_DBLE_CMPLX_ARG (h), n,
204 F77_DBLE_CMPLX_ARG (ptau), F77_DBLE_CMPLX_ARG (pwork), lwork, info));
205
206 m_unitary_hess_mat = m_hess_mat;
207 Complex *z = m_unitary_hess_mat.rwdata ();
208
209 F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, F77_DBLE_CMPLX_ARG (z), n,
210 F77_DBLE_CMPLX_ARG (ptau), F77_DBLE_CMPLX_ARG (pwork),
211 lwork, info));
212
213 F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
214 F77_CONST_CHAR_ARG2 (&side, 1),
215 n, ilo, ihi, pscale, n, F77_DBLE_CMPLX_ARG (z), n, info
216 F77_CHAR_ARG_LEN (1)
217 F77_CHAR_ARG_LEN (1)));
218
219 // If someone thinks of a more graceful way of
220 // doing this (or faster for that matter :-)),
221 // please let me know!
222
223 if (n > 2)
224 for (F77_INT j = 0; j < a_nc; j++)
225 for (F77_INT i = j+2; i < a_nr; i++)
226 m_hess_mat.elem (i, j) = 0;
227
228 return info;
229}
230
231template <>
232OCTAVE_API octave_idx_type
233hess<FloatComplexMatrix>::init (const FloatComplexMatrix& a)
234{
235 F77_INT a_nr = to_f77_int (a.rows ());
236 F77_INT a_nc = to_f77_int (a.cols ());
237
238 if (a_nr != a_nc)
239 {
240 (*current_liboctave_error_handler) ("hess: requires square matrix");
241 return -1;
242 }
243
244 char job = 'N';
245 char side = 'R';
246
247 F77_INT n = a_nc;
248 F77_INT lwork = 32 * n;
249 F77_INT info;
250 F77_INT ilo;
251 F77_INT ihi;
252
253 m_hess_mat = a;
254 FloatComplex *h = m_hess_mat.rwdata ();
255
256 Array<float> scale (dim_vector (n, 1));
257 float *pscale = scale.rwdata ();
258
259 F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
260 n, F77_CMPLX_ARG (h), n, ilo, ihi, pscale, info
261 F77_CHAR_ARG_LEN (1)));
262
263 Array<FloatComplex> tau (dim_vector (n-1, 1));
264 FloatComplex *ptau = tau.rwdata ();
265
266 Array<FloatComplex> work (dim_vector (lwork, 1));
267 FloatComplex *pwork = work.rwdata ();
268
269 F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, F77_CMPLX_ARG (h), n,
270 F77_CMPLX_ARG (ptau), F77_CMPLX_ARG (pwork), lwork, info));
271
272 m_unitary_hess_mat = m_hess_mat;
273 FloatComplex *z = m_unitary_hess_mat.rwdata ();
274
275 F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, F77_CMPLX_ARG (z), n,
276 F77_CMPLX_ARG (ptau), F77_CMPLX_ARG (pwork),
277 lwork, info));
278
279 F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
280 F77_CONST_CHAR_ARG2 (&side, 1),
281 n, ilo, ihi, pscale, n, F77_CMPLX_ARG (z), n, info
282 F77_CHAR_ARG_LEN (1)
283 F77_CHAR_ARG_LEN (1)));
284
285 // If someone thinks of a more graceful way of
286 // doing this (or faster for that matter :-)),
287 // please let me know!
288
289 if (n > 2)
290 for (F77_INT j = 0; j < a_nc; j++)
291 for (F77_INT i = j+2; i < a_nr; i++)
292 m_hess_mat.elem (i, j) = 0;
293
294 return info;
295}
296
297OCTAVE_END_NAMESPACE(math)
298OCTAVE_END_NAMESPACE(octave)
Array
N Dimensional Array with copy-on-write semantics.
Definition Array-base.h:130
Array::rows
octave_idx_type rows() const
Definition Array-base.h:485
Array::cols
octave_idx_type cols() const
Definition Array-base.h:495
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1790
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:92
hess
Definition hess.h:38
octave_idx_type
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
scale
void scale(Matrix &m, double x, double y, double z)
Definition graphics.cc:5741
OCTAVE_API
#define OCTAVE_API
Definition main.in.cc:55
Complex
std::complex< double > Complex
Definition oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition oct-cmplx.h:34