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 int main(int argc, char *argv[]) {
00016   Teuchos::oblackholestream blackhole;
00017   Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);
00018 
00019   //
00020   // Specify types used in this example
00021   // 
00022   typedef double Scalar;
00023   typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
00024   typedef int Ordinal;
00025   typedef Tpetra::DefaultPlatform::DefaultPlatformType           Platform;
00026   typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType Node;
00027   using Tpetra::global_size_t;
00028 
00029   // 
00030   // Get the default communicator and node
00031   //
00032   Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
00033   Teuchos::RCP<const Teuchos::Comm<int> > comm = platform.getComm();
00034   Teuchos::RCP<Node>             node = platform.getNode();
00035   const int myRank = comm->getRank();
00036 
00037   //
00038   // Get example parameters from command-line processor
00039   //  
00040   bool printMatrix = false;
00041   bool verbose = (myRank==0);
00042   int niters = 100;
00043   int numGlobalElements = 100;
00044   Magnitude tolerance = 1.0e-2;
00045   Teuchos::CommandLineProcessor cmdp(false,true);
00046   cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
00047   cmdp.setOption("numGlobalElements",&numGlobalElements,"Global problem size.");
00048   cmdp.setOption("tolerance",&tolerance,"Relative residual tolerance used for solver.");
00049   cmdp.setOption("iterations",&niters,"Maximum number of iterations.");
00050   cmdp.setOption("printMatrix","noPrintMatrix",&printMatrix,"Print the full matrix after reading it.");
00051   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00052     return -1;
00053   }
00054 
00055   // 
00056   // Say hello, print some communicator info
00057   //
00058   if (verbose) {
00059     std::cout << "\n" << Tpetra::version() << std::endl << std::endl;
00060   }
00061   std::cout << "Comm info: " << *comm;
00062 
00063   //
00064   // Construct the problem
00065   //
00066   // Construct a Map that puts approximately the same number of equations on each processor.
00067   Teuchos::RCP<const Tpetra::Map<Ordinal> > map = Tpetra::createUniformContigMap<Ordinal,Ordinal>(numGlobalElements, comm);
00068   // Get update list and number of local equations from newly created map.
00069   const size_t numMyElements = map->getNodeNumElements();
00070   Teuchos::ArrayView<const Ordinal> myGlobalElements = map->getNodeElementList();
00071   // Create an OTeger vector NumNz that is used to build the Petra Matrix.
00072   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
00073   // on this processor
00074   Teuchos::ArrayRCP<size_t> NumNz = Teuchos::arcp<size_t>(numMyElements);
00075   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00076   // So we need 2 off-diagonal terms (except for the first and last equation)
00077   for (size_t i=0; i < numMyElements; ++i) {
00078     if (myGlobalElements[i] == 0 || myGlobalElements[i] == numGlobalElements-1) {
00079       // boundary
00080       NumNz[i] = 2;
00081     }
00082     else {
00083       NumNz[i] = 3;
00084     }
00085   }
00086   // Create a Tpetra::Matrix using the Map, with a static allocation dictated by NumNz
00087   Teuchos::RCP< Tpetra::CrsMatrix<Scalar,Ordinal> > A;
00088   A = Teuchos::rcp( new Tpetra::CrsMatrix<Scalar,Ordinal>(map, NumNz, Tpetra::StaticProfile) );
00089   // We are done with NumNZ
00090   NumNz = Teuchos::null;
00091   // Add rows one-at-a-time
00092   // Off diagonal values will always be -1
00093   const Scalar two    = static_cast<Scalar>( 2.0);
00094   const Scalar negOne = static_cast<Scalar>(-1.0);
00095   for (size_t i=0; i<numMyElements; i++) {
00096     if (myGlobalElements[i] == 0) {
00097       A->insertGlobalValues( myGlobalElements[i],
00098                              Teuchos::tuple<Ordinal>( myGlobalElements[i], myGlobalElements[i]+1 ),
00099                              Teuchos::tuple<Scalar> ( two, negOne ) );
00100     }
00101     else if (myGlobalElements[i] == numGlobalElements-1) {
00102       A->insertGlobalValues( myGlobalElements[i],
00103                              Teuchos::tuple<Ordinal>( myGlobalElements[i]-1, myGlobalElements[i] ),
00104                              Teuchos::tuple<Scalar> ( negOne, two ) );
00105     }
00106     else {
00107       A->insertGlobalValues( myGlobalElements[i],
00108                              Teuchos::tuple<Ordinal>( myGlobalElements[i]-1, myGlobalElements[i], myGlobalElements[i]+1 ),
00109                              Teuchos::tuple<Scalar> ( negOne, two, negOne ) );
00110     }
00111   }
00112   // Finish up
00113   A->fillComplete(Tpetra::DoOptimizeStorage);
00114   if (printMatrix) {
00115     Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
00116     A->describe(*fos, Teuchos::VERB_EXTREME);
00117   }
00118   else if (verbose) {
00119     std::cout << std::endl << A->description() << std::endl << std::endl;
00120   }
00121 
00122   //
00123   // Iterate
00124   //
00125   Scalar lambda;
00126   lambda = TpetraExamples::powerMethod<Scalar,Ordinal>(A, niters, tolerance, verbose);
00127 
00128   // Increase diagonal dominance
00129   if (verbose) {
00130     std::cout << "\nIncreasing magnitude of first diagonal term, solving again\n"
00131               << std::endl;
00132   }
00133 
00134   //A->resumeFill();
00135   if (A->getRowMap()->isNodeGlobalElement(0)) {
00136     // get a copy of the row with with global index 0
00137     // modify the diagonal entry of that row
00138     // submit the modified values to the matrix
00139     const Ordinal ID = 0;
00140     size_t numVals = A->getNumEntriesInGlobalRow(ID);
00141     Teuchos::Array<Scalar>  rowvals(numVals);
00142     Teuchos::Array<Ordinal> rowinds(numVals);
00143     A->getGlobalRowCopy(ID, rowinds, rowvals, numVals);       // Get A(0,:)
00144     for (size_t i=0; i<numVals; i++) {
00145       if (rowinds[i] == ID) {
00146         // we have found the diagonal; modify it and break the loop
00147         rowvals[i] *= 10.0;
00148         break;
00149       }
00150     }
00151     A->replaceGlobalValues(ID, rowinds(), rowvals());
00152   }
00153   //A->fillComplete();
00154 
00155   // Iterate (again)
00156   lambda = TpetraExamples::powerMethod<Scalar,Ordinal>(A, niters, tolerance, verbose);  
00157 
00158   if (verbose) {
00159     std::cout << "\nEnd Result: TEST PASSED" << std::endl;
00160   }
00161   return 0;
00162 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines