Skip to content

Commit

Permalink
Add complex-valued tests for minres, bicgstab, and gcrodr
Browse files Browse the repository at this point in the history
  • Loading branch information
hkthorn committed Aug 7, 2024
1 parent d0b5e32 commit e15a45e
Show file tree
Hide file tree
Showing 11 changed files with 914 additions and 7 deletions.
29 changes: 29 additions & 0 deletions packages/belos/test/BiCGStab/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@


TRIBITS_INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING ${CMAKE_CURRENT_SOURCE_DIR}/../MVOPTester)

ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils)
IF (${PACKAGE_NAME}_ENABLE_Triutils)

ASSERT_DEFINED(Teuchos_ENABLE_COMPLEX)
IF(Teuchos_ENABLE_COMPLEX)

TRIBITS_INCLUDE_DIRECTORIES(../MVOPTester)

TRIBITS_ADD_EXECUTABLE_AND_TEST(
bicgstab_complex_hb
SOURCES test_bicgstab_complex_hb.cpp
ARGS
"--verbose --filename=mhd1280b.cua"
)

ASSERT_DEFINED(Anasazi_SOURCE_DIR)
TRIBITS_COPY_FILES_TO_BINARY_DIR(CopyTestBiCGStabComplexFiles
SOURCE_DIR ${Anasazi_SOURCE_DIR}/testmatrices
SOURCE_FILES mhd1280b.cua
EXEDEPS bicgstab_complex_hb
)

ENDIF(Teuchos_ENABLE_COMPLEX)

ENDIF(${PACKAGE_NAME}_ENABLE_Triutils)
215 changes: 215 additions & 0 deletions packages/belos/test/BiCGStab/test_bicgstab_complex_hb.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
// @HEADER
// *****************************************************************************
// Belos: Block Linear Solvers Package
//
// Copyright 2004-2016 NTESS and the Belos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
//
// This driver reads a problem from a Harwell-Boeing (HB) file.
// The right-hand-side from the HB file is used instead of random vectors.
// The initial guesses are all set to zero.
//
// NOTE: No preconditioner is used in this case.
//
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosBiCGStabSolMgr.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_StandardCatchMacros.hpp"

#ifdef HAVE_MPI
#include <mpi.h>
#endif

// I/O for Harwell-Boeing files
#ifdef HAVE_BELOS_TRIUTILS
#include "Trilinos_Util_iohb.h"
#endif

#include "MyMultiVec.hpp"
#include "MyBetterOperator.hpp"
#include "MyOperator.hpp"

using namespace Teuchos;

int main(int argc, char *argv[]) {
//
typedef std::complex<double> ST;
typedef ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Belos::MultiVec<ST> MV;
typedef Belos::Operator<ST> OP;
typedef Belos::MultiVecTraits<ST,MV> MVT;
typedef Belos::OperatorTraits<ST,MV,OP> OPT;
ST one = SCT::one();
ST zero = SCT::zero();

int info = 0;
int MyPID = 0;
bool norm_failure = false;

Teuchos::GlobalMPISession session(&argc, &argv, NULL);
//
using Teuchos::RCP;
using Teuchos::rcp;

bool success = false;
bool verbose = false;
try {
bool proc_verbose = false;
int frequency = -1; // how often residuals are printed by solver
int blocksize = 1;
int numrhs = 1;
std::string filename("mhd1280b.cua");
MT tol = 1.0e-5; // relative residual tolerance

CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
cmdp.setOption("tol",&tol,"Relative residual tolerance used by BiCGStab solver.");
cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}

proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
if (proc_verbose) {
std::cout << Belos::Belos_Version() << std::endl << std::endl;
}
if (!verbose)
frequency = -1; // reset frequency if test is not verbose


#ifndef HAVE_BELOS_TRIUTILS
std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
if (MyPID==0) {
std::cout << "End Result: TEST FAILED" << std::endl;
}
return -1;
#endif

// Get the data from the HB file
int dim,dim2,nnz;
MT *dvals;
int *colptr,*rowind;
ST *cvals;
nnz = -1;
info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
&colptr,&rowind,&dvals);
if (info == 0 || nnz < 0) {
if (MyPID==0) {
std::cout << "Error reading '" << filename << "'" << std::endl;
std::cout << "End Result: TEST FAILED" << std::endl;
}
return -1;
}
// Convert interleaved doubles to std::complex values
cvals = new ST[nnz];
for (int ii=0; ii<nnz; ii++) {
cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
}
// Build the problem matrix
RCP< MyBetterOperator<ST> > A
= rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
//
// ********Other information used by block solver***********
// *****************(can be user specified)******************
//
int maxits = dim/blocksize; // maximum number of iterations to run
//
ParameterList belosList;
belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
if (verbose) {
belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails );
if (frequency > 0)
belosList.set( "Output Frequency", frequency );
}
else
belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
//
// Construct the right-hand side and solution multivectors.
// NOTE: The right-hand side will be constructed such that the solution is
// a vectors of one.
//
RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
MVT::MvRandom( *soln );
OPT::Apply( *A, *soln, *rhs );
MVT::MvInit( *soln, zero );
//
// Construct an unpreconditioned linear problem instance.
//
RCP<Belos::LinearProblem<ST,MV,OP> > problem =
rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
bool set = problem->setProblem();
if (set == false) {
if (proc_verbose)
std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
return -1;
}
//
// *******************************************************************
// *************Start the BiCGStab iteration***********************
// *******************************************************************
//
Belos::BiCGStabSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );

//
// **********Print out information about problem*******************
//
if (proc_verbose) {
std::cout << std::endl << std::endl;
std::cout << "Dimension of matrix: " << dim << std::endl;
std::cout << "Number of right-hand sides: " << numrhs << std::endl;
std::cout << "Block size used by solver: " << blocksize << std::endl;
std::cout << "Max number of BiCGStab iterations: " << maxits << std::endl;
std::cout << "Relative residual tolerance: " << tol << std::endl;
std::cout << std::endl;
}
//
// Perform solve
//
Belos::ReturnType ret = solver.solve();
//
// Compute actual residuals.
//
RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
OPT::Apply( *A, *soln, *temp );
MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
MVT::MvNorm( *temp, norm_num );
MVT::MvNorm( *rhs, norm_denom );
for (int i=0; i<numrhs; ++i) {
if (proc_verbose)
std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
if ( norm_num[i] / norm_denom[i] > tol ) {
norm_failure = true;
}
}

// Clean up.
delete [] dvals;
delete [] colptr;
delete [] rowind;
delete [] cvals;

success = ret==Belos::Converged && !norm_failure;

if (success) {
if (proc_verbose)
std::cout << "End Result: TEST PASSED" << std::endl;
} else {
if (proc_verbose)
std::cout << "End Result: TEST FAILED" << std::endl;
}
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);

return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
} // end test_minres_complex_hb.cpp
3 changes: 2 additions & 1 deletion packages/belos/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@ ADD_SUBDIRECTORY(Util)
# mfh 10 Jan 2013: MINRES test is taking a long time and causing test timeouts.
# It's not obviously broken, so I'm disabling the test for now until we can find
# a different test problem that doesn't take so long.
#ADD_SUBDIRECTORY(MINRES)
ADD_SUBDIRECTORY(MINRES)
ADD_SUBDIRECTORY(BiCGStab)
ADD_SUBDIRECTORY(MVOPTester)

IF(${PACKAGE_NAME}_ENABLE_Experimental)
Expand Down
4 changes: 2 additions & 2 deletions packages/belos/test/MINRES/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@ IF (${PACKAGE_NAME}_ENABLE_Triutils)
)

ASSERT_DEFINED(Anasazi_SOURCE_DIR)
TRIBITS_COPY_FILES_TO_BINARY_DIR(CopyTestBlockCGComplexFiles
TRIBITS_COPY_FILES_TO_BINARY_DIR(CopyTestMINRESComplexFiles
SOURCE_DIR ${Anasazi_SOURCE_DIR}/testmatrices
SOURCE_FILES mhd1280b.cua
SOURCE_FILES mhd1280b.cua
EXEDEPS minres_complex_hb
)

Expand Down
1 change: 1 addition & 0 deletions packages/belos/test/MINRES/test_minres_complex_hb.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "BelosMinresSolMgr.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_StandardCatchMacros.hpp"

#ifdef HAVE_MPI
#include <mpi.h>
Expand Down
18 changes: 18 additions & 0 deletions packages/belos/tpetra/test/BiCGStab/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,24 @@ IF (${PACKAGE_NAME}_ENABLE_Triutils)
STANDARD_PASS_OUTPUT
)

IF (Tpetra_INST_COMPLEX_DOUBLE)

TRIBITS_ADD_EXECUTABLE_AND_TEST(
Tpetra_bicgstab_complex_hb
SOURCES test_bicgstab_complex_hb.cpp
ARGS
"--verbose --filename=mhd1280b.cua"
)

ASSERT_DEFINED(Anasazi_SOURCE_DIR)
TRIBITS_COPY_FILES_TO_BINARY_DIR(CopyTpetraBiCGStabComplexFiles
SOURCE_DIR ${Anasazi_SOURCE_DIR}/testmatrices
SOURCE_FILES mhd1280b.cua
EXEDEPS Tpetra_bicgstab_complex_hb
)

ENDIF()

TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestBiCGSTABFiles
SOURCE_DIR ${${PACKAGE_NAME}_SOURCE_DIR}/epetra/test/BiCGStab
SOURCE_FILES cage4.hb
Expand Down
Loading

0 comments on commit e15a45e

Please sign in to comment.