test/test_ml_laplace.cpp

Go to the documentation of this file.

//
// C++ Implementation: test_romberg
//
// Description: unit tests ML preconditioner
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <unistd.h>
#include <limits>

#include <tut.h>

#include "rlog/rlog.h"

#include "ml_include.h"
#include "Teuchos_ParameterList.hpp"
#include "Epetra_Comm.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "AztecOO.h"
#include "Epetra_Operator.h"
#include "Epetra_MultiVector.h"
#include "Epetra_LinearProblem.h"
#include "ml_epetra_preconditioner.h"

#include "OutputCapturer.h"
#include "balancedmtxreader.h"
#include "mtxreader.h"

#undef RLOG_TIME_TSC
#include "rlog/RLogTime.h"

namespace tut {
    Epetra_Comm& get_comm_world();
}

namespace
{

    struct testdata
    {
        testdata()
            : tol(100 * std::numeric_limits<double>::epsilon()),
              comm(tut::get_comm_world())
        {
            if (H == 0)
                read_matrices(tut::get_comm_world());
        }

        static void read_matrices(Epetra_Comm& comm) {
            rInfo("Reading test matrices...");
            string filename;
            BaseMtxReader *reader = 0;

            // read matrix Y_transp
            filename = "test/ml_laplace/H.mtx";
            if (access(filename.c_str(), R_OK)) {
                rError("Unable to access H matrix file.");
                exit(1);
            }
            reader = new BalancedMtxReader(filename, comm);
            H = reader->read();
            delete reader;
        }

        const double tol;
        Epetra_Comm& comm;
        static Epetra_CrsMatrix* H;
    };

    // Initialise static members
    Epetra_CrsMatrix* testdata::H = 0;

    typedef tut::test_group<testdata> testgroup;
    typedef testgroup::object testobject;
    testgroup group("ml_laplace");

    template<>
    template<>
    void testobject::test<1>()
    {
        // Construct ML preconditioner
        Teuchos::ParameterList MLList;
        ML_Epetra::SetDefaults("SA", MLList);
        MLList.set("aggregation: type", "Uncoupled-MIS");
        ostringstream buf;
        buf << "ML parameters:" << endl;
        MLList.print(buf, 8);
        if (comm.MyPID() == 0)
            rDebug(buf.str().c_str());

        OutputCapturer capt;
        capt.begin_capture();
        Epetra_Operator *prec = new ML_Epetra::MultiLevelPreconditioner(*H, MLList, true);
        capt.end_capture();
        rDebug("%s", capt.get_stdout().c_str());

        // Right hand side
        Epetra_Vector b(H->OperatorRangeMap());
        b.PutScalar(1.0);

        // Initial guess
        Epetra_Vector x(H->OperatorDomainMap());
        x.PutScalar(0.0);

        // Setup iterative solver
        Epetra_LinearProblem problem(H, &x, &b);
        AztecOO solver(problem);
        solver.SetPrecOperator(prec);

        // Solve
        capt.begin_capture();
        solver.Iterate(1000, 1.0e-10);
        capt.end_capture();
        rDebug("%s", capt.get_stdout().c_str());

        // Compute Residual
        Epetra_Vector bcomp(H->OperatorRangeMap());
        H->Apply(x, bcomp);
        Epetra_Vector resid(H->OperatorRangeMap());
        resid.Update(1.0, b, -1.0, bcomp, 0.0);
        double residual;
        resid.Norm2(&residual);
        rDebug("Residual=%e", residual);
        ensure(residual < 1.0e-7);
    }

};

---