[SciPy-user] Hessenberg form / Improving the wrapper for dgehrd

[Nils Wagner]

Hi all,

I was going to build a wrapper for dgehrd in order to have 
a function similar to hess in Matlab.

I am quite sure, that my solution is not optimal. The output 
contains not only the Hessenberg form of the input matrix but
the array tau. How can I circumvent this ?

>man dgehrd
...
the elements below the first subdiagonal, with the  array  TAU, 
represent
the orthogonal matrix Q as a product of elementary reflectors. 
See Further Details.  LDA     (input) INTEGER The leading
dimension of the array A.  LDA >= max(1,N).

Therefore I would appreciate any hint for improving the 
interface according to the rules for a possible implementation in scipy.

Thanks in advance.

Nils

!    -*- f90 -*-
python module hess ! in
    interface  ! in :hess
        subroutine dgehrd(n,ilo,ihi,a,lda,tau,work,lwork,info)
!in:hess:dgehrd.f

            callstatement
(*f2py_func)(&n,&ilo,&ihi,a,&lda,tau,work,&lwork,&info)

            callprotoargument
int*,int*,int*,double*,int*,double*,double*,int*,int*

            integer intent(hide):: n=shape(a,1)
            integer intent(hide):: ilo=1
            integer intent(hide):: ihi=n
            double precision dimension(lda,*),intent(in,out,copy,out=h)
:: a
            integer
optional,check(shape(a,0)==lda),depend(a)::lda=shape(a,0)
            double precision dimension(n),intent(out) :: tau
            double precision dimension(MAX(1,lwork)),intent(out) :: work
            integer intent(in):: lwork
            integer intent(out):: info
        end subroutine dgehrd
    end interface
end python module hess
-------------- next part --------------
from scipy import *
import hess
#
# First example 5.4.6
# Biswa Nath Datta
# Numerical linear algebra and applications
# Brooks/Cole Publishing Company
# (1995)
#
a=zeros((3,3),Float)
a[0,0] = 1
a[0,1] = 2
a[0,2] = 5
a[1,0] = 3
a[1,1] = 7
a[1,2] = 9
a[2,0] = 2
a[2,1] = 5
a[2,2] = 3
#
# test example 5.4.7
# page 152 
# Biswa Nath Datta
# Numerical linear algebra and applications
# Brooks/Cole Publishing Company
# (1995)
#
a[0,0] = 0
a[0,1] = 1
a[0,2] = 1
a[1,0] = 1
a[1,1] = 2
a[1,2] = 1
a[2,0] = 1
a[2,1] = 1
a[2,2] = 1
#
# solution to 5.4.7
#
h1 = zeros((3,3),Float)
h1[0,1] = -sqrt(2)
h1[1,0] = -sqrt(2)
h1[1,1] = 2.5
h1[1,2] = 0.5
h1[2,1] = 0.5
h1[2,2] = 0.5



# workspace query
a,tau,work,info=hess.dgehrd(a,lwork=-1)
lwork=int(work[0])
h,tau,work,info=hess.dgehrd(a,lwork=lwork)
print 'Hessenberg form'
print h


def do_hess(mat):
   h,tau,work,info=hess.dgehrd(a,lwork=-1)
   lwork=int(work[0])
   h,tau,work,info=hess.dgehrd(a,lwork=lwork)
   return h

h2 = do_hess(a)
print h2