Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BlockKrylovSchur/BlockKrylovSchurEpetraExSVD.cpp

Use Anasazi::BlockKrylovSchurSolMgr to compute a singular value decomposition (SVD), using Epetra data structures.

// @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 example shows how to use the block Krylov-Schur method to compute a few
// of the largest singular values (sigma) and corresponding right singular
// vectors (v) for the matrix A by solving the symmetric problem:
//
// (A'*A)*v = sigma*v
//
// where A is an m by n real matrix that is derived from the simplest finite
// difference discretization of the 2-dimensional kernel K(s,t)dt where
//
// K(s,t) = s(t-1) if 0 <= s <= t <= 1
// t(s-1) if 0 <= t <= s <= 1
//
// NOTE: This example came from the ARPACK SVD driver dsvd.f
//
#include "Epetra_CrsMatrix.h"
#include "Teuchos_StandardCatchMacros.hpp"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#include <mpi.h>
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
int main(int argc, char *argv[]) {
#ifdef EPETRA_MPI
// Initialize MPI
MPI_Init(&argc,&argv);
#endif
bool success = false;
try {
#ifdef EPETRA_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
#endif
int i, j, info;
const double one = 1.0;
const double zero = 0.0;
int MyPID = Comm.MyPID();
// Dimension of the matrix
int m = 500;
int n = 100;
// Construct a Map that puts approximately the same number of
// equations on each processor.
Epetra_Map RowMap(m, 0, Comm);
Epetra_Map ColMap(n, 0, Comm);
// Get update list and number of local equations from newly created Map.
int NumMyRowElements = RowMap.NumMyElements();
std::vector<int> MyGlobalRowElements(NumMyRowElements);
RowMap.MyGlobalElements(&MyGlobalRowElements[0]);
/* We are building an m by n matrix with entries
A(i,j) = k*(si)*(tj - 1) if i <= j
= k*(tj)*(si - 1) if i > j
where si = i/(m+1) and tj = j/(n+1) and k = 1/(n+1).
*/
// Create an Epetra_Matrix
Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy, RowMap, n) );
// Compute coefficients for discrete integral operator
std::vector<double> Values(n);
std::vector<int> Indices(n);
double inv_mp1 = one/(m+1);
double inv_np1 = one/(n+1);
for (i=0; i<n; i++) { Indices[i] = i; }
for (i=0; i<NumMyRowElements; i++) {
//
for (j=0; j<n; j++) {
//
if ( MyGlobalRowElements[i] <= j ) {
Values[j] = inv_np1 * ( (MyGlobalRowElements[i]+one)*inv_mp1 ) * ( (j+one)*inv_np1 - one ); // k*(si)*(tj-1)
}
else {
Values[j] = inv_np1 * ( (j+one)*inv_np1 ) * ( (MyGlobalRowElements[i]+one)*inv_mp1 - one ); // k*(tj)*(si-1)
}
}
info = A->InsertGlobalValues(MyGlobalRowElements[i], n, &Values[0], &Indices[0]);
assert( info==0 );
}
// Finish up
info = A->FillComplete(ColMap, RowMap);
assert( info==0 );
info = A->OptimizeStorage();
assert( info==0 );
A->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
//************************************
// Start the block Arnoldi iteration
//***********************************
//
// Variables used for the Block Arnoldi Method
//
int nev = 4;
int blockSize = 1;
int numBlocks = 10;
int maxRestarts = 20;
double tol = lapack.LAMCH('E');
std::string which = "LM";
//
// Create parameter list to pass into solver
//
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Which", which );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Num Blocks", numBlocks );
MyPL.set( "Maximum Restarts", maxRestarts );
MyPL.set( "Convergence Tolerance", tol );
// Create an Anasazi::EpetraMultiVec for an initial vector to start the solver.
// Note: This needs to have the same number of columns as the blocksize.
ivec->MvRandom();
// Call the constructor for the (A^T*A) operator
// Inform the eigenproblem that the matrix A is symmetric
MyProblem->setHermitian(true);
// Set the number of eigenvalues requested and the blocksize the solver should use
MyProblem->setNEV( nev );
// Inform the eigenproblem that you are finished passing it information
bool boolret = MyProblem->setProblem();
if (!boolret) {
throw "Anasazi::BasicEigenproblem::setProblem() returned with error.";
}
// Initialize the Block Arnoldi solver
// Solve the problem to the specified tolerances or length
Anasazi::ReturnType returnCode = MySolverMgr.solve();
if (returnCode != Anasazi::Converged && MyPID==0) {
std::cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
// Get the eigenvalues and eigenvectors from the eigenproblem
Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution();
std::vector<Anasazi::Value<double> > evals = sol.Evals;
int numev = sol.numVecs;
if (numev > 0) {
// Compute singular values/vectors and direct residuals.
//
// Compute singular values which are the square root of the eigenvalues
if (MyPID==0) {
std::cout<<"------------------------------------------------------"<<std::endl;
std::cout<<"Computed Singular Values: "<<std::endl;
std::cout<<"------------------------------------------------------"<<std::endl;
}
for (i=0; i<numev; i++) { evals[i].realpart = Teuchos::ScalarTraits<double>::squareroot( evals[i].realpart ); }
//
// Compute left singular vectors : u = Av/sigma
//
std::vector<double> tempnrm(numev), directnrm(numev);
std::vector<int> index(numev);
for (i=0; i<numev; i++) { index[i] = i; }
Anasazi::EpetraMultiVec Av(RowMap,numev), u(RowMap,numev);
Anasazi::EpetraMultiVec* evecs = dynamic_cast<Anasazi::EpetraMultiVec* >(sol.Evecs->CloneViewNonConst( index ));
A->Apply( *evecs, Av );
Av.MvNorm( tempnrm );
for (i=0; i<numev; i++) { S(i,i) = one/tempnrm[i]; };
u.MvTimesMatAddMv( one, Av, S, zero );
//
// Compute direct residuals : || Av - sigma*u ||
//
for (i=0; i<numev; i++) { S(i,i) = evals[i].realpart; }
Av.MvTimesMatAddMv( -one, u, S, one );
Av.MvNorm( directnrm );
if (MyPID==0) {
std::cout.setf(std::ios_base::right, std::ios_base::adjustfield);
std::cout<<std::setw(16)<<"Singular Value"
<<std::setw(20)<<"Direct Residual"
<<std::endl;
std::cout<<"------------------------------------------------------"<<std::endl;
for (i=0; i<numev; i++) {
std::cout<<std::setw(16)<<evals[i].realpart
<<std::setw(20)<<directnrm[i]
<<std::endl;
}
std::cout<<"------------------------------------------------------"<<std::endl;
}
}
success = true;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
#ifdef EPETRA_MPI
MPI_Finalize() ;
#endif
return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
}