Tpetra Matrix/Vector Services Version of the Day
Tpetra_Power_Method.cpp
00001 #include <Teuchos_GlobalMPISession.hpp>
00002 #include <Teuchos_oblackholestream.hpp>
00003 #include <Teuchos_CommandLineProcessor.hpp>
00004 #include <Teuchos_Array.hpp>
00005 
00006 #include "Tpetra_Power_Method.hpp"
00007 
00008 #include "Tpetra_DefaultPlatform.hpp"
00009 #include "Tpetra_Version.hpp"
00010 #include "Tpetra_Map.hpp"
00011 #include "Tpetra_MultiVector.hpp"
00012 #include "Tpetra_Vector.hpp"
00013 #include "Tpetra_CrsMatrix.hpp"
00014 
00015 #include <algorithm>
00016 #include <functional>
00017 
00018 int main(int argc, char *argv[]) {
00019   Teuchos::oblackholestream blackhole;
00020   Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);
00021 
00022   //
00023   // Specify types used in this example
00024   // 
00025   typedef double                                                  Scalar;
00026   typedef Teuchos::ScalarTraits<Scalar>::magnitudeType            Magnitude;
00027   typedef int                                                     Ordinal;
00028   typedef Tpetra::DefaultPlatform::DefaultPlatformType            Platform;
00029   typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType  Node;
00030   typedef Tpetra::Map<Ordinal,Ordinal,Node>                       Map;
00031   typedef Tpetra::CrsMatrix<Scalar,Ordinal,Ordinal,Node>          CrsMatrix;
00032   using Teuchos::RCP;
00033   using Teuchos::tuple;
00034 
00035   // 
00036   // Get the default communicator and node
00037   //
00038   Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
00039   RCP<const Teuchos::Comm<int> > comm = platform.getComm();
00040   RCP<Node>                      node = platform.getNode();
00041   const int myRank = comm->getRank();
00042 
00043   //
00044   // Get example parameters from command-line processor
00045   //  
00046   bool printMatrix = false;
00047   bool verbose = (myRank==0);
00048   int niters = 100;
00049   int numGlobalElements = 100;
00050   Magnitude tolerance = 1.0e-2;
00051   Teuchos::CommandLineProcessor cmdp(false,true);
00052   cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
00053   cmdp.setOption("numGlobalElements",&numGlobalElements,"Global problem size.");
00054   cmdp.setOption("tolerance",&tolerance,"Relative residual tolerance used for solver.");
00055   cmdp.setOption("iterations",&niters,"Maximum number of iterations.");
00056   cmdp.setOption("printMatrix","noPrintMatrix",&printMatrix,"Print the full matrix after reading it.");
00057   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00058     return -1;
00059   }
00060 
00061   // 
00062   // Say hello, print some communicator info
00063   //
00064   if (verbose) {
00065     std::cout << "\n" << Tpetra::version() << std::endl << std::endl;
00066   }
00067   std::cout << "Comm info: " << *comm;
00068 
00069   //
00070   // Construct the problem
00071   //
00072   // Construct a Map that puts approximately the same number of equations on each processor.
00073   RCP<const Map> map = Tpetra::createUniformContigMap<Ordinal,Ordinal>(numGlobalElements, comm);
00074   // Get update list and number of local equations from newly created map.
00075   const size_t numMyElements = map->getNodeNumElements();
00076   Teuchos::ArrayView<const Ordinal> myGlobalElements = map->getNodeElementList();
00077   // Create a CrsMatrix using the map, with a dynamic allocation of 3 entries per row
00078   RCP<CrsMatrix> A = Tpetra::createCrsMatrix<Scalar>(map,3);
00079   // Add rows one-at-a-time
00080   for (size_t i=0; i<numMyElements; i++) {
00081     if (myGlobalElements[i] == 0) {
00082       A->insertGlobalValues( myGlobalElements[i],
00083                              tuple<Ordinal>( myGlobalElements[i], myGlobalElements[i]+1 ),
00084                              tuple<Scalar> ( 2.0, -1.0 ) );
00085     }
00086     else if (myGlobalElements[i] == numGlobalElements-1) {
00087       A->insertGlobalValues( myGlobalElements[i],
00088                              tuple<Ordinal>( myGlobalElements[i]-1, myGlobalElements[i] ),
00089                              tuple<Scalar> ( -1.0, 2.0 ) );
00090     }
00091     else {
00092       A->insertGlobalValues( myGlobalElements[i],
00093                              tuple<Ordinal>( myGlobalElements[i]-1, myGlobalElements[i], myGlobalElements[i]+1 ),
00094                              tuple<Scalar> ( -1.0, 2.0, -1.0 ) );
00095     }
00096   }
00097   // Complete the fill, ask that storage be reallocated and optimized
00098   A->fillComplete(Tpetra::DoOptimizeStorage);
00099   if (printMatrix) {
00100     RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
00101     A->describe(*fos, Teuchos::VERB_EXTREME);
00102   }
00103   else if (verbose) {
00104     std::cout << std::endl << A->description() << std::endl << std::endl;
00105   }
00106 
00107   //
00108   // Iterate
00109   //
00110   TpetraExamples::powerMethod<Scalar,Ordinal>(A, niters, tolerance, verbose);
00111 
00112   // Increase diagonal dominance
00113   if (verbose) {
00114     std::cout << "\nIncreasing magnitude of first diagonal term, solving again\n"
00115               << std::endl;
00116   }
00117 
00118   A->resumeFill();
00119   if (A->getRowMap()->isNodeGlobalElement(0)) {
00120     // get a copy of the row with with global index 0
00121     // modify the diagonal entry of that row
00122     // submit the modified values to the matrix
00123     size_t numVals = A->getNumEntriesInGlobalRow(0);
00124     Teuchos::Array<Scalar>  rowvals(numVals);
00125     Teuchos::Array<Ordinal> rowinds(numVals);
00126     A->getGlobalRowCopy(0, rowinds, rowvals, numVals);
00127     // replace the diagonal with 10.0
00128     std::replace_if( rowvals.begin(), rowvals.end(), 
00129                      std::bind2nd(std::equal_to<Ordinal>(), 0), 
00130                      10.0 
00131                    );
00132     A->replaceGlobalValues(0, rowinds(), rowvals());
00133   }
00134   A->fillComplete();
00135 
00136   // Iterate again
00137   TpetraExamples::powerMethod<Scalar,Ordinal>(A, niters, tolerance, verbose);  
00138 
00139   if (verbose) {
00140     std::cout << "\nEnd Result: TEST PASSED" << std::endl;
00141   }
00142   return 0;
00143 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines