GNU Octave  4.2.1
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) 1994-2017 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 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "CMatrix.h"
28 #include "dMatrix.h"
29 #include "fCMatrix.h"
30 #include "fMatrix.h"
31 #include "hess.h"
32 #include "lo-error.h"
33 #include "lo-lapack-proto.h"
34 
35 namespace octave
36 {
37  namespace math
38  {
39  template <>
40  octave_idx_type
41  hess<Matrix>::init (const Matrix& a)
42  {
43  octave_idx_type a_nr = a.rows ();
44  octave_idx_type a_nc = a.cols ();
45 
46  if (a_nr != a_nc)
47  (*current_liboctave_error_handler) ("hess: requires square matrix");
48 
49  char job = 'N';
50  char side = 'R';
51 
52  octave_idx_type n = a_nc;
53  octave_idx_type lwork = 32 * n;
54  octave_idx_type info;
55  octave_idx_type ilo;
56  octave_idx_type ihi;
57 
58  hess_mat = a;
59  double *h = hess_mat.fortran_vec ();
60 
61  Array<double> scale (dim_vector (n, 1));
62  double *pscale = scale.fortran_vec ();
63 
64  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
65  n, h, n, ilo, ihi, pscale, info
66  F77_CHAR_ARG_LEN (1)));
67 
68  Array<double> tau (dim_vector (n-1, 1));
69  double *ptau = tau.fortran_vec ();
70 
71  Array<double> work (dim_vector (lwork, 1));
72  double *pwork = work.fortran_vec ();
73 
74  F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
75  lwork, info));
76 
77  unitary_hess_mat = hess_mat;
78  double *z = unitary_hess_mat.fortran_vec ();
79 
80  F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork,
81  lwork, info));
82 
83  F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
84  F77_CONST_CHAR_ARG2 (&side, 1),
85  n, ilo, ihi, pscale, n, z,
86  n, info
87  F77_CHAR_ARG_LEN (1)
88  F77_CHAR_ARG_LEN (1)));
89 
90  // If someone thinks of a more graceful way of doing
91  // this (or faster for that matter :-)), please let
92  // me know!
93 
94  if (n > 2)
95  for (octave_idx_type j = 0; j < a_nc; j++)
96  for (octave_idx_type i = j+2; i < a_nr; i++)
97  hess_mat.elem (i, j) = 0;
98 
99  return info;
100  }
101 
102  template <>
103  octave_idx_type
104  hess<FloatMatrix>::init (const FloatMatrix& a)
105  {
106  octave_idx_type a_nr = a.rows ();
107  octave_idx_type a_nc = a.cols ();
108 
109  if (a_nr != a_nc)
110  (*current_liboctave_error_handler) ("hess: requires square matrix");
111 
112  char job = 'N';
113  char side = 'R';
114 
115  octave_idx_type n = a_nc;
116  octave_idx_type lwork = 32 * n;
117  octave_idx_type info;
118  octave_idx_type ilo;
119  octave_idx_type ihi;
120 
121  hess_mat = a;
122  float *h = hess_mat.fortran_vec ();
123 
124  Array<float> scale (dim_vector (n, 1));
125  float *pscale = scale.fortran_vec ();
126 
127  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
128  n, h, n, ilo, ihi, pscale, info
129  F77_CHAR_ARG_LEN (1)));
130 
131  Array<float> tau (dim_vector (n-1, 1));
132  float *ptau = tau.fortran_vec ();
133 
134  Array<float> work (dim_vector (lwork, 1));
135  float *pwork = work.fortran_vec ();
136 
137  F77_XFCN (sgehrd, SGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
138  lwork, info));
139 
140  unitary_hess_mat = hess_mat;
141  float *z = unitary_hess_mat.fortran_vec ();
142 
143  F77_XFCN (sorghr, SORGHR, (n, ilo, ihi, z, n, ptau, pwork,
144  lwork, info));
145 
146  F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
147  F77_CONST_CHAR_ARG2 (&side, 1),
148  n, ilo, ihi, pscale, n, z,
149  n, info
150  F77_CHAR_ARG_LEN (1)
151  F77_CHAR_ARG_LEN (1)));
152 
153  // If someone thinks of a more graceful way of doing
154  // this (or faster for that matter :-)), please let
155  // me know!
156 
157  if (n > 2)
158  for (octave_idx_type j = 0; j < a_nc; j++)
159  for (octave_idx_type i = j+2; i < a_nr; i++)
160  hess_mat.elem (i, j) = 0;
161 
162  return info;
163  }
164 
165  template <>
166  octave_idx_type
167  hess<ComplexMatrix>::init (const ComplexMatrix& a)
168  {
169  octave_idx_type a_nr = a.rows ();
170  octave_idx_type a_nc = a.cols ();
171 
172  if (a_nr != a_nc)
173  (*current_liboctave_error_handler) ("hess: requires square matrix");
174 
175  char job = 'N';
176  char side = 'R';
177 
178  octave_idx_type n = a_nc;
179  octave_idx_type lwork = 32 * n;
180  octave_idx_type info;
181  octave_idx_type ilo;
182  octave_idx_type ihi;
183 
184  hess_mat = a;
185  Complex *h = hess_mat.fortran_vec ();
186 
187  Array<double> scale (dim_vector (n, 1));
188  double *pscale = scale.fortran_vec ();
189 
190  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
191  n, F77_DBLE_CMPLX_ARG (h), n, ilo, ihi, pscale, info
192  F77_CHAR_ARG_LEN (1)));
193 
194  Array<Complex> tau (dim_vector (n-1, 1));
195  Complex *ptau = tau.fortran_vec ();
196 
197  Array<Complex> work (dim_vector (lwork, 1));
198  Complex *pwork = work.fortran_vec ();
199 
200  F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, F77_DBLE_CMPLX_ARG (h), n,
201  F77_DBLE_CMPLX_ARG (ptau), F77_DBLE_CMPLX_ARG (pwork), lwork, info));
202 
203  unitary_hess_mat = hess_mat;
204  Complex *z = unitary_hess_mat.fortran_vec ();
205 
206  F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, F77_DBLE_CMPLX_ARG (z), n,
207  F77_DBLE_CMPLX_ARG (ptau), F77_DBLE_CMPLX_ARG (pwork),
208  lwork, info));
209 
210  F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
211  F77_CONST_CHAR_ARG2 (&side, 1),
212  n, ilo, ihi, pscale, n, F77_DBLE_CMPLX_ARG (z), n, info
213  F77_CHAR_ARG_LEN (1)
214  F77_CHAR_ARG_LEN (1)));
215 
216  // If someone thinks of a more graceful way of
217  // doing this (or faster for that matter :-)),
218  // please let me know!
219 
220  if (n > 2)
221  for (octave_idx_type j = 0; j < a_nc; j++)
222  for (octave_idx_type i = j+2; i < a_nr; i++)
223  hess_mat.elem (i, j) = 0;
224 
225  return info;
226  }
227 
228  template <>
229  octave_idx_type
230  hess<FloatComplexMatrix>::init (const FloatComplexMatrix& a)
231  {
232  octave_idx_type a_nr = a.rows ();
233  octave_idx_type a_nc = a.cols ();
234 
235  if (a_nr != a_nc)
236  {
237  (*current_liboctave_error_handler) ("hess: requires square matrix");
238  return -1;
239  }
240 
241  char job = 'N';
242  char side = 'R';
243 
244  octave_idx_type n = a_nc;
245  octave_idx_type lwork = 32 * n;
246  octave_idx_type info;
247  octave_idx_type ilo;
248  octave_idx_type ihi;
249 
250  hess_mat = a;
251  FloatComplex *h = hess_mat.fortran_vec ();
252 
253  Array<float> scale (dim_vector (n, 1));
254  float *pscale = scale.fortran_vec ();
255 
256  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
257  n, F77_CMPLX_ARG (h), n, ilo, ihi, pscale, info
258  F77_CHAR_ARG_LEN (1)));
259 
260  Array<FloatComplex> tau (dim_vector (n-1, 1));
261  FloatComplex *ptau = tau.fortran_vec ();
262 
263  Array<FloatComplex> work (dim_vector (lwork, 1));
264  FloatComplex *pwork = work.fortran_vec ();
265 
266  F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, F77_CMPLX_ARG (h), n,
267  F77_CMPLX_ARG (ptau), F77_CMPLX_ARG (pwork), lwork, info));
268 
269  unitary_hess_mat = hess_mat;
270  FloatComplex *z = unitary_hess_mat.fortran_vec ();
271 
272  F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, F77_CMPLX_ARG (z), n,
273  F77_CMPLX_ARG (ptau), F77_CMPLX_ARG (pwork),
274  lwork, info));
275 
276  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
277  F77_CONST_CHAR_ARG2 (&side, 1),
278  n, ilo, ihi, pscale, n, F77_CMPLX_ARG (z), n, info
279  F77_CHAR_ARG_LEN (1)
280  F77_CHAR_ARG_LEN (1)));
281 
282  // If someone thinks of a more graceful way of
283  // doing this (or faster for that matter :-)),
284  // please let me know!
285 
286  if (n > 2)
287  for (octave_idx_type j = 0; j < a_nc; j++)
288  for (octave_idx_type i = j+2; i < a_nr; i++)
289  hess_mat.elem (i, j) = 0;
290 
291  return info;
292  }
293  }
294 }
octave
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
a_nc
octave_idx_type a_nc
Definition: sylvester.cc:72
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:401
a
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
h
double h
Definition: graphics.cc:11205
a_nr
octave_idx_type a_nr
Definition: sylvester.cc:71
Matrix
Definition: dMatrix.h:37
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
i
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5150
octave_idx_type
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
octave::math::hess::init
octave_idx_type init(const T &a)
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
Array::fortran_vec
const T * fortran_vec(void) const
Definition: Array.h:584
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:409
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
Generated on Wed May 10 2017 15:42:58 for GNU Octave by   doxygen 1.8.8