CUSOLVER for linalg routines

[guillochon]

As the CULA library is effectively abandonware for the Mac (the free version links are broken and haven't been updated in 4 years), it is critical to have CUSOLVER working for all functions in order for skcuda to be useful on OS X.

Here's a running list of functions that should have the CUSOLVER option available within linalg that I believe currently don't:

• [x] cho_solve
• [x] inv
• [x] qr
• [ ] eig (partially implemented)
• [ ] dmd

[lebedov]

svd and pinv already have cusolver support in the latest GitHub revision. eig has cusolver support when the input is symmetric/Hermitian.

[guillochon]

Ah, should the svd and pinv tests be enabled then for cusolver?

[lebedov]

There already are unit tests in place for svd/pinv cusolver support, e.g., here and here.

[lebedov]

cho_factor contains cusolver support.

[guillochon]

Oh great! For me at least, the one remaining function that would be great to have is cho_solve. Also inv seems like it should be higher priority.

[lebedov]

cho_solve now contains cusolver support.

[lebedov]

inv now contains cusolver support.

[lebedov]

qr now contains cusolver support.

[botev]

@lebedov Hi, does there exist any interface of cusolver for solving the linear system Ax=b where A is already upper or lower triangular (e.g. does not need the be factored). That is if we want to solve the system A^{-1/2} x = b and we already have L=A^{1/2} is there any interface to directly solve the system without going to factor L further?

[lebedov]

@botev, the cusolver*potrs and cusolver*getrs functions can be used to solve systems that have already been factorized with Cholesky or LU decomposition, respectively (the former is utilized in the cho_solve function, while the latter is used in the inv function). LAPACK does have a function for solving triangular systems (*tptrs), but I'm not aware of an equivalent in CUSOLVER.

[botev]

@lebedov So if I follow correctly the cusovler code for solving A^{-1} x = b where A is symmetric it proceeds like this:

1. Use DnSpotrf to find L such that LL^T = A
2. Use DnSpotrs to solve (LL^T)^{-1}x = b

In my case, however, I want to solve L^{-1}x = b where L is lower triangular, however, note that this is not the factorization of the matrix I want to inverse, but the actual matrix. If I understand correctly this currently does not have a specialised procedure but one has to apply the same two steps?

[lebedov]

@botev Right - I'm not aware of a specialized CUSOLVER procedure when the actual matrix is triangular.

[botev]

Thanks a lot for the clarification.

[botev]

@lebedov What about this - http://scikit-cuda.readthedocs.io/en/latest/generated/skcuda.cublas.cublasStrsm.html isn't that exactly for solving L^{-1} x = b?

[lebedov]

@botev You're right - I haven't tried it, but it looks like it will do the trick.

[Randl]

@lebedov Seems like eig doesn't work for complex (Hermitian) input

[lebedov]

@Randl can you provide an example? The unit test for Hermitian input using the CUSOLVER backend seems to pass in the latest revision.