Owijanie funkcji LAPACKE za pomocą Cython

Question

Próbuję zawinąć funkcję LAPACK dgtsv (solver dla tridiagonal systemów równań) przy użyciu Cython.

Natknąłem się na this previous answer, ale ponieważ dgtsv nie jest jedną z funkcji LAPACK owiniętych w scipy.linalg, nie sądzę, żebym mógł użyć tego konkretnego podejścia. Zamiast tego próbuję podążać za this example.

Oto treść mojego lapacke.pxd pliku:

ctypedef int lapack_int 

cdef extern from "lapacke.h" nogil: 

    int LAPACK_ROW_MAJOR 
    int LAPACK_COL_MAJOR 

    lapack_int LAPACKE_dgtsv(int matrix_order, 
          lapack_int n, 
          lapack_int nrhs, 
          double * dl, 
          double * d, 
          double * du, 
          double * b, 
          lapack_int ldb) 

... oto mój cienki Cython wrapper w _solvers.pyx:

#!python 

cimport cython 
from lapacke cimport * 

cpdef TDMA_lapacke(double[::1] DL, double[::1] D, double[::1] DU, 
        double[:, ::1] B): 

    cdef: 
     lapack_int n = D.shape[0] 
     lapack_int nrhs = B.shape[1] 
     lapack_int ldb = B.shape[0] 
     double * dl = &DL[0] 
     double * d = &D[0] 
     double * du = &DU[0] 
     double * b = &B[0, 0] 
     lapack_int info 

    info = LAPACKE_dgtsv(LAPACK_ROW_MAJOR, n, nrhs, dl, d, du, b, ldb) 

    return info 

... i oto skrypt otoki Python i testy:

import numpy as np 
from scipy import sparse 
from cymodules import _solvers 


def trisolve_lapacke(dl, d, du, b, inplace=False): 

    if (dl.shape[0] != du.shape[0] or dl.shape[0] != d.shape[0] - 1 
      or b.shape != d.shape): 
     raise ValueError('Invalid diagonal shapes') 

    if b.ndim == 1: 
     # b is (LDB, NRHS) 
     b = b[:, None] 

    # be sure to force a copy of d and b if we're not solving in place 
    if not inplace: 
     d = d.copy() 
     b = b.copy() 

    # this may also force copies if arrays are improperly typed/noncontiguous 
    dl, d, du, b = (np.ascontiguousarray(v, dtype=np.float64) 
        for v in (dl, d, du, b)) 

    # b will now be modified in place to contain the solution 
    info = _solvers.TDMA_lapacke(dl, d, du, b) 
    print info 

    return b.ravel() 


def test_trisolve(n=20000): 

    dl = np.random.randn(n - 1) 
    d = np.random.randn(n) 
    du = np.random.randn(n - 1) 

    M = sparse.diags((dl, d, du), (-1, 0, 1), format='csc') 
    x = np.random.randn(n) 
    b = M.dot(x) 

    x_hat = trisolve_lapacke(dl, d, du, b) 

    print "||x - x_hat|| = ", np.linalg.norm(x - x_hat) 

Niestety, test_trisolve tylko se gfaults na połączenie z _solvers.TDMA_lapacke. Jestem prawie pewien, że moja setup.py jest poprawna - pokazuje, że _solvers.so jest połączone z poprawnymi bibliotekami współdzielonymi w środowisku wykonawczym.

Nie jestem pewien, jak przejść dalej - jakiekolwiek pomysły?

---

Krótka aktualizacja:

dla mniejszych wartości n ja zwykle nie od razu się naruszenia ochrony pamięci, ale nie dostać wyniki nonsensowne (|| x - x_hat || powinien być bardzo bliski 0):

In [28]: test_trisolve2.test_trisolve(10) 
0 
||x - x_hat|| = 6.23202576396 

In [29]: test_trisolve2.test_trisolve(10) 
-7 
||x - x_hat|| = 3.88623414288 

In [30]: test_trisolve2.test_trisolve(10) 
0 
||x - x_hat|| = 2.60190676562 

In [31]: test_trisolve2.test_trisolve(10) 
0 
||x - x_hat|| = 3.86631743386 

In [32]: test_trisolve2.test_trisolve(10) 
Segmentation fault 

Zazwyczaj LAPACKE_dgtsv powraca z kodem 0 (który powinien wskazywać sukces), ale od czasu do czasu otrzymuję -7, co oznacza, że ​​argument 7 (b) miał nieprawidłową wartość. To, co się dzieje, to faktyczna modyfikacja tylko pierwszej wartości b. Jeśli nadal dzwonię pod numer test_trisolve, w końcu trafię na segment, nawet jeśli n jest mały.

Answer 1

OK, doszedłem do tego w końcu - wydaje mi się, że źle zrozumiałem, o czym w tym przypadku mówi wiersz i kolumna-major.

Ponieważ tablice sąsiadujące z C podążają za rzędem-głównym porządkiem, założyłem, że powinienem podać LAPACK_ROW_MAJOR jako pierwszy argument dla LAPACKE_dgtsv.

W rzeczywistości, jeśli zmienię

info = LAPACKE_dgtsv(LAPACK_ROW_MAJOR, ...) 

do

info = LAPACKE_dgtsv(LAPACK_COL_MAJOR, ...) 

wtedy moja funkcja działa:

test_trisolve2.test_trisolve() 
0 
||x - x_hat|| = 6.67064747632e-12 

To wydaje się dość intuicyjne do mnie - może ktoś wyjaśnić, dlaczego to jest przypadek?

Answer 2

Mimo że pytanie jest dość stare, wydaje się, że nadal jest ono istotne. Obserwowane zachowanie się jest wynikiem błędnej parametru LDB:

• Fortran tablice są kol głównym a głównym wymiar tablicy B odpowiada N. Zatem LDB> = max (1, Y).
• Z głównym wierszem LDB odpowiada NRHS i dlatego warunek LDB> = max (1, NRHS) musi być spełniony.

Komentarz # B oznacza (LDB, NRHS) nie jest prawidłowe, ponieważ b ma wymiar (LDB, N) i LDB powinien wynosić 1 w tym przypadku.

Przejście z LAPACK_ROW_MAJOR na LAPACK_COL_MAJOR rozwiązuje problem, o ile NRHS jest równe 1. Układ pamięci col-dur (N, 1) jest taki sam, jak wiersz główny (1, N). Jednak zawiedzie, jeśli NRHS jest większa niż 1.