[happy_clp]

Recently, I utilize PETSc to solve Ax=b,but A is pretty ill-conditioned, all kinds of iterative methods were diverged, Only LU has good performance.But LU will encounter memory bottleneck. Because PETSc's precontioner is very poor when I use iterative methods, I want to use trilinos package!Do you think this is a right way?Thank you a lot! Or use hypre + petsc, but I don't konw how call pc of hypre in petsc. I have not these knowledge and experience.Whether I need to modify my original program? As for trilinos, I cannot download it from web,because the web is very very slow.I am puzzled!

[jed]

Originally Posted by happy_clp

Recently, I utilize PETSc to solve Ax=b,but A is pretty ill-conditioned, all kinds of iterative methods were diverged, Only LU has good performance.But LU will encounter memory bottleneck. Because PETSc's precontioner is very poor when I use iterative methods, I want to use trilinos package!Do you think this is a right way?Thank you a lot! Or use hypre + petsc, but I don't konw how call pc of hypre in petsc. I have not these knowledge and experience.Whether I need to modify my original program? As for trilinos, I cannot download it from web,because the web is very very slow.I am puzzled!

There are no good black-box preconditioners for indefinite problems. That is pretty fundamental and using Trilinos will not fix that. In my experience, it is easier to create first-class matrix and preconditioner objects using PETSc. PETSc provides an interface to more third-party preconditioners. The best preconditioning package available in Trilinos is ML which is also available from PETSc.

Trilinos is relatively weaker on flexible object creation. That is, you end up having to do a lot of work to avoid hard-coding the types. Since you want to be able to try out many combinations, especially when using multilevel and/or domain decomposition preconditioners, this is pretty serious. PETSc is more dynamic in this regard, you should never need to determine the matrix type, preconditioner, etc at compile time.

[jugghead]

Is there any commercial alternative to petsc?

[maherou]

Originally Posted by jed

There are no good black-box preconditioners for indefinite problems. That is pretty fundamental and using Trilinos will not fix that. In my experience, it is easier to create first-class matrix and preconditioner objects using PETSc. PETSc provides an interface to more third-party preconditioners. The best preconditioning package available in Trilinos is ML which is also available from PETSc.

Trilinos is relatively weaker on flexible object creation. That is, you end up having to do a lot of work to avoid hard-coding the types. Since you want to be able to try out many combinations, especially when using multilevel and/or domain decomposition preconditioners, this is pretty serious. PETSc is more dynamic in this regard, you should never need to determine the matrix type, preconditioner, etc at compile time.

Although it is true that ML is available via PETSc, the more advanced features are not accessible, and the version of ML in PETSc is not the latest. The latest Trilinos release (9.0.x) supports PETSc data structures, allowing you to use PETSc matrices and vectors with Trilinos, including ML.

I am not sure what is meant by the inflexible object creation comment. There are several native data classes that are only restrictive for the purposes of performance. The FECrsMatrix classes allows any processor to register any data needed by any other processor. Trilinos ultimately provides complete flexibility because all of its solvers and preconditioners access matrix operations and data via abstract interfaces. If you don't like the native data classes, you can use PETSc's or roll your own. You can also easily compose advanced preconditioners using any combination of other preconditioners, something that is essential for advanced multiphysics. In addition, the next version of Trilinos will have a new data package called Tpetra that will support templated data types, allowing mixing of float and double, and complex types all at the same time with a simple parameter redefinition.

Finally, the right package in Trilinos to use for accessing solvers and preconditioners is Stratimikos. It provides a common access to all packages driven by a parameter list class that can be constructed in code or read from an XML input file.

PETSc is a great package, so I am not discouraging people from using it. I just don't want misconceptions to exist about Trilinos.

[jed]

Originally Posted by maherou

Although it is true that ML is available via PETSc, the more advanced features are not accessible, and the version of ML in PETSc is not the latest. The latest Trilinos release (9.0.x) supports PETSc data structures, allowing you to use PETSc matrices and vectors with Trilinos, including ML.

PETSc has supported version of ML in Trilinos 9 (ML-6.2) since a couple weeks after the release. You are right that not all features of ML are available, but it is more flexible in some ways. In particular, the interface actually only uses ML to produce the restriction/interpolation operators. The rest of the multigrid mechanics employs PETSc's multigrid framework so you can control every detail of smoother and solve on every level from the command line (or from your code if you prefer).

Originally Posted by maherou

I am not sure what is meant by the inflexible object creation comment. There are several native data classes that are only restrictive for the purposes of performance. The FECrsMatrix classes allows any processor to register any data needed by any other processor. Trilinos ultimately provides complete flexibility because all of its solvers and preconditioners access matrix operations and data via abstract interfaces. If you don't like the native data classes, you can use PETSc's or roll your own. You can also easily compose advanced preconditioners using any combination of other preconditioners, something that is essential for advanced multiphysics. In addition, the next version of Trilinos will have a new data package called Tpetra that will support templated data types, allowing mixing of float and double, and complex types all at the same time with a simple parameter redefinition.

Finally, the right package in Trilinos to use for accessing solvers and preconditioners is Stratimikos. It provides a common access to all packages driven by a parameter list class that can be constructed in code or read from an XML input file.

I was not aware of Stratimikos, it looks like it offers a rough equivalent to the "flexible object creation" I was referring to (at the expense of yet another package...) and gives similar functionality to the PETSc options database. Depending on precisely what is meant by "composition" of preconditioners, this is either available on the command line or by writing shell preconditioners. Yes, being able to write such preconditioners is essential for multiphysics, but I haven't seen Trilinos code that makes this easier than it is in PETSc.

The functionality you describe in FECrsMatrix is present in all PETSc parallel matrix formats. Likewise for data structure neutrality. Incidentally, I think that C++ inheritance is a bad model for matrices since there are so many overlapping sets of possible functionality. You either end up forcing people to implement pure virtual functions that would be unreasonably expensive (so they just error) or you give everything a default implementation. What ends up happeningg is an attempt to partition the interface into orthogonal subsets, but this either means massively multiple inheritance or a class hiererchy in which the users' matrix doesn't fit. With OO in C using pimpl (PETSc's model), you can have a flat hierarchy and just check whether function pointers are NULL if you want to find out whether a given type implements an interface.

Being able to try fairly exotic preconditioning strategies on the command line is a strength of PETSc that I haven't seen available from Trilinos. Experimentation is essential and the overhead of recompiling or even editing an XML file is significant.

Tpetra has been around for a long time so I'm skeptical that it will be usable from all of Trilinos at the next release. It will offer a significant advantage for people who need complex-valued matrices for a small part of a multiphysics problem where the bulk of the computation is real-valued (if the complex-valued part makes up most of the computation, there is little hit to just using complex-valued matrices throughout). PETSc has had support for lower precision matrices (usually used for the preconditioning matrix with a matrix-free Jacobian) for ages, but you cannot use real- and complex-valued matrices in the same code.

Originally Posted by maherou

PETSc is a great package, so I am not discouraging people from using it. I just don't want misconceptions to exist about Trilinos.

Trilinos is also quite good, but I have frequently run into the attitude that "it's C++ therefore more OO/flexible/better than C" (maybe it's just particularly strong at my university). On multiple occasions, I have talked with people at conferences who are using Trilinos (or even Matlab) and been told that they "are working on implementing this feature" which turns out to be possible on the command line with PETSc. I'm sure that a sufficiently experienced Trilinos user could have implemented these features in very few lines of code, so don't take this anecdote too seriously.

[maherou]

Originally Posted by jed

I was not aware of Stratimikos, it looks like it offers a rough equivalent to the "flexible object creation" I was referring to (at the expense of yet another package...)

The nice thing about this "yet another package" is that it encapsulates a bunch of others

Originally Posted by jed

The functionality you describe in FECrsMatrix is present in all PETSc parallel matrix formats. Likewise for data structure neutrality.

Outside of the FECrsMatrix, Trilinos specifically prohibits implicit data communication. This is a philosophy that some find restrictive, but it does reduce the likelihood of performance problems.

Originally Posted by jed

Incidentally, I think that C++ inheritance is a bad model for matrices since there are so many overlapping sets of possible functionality. You either end up forcing people to implement pure virtual functions that would be unreasonably expensive (so they just error) or you give everything a default implementation. What ends up happeningg is an attempt to partition the interface into orthogonal subsets, but this either means massively multiple inheritance or a class hiererchy in which the users' matrix doesn't fit. With OO in C using pimpl (PETSc's model), you can have a flat hierarchy and just check whether function pointers are NULL if you want to find out whether a given type implements an interface.

Use of C++ by Trilinos has never been second-guessed. I have never experienced the problems you mention above. In fact, it is very nice to have language support for OO programming.

Originally Posted by jed

Being able to try fairly exotic preconditioning strategies on the command line is a strength of PETSc that I haven't seen available from Trilinos. Experimentation is essential and the overhead of recompiling or even editing an XML file is significant.

Trilinos also has a command line parser, which I agree is convenient to have.

Originally Posted by jed

Tpetra has been around for a long time so I'm skeptical that it will be usable from all of Trilinos at the next release.

Tpetra is for real now. It is intended to be the workhorse for Trilinos starting with the 10.0 release in September. In addition to being able to support all common precisions, it will allow mixed-precision computations and extended precision. This will be very useful in many settings.

[jed]

Originally Posted by maherou

Use of C++ by Trilinos has never been second-guessed. I have never experienced the problems you mention above. In fact, it is very nice to have language support for OO programming.

C++ is certainly the natural choice and I don't blame you for choosing it. I have encountered matrix hierarchy issues when dealing with matrices such as low-rank corrections, Schur complements, and matrices with block tensor-product structure. Some such matrices have various types of incomplete solves or may even be able to provide certain entries (a diagonal, for instance) at modest costs, but certainly can't provide all entries.

A specific example in Trilinos is Epetra_Operator::ApplyInverse(). A user must implement this function, but it doesn't make sense for many matrix types (in principle, the operation exists but it's not affordable) so have to implement it and throw an error. But now, a caller cannot determine if ApplyInverse is implemented without calling it. The same applies for applying a transpose. You can add flags to determine whether the matrix has this operation, but that's working outside the type system and someone is responsible for keeping it consistent. You could split these operations into separate base classes Epetra_InvertableOperator, Epetra_TransposableOperator, multiply inherit from these, and dynamic_cast when you need that operation. This might be the most "OO" operation, but there are reasons everyone doesn't do it.

The choice in Trilinos is to have a small number of matrix interfaces and require the user to provide dummy implementations of all the methods that don't make sense for their specific implementation. If a method is added to the interface, the user's code won't compile even though that method may not make any sense for them. The decision to structure the hierarchy this way is certainly reasonable, but I wouldn't claim that the type system is helping or that it's nicer to work with than a v-table supporting reflection.

Originally Posted by maherou

Trilinos also has a command line parser, which I agree is convenient to have.

Yes, but this is quite different from assigning prefixes to objects so that you can change their type and behavior from the command line without the user having to write code. As an example, consider

mpirun -n 80 ./app -snes_mf_operator -snes_ksp_ew -pc_type mg -pc_mg_type kaskade -mg_levels_2_pc_type hypre -mg_levels_2_pc_hypre_bj TRUE -mg_levels_2_pc_hypre_type euclid -mg_coarse_pc_redundant_number 10 -mg_coarse_redundant_pc_factor_mat_solver_package mumps

This does Jacobian-free Newton-Krylov with linear tolerances adjusted according to the Eisenstat-Walker heuristic, using a multigrid preconditioner. The preconditioner is using a third-party smoother on level 2 and solving the coarse grid component semi-redundantly (10 groups of 8 processors) using a third-party parallel direct solver. To have these options available takes no user code. My understanding is that the user would have to write a bit of code to support this sort of command-line manipulation using Trilinos. If one is using Stratimikos, can this be done from XML (without the user writing any extra code)? If so, it should be possible to adopt a name-mangling convention and make it available on the command line as well.

Originally Posted by maherou

Tpetra is for real now. It is intended to be the workhorse for Trilinos starting with the 10.0 release in September. In addition to being able to support all common precisions, it will allow mixed-precision computations and extended precision. This will be very useful in many settings.

Excellent, I look forward to checking it out. Will AztecOO still be the workhorse, or is Belos "coming of age"?

[maherou]

Originally Posted by jed

A specific example in Trilinos is Epetra_Operator::ApplyInverse(). A user must implement this function, but it doesn't make sense for many matrix types (in principle, the operation exists but it's not affordable) so have to implement it and throw an error. But now, a caller cannot determine if ApplyInverse is implemented without calling it. The same applies for applying a transpose. You can add flags to determine whether the matrix has this operation, but that's working outside the type system and someone is responsible for keeping it consistent. You could split these operations into separate base classes Epetra_InvertableOperator, Epetra_TransposableOperator, multiply inherit from these, and dynamic_cast when you need that operation. This might be the most "OO" operation, but there are reasons everyone doesn't do it.

In most realistic implementations of Epetra_Operator, the Apply and ApplyInverse methods are implemented, the former as the standard "matrix" operator and the latter as the preconditioner. ApplyInverse is the method used by the iterative solvers for preconditioning purposes. It is very convenient to have both methods in the same class since the data of the matrix and the preconditioner often have common elements. If you don't want to implement both methods, having a no-op Apply or ApplyInverse is not much different that a null pointer that is used by PETSc. You just use the C++ try-and-catch mechanism.

Originally Posted by jed

Yes, but this is quite different from assigning prefixes to objects so that you can change their type and behavior from the command line without the user having to write code. As an example, consider

Stratimikos can provide this similar polymorphic behavior from the command line, with a simple driver program. But Trilinos does not put as much emphasis on command line parameterization since the primary focus is integration into large-scale applications where the ability to read an XML file is more effective.

Originally Posted by jed

Excellent, I look forward to checking it out. Will AztecOO still be the workhorse, or is Belos "coming of age"?

Belos is coming of age and is now the recommended solver for new usage of Trilinos. It is compatible already with Tpetra and available via Stratimikos. Belos contains a variety of block iterative methods and recycling methods in addition to the standard CG and GMRES. It is especially useful for optimization, sensitivity analysis, stochastic simulations, etc. (which are also part of Trilinos now), where multiple RHS are available.

Thanks for the good exchange.

[jed]

Originally Posted by maherou

In most realistic implementations of Epetra_Operator, the Apply and ApplyInverse methods are implemented, the former as the standard "matrix" operator and the latter as the preconditioner. ApplyInverse is the method used by the iterative solvers for preconditioning purposes. It is very convenient to have both methods in the same class since the data of the matrix and the preconditioner often have common elements. If you don't want to implement both methods, having a no-op Apply or ApplyInverse is not much different that a null pointer that is used by PETSc. You just use the C++ try-and-catch mechanism.

Yes, the relationship between matrix and preconditioner is tricky, there is always a trade-off of how much to put in the matrix object and how much to put in the preconditioner object. I find it clumsy to be required to provide no-op implementations of all the functions that don't make sense, consider implementing Epetra_RowMatrix when all you really want to implement is ExtractDiagonalCopy(). In addition, the C++ try-and-catch is only equivalent if you never want to check whether that method is implemented without calling it. The PETSc approach is more dynamic in this regard, and I prefer it, but there are good reasons for the Trilinos design too.

Originally Posted by maherou

Stratimikos can provide this similar polymorphic behavior from the command line, with a simple driver program. But Trilinos does not put as much emphasis on command line parameterization since the primary focus is integration into large-scale applications where the ability to read an XML file is more effective.

It's good to hear this sort of runtime polymorphism is being made easier. I understand the rationale for XML, but I still like my command line

Originally Posted by maherou

Belos is coming of age and is now the recommended solver for new usage of Trilinos. It is compatible already with Tpetra and available via Stratimikos. Belos contains a variety of block iterative methods and recycling methods in addition to the standard CG and GMRES. It is especially useful for optimization, sensitivity analysis, stochastic simulations, etc. (which are also part of Trilinos now), where multiple RHS are available.

Cool stuff.

Originally Posted by maherou

Thanks for the good exchange.

Thanks you, it's been informative.