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