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
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
