[Trilinos-Users] Problems with AztecOO-NOX
andrea3.villa at mail.polimi.it
andrea3.villa at mail.polimi.it
Thu Mar 10 14:04:29 MST 2011
I'm experiencing a tough problem using at the same time NOX and
AztecOO. I attach below a small code that reproduces the error. The
code is made of two parts: the part A and the part B. If run alone
each part goes flawlessly and I have also checked with valgrind that
no memory error occurs. However if the two parts are run together a
strange segmentation fault takes place when AztecOO constructs the
preconditioner. In this case the valgrind output suggests a conflict
between some AztecOO functions and the RCP of
NOX::Epetra::LinearSystemAztecOO. I would appreciate if you had any
work arount. I attach below the valgrind output and the code.
Thankyou very much for any suggestion.
Invalid read of size 1
==13508== at 0x5D4651C: std::basic_ostream<char,
std::char_traits<char> >& std::__ostream_insert<char,
std::char_traits<char> >(std::basic_ostream<char,
std::char_traits<char> >&, char const*, long) (in
/usr/lib/libstdc++.so.6.0.13)
==13508== by 0x746A07: AZOO_printf_out(char const*, __va_list_tag*)
(in /home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x75C6C3: AZ_printf_out (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x768230: AZ_iterate (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x73AF03: AztecOO::Iterate(int, double) (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x4095F3: main (in /home/andrea/Documenti/spark2/bin/spark2)
==13508== Address 0x6ad10f1 is 33 bytes inside a block of size 272 free'd
==13508== at 0x4C26DCF: operator delete(void*) (vg_replace_malloc.c:387)
==13508== by 0x82655C: Teuchos::RCPNodeHandle::unbindOne() (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x476349:
NOX::Epetra::LinearSystemAztecOO::~LinearSystemAztecOO() (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x82655C: Teuchos::RCPNodeHandle::unbindOne() (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x40E9D4: solverNL::solve() (in
/home/andrea/Documenti/spark2/bin/spark2)
==13508== by 0x408445: main (in /home/andrea/Documenti/spark2/bin/spark2)
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include "Epetra_ConfigDefs.h"
#include "Epetra_SerialComm.h"
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_FEVector.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_FECrsMatrix.h"
#include "NOX.H"
#include "NOX_Epetra.H"
#include "NOX_Epetra_Interface_Required.H"
#include "NOX_Epetra_Interface_Jacobian.H"
#include "NOX_Epetra_LinearSystem_AztecOO.H"
#include "NOX_Epetra_Group.H"
#include <Teuchos_RCPDecl.hpp>
#include <NOX_Epetra_Interface_Required.H>
#include "AztecOO.h"
using namespace std;
using namespace Teuchos;
class SimpleProblemInterface : public NOX::Epetra::Interface::Required,
public NOX::Epetra::Interface::Jacobian
{
public:
//! Constructor
SimpleProblemInterface() { };
//! Destructor
~SimpleProblemInterface()
{
};
bool computeF(const Epetra_Vector & x, Epetra_Vector & f,
NOX::Epetra::Interface::Required::FillType F )
{
f[0] = x[0] + 1.0;
f[1] = x[1] + 1.0;
return true;
};
bool computeJacobian(const Epetra_Vector & x, Epetra_Operator & Jac)
{
Epetra_CrsMatrix * J;
J = dynamic_cast<Epetra_CrsMatrix*>(&Jac);
if (J == NULL) {
cout << "*ERR* Problem_Interface::computeJacobian() - The
supplied" << endl;
cout << "*ERR* Epetra_Operator is NOT an Epetra_CrsMatrix!" << endl;
throw;
}
std::vector<int> indices(2);
std::vector<double> values(2);
indices[0] = 0;
indices[1] = 1;
// Row 0
values[0] = 1.0;
values[1] = 0.0;
J->ReplaceGlobalValues(0, 2, &values[0], &indices[0]);
// Row 1
values[0] = 0.0;
values[1] = 1.0;
J->ReplaceGlobalValues(1, 2, &values[0], &indices[0]);
return true;
}
};
class solverNL
{
public:
solverNL() { }
void solve()
{
typedef NOX::Epetra::LinearSystemAztecOO LINSOLVER;
typedef NOX::Epetra::Vector NOXVECTOR;
Epetra_SerialComm Comm;
Epetra_Map Map(2,0,Comm);
Epetra_Vector x(Map);
//Jacobian declaration
RCP<Epetra_CrsMatrix> Jac = rcp(new Epetra_CrsMatrix(Copy,Map,1));
std::vector<int> indices(2);
std::vector<double> values(2);
indices[0]=0;
indices[1]=1;
values[0] = 2.0;
values[1] = 2.0;
Jac->InsertGlobalValues(0, 2, &values[0], &indices[0]);
values[0] = - 2.0;
values[1] = 1.0;
Jac->InsertGlobalValues(1, 2, &values[0], &indices[0]);
Jac->FillComplete();
//Solver List
Teuchos::RCP<Teuchos::ParameterList> listaPtr =
Teuchos::rcp(new Teuchos::ParameterList);
Teuchos::ParameterList& lista = *(listaPtr.get());
lista.set("Nonlinear Solver", "Line Search Based");
Teuchos::ParameterList& listaPrint =
lista.sublist("Printing");
listaPrint.set("MyPID", Comm.MyPID());
listaPrint.set("Output Precision", 3);
listaPrint.set("Output Processor", 0);
listaPrint.set("Output Information", NOX::Utils::OuterIteration
+ NOX::Utils::InnerIteration + NOX::Utils::Warning);
Teuchos::ParameterList & listaSearch = lista.sublist("Line Search");
listaSearch.set("Method", "Full Step");
Teuchos::ParameterList & listaDir = lista.sublist("Direction");
listaDir.set("Method", "Newton");
Teuchos::ParameterList & listaNewton = listaDir.sublist("Newton");
listaNewton.set("Forcing Term Method", "Constant");
Teuchos::ParameterList & listaSolver =
listaNewton.sublist("Linear Solver");
listaSolver.set("Aztec Solver", "GMRES");
listaSolver.set("Max Iterations", 1000);
listaSolver.set("Tolerance", 1e-7);
listaSolver.set("Output Frequency", 50);
listaSolver.set("Aztec Preconditioner", "ilu");
//Declaration of the problem
RCP<SimpleProblemInterface> interface = rcp(new SimpleProblemInterface);
//Nonlinear solver solution
RCP<NOX::Epetra::Interface::Required> iReq = interface;
RCP<NOX::Epetra::Interface::Jacobian> iJac = interface;
RCP<LINSOLVER> linSys = Teuchos::rcp(new
LINSOLVER(listaPrint,listaSolver,iReq,iJac,Jac,x));
NOXVECTOR noxInitGuess(x,NOX::DeepCopy);
RCP<NOX::Epetra::Group> gruppo = Teuchos::rcp(new
NOX::Epetra::Group(listaPrint,iReq,noxInitGuess,linSys));
RCP<NOX::StatusTest::NormF> testNormF = Teuchos::rcp(new
NOX::StatusTest::NormF(1e-6));
RCP<NOX::StatusTest::MaxIters> testMaxIters = Teuchos::rcp(new
NOX::StatusTest::MaxIters(100));
RCP<NOX::StatusTest::Combo> combo = Teuchos::rcp(new
NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR,testNormF,
testMaxIters));
Teuchos::RCP<NOX::Solver::Generic> solver =
NOX::Solver::buildSolver(gruppo, combo, listaPtr);
NOX::StatusTest::StatusType status = solver->solve();
const NOX::Epetra::Group & finalGroup = dynamic_cast<const
NOX::Epetra::Group&>(solver->getSolutionGroup());
const Epetra_Vector & finalSolution = (dynamic_cast<const
NOX::Epetra::Vector&>(finalGroup.getX())).getEpetraVector();
}
};
int main(int argc, char *argv[])
{
//--------//
// PART A //
//--------//
//NonLinear
Solver_______________________________________________________________________________
solverNL nonlinearSolution;
nonlinearSolution.solve();
//--------//
// PART B //
//--------//
//PARAMETERS
LIST________________________________________________________________________________
Teuchos::ParameterList List("Solver");
List.set("solver", "AZ_gmres");
List.set("conv", "r0");
List.set("tol", 1e-6);
List.set("max_iter", 1000);
List.set("kspace", 100);
List.set("output", "last");
List.set("precond", "dom_decomp");
List.set("overlap", 0);
List.set("subdomain_solve", "ilut");
List.set("drop", 1e-12);
//Linear
Solver__________________________________________________________________________________
Epetra_SerialComm Comm;
Epetra_Map Map(2,0,Comm);
Epetra_FEVector Y(Map);
Epetra_FEVector T(Map);
Epetra_FECrsMatrix M(Copy,Map,1);
//Matrix assembling
std::vector<int> indices(2);
std::vector<double> values(2);
indices[0] = 0;
indices[1] = 1;
values[0] = 1.0;
values[1] = 0.0;
M.InsertGlobalValues(0, 2, &values[0], &indices[0]);
values[0] = 0.0;
values[1] = 1.0;
M.InsertGlobalValues(1, 2, &values[0], &indices[0]);
M.FillComplete();
//Linear solution
AztecOO linSolver;
double Condest;
linSolver.SetUserMatrix(&M);
linSolver.SetLHS(&Y);
linSolver.SetRHS(&T);
linSolver.SetParameters(List);
cout << "Construct preconditioner" << endl;
linSolver.ConstructPreconditioner(Condest);
cout << "Iterate" << endl;
linSolver.Iterate(1000,1e-6);
cout << "End Iterate" << endl;
}
----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.
More information about the Trilinos-Users
mailing list