test_block_op.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //           TSFExtended: Trilinos Solver Framework Extended
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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 
00030 
00031 #include "Teuchos_GlobalMPISession.hpp"
00032 #include "Teuchos_DefaultComm.hpp"
00033 #include "Thyra_VectorImpl.hpp"
00034 #include "Thyra_VectorSpaceImpl.hpp"
00035 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00036 #include "Thyra_LinearOperatorImpl.hpp"
00037 #include "Thyra_LinearCombinationImpl.hpp"
00038 #include "Thyra_DefaultBlockedLinearOp.hpp"
00039 #include "Teuchos_CommandLineProcessor.hpp"
00040 #include "Teuchos_ScalarTraits.hpp"
00041 #include "Teuchos_StandardCatchMacros.hpp"
00042 #include "Teuchos_VerboseObject.hpp"
00043 
00044 
00045 using namespace Teuchos;
00046 using namespace Thyra;
00047 
00048 
00049 template <class Scalar> inline
00050 LinearOperator<Scalar> makeRandomDenseOperator(int nc, const VectorSpace<Scalar>& rowSp)
00051 {
00052   typedef typename Teuchos::ScalarTraits<Scalar> ST;
00053   RefCountPtr<Thyra::MultiVectorBase<Scalar> > mv = rowSp.createMembers(nc);
00054   Thyra::randomize(-ST::one(), ST::one(), &*mv);
00055   RefCountPtr<Thyra::LinearOpBase<Scalar> > rtn = mv;
00056   return rtn;
00057 }
00058 
00059 template <class Scalar>
00060 bool runTest(Teuchos::RefCountPtr<Teuchos::FancyOStream>& out);
00061 
00062 int main(int argc, char *argv[]) 
00063 {
00064   bool success = false;
00065   
00066   GlobalMPISession mpiSession(&argc, &argv);
00067   typedef Teuchos::ScalarTraits<double> ST;
00068   
00069   // Get stream that can print to just root or all streams!
00070   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00071     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00072   
00073   try
00074     {
00075       CommandLineProcessor  clp;
00076       clp.throwExceptions(false);
00077       clp.addOutputSetupOptions(true);
00078       bool verbose = false;
00079       clp.setOption( "verbose", "quiet", &verbose, 
00080                      "Determines if any output is printed or not." );
00081 
00082       
00083       CommandLineProcessor::EParseCommandLineReturn parse_return 
00084         = clp.parse(argc,argv);
00085       if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00086 
00087       if (!verbose) out = rcp(new FancyOStream(rcp(new oblackholestream())));
00088 
00089 
00090       success = runTest<double>(out);
00091 
00092       success = runTest<float>(out) && success;
00093 
00094 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00095       success = runTest<std::complex<double> >(out) && success;
00096 #endif
00097 
00098     }
00099 
00100   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,success)
00101 
00102     if (success)
00103       {
00104         *out << "all tests PASSED!" << std::endl;
00105         return 0;
00106       }
00107     else
00108       {
00109         *out << "at least one test FAILED!" << std::endl;
00110         return 1;
00111       }
00112 }
00113 
00114 
00115 
00116 
00117 
00118 
00119 template <class Scalar> inline
00120 bool runTest(Teuchos::RefCountPtr<Teuchos::FancyOStream>& out)
00121 {
00122   typedef typename Teuchos::ScalarTraits<Scalar> ST;
00123   typedef typename ST::magnitudeType ScalarMag;
00124   *out << "==========================================================================="
00125        << endl;
00126   *out << "running block operator test for " << ST::name() << endl;
00127   *out << "==========================================================================="
00128        << endl;
00129   Array<int> rangeSpaceSizes = tuple(5, 10);
00130   Array<int> domainSpaceSizes = tuple(5, 8);
00131   
00132   Array<VectorSpace<Scalar> > rangeBlocks(rangeSpaceSizes.size());
00133   
00134   Array<Array<LinearOperator<Scalar> > > blocks(rangeSpaceSizes.size());
00135   
00136   for (unsigned int br=0; br<rangeSpaceSizes.size(); br++)
00137     {
00138       int n = rangeSpaceSizes[br];
00139       rangeBlocks[br] 
00140         = new DefaultSpmdVectorSpace<Scalar>(DefaultComm<Index>::getComm(),n,-1);
00141       
00142       blocks[br].resize(domainSpaceSizes.size());
00143       
00144       for (unsigned int bc=0; bc<domainSpaceSizes.size(); bc++)
00145         {
00146           blocks[br][bc] = makeRandomDenseOperator<Scalar>(domainSpaceSizes[bc],
00147                                                            rangeBlocks[br]);
00148         }
00149     }
00150   
00151   
00152   LinearOperator<Scalar> A = block2x2(blocks[0][0], blocks[0][1],
00153                                       blocks[1][0], blocks[1][1]);
00154   
00155   VectorSpace<Scalar> domain = A.domain();
00156   VectorSpace<Scalar> range = A.range();
00157   
00158   *out << "A num block rows = " << A.numBlockRows() << endl;
00159   *out << "A num block cols = " << A.numBlockCols() << endl;
00160   
00161   *out << "A domain size = " << domain.dim() << endl;
00162   *out << "A range size = " << range.dim() << endl;
00163   
00164   Vector<Scalar> x = domain.createMember();
00165   *out << "randomizing trial vector" << endl;
00166   Thyra::randomize(-ST::one(), ST::one(), x.ptr().get());
00167   
00168   Array<Vector<Scalar> > xBlock(domain.numBlocks());
00169   for (unsigned int i=0; i<xBlock.size(); i++)
00170     {
00171       xBlock[i] = x.getBlock(i);
00172     }
00173   
00174   *out << "x size = " << space(x).dim() << endl;
00175   
00176   *out << "------------------------------------------------------------" << endl;
00177   *out << "computing A*x..." << endl;
00178   Vector<Scalar> y0;
00179   y0 = A * x;
00180   
00181   
00182   Vector<Scalar> y1 = range.createMember();
00183   *out << "------------------------------------------------------------" << endl;
00184   *out << "computing A*x block-by-block..." << endl;
00185   Array<Vector<Scalar> > yBlock(range.numBlocks());
00186   for (unsigned int i=0; i<yBlock.size(); i++)
00187     {
00188       yBlock[i] = range.getBlock(i).createMember();
00189       zeroOut(yBlock[i]);
00190       for (unsigned int j=0; j<xBlock.size(); j++)
00191         {
00192           LinearOperator<Scalar> Aij = A.getBlock(i,j);
00193           if (Aij.ptr().get()==0) continue;
00194           yBlock[i] = yBlock[i] + Aij * xBlock[j];
00195         }
00196       y1.setBlock(i, yBlock[i]);
00197     }
00198   
00199   ScalarMag err = norm2(y1 - y0);
00200   *out << "error = " << err << endl;
00201   
00202   ScalarMag tol = 1.0e2 * ST::prec();
00203 
00204   bool ok = err < tol;
00205   if (ok)
00206     {
00207       *out << "err=" << err << " tol=" << tol << ", test PASSED!" << endl;
00208     }
00209   else
00210     {
00211       *out << "err=" << err << " tol=" << tol << ", test FAILED!" << endl;
00212     }
00213 
00214   return err < tol;
00215 }
00216 

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1