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