 #include <cusp/array2d.h>
 #include <cusp/multiply.h>
 
-template <class MatrixType>
-void _TestJacobiRelaxation(void)
+template <typename Matrix>
+void TestJacobiRelaxation(void)
 {
-    typedef typename MatrixType::memory_space Space;
+    typedef typename Matrix::memory_space Space;
 
     cusp::array2d<float, Space> M(5,5);
     M(0,0) = 1.0;  M(0,1) = 1.0;  M(0,2) = 2.0;  M(0,3) = 0.0;  M(0,4) = 0.0; 

    expected[4] =  1.625;  // (5 + 8    ) / 8   = 1.625


    MatrixType A(M);
    cusp::relaxation::jacobi<float, Space> relax(A);

    relax(A, b, x);

    ASSERT_EQUAL(x, expected);
}

template <class SparseMatrix>
void TestJacobiRelaxation(void)
{
    _TestJacobiRelaxation<SparseMatrix>();
}
DECLARE_SPARSE_MATRIX_UNITTEST(TestJacobiRelaxation);


template <class MatrixType>
void _TestJacobiRelaxationWithWeighting(void)
{
    typedef typename MatrixType::memory_space Space;

    cusp::array2d<float, Space> M(2,2);
    M(0,0) = 2.0;  M(0,1) = 1.0;
    M(1,0) = 1.0;  M(1,1) = 3.0;

    MatrixType A(M);

    // use default omega
    {

        cusp::array1d<float, Space> x(2, -1.0);
        cusp::relaxation::jacobi<float, Space> relax(A, 0.5);
        relax(A, b, x);
        ASSERT_ALMOST_EQUAL(x[0], 1.0);
        ASSERT_ALMOST_EQUAL(x[1], 0.5);
    }
    
    // override default omega
