crash while building a sym tridiagonal eigmax

[alanedelman]

Writing my first lapack interface. Not sure if I have all the pointers right, etc., but
anyway as long as code just crashes, it's hard to debug. I may try this on a
mac later today, but if anyone can just help me sooner, that would be great.
Here's the exact code:

const liblapack = Base.liblapack_name
import Base.BlasInt

## (ST) Symmetric tridiagonal - maxeig
for (stebz, elty) in
    ((:dstebz_,:Float64),
     (:sstebz_,:Float32))
    @eval begin
        function stebz!(dv::Vector{$elty}, ev::Vector{$elty})
            n = length(dv)
            if length(ev) != (n-1) throw(LapackDimMisMatch("stebz!")) end
            z=Array($elty,1)
            z[1]=0.0    #just to be sure
            m=Array(BlasInt,1)
            nsplit =Array(BlasInt,1)
            w  = Array($elty, n)
            iblock = Array(BlasInt,n)
            isplit = Array(BlasInt,n)
            work = Array($elty, 4*n)
            iwork = Array(BlasInt,3*n)
            info = Array(BlasInt, 1)
            ccall(($(string(stebz)),liblapack), Void,
                  (Ptr{Uint8},Ptr{Uint8},
                   Ptr{BlasInt},Ptr{BlasInt},Ptr{BlasInt},Ptr{BlasInt},Ptr{BlasInt},
                   Ptr{$elty},Ptr{$elty},Ptr{$elty},
                   Ptr{BlasInt},Ptr{BlasInt},
                   Ptr{$elty},Ptr{BlasInt},Ptr{$elty},Ptr{BlasInt},Ptr{BlasInt}),
                   'I','B',
                   &n,0,0,&n,&n,
                   z,dv,ev,
                   m, nsplit,
                   w, iblock,work,iwork,info)
            if info[1] != 0 throw(LapackException(info[1])) end
            w
        end
    end
end

[JeffBezanson]

According to http://www.netlib.org/lapack/double/dstebz.f, the fourth and fifth arguments are float, not integer, and all scalars passed to fortran functions need &, including constants.

[andreasnoack]

Also, the isplit argument was missing in ccall. This one works for me

const liblapack = Base.liblapack_name
import Base.BlasInt

for (stebz, elty) in ((:dstebz_,:Float64),
                      (:sstebz_,:Float32))
    @eval begin
        function stebz!(dv::Vector{$elty}, ev::Vector{$elty})
            n = length(dv)
            if length(ev) != (n-1) throw(LapackDimMisMatch("stebz!")) end
            z = 0.0 #just to be sure
            m = Array(BlasInt,1)
            nsplit = Array(BlasInt,1)
            w = Array($elty, n)
            iblock = Array(BlasInt,n)
            isplit = Array(BlasInt,n)
            work = Array($elty, 4*n)
            iwork = Array(BlasInt,3*n)
            info = Array(BlasInt, 1)
            ccall(($(string(stebz)),liblapack), Void,
                (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty},
                Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt},
                Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                Ptr{BlasInt}, Ptr{BlasInt}),
                &'I', &'B', &n, &0.0, 
                &0.0, &n, &n, &z, 
                dv, ev, m, nsplit,
                w, iblock, isplit, work, 
                iwork, info)
                if info[1] != 0 throw(LapackException(info[1])) end
            w
        end
    end
end

[alanedelman]

cool! guess i was close.

[alanedelman]

THANK YOU BOTH SO MUCH! I can do a 100,000 eigmax in 100 msecs
while a 20,000 eigvals already rakes 18 seconds (factor of 200 more time!)

Something to look into is the different answers slightly numerically. The eigenvalue I'm testing
is ill-conditioned, and not the one I'll need for my real application, so perhaps it falls into the
usual backward error analysis situations.

[ViralBShah]

We should probably implement fcall as discussed in #2167. Is there anything we can do to avoid a crash in such situations?

Also, after the 0.1 freeze is lifted, let's get stebz in. Closing, since this immediate issue is fixed.