void MultiGridSolver::RelaxError(unsigned int nIterations, ublas::matrix<double> & input_f, ublas::matrix<double> & inout_u, ublas::matrix<double> & output_r) {

    //ensure consistent sizes
    assert(input_f.size1() == output_r.size1() && output_r.size1() == inout_u.size1());
    assert(input_f.size2() == output_r.size2() && output_r.size2() == inout_u.size2());

    //size
    unsigned int n = inout_u.size1();

    //spacing
    double spacing = 1.0/n;

    for (unsigned int n_iter = 0; n_iter < nIterations; n_iter++) {

        //red halfiteration
        for (unsigned int i = 1; i < n-1; i++) {
            for (unsigned int j = 1; j < n-1; j++) {
                if ( (i+j) %2 == 0) {   //optimize!
                    inout_u(i,j) = (inout_u(i-1,j) + inout_u(i+1,j)+ inout_u(i,j-1) + inout_u(i,j+1) - spacing*spacing*input_f(i,j))/4.0;
                }
            }
        }

        //black halfiteration
        for (unsigned int i = 1; i < n-1; i++) {
            for (unsigned int j = 1; j < n-1; j++) {
                if ( (i+j) %2 == 1) {   //optimize!
                    inout_u(i,j) = (inout_u(i-1,j) + inout_u(i+1,j)+ inout_u(i,j-1) + inout_u(i,j+1) - spacing*spacing*input_f(i,j))/4.0;
                }
            }
        }

    }

    //residual
    for (unsigned int i = 1; i < n-1; i++) {
        for (unsigned int j = 1; j < n-1; j++) {
            double laplacian = (-inout_u(i-1,j) - inout_u(i+1,j) - inout_u(i,j-1) - inout_u(i,j+1) + 4*inout_u(i,j)) / (spacing*spacing);
            output_r(i,j) = (laplacian - input_f(i,j))*(laplacian-input_f(i,j));
        }
    }
}