Tpetra Matrix/Vector Services Version of the Day
Tpetra_Power_Method.cpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 //
00005 //          Tpetra: Templated Linear Algebra Services Package
00006 //                 Copyright (2008) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 // @HEADER
00042 */
00043 
00044 #include <Teuchos_GlobalMPISession.hpp>
00045 #include <Teuchos_oblackholestream.hpp>
00046 #include <Teuchos_CommandLineProcessor.hpp>
00047 #include <Teuchos_Array.hpp>
00048 
00049 #include "Tpetra_Power_Method.hpp"
00050 
00051 #include "Tpetra_DefaultPlatform.hpp"
00052 #include "Tpetra_Version.hpp"
00053 #include "Tpetra_Map.hpp"
00054 #include "Tpetra_MultiVector.hpp"
00055 #include "Tpetra_Vector.hpp"
00056 #include "Tpetra_CrsMatrix.hpp"
00057 
00058 #include <algorithm>
00059 #include <functional>
00060 
00061 int main(int argc, char *argv[]) {
00062   Teuchos::oblackholestream blackhole;
00063   Teuchos::GlobalMPISession mpiSession (&argc, &argv, &blackhole);
00064 
00065   //
00066   // Specify types used in this example
00067   //
00068   typedef double                                                  Scalar;
00069   typedef Teuchos::ScalarTraits<Scalar>::magnitudeType            Magnitude;
00070   typedef int                                                     Ordinal;
00071   typedef Tpetra::DefaultPlatform::DefaultPlatformType            Platform;
00072   typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType  Node;
00073   typedef Tpetra::Map<Ordinal,Ordinal,Node>                       Map;
00074   typedef Tpetra::CrsMatrix<Scalar,Ordinal,Ordinal,Node>          CrsMatrix;
00075   using Teuchos::Array;
00076   using Teuchos::ArrayView;
00077   using Teuchos::as;
00078   using Teuchos::getFancyOStream;
00079   using Teuchos::RCP;
00080   using Teuchos::rcpFromRef;
00081   using Teuchos::tuple;
00082   using std::cout;
00083   using std::endl;
00084   typedef Array<Scalar>::size_type size_type;
00085 
00086   //
00087   // Get the default communicator and node
00088   //
00089   Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
00090   RCP<const Teuchos::Comm<int> > comm = platform.getComm();
00091   RCP<Node>                      node = platform.getNode();
00092   const int myRank = comm->getRank();
00093 
00094   //
00095   // Get example parameters from command-line processor
00096   //
00097   bool printMatrix = false;
00098   bool verbose = (myRank==0);
00099   int niters = 100;
00100   int numGlobalElements = 100;
00101   Magnitude tolerance = 1.0e-2;
00102   Teuchos::CommandLineProcessor cmdp (false, true);
00103   cmdp.setOption ("verbose", "quiet", &verbose, "Print messages and results.");
00104   cmdp.setOption ("numGlobalElements", &numGlobalElements,
00105                   "Global problem size.");
00106   cmdp.setOption ("tolerance", &tolerance,
00107                   "Relative residual tolerance used for solver.");
00108   cmdp.setOption ("iterations", &niters, "Maximum number of iterations.");
00109   cmdp.setOption ("printMatrix", "noPrintMatrix", &printMatrix,
00110                   "Print the full matrix after reading it.");
00111   if (cmdp.parse (argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00112     return -1;
00113   }
00114 
00115   //
00116   // Say hello, print some communicator info
00117   //
00118   if (verbose) {
00119     cout << endl << Tpetra::version() << endl << endl;
00120   }
00121   cout << "Comm info: " << *comm;
00122 
00123   //
00124   // Construct the problem
00125   //
00126   // Construct a Map that puts approximately the same number of equations on each processor.
00127   RCP<const Map> map = Tpetra::createUniformContigMap<Ordinal,Ordinal> (numGlobalElements, comm);
00128   // Get update list and number of local equations from newly created map.
00129   const size_type numMyElements = as<size_type> (map->getNodeNumElements ());
00130   ArrayView<const Ordinal> myGlobalElements = map->getNodeElementList ();
00131   // Create a CrsMatrix using the map, with a dynamic allocation of 3 entries per row.
00132   RCP<CrsMatrix> A = Tpetra::createCrsMatrix<Scalar> (map, 3);
00133   // Add rows one at a time.
00134   for (size_type i = 0; i < numMyElements; ++i) {
00135     // Teuchos::tuple<T> constructs an in-line fixed-length array with
00136     // entries of type T.  It's convenient for examples like this.
00137     if (myGlobalElements[i] == 0) {
00138       A->insertGlobalValues (myGlobalElements[i],
00139                              tuple<Ordinal> (myGlobalElements[i], myGlobalElements[i]+1),
00140                              tuple<Scalar> (2.0, -1.0));
00141     }
00142     else if (myGlobalElements[i] == numGlobalElements-1) {
00143       A->insertGlobalValues (myGlobalElements[i],
00144                              tuple<Ordinal> (myGlobalElements[i]-1, myGlobalElements[i]),
00145                              tuple<Scalar> (-1.0, 2.0));
00146     }
00147     else {
00148       A->insertGlobalValues (myGlobalElements[i],
00149                              tuple<Ordinal> (myGlobalElements[i]-1,
00150                                              myGlobalElements[i],
00151                                              myGlobalElements[i]+1),
00152                              tuple<Scalar> (-1.0, 2.0, -1.0));
00153     }
00154   }
00155   // Complete the fill.  Note that if the row Map, domain Map, and
00156   // range Map are not all the same, you will need to supply the
00157   // domain and range Maps to fillComplete().
00158   A->fillComplete ();
00159   if (printMatrix) {
00160     RCP<Teuchos::FancyOStream> fos = fancyOStream (rcpFromRef (cout));
00161     A->describe (*fos, Teuchos::VERB_EXTREME);
00162   }
00163   else if (verbose) {
00164     cout << endl << A->description() << endl << endl;
00165   }
00166 
00167   //
00168   // Iterate
00169   //
00170   TpetraExamples::powerMethod<Scalar,Ordinal> (A, niters, tolerance, verbose);
00171 
00172   // Increase diagonal dominance
00173   if (verbose) {
00174     cout << endl << "Increasing magnitude of first diagonal term, solving again"
00175          << endl << endl;
00176   }
00177 
00178   A->resumeFill();
00179 
00180   // Get a copy of the row with with global index 0.  Modify the
00181   // diagonal entry of that row, and submit the modified values to the
00182   // matrix.  This code assumes that row 0 is uniquely owned.
00183   if (A->getRowMap ()->isNodeGlobalElement (0)) {
00184     // Row 0 belongs to the calling process.  Get the number of
00185     // entries in row 0 on my process.  Teuchos::as casts between
00186     // types like static_cast, but also checks in a debug build
00187     // whether the cast will overflow.
00188     size_t numVals = A->getNumEntriesInGlobalRow (0);
00189     Array<Scalar>  rowVals (as<size_type> (numVals));
00190     Array<Ordinal> rowInds (as<size_type> (numVals));
00191     // The parentheses indicate that we are passing in a view of the
00192     // array.  The callee can't modify the length of the view.  If we
00193     // didn't supply enough space, the method will throw an exception.
00194     // That's why we ask above for the number of entries in the row.
00195     // numVals on output is set to the number of entries in the row,
00196     // which is why it was declared nonconst above.
00197     A->getGlobalRowCopy (0, rowInds (), rowVals (), numVals);
00198     // Replace the diagonal entry with 10.0.
00199     for (size_type k = 0; k < as<size_type> (numVals); ++k) {
00200       if (rowInds[k] == 0) {
00201         rowVals[k] = as<Scalar> (10.0);
00202       }
00203     }
00204     A->replaceGlobalValues (0, rowInds (), rowVals ());
00205   }
00206   // This call assumes that the row Map, domain Map, and range Map are
00207   // all the same.  If this is not the case, you should supply the
00208   // domain Map and range Map to fillComplete().
00209   A->fillComplete ();
00210 
00211   // Iterate again
00212   TpetraExamples::powerMethod<Scalar,Ordinal> (A, niters, tolerance, verbose);
00213 
00214   if (verbose) {
00215     cout << endl << "End Result: TEST PASSED" << endl;
00216   }
00217   return 0;
00218 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines