Sparse matrix support in Armadillo is not yet complete.
You can use ARPACK for (eigen-)decomposition of sparse matrix. Sparse matrix solvers will probably comes in the next release which may use the CXSparse library from the SuiteSparse project.
Question
I'm writing a basic FEM program using Armadillo. I use sp_mat
and vec
as matrix and vector type. The problem is, that when I do solve(X, b)
i get an error. Could it be that solve
does not support sp_mat
. Any alternatives except of using dense matrices? Below is the code, where the sp_mat
doesn't compile. If I use the the commented line mat A
instead it works fine.
int N = 3;
double h = 1./N;
//mat A = zeros<mat>(N+1,N+1);
sp_mat A(N+1,N+1);
for(int i=0;i<=N;i++) {
if(i>0) {A(i,i-1)=-1.;}
A(i,i)=2.;
if(i<N) {A(i,i+1)=-1.;}
}
A(N,N)=1;
vec b = zeros(N+1);
for(int i=0;i<=N;i++) {
b(i)=h;
}
vec zeta = solve(A,b);
cout << zeta;
The error:
make all
Building file: ../src/FEM.cpp
Invoking: GCC C++ Compiler
g++ -O0 -g3 -Wall -c -fmessage-length=0 -MMD -MP -MF"src/FEM.d" -MT"src/FEM.d" -o "src/FEM.o" "../src/FEM.cpp"
../src/FEM.cpp: In function ‘int main()’:
../src/FEM.cpp:37:22: error: no matching function for call to ‘solve(arma::sp_mat&, arma::vec&)’
vec zeta = solve(A,b);
^
../src/FEM.cpp:37:22: note: candidates are:
In file included from /usr/include/armadillo:397:0,
from ../src/FEM.cpp:10:
/usr/include/armadillo_bits/fn_solve.hpp:25:1: note: template<class T1, class T2> const arma::Glue<T1, T2, arma::glue_solve> arma::solve(const arma::Base<typename T1::elem_type, T1>&, const arma::Base<typename T1::elem_type, T2>&, bool, const typename arma::arma_blas_type_only<typename T1::elem_type>::result*)
solve
^
/usr/include/armadillo_bits/fn_solve.hpp:25:1: note: template argument deduction/substitution failed:
../src/FEM.cpp:37:22: note: ‘arma::sp_mat {aka arma::SpMat<double>}’ is not derived from ‘const arma::Base<typename T1::elem_type, T1>’
vec zeta = solve(A,b);
^
In file included from /usr/include/armadillo:397:0,
from ../src/FEM.cpp:10:
/usr/include/armadillo_bits/fn_solve.hpp:44:1: note: template<class T1, class T2> const arma::Glue<T1, T2, arma::glue_solve> arma::solve(const arma::Base<typename T1::elem_type, T1>&, const arma::Base<typename T1::elem_type, T2>&, const char*, const typename arma::arma_blas_type_only<typename T1::elem_type>::result*)
solve
^
/usr/include/armadillo_bits/fn_solve.hpp:44:1: note: template argument deduction/substitution failed:
../src/FEM.cpp:37:22: note: ‘arma::sp_mat {aka arma::SpMat<double>}’ is not derived from ‘const arma::Base<typename T1::elem_type, T1>’
vec zeta = solve(A,b);
^
In file included from /usr/include/armadillo:397:0,
from ../src/FEM.cpp:10:
/usr/include/armadillo_bits/fn_solve.hpp:67:1: note: template<class T1, class T2> const arma::Glue<T1, T2, arma::glue_solve_tr> arma::solve(const arma::Op<T1, arma::op_trimat>&, const arma::Base<typename T1::elem_type, T2>&, bool, const typename arma::arma_blas_type_only<typename T1::elem_type>::result*)
solve
^
/usr/include/armadillo_bits/fn_solve.hpp:67:1: note: template argument deduction/substitution failed:
../src/FEM.cpp:37:22: note: ‘arma::sp_mat {aka arma::SpMat<double>}’ is not derived from ‘const arma::Op<T1, arma::op_trimat>’
vec zeta = solve(A,b);
^
In file included from /usr/include/armadillo:397:0,
from ../src/FEM.cpp:10:
/usr/include/armadillo_bits/fn_solve.hpp:87:1: note: template<class T1, class T2> const arma::Glue<T1, T2, arma::glue_solve_tr> arma::solve(const arma::Op<T1, arma::op_trimat>&, const arma::Base<typename T1::elem_type, T2>&, const char*, const typename arma::arma_blas_type_only<typename T1::elem_type>::result*)
solve
^
/usr/include/armadillo_bits/fn_solve.hpp:87:1: note: template argument deduction/substitution failed:
../src/FEM.cpp:37:22: note: ‘arma::sp_mat {aka arma::SpMat<double>}’ is not derived from ‘const arma::Op<T1, arma::op_trimat>’
vec zeta = solve(A,b);
^
In file included from /usr/include/armadillo:397:0,
from ../src/FEM.cpp:10:
/usr/include/armadillo_bits/fn_solve.hpp:110:1: note: template<class T1, class T2> bool arma::solve(arma::Mat<typename T1::elem_type>&, const arma::Base<typename T1::elem_type, T1>&, const arma::Base<typename T1::elem_type, T2>&, bool, const typename arma::arma_blas_type_only<typename T1::elem_type>::result*)
solve
^
/usr/include/armadillo_bits/fn_solve.hpp:110:1: note: template argument deduction/substitution failed:
../src/FEM.cpp:37:22: note: ‘arma::sp_mat {aka arma::SpMat<double>}’ is not derived from ‘arma::Mat<typename T1::elem_type>’
vec zeta = solve(A,b);
^
In file included from /usr/include/armadillo:397:0,
from ../src/FEM.cpp:10:
/usr/include/armadillo_bits/fn_solve.hpp:139:1: note: template<class T1, class T2> bool arma::solve(arma::Mat<typename T1::elem_type>&, const arma::Base<typename T1::elem_type, T1>&, const arma::Base<typename T1::elem_type, T2>&, const char*, const typename arma::arma_blas_type_only<typename T1::elem_type>::result*)
solve
^
/usr/include/armadillo_bits/fn_solve.hpp:139:1: note: template argument deduction/substitution failed:
../src/FEM.cpp:37:22: note: ‘arma::sp_mat {aka arma::SpMat<double>}’ is not derived from ‘arma::Mat<typename T1::elem_type>’
vec zeta = solve(A,b);
^
make: *** [src/FEM.o] Error 1
Solution
Sparse matrix support in Armadillo is not yet complete.
You can use ARPACK for (eigen-)decomposition of sparse matrix. Sparse matrix solvers will probably comes in the next release which may use the CXSparse library from the SuiteSparse project.
OTHER TIPS
As of version 5.0, Armadillo has the spsolve() function for solving sparse systems.