[Trilinos-Users] amesos + superlu_dist weirdness

[Nico Schlömer]

Hi all,

I'm trying to get Amesos to work with SuperLU_DIST here, and things don't work 
out so well. What I have is a C++ source that contructs a simple test problem, 
and I can solve this with Amesos_Klu and Amesos_Superlu.

When using Amesos_Superludist, though, the solver runs fine, and just returns 
wrong results. :/

Using Superlu_DIST 2.0 here, and observed the same thing for 2.3. This is 
tested with trilinos 9.0.3 as well as the latest pre-10.0 version.

It may very well be that I'm doing something wrong in the code (I appended the 
C++ source), but as mentioned, SuperLU (no dist) and the built-in Klu work 
fine.

Any clues at all? Cheers,
Nico



-------------- next part --------------
#include "Galeri_Maps.h"
#include "Galeri_CrsMatrices.h"
#include "Galeri_Utils.h" // for ComputeNorm

#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#include "mpi.h"
#else
#include "Epetra_SerialComm.h"
#endif

#include "Teuchos_ParameterList.hpp"

#include "Epetra_Vector.h"

#include "Epetra_Time.h"

// Amesos headers
#include "Amesos.h"
#include "Amesos_BaseSolver.h"
#include "Amesos_ConfigDefs.h"


int main(int argv, char* argc[])
{
#ifdef HAVE_MPI
  MPI_Init(&argv, &argc);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  Teuchos::ParameterList GaleriList;
  GaleriList.set("nx", 100);
  GaleriList.set("ny", 100);

  GaleriList.set("mx", Comm.NumProc());
  GaleriList.set("my", 1);

  // create domain map out of this
  Epetra_Map*  Map = 0;
  Map = Galeri::CreateMap("Cartesian2D", Comm, GaleriList);

  // create matrix layout
  Epetra_RowMatrix* Matrix = Galeri::CreateCrsMatrix("Laplace2D", Map, GaleriList);

  Epetra_Vector Sol(*Map,true);    // the solution
  Sol.Random();               // set random

  Epetra_Vector RHS(*Map);        // right-hand side

  Matrix->Multiply(false,Sol,RHS);

  Epetra_Vector U(*Map);  // vector that will contain the solution
  U.PutScalar(0.0);       // set back to 0

  // define problem
  Epetra_LinearProblem Problem(Matrix, &U, &RHS);

  // define solver properties
  Amesos_BaseSolver * Solver;
  Amesos Factory;
//   string SolverType = "Amesos_Klu";
//   string SolverType = "Amesos_Scalapack";
//  string SolverType = "Amesos_Superlu";
  string SolverType = "Amesos_Superludist";
  Solver = Factory.Create(SolverType, Problem);
  assert( Solver!=0 );

  Teuchos::ParameterList List;
  List.set("PrintTiming",true);
  List.set("PrintStatus",true);
//   List.set("MaxProcs",Comm.NumProc());
  Solver->SetParameters(List);

  // setup the timer
  Epetra_Time Time(Comm);

  // -------------------------------------------------------------------------
//   Time.ResetStartTime();
//   AMESOS_CHK_ERR(Solver->SymbolicFactorization());
//     std::cout << std::endl
//           << "symbolic factorization time = " 
//           << Time.ElapsedTime() << std::endl;
//   Comm.Barrier() ; 
// 
//   Time.ResetStartTime();
//   AMESOS_CHK_ERR(Solver->NumericFactorization());
//     std::cout << "numeric factorization time = " 
//           << Time.ElapsedTime() << std::endl;
//   Comm.Barrier() ;

  Time.ResetStartTime();
  AMESOS_CHK_ERR(Solver->Solve());
  if ( Comm.MyPID()==0 ){
    std::cout << "solve time = " 
          << Time.ElapsedTime() << std::endl;
  }
  Comm.Barrier() ;
  // -------------------------------------------------------------------------

  // get the residual norm
  double ResNorm = Galeri::ComputeNorm( Matrix, &U, &RHS);

  // get the error norm
  // -- there might actually be easier ways in trilinos to achieve this..
  double ErrNormComponent = 0.0;
  for (int j = 0 ; j< U.Map().NumMyElements() ; ++j)
        ErrNormComponent += (U[j] - Sol[j]) * (U[j] - Sol[j]);

  // sum it all up
  double ErrNorm = 0.0;
  Comm.SumAll( &ErrNormComponent, &ErrNorm, 1 );

  if ( Comm.MyPID()==0 ){
      std::cout << "||r|| = " << ResNorm << std::endl;
      ErrNorm = sqrt(ErrNorm);
      std::cout << "||e|| = " << ErrNorm << std::endl;
  }

  // clean up
  delete Map;
  delete Matrix;

#ifdef HAVE_MPI
  MPI_Finalize();
#endif
}