2016-01-04 5 views
3

В настоящее время я пытаюсь перевести мой существующий код Python в Julia, и мне нужно вычислить Cholesky Decomposition сложной комплексной матрицы. Правильной процедурой LAPACK является cpbtrf (тот, который в настоящее время называется SciPy), и я изо всех сил пытаюсь заставить его работать в Julia.Как вызвать код LAPACK (cpbtrf) в Julia

Я не уверен, какие дополнительные детали дать, я довольно новичок в Джулии, и я уверен, что делаю что-то глупое. Вызов LAPACK возвращает 1 в информационной переменной, что указывает на то, что что-то не является положительным, но я знаю, что это (SciPy с радостью разлагает одну и ту же матрицу).

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) 

Редактировать: Просто пояснить, это пример скелета. Код, который я пишу, касается матриц порядка 10^5 или больше и может нуждаться в пента-, гекса-, гепта-диагональных матрицах и т. Д. Мне нужен алгоритм с привязкой к диапазону.

ответ

2

Это все правильно, за исключением подпрограммы LAPACK. Вы используете 128-битные комплексные номера, поэтому вы должны использовать :zpbtrf_ вместо :cpbtrf_.

+0

Подчеркивая также отсутствие моего Фортрана! Большое спасибо, это работает отлично! – jimbarrett27

 Смежные вопросы

  • Нет связанных вопросов^_^