Thyra Version of the Day
sillyCgSolve.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_SILLY_CG_SOLVE_HPP
00043 #define THYRA_SILLY_CG_SOLVE_HPP
00044 
00045 #include "Thyra_LinearOpBase.hpp"
00046 #include "Thyra_VectorStdOps.hpp"
00047 #include "Thyra_AssertOp.hpp"
00048 
00049 
00061 template<class Scalar>
00062 bool sillyCgSolve(
00063   const Thyra::LinearOpBase<Scalar> &A,
00064   const Thyra::VectorBase<Scalar> &b,
00065   const int maxNumIters,
00066   const typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance,
00067   const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x,
00068   std::ostream &out
00069   )
00070 {
00071 
00072   // Create some typedefs and some other stuff to make the code cleaner
00073   typedef Teuchos::ScalarTraits<Scalar> ST; typedef typename ST::magnitudeType ScalarMag;
00074   const Scalar one = ST::one(), zero = ST::zero();  using Teuchos::as;
00075   using Teuchos::RCP; using Thyra::VectorSpaceBase; using Thyra::VectorBase;
00076   using Thyra::NOTRANS; using Thyra::V_V; using Thyra::apply;
00077   
00078 
00079   // Validate input
00080   THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES("sillyCgSolve()", A, Thyra::NOTRANS, *x, &b);
00081   Teuchos::EVerbosityLevel vl = Teuchos::VERB_MEDIUM;
00082 
00083   out << "\nStarting CG solver ...\n" << std::scientific << "\ndescribe A:\n"<<describe(A, vl)
00084       << "\ndescribe b:\n"<<describe(b, vl)<<"\ndescribe x:\n"<<describe(*x, vl)<<"\n";
00085 
00086   // Initialization
00087   const RCP<const VectorSpaceBase<Scalar> > space = A.domain();
00088   const RCP<VectorBase<Scalar> > r = createMember(space);
00089   // r = -A*x + b
00090   V_V(r.ptr(), b); apply<Scalar>(A, NOTRANS, *x, r.ptr(), -one, one);
00091   const ScalarMag r0_nrm = norm(*r);
00092   if (r0_nrm==zero) return true;
00093   const RCP<VectorBase<Scalar> > p = createMember(space), q = createMember(space);
00094   Scalar rho_old = -one;
00095 
00096   // Perform the iterations
00097   for( int iter = 0; iter <= maxNumIters; ++iter ) {
00098 
00099     // Check convergence and output iteration
00100     const ScalarMag r_nrm = norm(*r);
00101     const bool isConverged = r_nrm/r0_nrm <= tolerance;
00102     if( iter%(maxNumIters/10+1) == 0 || iter == maxNumIters || isConverged ) {
00103       out << "Iter = " << iter << ", ||b-A*x||/||b-A*x0|| = " << (r_nrm/r0_nrm) << std::endl;
00104       if( r_nrm/r0_nrm < tolerance ) return true; // Success!
00105     }
00106 
00107     // Compute iteration
00108     const Scalar rho = inner(*r, *r);        // <r,r>              -> rho
00109     if (iter==0) V_V(p.ptr(), *r);           // r                  -> p   (iter == 0)
00110     else Vp_V( p.ptr(), *r, rho/rho_old );   // r+(rho/rho_old)*p  -> p   (iter  > 0)
00111     apply<Scalar>(A, NOTRANS, *p, q.ptr());  // A*p                -> q
00112     const Scalar alpha = rho/inner(*p, *q);  // rho/<p,q>          -> alpha
00113     Vp_StV( x, +alpha, *p );                 // +alpha*p + x       -> x
00114     Vp_StV( r.ptr(), -alpha, *q );           // -alpha*q + r       -> r
00115     rho_old = rho;                           // rho                -> rho_old (for next iter)
00116 
00117   }
00118 
00119   return false; // Failure
00120 
00121 } // end sillyCgSolve
00122 
00123 #endif // THYRA_SILLY_CG_SOLVE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines