Problems setting up matrix = solve(matrix, matrix, tag)

[Mr Dan]

Hi,

I am struggling with setting up a Matrix-Matrix solver for matrix inversion. My problem occurs at compile time with Visual Studio 2012.

The code snippet looks as follows:

// Initialize Matrix 1
viennacl::matrix<ScalarType> gpu_covmat(nPts, nPts); // Initialize covariance matrix
viennacl::copy(cpu_covmat, gpu_covmat);

// Do some computation
....

// Initialize Matrix 2
viennacl::matrix<ScalarType> gpu_invCovmat(nPts, nPts);

// Set up solver and invert
viennacl::linalg::cg_tag solver_tag(1e-8, 300);
viennacl::identity_matrix<ScalarType> iMatrix(nPts);
viennacl::matrix<ScalarType> sMatrix(iMatrix);
gpu_invCovmat = viennacl::linalg::solve(gpu_covmat, sMatrix, solver_tag);

I get a bunch of errors related to the last line which all look somehow like this:

Error 4 error C2784: 'viennacl::vector_expression<const viennacl::vector_base<T="">,const viennacl::vector_tuple<T>,viennacl::op_inner_prod> viennacl::linalg::inner_prod(const viennacl::vector_base<T> &,const viennacl::vector_tuple<T> &)' : could not deduce template argument for 'const viennacl::vector_base<T> &' from 'const viennacl::matrix<SCALARTYPE,F>' ...\ViennaCL-1.5.1\viennacl\linalg\cg.hpp 103 1 custom-kernels

Which is odd, since I want to solve matrix/matrix and not matrix/vector.

Can anybody help?

Thanks,

Daniel

[Karl Rupp]

Hi,
the iterative solvers are only available for the standard system form
Ax = b.
So, if your matrices are sparse, then you need to call the solver separately for each column of the matrix. On the other hand, if your matrices are dense, then you need to LU-factorize the system matrix instead, followed by a forward lower triangular and a backward upper triangular solve. Have a look at examples/tutorial/blas2.cpp on how to do this.

(Note that the performance of LU factorization in ViennaCL is not great at the moment)

Best regards,
Karli

[Mr Dan]

Hi Karli,

thanks for the quick reply! First I have to say, that I like the concept of ViennaCL and ease of access to OpenCL devices. By using the provided base functionality, I can go pretty far without any need to dig to deep into custom kernel creation. It is pretty cool!

Anyway, I am a little bit confused, regarding the documentation/ manual. Both declare many overloads of my required solve, including the one for matrix/ matrix:

http://viennacl.sourceforge.net/doc/namespaceviennacl_1_1linalg.html#ae4e9abff733ebe1feda607b4ccea8fc1

or respectively:

http://viennacl.sourceforge.net/viennacl-manual-current.pdf page 25).

As for my problem, I could also use cholesky-decomposition in this particular case, since my mytrix is positive definit and symmetric. The documentation indicates again a method for this:

http://viennacl.sourceforge.net/doc/namespaceviennacl_1_1linalg_1_1detail_1_1spai.html.

Later I'll also need a determinant which is declared here, though only for ublas types:

http://viennacl.sourceforge.net/doc/namespaceviennacl_1_1linalg_1_1detail_1_1spai.html#af20412963cd669e1fe798c45e07bdf74

Could you briefly give a hint for proper usage of the documenetation. When designing my algorithm, how can I figure out which methods are implemented and which are not?

Thanks,

Daniel

[Karl Rupp]

Hi Daniel,

currently we don't have a Cholesky factorization for symmetric dense matrices implemented, the one you found is used within the sparse approximate inverse (SPAI) preconditioner only and factorizes many small matrix blocks simultaneously. So I can only advise you to use the LU factorization for now. On parallel hardware it's fairly hard to really make use of symmetry, so you don't lose much in terms of performance.

As for the determinant, it's best to compute it from the LU-factors:
det(A) = det(LU) = det(L)det(U)
where det(L) and det(U) are just the products of the diagonal entries, since L and U are triangular matrices.

With respect to documentation the preferred source for picking the respective functionality is the PDF manual. Once you know the functions you want to use, browse the Doxygen-pages to find the details. We want to make this less cumbersome and will provide an everything-HTML solution with the 1.6.0 release soon (see ticket here: https://github.com/viennacl/viennacl-dev/issues/18 for some explanations and here: http://viennashe.sourceforge.net/doc/ for another project where this is all integrated already)

Best regards,
Karli