void VCycle::Solve(ublas::matrix<double> & inout_u, ublas::matrix<double> & output_r, ublas::matrix<double> & input_f) {

    unsigned int size = inout_u.size1();

    //printMatrix(inout_u);

    //relax first
    RelaxError(m_nInitIterations, input_f, inout_u, output_r);

    //printMatrix(inout_u);

    //restrict
    ublas::matrix<double> restrictedResid = Restrict(output_r);

    //new size
    unsigned int new_size = restrictedResid.size1();

    ublas::matrix<double> error;
    //if further restriction makes sense, go to coarser grid
    if (new_size >= 3) {

        //new (zero) initial guess
        ublas::matrix<double> * new_U = new ublas::matrix<double>(new_size,new_size,0);

        //new residuals of correct size
        ublas::matrix<double> * new_residuals = new ublas::matrix<double>(new_size,new_size,0);

        //go to coarser grid
        Solve(*new_U,*new_residuals,restrictedResid);

        //prolong error afterwards
        error = Prolong(*new_U);
    } else {
        error = ublas::matrix<double>(size,size,0.0);
    }

    //correct
    inout_u = inout_u + error;

    //relax again
    RelaxError(m_nEndIterations,input_f,inout_u,output_r);
}