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