GNU Octave  3.8.0
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
fCmplxHESS.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2013 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 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "fCmplxHESS.h"
28 #include "f77-fcn.h"
29 #include "lo-error.h"
30 
31 extern "C"
32 {
33  F77_RET_T
34  F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL,
36  const octave_idx_type&, octave_idx_type&,
37  octave_idx_type&, float*, octave_idx_type&
39 
40  F77_RET_T
41  F77_FUNC (cgehrd, CGEHRD) (const octave_idx_type&, const octave_idx_type&,
42  const octave_idx_type&, FloatComplex*,
43  const octave_idx_type&, FloatComplex*,
44  FloatComplex*, const octave_idx_type&,
45  octave_idx_type&);
46 
47  F77_RET_T
48  F77_FUNC (cunghr, CUNGHR) (const octave_idx_type&, const octave_idx_type&,
49  const octave_idx_type&, FloatComplex*,
50  const octave_idx_type&, FloatComplex*,
51  FloatComplex*, const octave_idx_type&,
52  octave_idx_type&);
53 
54  F77_RET_T
55  F77_FUNC (cgebak, CGEBAK) (F77_CONST_CHAR_ARG_DECL,
57  const octave_idx_type&, const octave_idx_type&,
58  const octave_idx_type&, float*,
59  const octave_idx_type&, FloatComplex*,
60  const octave_idx_type&, octave_idx_type&
63 }
64 
67 {
68  octave_idx_type a_nr = a.rows ();
69  octave_idx_type a_nc = a.cols ();
70 
71  if (a_nr != a_nc)
72  {
73  (*current_liboctave_error_handler)
74  ("FloatComplexHESS requires square matrix");
75  return -1;
76  }
77 
78  char job = 'N';
79  char side = 'R';
80 
81  octave_idx_type n = a_nc;
82  octave_idx_type lwork = 32 * n;
83  octave_idx_type info;
84  octave_idx_type ilo;
85  octave_idx_type ihi;
86 
87  hess_mat = a;
89 
90  Array<float> scale (dim_vector (n, 1));
91  float *pscale = scale.fortran_vec ();
92 
93  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
94  n, h, n, ilo, ihi, pscale, info
95  F77_CHAR_ARG_LEN (1)));
96 
97  Array<FloatComplex> tau (dim_vector (n-1, 1));
98  FloatComplex *ptau = tau.fortran_vec ();
99 
100  Array<FloatComplex> work (dim_vector (lwork, 1));
101  FloatComplex *pwork = work.fortran_vec ();
102 
103  F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info));
104 
107 
108  F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
109  lwork, info));
110 
111  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
112  F77_CONST_CHAR_ARG2 (&side, 1),
113  n, ilo, ihi, pscale, n, z, n, info
114  F77_CHAR_ARG_LEN (1)
115  F77_CHAR_ARG_LEN (1)));
116 
117  // If someone thinks of a more graceful way of
118  // doing this (or faster for that matter :-)),
119  // please let me know!
120 
121  if (n > 2)
122  for (octave_idx_type j = 0; j < a_nc; j++)
123  for (octave_idx_type i = j+2; i < a_nr; i++)
124  hess_mat.elem (i, j) = 0;
125 
126  return info;
127 }