tminres
SerialExample/ex1.cpp
Go to the documentation of this file.
00001 // tminres is free software; you can redistribute it and/or modify it under the
00002 // terms of the GNU Lesser General Public License (as published by the Free
00003 // Software Foundation) version 2.1 dated February 1999.
00004 //
00005 // Authors:
00006 // - Umberto Villa, Emory University - uvilla@emory.edu
00007 // - Michael Saunders, Stanford University
00008 // - Santiago Akle, Stanford University
00009 
00016 #include <tminres.hpp>
00017 #include "SimpleVector.hpp"
00018 #include <cmath>
00019 #include <fstream>
00020 
00022 /*
00023  * @brief An abstract preconditioner class.
00024  * It should be used to call minres without preconditioner
00025  */
00026 class Preconditioner
00027 {
00028 public:
00030         virtual void Apply(const SimpleVector & X, SimpleVector & Y) const = 0;
00031 };
00032 
00034 /*
00035  * @brief A simple diagonal linear operator.
00036  * The first half of the diagonal entries are positive, the second half are negative
00037  */
00038 class Operator
00039 {
00040 public:
00041 
00043         /*
00044          * @param size_ problem size
00045          */
00046         Operator(int size_) :
00047                 size(size_)
00048         {
00049                 vals = new double[size];
00050                 int i(0);
00051                 for( ; i < size/2; ++i)
00052                         vals[i] = i+1;
00053 
00054                 for( ; i < size; ++i)
00055                         vals[i] = -i;
00056         }
00057 
00058         ~Operator()
00059         {
00060                 delete[] vals;
00061         }
00062 
00064         void Apply(const SimpleVector & X, SimpleVector & Y) const
00065         {
00066                 for(int i(0); i < size; ++i)
00067                         Y[i] = vals[i] * X[i];
00068         }
00069 
00070         void Print(std::ostream & os)
00071         {
00072                 os<< "d = [";
00073                 int i(0);
00074                 for( ; i<size-1; ++i)
00075                         os<<vals[i]<<"; ";
00076                 os<<vals[i]<<"];\n";
00077                 os<<"Op = spdiags(d, 0,"<<size<<", "<<size<<");\n";
00078         }
00079 
00080 private:
00081         int size;
00082         double *vals;
00083 };
00084 
00085 
00091 int main()
00092 {
00093         //(1) Define the size of the problem we want to solve
00094         int size(10);
00095         //(2) Define the linear operator "op" we want to solve.
00096         Operator op(size);
00097         //(3) Define the exact solution (at random)
00098         SimpleVector sol(size);
00099         sol.Randomize( 1 );
00100 
00101         //(4) Define the "rhs" as "rhs = op*sol"
00102         SimpleVector rhs(size);
00103         op.Apply(sol, rhs);
00104         double rhsNorm( sqrt(InnerProduct(rhs,rhs)) );
00105         std::cout << "|| rhs || = " << rhsNorm << "\n";
00106 
00107         //(5) We don't use any preconditioner. Let prec be a null pointer.
00108         Preconditioner * prec = NULL;
00109 
00110         //(6) Use an identically zero initial guess
00111         SimpleVector x(size);
00112         x = 0;
00113 
00114         //(7) Set the minres parameters
00115         double shift(0);
00116         int max_iter(100);
00117         double tol(1e-6);
00118         bool show(true);
00119 
00120         //(8) Solve the problem with minres
00121         MINRES(op, x, rhs, prec, shift, max_iter, tol, show);
00122 
00123         //(9) Compute the error || x_ex - x_minres ||_2
00124         subtract(x, sol, x);
00125         double err2 = InnerProduct(x,x);
00126         std::cout<< "|| x_ex - x_n || = " << sqrt(err2) << "\n";
00127 
00128         std::ofstream fid("ex1.m");
00129         op.Print(fid);
00130         fid<< "rhs = [";
00131         for(int i(0); i<size-1; ++i)
00132                 fid<<rhs[i] <<"; ";
00133         fid<<rhs[size-1] <<"]; \n";
00134         fid<<"[ x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm ] = minres(Op, rhs, [], 0, true, false, 100, 1e-6);\n";
00135 
00136 
00137         return 0;
00138 }
 All Classes Files Functions Friends