4

I'm currently trying to translate my existing Python code into Julia, and I need to compute a Cholesky Decomposition of a banded, complex matrix. The correct LAPACK routine is cpbtrf (the one currently called by SciPy), and I'm struggling to get it to work in Julia.

I'm not sure what extra details to give, I'm pretty new to Julia and I'm sure I'm doing something stupid. The LAPACK call returns a 1 in the info variable, indicating that something isn't positive definite, but I know it is (SciPy happily decomposes the same matrix).

BlasInt = Base.LinAlg.BlasInt
chk = Base.LinAlg.chkstride1

function cholesky_banded!(ab::StridedMatrix{Complex128}, uplo::Char, n::Integer, kd::Integer)

    chk(ab)

    ldab = size(ab,1)

    info = Ref{BlasInt}()

    ccall((:cpbtrf_,Base.liblapack_name),Void,(Ptr{UInt8},Ptr{BlasInt},Ptr{BlasInt},
    Ptr{Complex128},Ptr{BlasInt},Ptr{BlasInt}),&uplo,&n,&kd,ab,&ldab,info)

    ab, info[]

end

mat = zeros(Complex128,2,3)
mat[1,1:end] = 2
mat[2,1:end-1] = -1

cholesky_banded!(mat,'L',3,1)

edit: Just to clarify, this is a skeleton example. The code I'm writing deals with matrices of order 10^5 or bigger, and can need penta-, hexa-, hepta-diagonal matrices and so on. I need a banded-specific algorithm.

jimbarrett27
  • 157
  • 1
  • 2
  • 9

1 Answers1

3

It's all correct except for the LAPACK subroutine. You are using 128 bit complex numbers so you should use :zpbtrf_ instead of :cpbtrf_.

Andreas Noack
  • 1,370
  • 6
  • 11