Combining Trilinos Packages To Solve Linear Systems

1. Combining Trilinos PackagesTo Solve Linear Systems Michael A. Heroux Sandia National Labs

2. Epetra User Class Categories • Sparse Matrices: RowMatrix, (CrsMatrix, VbrMatrix, FECrsMatrix, FEVbrMatrix) • Linear Operator: Operator: (AztecOO, ML, Ifpack) • Dense Matrices: DenseMatrix, DenseVector, BLAS, LAPACK, SerialDenseSolver • Vectors: Vector, MultiVector • Graphs: CrsGraph • Data Layout: Map, BlockMap, LocalMap • Redistribution: Import, Export, LbGraph, LbMatrix • Aggregates: LinearProblem • Parallel Machine: Comm, (SerialComm, MpiComm, MpiSmpComm) • Utilities: Time, Flops

3. Matrix Class Inheritance Details • CrsMatrix and VbrMatrix inherit from: • Distributed Object: How data is spread across machine. • Computational Object: Performs FLOPS. • BLAS: Use BLAS kernels. • RowMatrix: An object from either class has a common row access interface (used by AztecOO).

4. LinearProblem Class • A linear problem is defined by: • Matrix A : An Epetra_RowMatrix or Epetra_Operator object. (often a CrsMatrix or VbrMatrix object.) • Vectors x, b : Vector objects. • To call AztecOO, first define a LinearProblem: • Constructed from A, x and b. • Once defined, can: • Scale the problem (explicit preconditioning). • Precondition it (implicitly). • Change x and b.

5. AztecOO • Aztec is the workhorse solver at Sandia: • Extracted from the MPSalsa reacting flow code. • Installed in dozens of Sandia apps. • 800+ external licenses. • AztecOO leverages the investment in Aztec: • Uses Aztec iterative methods and preconditioners. • AztecOO improves on Aztec by: • Using Epetra objects for defining matrix and RHS. • Providing more preconditioners/scalings. • Using C++ class design to enable more sophisticated use. • AztecOO interfaces allows: • Continued use of Aztec for functionality. • Introduction of new solver capabilities outside of Aztec.

6. AztecOO Extensibility • AztecOO is designed to accept externally defined: • Operators (both A and M): • The linear operator A is accessed as an Epetra_Operator. • Users can register a preconstructed preconditioner as an Epetra_Operator. • RowMatrix: • If A is registered as a RowMatrix, Aztec’s preconditioners are accessible. • Alternatively M can be registered separately as an Epetra_RowMatrix, and Aztec’s preconditioners are accessible. • StatusTests: • Aztec’s standard stopping criteria are accessible. • Can override these mechanisms by registering a StatusTest Object.

7. AztecOO understands Epetra_Operator • AztecOO is designed to accept externally defined: • Operators (both A and M). • RowMatrix (Facilitates use of AztecOO preconditioners with external A). • StatusTests (externally-defined stopping criteria). Epetra_Operator Methods Documentation

8. A Simple Epetra/AztecOO Program // Header files omitted… int main(int argc, char *argv[]) { MPI_Init(&argc,&argv); // Initialize MPI, MpiComm Epetra_MpiComm Comm( MPI_COMM_WORLD ); // ***** Create x and b vectors ***** Epetra_Vector x(Map); Epetra_Vector b(Map); b.Random(); // Fill RHS with random #s // ***** Map puts same number of equations on each pe ***** int NumMyElements = 1000 ; Epetra_Map Map(-1, NumMyElements, 0, Comm); int NumGlobalElements = Map.NumGlobalElements(); // ***** Create Linear Problem ***** Epetra_LinearProblem problem(&A, &x, &b); // ***** Create/define AztecOO instance, solve ***** AztecOO solver(problem); solver.SetAztecOption(AZ_precond, AZ_Jacobi); solver.Iterate(1000, 1.0E-8); // ***** Create an Epetra_Matrix tridiag(-1,2,-1) ***** Epetra_CrsMatrix A(Copy, Map, 3); double negOne = -1.0; double posTwo = 2.0; for (int i=0; i<NumMyElements; i++) { int GlobalRow = A.GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1; if (RowLess1!=-1) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1); if (RowPlus1!=NumGlobalElements) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1); A.InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow); } A.TransformToLocal(); // Transform from GIDs to LIDs // ***** Report results, finish *********************** cout << "Solver performed " << solver.NumIters() << " iterations." << endl << "Norm of true residual = " << solver.TrueResidual() << endl; MPI_Finalize() ; return 0; }

9. AztecOO with MLTrilinos/packages/aztecoo/example/MLAztecOO // Header files omitted… int main(int argc, char *argv[]) { MPI_Init(&argc,&argv); // Initialize MPI, MpiComm Epetra_MpiComm Comm( MPI_COMM_WORLD ); // ***** Create x and b vectors ***** Epetra_Vector x(Map); Epetra_Vector b(Map); b.Random(); // Fill RHS with random #s // ***** Map puts same number of equations on each pe ***** int NumMyElements = 1000 ; Epetra_Map Map(-1, NumMyElements, 0, Comm); int NumGlobalElements = Map.NumGlobalElements(); // ***** Create Linear Problem ***** Epetra_LinearProblem problem(&A, &x, &b); // ***** Create/define AztecOO instance, solve ***** AztecOO solver(problem); // Code from next slide gets inserted here. solver.Iterate(1000, 1.0E-8); // ***** Create an Epetra_Matrix tridiag(-1,2,-1) ***** Epetra_CrsMatrix A(Copy, Map, 3); double negOne = -1.0; double posTwo = 2.0; for (int i=0; i<NumMyElements; i++) { int GlobalRow = A.GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1; if (RowLess1!=-1) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1); if (RowPlus1!=NumGlobalElements) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1); A.InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow); } A.TransformToLocal(); // Transform from GIDs to LIDs // ***** Report results, finish *********************** cout << "Solver performed " << solver.NumIters() << " iterations." << endl << "Norm of true residual = " << solver.TrueResidual() << endl; MPI_Finalize() ; return 0; }

10. ML Preconditioner Setup cout << "Creating multigrid hierarchy" << endl; ML *ml_handle; int N_levels = 10; ML_Set_PrintLevel(3); ML_Create(&ml_handle,N_levels); EpetraMatrix2MLMatrix(ml_handle, 0, &A); ML_Aggregate *agg_object; ML_Aggregate_Create(&agg_object); ML_Aggregate_Set_MaxCoarseSize(agg_object,1); N_levels = ML_Gen_MGHierarchy_UsingAggregation(ml_handle, 0, ML_INCREASING, agg_object); // Set a symmetric Gauss-Seidel smoother for the MG method ML_Gen_Smoother_SymGaussSeidel(ml_handle, ML_ALL_LEVELS, ML_BOTH, 1, ML_DEFAULT); ML_Gen_Solver (ml_handle, ML_MGV, 0, N_levels-1); cout << "Creating Epetra_ML_Operator wrapper" << endl; Epetra_ML_Operator MLop(ml_handle,comm,map,map); //note: SetPrecOperator() also sets AZ_user_precond solver.SetPrecOperator(&MLop);

11. AztecOO UserOp/UserMat Recursive Call ExampleTrilinos/packages/aztecoo/example/AztecOO_RecursiveCall • Poisson2dOperator A(nx, ny, comm); // Generate nx by ny Poisson operator • Epetra_CrsMatrix * precMatrix = A.GeneratePrecMatrix(); // Build tridiagonal approximate Poisson • Epetra_Vector xx(A.OperatorDomainMap()); // Generate vectors (xx will be used to generate RHS b) • Epetra_Vector x(A.OperatorDomainMap()); • Epetra_Vector b(A.OperatorRangeMap()); • xx.Random(); // Generate exact x and then rhs b • A.Apply(xx, b); • // Build AztecOO solver that will be used as a preconditioner • Epetra_LinearProblem precProblem; • precProblem.SetOperator(precMatrix); • AztecOO precSolver(precProblem); • precSolver.SetAztecOption(AZ_precond, AZ_ls); • precSolver.SetAztecOption(AZ_output, AZ_none); • precSolver.SetAztecOption(AZ_solver, AZ_cg); • AztecOO_Operator precOperator(&precSolver, 20); • Epetra_LinearProblem problem(&A, &x, &b); // Construct linear problem • AztecOO solver(problem); // Construct solver • solver.SetPrecOperator(&precOperator); // Register Preconditioner operator • solver.SetAztecOption(AZ_solver, AZ_cg); • solver.Iterate(Niters, 1.0E-12);

12. AztecOO Understands Epetra_RowMatrix Epetra_RowMatrix Methods

13. Ifpack/AztecOO Example Trilinos/packages/aztecoo/example/IfpackAztecOO • // Assume A, x, b are define, LevelFill and Overlap are specified • Ifpack_IlukGraph IlukGraph(A.Graph(), LevelFill, Overlap); • IlukGraph.ConstructFilledGraph(); • Ifpack_CrsRiluk ILUK (IlukGraph); • ILUK.InitValues(A); • assert(ILUK->Factor()==0); // Note: All Epetra/Ifpack/AztecOO method return int err codes • double Condest; • ILUK.Condest(false, Condest); // Get condition estimate • if (Condest > tooBig) { • ILUK.SetAbsoluteThreshold(Athresh); • ILUK.SetRelativeThreshold(Rthresh); • Go back to line 4 and try again • } • Epetra_LinearProblem problem(&A, &x, &b); // Construct linear problem • AztecOO solver(problem); // Construct solver • solver.SetPrecOperator(&ILUK); // Register Preconditioner operator • solver.SetAztecOption(AZ_solver, AZ_cg); • solver.Iterate(Niters, 1.0E-12); • // Once this linear solutions complete and the next nonlinear step is advanced, • // we will return to the solver, but only need to execute steps 4 on down…

14. Multiple Stopping Criteria • Possible scenario for stopping an iterative solver: • Test 1: Make sure residual is decreased by 6 orders of magnitude. And • Test 2: Make sure that the inf-norm of true residual is no more 1.0E-8. But • Test 3: do no more than 200 iterations. • Note: Test 1 is cheap. Do it before Test 2.

15. AztecOO StatusTest classes • AztecOO_StatusTest: • Abstract base class for defining stopping criteria. • Combo class: OR, AND, SEQ AztecOO_StatusTest Methods

16. AztecOO/StatusTest Example Trilinos/packages/aztecoo/example/AztecOO • // Assume A, x, b are define • Epetra_LinearProblem problem(&A, &x, &b); // Construct linear problem • AztecOO solver(problem); // Construct solver • AztecOO_StatusTestResNorm restest1(A, x, bb, 1.0E-6); • restest1.DefineResForm(AztecOO_StatusTestResNorm::Implicit, AztecOO_StatusTestResNorm::TwoNorm); • restest1.DefineScaleForm(AztecOO_StatusTestResNorm::NormOfInitRes, AztecOO_StatusTestResNorm::TwoNorm); • AztecOO_StatusTestResNorm restest2(A, x, bb, 1.0E-8); • restest2.DefineResForm(AztecOO_StatusTestResNorm::Explicit, AztecOO_StatusTestResNorm::InfNorm); • restest2.DefineScaleForm(AztecOO_StatusTestResNorm::NormOfRHS, AztecOO_StatusTestResNorm::InfNorm); • AztecOO_StatusTestCombo comboTest1(AztecOO_StatusTestCombo::SEQ, restest1, restest2); • AztecOO_StatusTestMaxIters maxItersTest(200); • AztecOO_StatusTestCombo comboTest2(AztecOO_StatusTestCombo::OR, maxItersTest1, comboTest1); • solver.SetStatusTest(&comboTest2); • solver.SetAztecOption(AZ_solver, AZ_cg); • solver.Iterate(Niters, 1.0E-12);

17. Summary • Trilinos packages are designed to interoperate. • Next generation linear solver packages are under development: Tpetra, TSF packages, Belos, Ifpack 3.0 • AztecOO is the production “manager” class. • Flexibility comes from abstract base classes: • Epetra_Operator: • All Epetra matrix classes implement. • Best way to define A and M when coefficient info not needed. • Epetra_RowMatrix: • All Epetra matrix classes implement. • Best way to define A and M when coefficient info is needed. • AztecOO_StatusTest: • A suite of parametrized status tests. • An abstract interface for users to define their own. • Ability to combine tests for sophisticated control of stopping.