tminres
SerialExample/ex1.cpp

A simple example of MINRES() usage without preconditioner.

// tminres is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
// Authors:
// - Umberto Villa, Emory University - uvilla@emory.edu
// - Michael Saunders, Stanford University
// - Santiago Akle, Stanford University

#include <tminres.hpp>
#include "SimpleVector.hpp"
#include <cmath>
#include <fstream>

/*
 * @brief An abstract preconditioner class.
 * It should be used to call minres without preconditioner
 */
class Preconditioner
{
public:
        virtual void Apply(const SimpleVector & X, SimpleVector & Y) const = 0;
};

/*
 * @brief A simple diagonal linear operator.
 * The first half of the diagonal entries are positive, the second half are negative
 */
class Operator
{
public:

        /*
         * @param size_ problem size
         */
        Operator(int size_) :
                size(size_)
        {
                vals = new double[size];
                int i(0);
                for( ; i < size/2; ++i)
                        vals[i] = i+1;

                for( ; i < size; ++i)
                        vals[i] = -i;
        }

        ~Operator()
        {
                delete[] vals;
        }

        void Apply(const SimpleVector & X, SimpleVector & Y) const
        {
                for(int i(0); i < size; ++i)
                        Y[i] = vals[i] * X[i];
        }

        void Print(std::ostream & os)
        {
                os<< "d = [";
                int i(0);
                for( ; i<size-1; ++i)
                        os<<vals[i]<<"; ";
                os<<vals[i]<<"];\n";
                os<<"Op = spdiags(d, 0,"<<size<<", "<<size<<");\n";
        }

private:
        int size;
        double *vals;
};


int main()
{
        //(1) Define the size of the problem we want to solve
        int size(10);
        //(2) Define the linear operator "op" we want to solve.
        Operator op(size);
        //(3) Define the exact solution (at random)
        SimpleVector sol(size);
        sol.Randomize( 1 );

        //(4) Define the "rhs" as "rhs = op*sol"
        SimpleVector rhs(size);
        op.Apply(sol, rhs);
        double rhsNorm( sqrt(InnerProduct(rhs,rhs)) );
        std::cout << "|| rhs || = " << rhsNorm << "\n";

        //(5) We don't use any preconditioner. Let prec be a null pointer.
        Preconditioner * prec = NULL;

        //(6) Use an identically zero initial guess
        SimpleVector x(size);
        x = 0;

        //(7) Set the minres parameters
        double shift(0);
        int max_iter(100);
        double tol(1e-6);
        bool show(true);

        //(8) Solve the problem with minres
        MINRES(op, x, rhs, prec, shift, max_iter, tol, show);

        //(9) Compute the error || x_ex - x_minres ||_2
        subtract(x, sol, x);
        double err2 = InnerProduct(x,x);
        std::cout<< "|| x_ex - x_n || = " << sqrt(err2) << "\n";

        std::ofstream fid("ex1.m");
        op.Print(fid);
        fid<< "rhs = [";
        for(int i(0); i<size-1; ++i)
                fid<<rhs[i] <<"; ";
        fid<<rhs[size-1] <<"]; \n";
        fid<<"[ x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm ] = minres(Op, rhs, [], 0, true, false, 100, 1e-6);\n";


        return 0;
}
 All Classes Files Functions Friends