GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
hess.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2021 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 
39 namespace octave
40 {
41  namespace math
42  {
43  template <>
44  octave_idx_type
45  hess<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  hess_mat = a;
63  double *h = hess_mat.fortran_vec ();
64 
65  Array<double> scale (dim_vector (n, 1));
66  double *pscale = scale.fortran_vec ();
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.fortran_vec ();
74 
75  Array<double> work (dim_vector (lwork, 1));
76  double *pwork = work.fortran_vec ();
77 
78  F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
79  lwork, info));
80 
81  unitary_hess_mat = hess_mat;
82  double *z = unitary_hess_mat.fortran_vec ();
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  hess_mat.elem (i, j) = 0;
102 
103  return info;
104  }
105 
106  template <>
107  octave_idx_type
108  hess<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  hess_mat = a;
126  float *h = hess_mat.fortran_vec ();
127 
128  Array<float> scale (dim_vector (n, 1));
129  float *pscale = scale.fortran_vec ();
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.fortran_vec ();
137 
138  Array<float> work (dim_vector (lwork, 1));
139  float *pwork = work.fortran_vec ();
140 
141  F77_XFCN (sgehrd, SGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
142  lwork, info));
143 
144  unitary_hess_mat = hess_mat;
145  float *z = unitary_hess_mat.fortran_vec ();
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  hess_mat.elem (i, j) = 0;
165 
166  return info;
167  }
168 
169  template <>
170  octave_idx_type
171  hess<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  hess_mat = a;
189  Complex *h = hess_mat.fortran_vec ();
190 
191  Array<double> scale (dim_vector (n, 1));
192  double *pscale = scale.fortran_vec ();
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.fortran_vec ();
200 
201  Array<Complex> work (dim_vector (lwork, 1));
202  Complex *pwork = work.fortran_vec ();
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  unitary_hess_mat = hess_mat;
208  Complex *z = unitary_hess_mat.fortran_vec ();
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  hess_mat.elem (i, j) = 0;
228 
229  return info;
230  }
231 
232  template <>
233  octave_idx_type
234  hess<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  hess_mat = a;
255  FloatComplex *h = hess_mat.fortran_vec ();
256 
257  Array<float> scale (dim_vector (n, 1));
258  float *pscale = scale.fortran_vec ();
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.fortran_vec ();
266 
267  Array<FloatComplex> work (dim_vector (lwork, 1));
268  FloatComplex *pwork = work.fortran_vec ();
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  unitary_hess_mat = hess_mat;
274  FloatComplex *z = unitary_hess_mat.fortran_vec ();
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  hess_mat.elem (i, j) = 0;
294 
295  return info;
296  }
297  }
298 }
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:423
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:415
Array::fortran_vec
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
Matrix
Definition: dMatrix.h:42
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:95
octave::math::hess::init
octave_idx_type init(const T &a)
octave_idx_type
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:315
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:44
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5846
n
octave_idx_type n
Definition: mx-inlines.cc:753
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34