tminres
|
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 }