[Trilinos-Users] Epetra local/global discrepancy

[Sensei]

Dear all,

I'm playing with Belos to solve a multi-vector problem following 
Trilinos10.12Tutorial.pdf page 73 (8.4 Using the Belos adapter to Epetra).



First of all, Belos works and I owe you guys! :)



Second, the question. Why don't setting values with local IDs work in my 
code?

If I use local IDs to set the matrix values (both InsertMyValues and 
SumIntoMyValues), I get the column values equal to -1, but if I use the 
InsertGlobalValues method then it works just fine.

The code I'm using is this:

     int    col[1];
     double val[1];

     bool local = false;

     for (int i = 0; i <= map->MaxLID(); i++)
     {
         std::cout << rank << " set LID " << i << " GID " << map->GID(i) 
<< std::endl;
         val[0] = 1.0; //(i + 1) * 1.1;

         if (local)
         {
             col[0] = i;
             A->InsertMyValues(i, 1, val, col); // Column == -1
             //A->SumIntoMyValues(i, 1, val, col); // Column == -1
         }
         else
         {
             col[0] = map->GID(i);
             A->InsertGlobalValues(map->GID(i), 1, val, col); // Works
         }
     }

     A->FillComplete();


As you can see, I am using a useless array because I thought the problem 
arised when passing the pointer to index i directly, but it doesn't make 
any difference.

The full code is at the bottom of this email.

I don't understand where I inserted a bug here. The matrix should be an 
identity matrix, but using local==true, being with InsertMyValues or 
SumIntoMyValues, I get -1 as column.

Can you see my error?

Thanks!






#include <iostream>

#include "mpi.h"

#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosEpetraAdapter.hpp"
#include "BelosBlockGmresSolMgr.hpp"

#include "Epetra_Map.h"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_CrsMatrix.h"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_StandardCatchMacros.hpp"

int main(int argc, char * argv[])
{
     int rank = 0;

#ifdef EPETRA_MPI
     // Initialize MPI
     MPI_Init(&argc,&argv);
     Epetra_MpiComm Comm(MPI_COMM_WORLD);
     rank = Comm.MyPID();
#else
     Epetra_SerialComm Comm;
#endif

     typedef double                            ST;
     typedef Teuchos::ScalarTraits<ST>        SCT;
     typedef SCT::magnitudeType                MT;
     typedef Epetra_MultiVector                MV;
     typedef Epetra_Operator                   OP;
     typedef Belos::MultiVecTraits<ST,MV>     MVT;
     typedef Belos::OperatorTraits<ST,MV,OP>  OPT;

     using Teuchos::ParameterList;
     using Teuchos::RCP;
     using Teuchos::rcp;

     RCP<Epetra_Map> map(new Epetra_Map(20, 0, Comm));
     RCP<Epetra_CrsMatrix> A(new Epetra_CrsMatrix(Copy, *map, 1));
     RCP<Epetra_MultiVector> X(new Epetra_MultiVector(*map, 2)), B(new 
Epetra_MultiVector(*map, 2));

     std::cout << "RCP " << rank << " #global " << 
map->NumGlobalElements() << " #local " << map->NumMyElements()
               << "LID " << map->MinLID() << "-" << map->MaxLID() << " "
               << "GID " << map->MinMyGID() << "-" << map->MaxMyGID() << 
std::endl;

     Comm.Barrier();

     int    col[1];
     double val[1];

     bool local = false;

     for (int i = 0; i <= map->MaxLID(); i++)
     {
         std::cout << rank << " set LID " << i << " GID " << map->GID(i) 
<< std::endl;
         val[0] = 1.0; //(i + 1) * 1.1;

         if (local)
         {
             col[0] = i;
             A->InsertMyValues(i, 1, val, col); // Column == -1
             //A->SumIntoMyValues(i, 1, val, col); // Column == -1
         }
         else
         {
             col[0] = map->GID(i);
             A->InsertGlobalValues(map->GID(i), 1, val, col); // Works
         }
     }

     A->FillComplete();

     Comm.Barrier();
     std::cout << "A" << std::endl;
     A->Print(std::cout);
     Comm.Barrier();
     std::cout << "========" << std::endl;
     Comm.Barrier();

     if (rank == 0)
     {
         std::cout << B->ReplaceMyValue(0, 0, 3.1415);
         std::cout << B->ReplaceMyValue(0, 1, -3.1415);
         std::cout << B->ReplaceMyValue(4, 0, 123.456);

         std::cout << B->ReplaceMyValue(9, 0, 3.1415);
     }
     std::cout << std::endl;

     std::cout << "X INITIAL GUESS" << std::endl;
     X->Print(std::cout);

     std::cout << "B RIGHT HAND SIDE" << std::endl;
     B->Print(std::cout);

     Comm.Barrier();
     std::cout << "========" << std::endl;
     Comm.Barrier();

     Teuchos::RCP<Epetra_Operator> op(A);

     std::cout << A.total_count();

     Comm.Barrier();
     std::cout << "========" << std::endl;
     Comm.Barrier();


     Teuchos::RCP<Belos::LinearProblem<ST, MV, OP>> linproblem(new 
Belos::LinearProblem<ST, MV, OP>(A, X, B));

     std::cout << linproblem->setProblem() << std::endl;

     int verb = Belos::Warnings + Belos::Errors + Belos::FinalSummary + 
Belos::TimingDetails;

     Teuchos::RCP<Teuchos::ParameterList> params(new 
Teuchos::ParameterList());

     params->set("Verbosity", verb);
     params->set("Maximum Iterations", 100);
     params->set("Convergence Tolerance", 1.0e-8);

     Teuchos::RCP<Belos::BlockGmresSolMgr<ST, MV, OP>> solver(new 
Belos::BlockGmresSolMgr<ST, MV, OP>(linproblem, params));

     Belos::ReturnType solverRet = solver->solve();

     Comm.Barrier();
     std::cout << "========" << std::endl;
     Comm.Barrier();

     X = linproblem->getLHS();

     std::cout << "X SOLUTION (INITIAL GUESS OVERWRITTEN)" << std::endl;
     X->Print(std::cout);
     Comm.Barrier();

     std::cout << "B RIGHT HAND SIDE" << std::endl;
     B->Print(std::cout);
     Comm.Barrier();

#ifdef EPETRA_MPI
     MPI_Finalize();
#endif
     return 0;
}