Thanks for the reply.
There is a lot of staff on OpenCurrent, even a complete implementation of finite difference in 3D. However, I only have a simple problem solved with a fast coefficient computation and a sparse matrices multiplication (I was thinking at the gradient descent in the Krylov method (since I have sparce matrices with a bunch of zeros and just few coefficients).
However It is not so easy to interpret what is going on inside Cusp. How can I display (eg. with printf) the content of a cusp-array1D (the last line in my code gives an error).
Assume I have this code with Ax=b, where A is a matrix of coefficient, x is a discretization of point of a function and b is the derivative as in the picture (solving the equations I obtain the derivatives as b):
// create an empty sparse matrix structure (CSR format)
cusp::csr_matrix<int, float, cusp::device_memory> A;
// initialize matrix
cusp::gallery::poisson5pt(A, 10, 10);
// allocate storage for solution (x) and right hand side (b)
cusp::array1d<float, cusp::device_memory> x(A.num_rows, 0);
cusp::array1d<float, cusp::device_memory> b(A.num_rows, 1);
// set stopping criteria:
// iteration_limit = 100
// relative_tolerance = 1e-6
cusp::verbose_monitor<float> monitor(b, 100, 1e-6);
// set preconditioner (identity)
cusp::identity_operator<float, cusp::device_memory> M(A.num_rows, A.num_rows);
// solve the linear system A x = b
cusp::krylov::cg(A, x, b, monitor, M);
printf("Val of A = %i \n", x.begin() );