man dgehrd ...

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