From c12666e22321accda08da0fc8821300edfbf6782 Mon Sep 17 00:00:00 2001 From: Caleb Schilly Date: Tue, 26 Sep 2023 16:44:25 -0400 Subject: [PATCH] #215: add cxx_main_lap.cpp --- .../test/GeneralizedDavidson/CMakeLists.txt | 7 + .../test/GeneralizedDavidson/cxx_main.cpp | 4 +- .../test/GeneralizedDavidson/cxx_main_lap.cpp | 257 ++++++++++++++++++ 3 files changed, 267 insertions(+), 1 deletion(-) create mode 100644 packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main_lap.cpp diff --git a/packages/anasazi/tpetra/test/GeneralizedDavidson/CMakeLists.txt b/packages/anasazi/tpetra/test/GeneralizedDavidson/CMakeLists.txt index 0835f99a8add..d7836169afbd 100644 --- a/packages/anasazi/tpetra/test/GeneralizedDavidson/CMakeLists.txt +++ b/packages/anasazi/tpetra/test/GeneralizedDavidson/CMakeLists.txt @@ -1,4 +1,11 @@ +TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_GeneralizedDavidson_Lap_test + SOURCES cxx_main_lap.cpp + ARGS + COMM serial mpi + ) + ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) IF (${PACKAGE_NAME}_ENABLE_Triutils) TRIBITS_ADD_EXECUTABLE_AND_TEST( diff --git a/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main.cpp b/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main.cpp index 779c6f60e0aa..4da53661263a 100644 --- a/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main.cpp +++ b/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main.cpp @@ -101,7 +101,7 @@ int run(int argc, char *argv[]) { const ST ONE = SCT::one(); int info = 0; - const auto comm = Tpetra::getDefaultComm (); + RCP > comm = Tpetra::getDefaultComm (); const int MyPID = comm->getRank(); bool testFailed; @@ -179,7 +179,9 @@ int run(int argc, char *argv[]) { free( colptr ); free( rowind ); } + std::cout << "Rank " << MyPID << " before fillComplete(): " << std::endl; K->fillComplete(); + std::cout << "Rank " << MyPID << " after fillComplete(): " << std::endl; // Create initial vectors RCP ivec = rcp( new MV(map,blockSize) ); diff --git a/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main_lap.cpp b/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main_lap.cpp new file mode 100644 index 000000000000..bce5540191ca --- /dev/null +++ b/packages/anasazi/tpetra/test/GeneralizedDavidson/cxx_main_lap.cpp @@ -0,0 +1,257 @@ +// @HEADER +// *********************************************************************** +// +// Anasazi: Block Eigensolvers Package +// Copyright 2004 Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Michael A. Heroux (maherou@sandia.gov) +// +// *********************************************************************** +// @HEADER +// +// This test is for GeneralizedDavidson solving a standard (Ax=xl) Hermitian +// eigenvalue problem where the operator (A) is the 1D finite-differenced Laplacian +// operator. + +#include "AnasaziConfigDefs.hpp" +#include "AnasaziTypes.hpp" + +#include "AnasaziTpetraAdapter.hpp" +#include "AnasaziBasicEigenproblem.hpp" +#include "AnasaziGeneralizedDavidsonSolMgr.hpp" +#include + +#include +#include + + +template +int run(int argc, char *argv[]) { + + using ST = typename Tpetra::MultiVector::scalar_type; + using LO = typename Tpetra::MultiVector<>::local_ordinal_type; + using GO = typename Tpetra::MultiVector<>::global_ordinal_type; + using NT = typename Tpetra::MultiVector<>::node_type; + + using SCT = typename Teuchos::ScalarTraits; + using MT = typename SCT::magnitudeType; + + using OP = Tpetra::Operator; + using MV = Tpetra::MultiVector; + using OPT = Anasazi::OperatorTraits; + using MVT = Anasazi::MultiVecTraits; + + using tmap_t = Tpetra::Map; + using tcrsmatrix_t = Tpetra::CrsMatrix; + + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::tuple; + using std::cout; + using std::endl; + + const ST ONE = SCT::one(); + + Tpetra::ScopeGuard tpetraScope (&argc,&argv); + RCP > comm = Tpetra::getDefaultComm (); + + const int MyPID = comm->getRank (); + const int NumImages = comm->getSize (); + + bool testFailed; + bool verbose = false; + bool debug = false; + std::string which("LM"); + int nev = 5; + int blockSize = 3; + MT tol = 1.0e-6; + int maxRestarts = 25; + int maxDim = 50; + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); + cmdp.setOption("debug","nodebug",&debug,"Print debugging information."); + cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM)."); + cmdp.setOption("nev",&nev,"Number of eigenvalues to compute."); + cmdp.setOption("blockSize",&blockSize,"Block size for the algorithm."); + cmdp.setOption("maxRestarts",&maxRestarts,"Number of restarts allowed."); + + if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (debug) verbose = true; + + if (MyPID == 0) { + cout << Anasazi::Anasazi_Version() << endl << endl; + } + + // -- Set finite difference grid + const int nx = 10; + int dim = nx * nx; + + // create map + RCP map = rcp (new tmap_t(dim,0,comm)); + RCP K = rcp (new tcrsmatrix_t(map, 4)); + int base = MyPID*nx; + if (MyPID != NumImages-1) { + for (int i=0; iinsertGlobalValues(static_cast(base+i ), tuple(base+i ), tuple( 2)); + K->insertGlobalValues(static_cast(base+i ), tuple(base+i+1), tuple(-1)); + K->insertGlobalValues(static_cast(base+i+1), tuple(base+i ), tuple(-1)); + K->insertGlobalValues(static_cast(base+i+1), tuple(base+i+1), tuple( 2)); + } + } + else { + for (int i=0; iinsertGlobalValues(static_cast(base+i ), tuple(base+i ), tuple( 2)); + K->insertGlobalValues(static_cast(base+i ), tuple(base+i+1), tuple(-1)); + K->insertGlobalValues(static_cast(base+i+1), tuple(base+i ), tuple(-1)); + K->insertGlobalValues(static_cast(base+i+1), tuple(base+i+1), tuple( 2)); + } + } + K->fillComplete(); + + // Create initial vectors + RCP ivec = rcp (new MV (map,blockSize)); + ivec->randomize (); + + // Create eigenproblem + RCP > problem = + rcp (new Anasazi::BasicEigenproblem (K, ivec)); + // + // Inform the eigenproblem that the operator K is symmetric + problem->setHermitian (true); + // + // Set the number of eigenvalues requested + problem->setNEV (nev); + // + // Inform the eigenproblem that you are done passing it information + bool boolret = problem->setProblem (); + if (! boolret) { + if (MyPID == 0) { + cout << "Anasazi::BasicEigenproblem::SetProblem() returned with error." << endl + << "End Result: TEST FAILED" << endl; + } + return -1; + } + + // Set verbosity level + int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails; + if (verbose) { + verbosity += Anasazi::IterationDetails; + } + if (debug) { + verbosity += Anasazi::Debug; + } + + // Eigensolver parameters + // + // Create parameter list to pass into the solver manager + Teuchos::ParameterList MyPL; + MyPL.set( "Verbosity", verbosity ); + MyPL.set( "Which", which ); + MyPL.set( "Maximum Subspace Dimension", maxDim ); + MyPL.set( "Block Size", blockSize ); + MyPL.set( "Maximum Restarts", maxRestarts ); + MyPL.set( "Convergence Tolerance", tol ); + // + // Create the solver manager + Anasazi::GeneralizedDavidsonSolMgr MySolverMgr(problem, MyPL); + + // Solve the problem to the specified tolerances or length + Anasazi::ReturnType returnCode = MySolverMgr.solve(); + testFailed = false; + if (returnCode != Anasazi::Converged) { + testFailed = true; + } + + // Get the eigenvalues and eigenvectors from the eigenproblem + Anasazi::Eigensolution sol = problem->getSolution(); + RCP evecs = sol.Evecs; + int numev = sol.numVecs; + + if (numev > 0) { + std::ostringstream os; + os.setf(std::ios::scientific, std::ios::floatfield); + os.precision(6); + + // Compute the direct residual + std::vector normV( numev ); + Teuchos::SerialDenseMatrix T (numev, numev); + for (int i = 0; i < numev; ++i) { + T(i,i) = sol.Evals[i].realpart; + } + RCP Kvecs = MVT::Clone( *evecs, numev ); + + OPT::Apply( *K, *evecs, *Kvecs ); + + MVT::MvTimesMatAddMv( -ONE, *evecs, T, ONE, *Kvecs ); + MVT::MvNorm( *Kvecs, normV ); + + os << "Direct residual norms computed in Tpetra_GeneralizedDavidson_lap_test.exe" << endl + << std::setw(20) << "Eigenvalue" << std::setw(20) << "Residual " << endl + << "----------------------------------------" << endl; + for (int i=0; i tol ) { + testFailed = true; + } + } + if (MyPID==0) { + cout << endl << os.str() << endl; + } + } + + if (testFailed) { + if (MyPID==0) { + cout << "End Result: TEST FAILED" << endl; + } + return -1; + } + // + // Default return value + // + if (MyPID==0) { + cout << "End Result: TEST PASSED" << endl; + } + return 0; + +} + +int main(int argc, char *argv[]) { + return run(argc,argv); + // run(argc,argv); +}