Tpetra Matrix/Vector Services Version of the Day
Tpetra_Power_Method_VbrMatrix.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #include <Teuchos_Array.hpp>
00043 #include <Teuchos_ScalarTraits.hpp>
00044 #include <Teuchos_OrdinalTraits.hpp>
00045 #include <Teuchos_RCP.hpp>
00046 #include <Teuchos_GlobalMPISession.hpp>
00047 #include <Teuchos_oblackholestream.hpp>
00048 
00049 #include "Tpetra_DefaultPlatform.hpp"
00050 #include "Tpetra_Version.hpp"
00051 #include "Tpetra_BlockMap.hpp"
00052 #include "Tpetra_BlockMultiVector.hpp"
00053 #include "Tpetra_Vector.hpp"
00054 #include "Tpetra_VbrMatrix.hpp"
00055 
00056 
00057 // prototype
00058 template <class Scalar, class Ordinal>
00059 Scalar power_method(const Teuchos::RCP<const Tpetra::Operator<Scalar,Ordinal> > &A, size_t niters, typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance, bool verbose);
00060 
00061 
00062 int main(int argc, char *argv[]) {
00063   Teuchos::oblackholestream blackhole;
00064   Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);
00065   typedef double Scalar;
00066   typedef int Ordinal;
00067   using Tpetra::global_size_t;
00068 
00069   Teuchos::RCP<const Teuchos::Comm<int> > comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
00070 
00071   size_t myRank = comm->getRank();
00072   size_t numProc = comm->getSize();
00073   bool verbose = (myRank==0);
00074 
00075   if (verbose) {
00076     std::cout << Tpetra::version() << std::endl << std::endl;
00077   }
00078   std::cout << *comm;
00079 
00080   // Get the number of global equations from the command line
00081   if (argc != 2) {
00082     if (verbose) {
00083       std::cout << "Usage: " << argv[0] << " number_of_equations" << std::endl;
00084     }
00085     return -1;
00086   }
00087   const int blockSize = 3;
00088   global_size_t numGlobalElements = std::atoi(argv[1]);
00089 
00090   //make sure numGlobalElements is an integer multiple of blockSize.
00091   numGlobalElements += blockSize - numGlobalElements%blockSize;
00092   global_size_t numGlobalBlocks = numGlobalElements/blockSize;
00093 
00094   if (numGlobalBlocks < numProc) {
00095     if (verbose) {
00096       std::cout << "numGlobalBlocks = " << numGlobalBlocks
00097                 << " cannot be less than the number of processors = " << numProc << std::endl;
00098     }
00099     return -1;
00100   }
00101 
00102   if (verbose) {
00103     std::cout << "numGlobalBlocks = " << numGlobalBlocks << ", numGlobalElements (point-rows): " << numGlobalElements << std::endl;
00104   }
00105 
00106   // Construct a Map that puts approximately the same number of equations on each processor.
00107 
00108   Ordinal indexBase = 0;
00109   Teuchos::RCP<const Tpetra::BlockMap<Ordinal> > blkmap = Teuchos::rcp(new Tpetra::BlockMap<Ordinal>(numGlobalBlocks, blockSize, indexBase, comm));
00110 
00111   // Get update list and number of local equations from newly created map.
00112 
00113   const size_t numMyBlocks = blkmap->getNodeNumBlocks();
00114 
00115   Teuchos::ArrayView<const Ordinal> myGlobalBlocks = blkmap->getNodeBlockIDs();
00116 
00117   // Create an OTeger vector NumNz that is used to build the Petra Matrix.
00118   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
00119   // on this processor
00120 
00121   // Create a Tpetra::VbrMatrix using the BlockMap.
00122   Teuchos::RCP< Tpetra::VbrMatrix<Scalar,Ordinal> > A;
00123   A = Teuchos::rcp( new Tpetra::VbrMatrix<Scalar,Ordinal>(blkmap, 3) );
00124 
00125   // Add  rows one-at-a-time.
00126   // We will build a block-tridiagonal matrix where the diagonal block has 2
00127   // on the diagonal and the off-diagonal blocks have -1 on the diagonal.
00128   Teuchos::Array<Scalar> two(blockSize*blockSize, 0);
00129   two[0] = 2.0; two[blockSize+1] = 2.0; two[blockSize*2+2] = 2.0;
00130   Teuchos::Array<Scalar> negOne(blockSize*blockSize, 0);
00131   negOne[0] = -1.0; negOne[blockSize+1] = -1.0; negOne[blockSize*2+2] = -1.0;
00132 
00133   for (size_t i=0; i<numMyBlocks; i++) {
00134     if (myGlobalBlocks[i] != 0) {
00135       A->setGlobalBlockEntry( myGlobalBlocks[i], myGlobalBlocks[i]-1,
00136                               blockSize, blockSize, blockSize, negOne() );
00137     }
00138 
00139     //always set the diagonal block:
00140     A->setGlobalBlockEntry( myGlobalBlocks[i], myGlobalBlocks[i],
00141                             blockSize, blockSize, blockSize, two() );
00142 
00143     if (static_cast<global_size_t>(myGlobalBlocks[i]) != numGlobalBlocks-1) {
00144       A->setGlobalBlockEntry( myGlobalBlocks[i], myGlobalBlocks[i]+1,
00145                               blockSize, blockSize, blockSize, negOne() );
00146     }
00147   }
00148 
00149   // Finish up matrix fill
00150   A->fillComplete();
00151   if (verbose) std::cout << std::endl << A->description() << std::endl << std::endl;
00152 
00153   // variable needed for iteration
00154   Scalar lambda; (void)lambda;
00155   const size_t niters = static_cast<size_t>(numGlobalElements*10);
00156   const Scalar tolerance = 1.0e-2;
00157 
00158   // Iterate
00159   lambda = power_method<Scalar,Ordinal>(A, niters, tolerance, verbose);
00160 
00161   // Increase diagonal dominance
00162   if (verbose) {
00163     std::cout << "\nIncreasing magnitude of first diagonal term, solving again\n"
00164               << std::endl;
00165   }
00166 
00167   if (A->getBlockRowMap()->getPointMap()->isNodeGlobalElement(0)) {
00168     // modify the diagonal entry of the row with global-id 0
00169     const Ordinal ID = 0;
00170     Ordinal numPtRows, numPtCols;
00171     Teuchos::ArrayRCP<Scalar> blockEntry;
00172     A->getGlobalBlockEntryViewNonConst(ID,ID, numPtRows,numPtCols, blockEntry);
00173     blockEntry[0] *= 10.0;
00174   }
00175 
00176   // Iterate (again)
00177   lambda = power_method<Scalar,Ordinal>(A, niters, tolerance, verbose);
00178 
00179   /* end main
00180    */
00181   if (verbose) {
00182     std::cout << "\nEnd Result: TEST PASSED" << std::endl;
00183   }
00184   return 0;
00185 }
00186 
00187 
00188 template <class Scalar, class Ordinal>
00189 Scalar power_method(const Teuchos::RCP<const Tpetra::Operator<Scalar,Ordinal> > &A, size_t niters, typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance, bool verbose) {
00190   typedef Scalar scalar_type;
00191   typedef Teuchos::ScalarTraits<scalar_type> STS;
00192   typedef typename STS::magnitudeType magnitude_type;
00193   typedef Teuchos::ScalarTraits<magnitude_type> STM;
00194 
00195   const bool NO_INITIALIZE_TO_ZERO = false;
00196   // create three vectors; do not bother initializing q to zero, as we will fill it with random below
00197   Tpetra::Vector<scalar_type, Ordinal> z (A->getRangeMap (), NO_INITIALIZE_TO_ZERO);
00198 
00199   // mfh 26 Nov 2012: For some reason, not initializing r to zero
00200   // causes the norm of r to be NaN in the loop below.  This fixes
00201   // a test failure on my workstation.
00202   Tpetra::Vector<scalar_type, Ordinal> q (A->getRangeMap ()); //, NO_INITIALIZE_TO_ZERO);
00203   Tpetra::Vector<scalar_type, Ordinal> r (A->getRangeMap ()); //, NO_INITIALIZE_TO_ZERO);
00204 
00205   // Fill z with random numbers
00206   z.randomize ();
00207   // Variables needed for iteration
00208   scalar_type lambda = STS::zero ();
00209   magnitude_type normz, residual = STM::zero ();
00210 
00211   // power iteration
00212   for (size_t iter = 0; iter < niters; ++iter) {
00213     normz = z.norm2 ();               // Compute 2-norm of z
00214     q.scale (STS::one () / normz, z); // Compute q = z / normz
00215     A->apply (q, z);                  // Compute z = A*q
00216     lambda = q.dot (z);               // Approximate maximum eigenvalue: lambda = dot(q, z)
00217     if ( iter % 100 == 0 || iter + 1 == niters ) {
00218       r.update (STS::one (), z, -lambda, q, STS::zero ());     // Compute A*q - lambda*q
00219 
00220       magnitude_type r_norm = r.norm2 ();
00221       residual = STS::magnitude (r_norm / lambda);
00222       if (verbose) {
00223         std::cout << "Iter = " << iter << ": lambda = " << lambda
00224                   << ", ||A*q - lambda*q||_2 = " << r_norm
00225                   << ", ||A*q - lambda*q||_2 / lambda = " << residual
00226                   << std::endl;
00227       }
00228     }
00229     if (residual < tolerance) {
00230       break;
00231     }
00232   }
00233   return lambda;
00234 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines