[Trilinos-Users] Inverse of Serial Dense Matrix

Nico Schlömer nico.schloemer at gmail.com
Fri Feb 17 08:21:36 MST 2012


Hi Paul,

there are specific classes for it; check out

http://trilinos.sandia.gov/packages/docs/dev/packages/teuchos/doc/html/classTeuchos_1_1SerialDenseSolver.html
http://trilinos.sandia.gov/packages/docs/dev/packages/teuchos/doc/html/classTeuchos_1_1SerialSpdDenseSolver.html

Here's how I used SerialSpdDenseSolver:

Quick & dirty:
================= *snip* =================
  Teuchos::SerialSpdDenseSolver<int,double> solver;
  TEUCHOS_ASSERT_EQUALITY( 0, solver.setMatrix( A ) );
  TEUCHOS_ASSERT_EQUALITY( 0, solver.setVectors( alpha, rhs ) );
  TEUCHOS_ASSERT_EQUALITY( 0, solver.solve() );
================= *snap* =================

Advanced:
================= *snip* =================
  Teuchos::SerialSpdDenseSolver<int,double> solver;
  TEUCHOS_ASSERT_EQUALITY( 0, solver.setMatrix( A ) );
  TEUCHOS_ASSERT_EQUALITY( 0, solver.setVectors( alpha, rhs ) );
  if ( solver.shouldEquilibrate() )
  {
    TEUCHOS_ASSERT_EQUALITY( 0, solver.equilibrateRHS() );
#if TRILINOS_MAJOR_MINOR_VERSION > 100803
    TEUCHOS_ASSERT_EQUALITY( 0, solver.equilibrateMatrix() );
    TEUCHOS_ASSERT_EQUALITY( 0, solver.solve() );
    TEUCHOS_ASSERT_EQUALITY( 0, solver.unequilibrateLHS() );
#else
    // A bug in Trilinos<=10.8.3 makes unequilibrateLHS() fail.
    // Work around by doing the unequilibration manually.
    // Note that this relies on the scaling being
    // s(i) = 1 / sqrt(A(i,i)).
    Teuchos::RCP<DoubleVector> diagA =
      Teuchos::rcp( new DoubleVector( numEdges ) );
    for ( int k=0; k<numEdges; k++ )
      (*diagA)(k) = (*A)(k,k);
    TEUCHOS_ASSERT_EQUALITY( 0, solver.equilibrateMatrix() );
    TEUCHOS_ASSERT_EQUALITY( 0, solver.solve() );
    for ( int k=0; k<numEdges; k++ )
      (*alpha)(k) = (*alpha)(k) / sqrt( (*diagA)(k) );
#endif
  }
  else
  {
    TEUCHOS_ASSERT_EQUALITY( 0, solver.solve() );
  }
================= *snap* =================

If you're running on a recent Trilinos version, you can forget about
the #else part of course.

Cheers,
Nico



On Fri, Feb 17, 2012 at 2:54 PM, Paul Dionne <pjd at cfdrc.com> wrote:
> What is the best method to obtain the inverse of a serial dense matrix?
> I think I'm misunderstanding the ApplyInverse() function.  I also tried
> setting up a SerialDenseSolver problem and calling that invert()
> function, but that left the matrix unchanged.
>
> Paul Dionne
>
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> http://software.sandia.gov/mailman/listinfo/trilinos-users
>




More information about the Trilinos-Users mailing list