]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add UMFPACK. Also add SuperLU for the moment, though I will remove this again since...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 20 Aug 2004 20:46:18 +0000 (20:46 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 20 Aug 2004 20:46:18 +0000 (20:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@9566 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_direct.h
deal.II/lac/source/sparse_direct.cc

index dcd119d01d54762b34984778adf1a60b0a0f91b9..1e56ead71c5e3031354d2f28d4393638b6615580 100644 (file)
  * indicating error codes, etc. It also manages allocation of the
  * right amount of temporary storage required by these functions.
  *
- * For a description of the steps necessary for the installation of
- * HSL subroutines, read the section on external libraries in the
- * <tt>deal.II</tt> ReadMe file.
+ * Note that this class only works if configuration of the deal.II library has
+ * detected the presence of this solver. Please read the README file on what
+ * the configure script is looking for and how to provide it.
+ *
  *
  * @section SPDMA1 Interface and Method
  *
@@ -610,7 +611,11 @@ class SparseDirectMA27 : public Subscriptor
  * indicating error codes, etc. It also manages allocation of the
  * right amount of temporary storage required by these functions.
  *
+ * Note that this class only works if configuration of the deal.II library has
+ * detected the presence of this solver. Please read the README file on what
+ * the configure script is looking for and how to provide it.
  *
+ * 
  * @section SPDMA47a Interface and Method
  *
  * For the meaning of the three functions initialize(), factorize(),
@@ -983,6 +988,289 @@ class SparseDirectMA47 : public Subscriptor
                       const unsigned int *ICNTL);
 };
 
+
+
+
+/**
+ * This class provides an interface to the sparse direct solver UMFPACK (see
+ * <a href="http://www.cise.ufl.edu/research/sparse/umfpack">this
+ * link</a>). UMFPACK is a set of routines for solving unsymmetric sparse
+ * linear systems, Ax=b, using the Unsymmetric-pattern MultiFrontal method and
+ * direct sparse LU factorization. Matrices may have symmetric or unsymmetrix
+ * sparsity patterns, and may have unsymmetric entries.
+ *
+ * Note that this class only works if configuration of the deal.II library has
+ * detected the presence of this solver. Please read the README file on
+ * what the configure script is looking for and how to provide it.
+ *
+ * @author Wolfgang Bangerth, 2004
+ */
+class SparseDirectUMFPACK : public Subscriptor
+{
+  public:
+                                    /**
+                                     * Constructor. See the
+                                     * documentation of this class
+                                     * for the meaning of the
+                                     * parameters to this function.
+                                     */
+    SparseDirectUMFPACK ();
+
+                                     /**
+                                      * Destructor.
+                                      */
+    ~SparseDirectUMFPACK ();    
+    
+                                    /**
+                                     * This function does nothing. It is only
+                                     * here to provide an interface that is
+                                     * consistent with that of the HSL MA27
+                                     * and MA47 solver classes.
+                                     */
+    void initialize (const SparsityPattern &sparsity_pattern);
+
+                                    /**
+                                     * Factorize the matrix. This function
+                                     * may be called multiple times for
+                                     * different matrices, after the object
+                                     * of this class has been initialized for
+                                     * a certain sparsity pattern. You may
+                                     * therefore save some computing time if
+                                     * you want to invert several matrices
+                                     * with the same sparsity
+                                     * pattern. However, note that the bulk
+                                     * of the computing time is actually
+                                     * spent in the factorization, so this
+                                     * functionality may not always be of
+                                     * large benefit.
+                                     *
+                                     * If the initialization step has
+                                     * not been performed yet, then
+                                     * the initialize() function is
+                                     * called at the beginning of
+                                     * this function.
+                                     *
+                                     * This function copies the contents of
+                                     * the matrix into its own storage; the
+                                     * matrix can therefore be deleted after
+                                     * this operation, even if subsequent
+                                     * solves are required.
+                                     */
+    void factorize (const SparseMatrix<double> &matrix);
+
+                                    /**
+                                     * Solve for a certain right hand
+                                     * side vector. This function may
+                                     * be called multiple times for
+                                     * different right hand side
+                                     * vectors after the matrix has
+                                     * been factorized. This yields a
+                                     * big saving in computing time,
+                                     * since the actual solution is
+                                     * fast, compared to the
+                                     * factorization of the matrix.
+                                     *
+                                     * The solution will be returned
+                                     * in place of the right hand
+                                     * side vector.
+                                     *
+                                     * If the factorization has not
+                                     * happened before, strange
+                                     * things will happen. Note that
+                                     * we can't actually call the
+                                     * factorize() function from
+                                     * here if it has not yet been
+                                     * called, since we have no
+                                     * access to the actual matrix.
+                                     */
+    void solve (Vector<double> &rhs_and_solution) const;
+
+                                    /**
+                                     * Call the three functions above
+                                     * in that order, i.e. perform
+                                     * the whole solution process for
+                                     * the given right hand side
+                                     * vector.
+                                     *
+                                     * The solution will be returned
+                                     * in place of the right hand
+                                     * side vector.
+                                     */
+    void solve (const SparseMatrix<double> &matrix,
+               Vector<double>             &rhs_and_solution);
+
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcMatrixNotSquare);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcUMFPACKError);
+    
+  private:
+                                     /**
+                                      * The UMFPACK routines allocate objects
+                                      * in which they store information about
+                                      * symbolic and numeric values of the
+                                      * decomposition. The actual data type of
+                                      * these objects is opaque, and only
+                                      * passed around as void pointers.
+                                      */
+    void *symbolic_decomposition;
+    void *numeric_decomposition;
+
+                                     /**
+                                      * Free all memory that hasn't been freed
+                                      * yet.
+                                      */
+    void clear ();
+
+                                     /**
+                                      * The arrays in which we store the data
+                                      * for the solver.
+                                      */
+    std::vector<int> Ap;
+    std::vector<int> Ai;
+    std::vector<double> Ax;
+
+                                     /**
+                                      * Control and info arrays for the solver
+                                      * routines.
+                                      */
+    std::vector<double> control;
+};
+
+
+
+class SparseDirectSuperLU : public Subscriptor
+{
+  public:
+                                    /**
+                                     * Constructor. See the
+                                     * documentation of this class
+                                     * for the meaning of the
+                                     * parameters to this function.
+                                     */
+    SparseDirectSuperLU ();
+
+                                     /**
+                                      * Destructor.
+                                      */
+    ~SparseDirectSuperLU ();    
+    
+                                    /**
+                                     * This function does nothing. It is only
+                                     * here to provide an interface that is
+                                     * consistent with that of the HSL MA27
+                                     * and MA47 solver classes.
+                                     */
+    void initialize (const SparsityPattern &sparsity_pattern);
+
+                                    /**
+                                     * Factorize the matrix. This function
+                                     * may be called multiple times for
+                                     * different matrices, after the object
+                                     * of this class has been initialized for
+                                     * a certain sparsity pattern. You may
+                                     * therefore save some computing time if
+                                     * you want to invert several matrices
+                                     * with the same sparsity
+                                     * pattern. However, note that the bulk
+                                     * of the computing time is actually
+                                     * spent in the factorization, so this
+                                     * functionality may not always be of
+                                     * large benefit.
+                                     *
+                                     * If the initialization step has
+                                     * not been performed yet, then
+                                     * the initialize() function is
+                                     * called at the beginning of
+                                     * this function.
+                                     *
+                                     * This function copies the contents of
+                                     * the matrix into its own storage; the
+                                     * matrix can therefore be deleted after
+                                     * this operation, even if subsequent
+                                     * solves are required.
+                                     */
+    void factorize (const SparseMatrix<double> &matrix);
+
+                                    /**
+                                     * Solve for a certain right hand
+                                     * side vector. This function may
+                                     * be called multiple times for
+                                     * different right hand side
+                                     * vectors after the matrix has
+                                     * been factorized. This yields a
+                                     * big saving in computing time,
+                                     * since the actual solution is
+                                     * fast, compared to the
+                                     * factorization of the matrix.
+                                     *
+                                     * The solution will be returned
+                                     * in place of the right hand
+                                     * side vector.
+                                     *
+                                     * If the factorization has not
+                                     * happened before, strange
+                                     * things will happen. Note that
+                                     * we can't actually call the
+                                     * factorize() function from
+                                     * here if it has not yet been
+                                     * called, since we have no
+                                     * access to the actual matrix.
+                                     */
+    void solve (Vector<double> &rhs_and_solution) const;
+
+                                    /**
+                                     * Call the three functions above
+                                     * in that order, i.e. perform
+                                     * the whole solution process for
+                                     * the given right hand side
+                                     * vector.
+                                     *
+                                     * The solution will be returned
+                                     * in place of the right hand
+                                     * side vector.
+                                     */
+    void solve (const SparseMatrix<double> &matrix,
+               Vector<double>             &rhs_and_solution);
+
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcMatrixNotSquare);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcSuperLUError);
+    
+  private:
+                                     /**
+                                      * A data type that holds all the data we
+                                      * need to preserve between calls to
+                                      * factorize() and solve(). The actual
+                                      * definition of this structure is in the
+                                      * source file since it depends on
+                                      * SuperLU's data types and we don't want
+                                      * to include their header file into this
+                                      * one.
+                                      */
+    struct Data;
+
+                                     /**
+                                      * One such object.
+                                      */
+    Data *data;
+    
+                                     /**
+                                      * Free all memory that hasn't been freed
+                                      * yet.
+                                      */
+    void clear ();
+};
+
 /*@}*/
 
 
index ac3ff6d620c8b4c3626ac3f317cc979651709620..cfe1eae35024147a3a3497e08d2b898b9c38e651 100644 (file)
 #  include <hsl/hsl.h>
 #endif
 
+
+// if configured for UMFPACK, include respective file. annoyingly the UMFPACK
+// files don't seem to have extern "C" wrapped around their headers...
+#ifdef DEAL_II_USE_UMFPACK
+extern "C" {
+#  include <umfpack.h>
+}
+#endif
+
+#include "/home/bangerth/tmp/superlu/SuperLU_3.0/SRC/dsp_defs.h"
+
 // if the HSL functions are not there, define them empty and throw an
 // exception
 #ifndef HAVE_HSL_MA27
@@ -1574,6 +1585,539 @@ call_ma47cd (const unsigned int *n_rows,           //scalar
 
 
 
+
+SparseDirectUMFPACK::SparseDirectUMFPACK ()
+                :
+                symbolic_decomposition (0),
+                numeric_decomposition (0),
+                control (UMFPACK_CONTROL)
+{
+  umfpack_di_defaults (&control[0]);
+}
+
+
+
+SparseDirectUMFPACK::~SparseDirectUMFPACK ()
+{
+  clear ();
+}
+
+
+void
+SparseDirectUMFPACK::clear ()
+{
+                                   // delete objects that haven't been deleted
+                                   // yet
+  if (symbolic_decomposition != 0)
+    {
+      umfpack_di_free_symbolic (&symbolic_decomposition);
+      symbolic_decomposition = 0;
+    }
+
+  if (numeric_decomposition != 0)
+    {
+      umfpack_di_free_numeric (&numeric_decomposition);
+      numeric_decomposition = 0;
+    }
+  
+  {
+    std::vector<int> tmp;
+    tmp.swap (Ap);
+  }
+
+  {
+    std::vector<int> tmp;
+    tmp.swap (Ai);
+  }
+  
+  {
+    std::vector<double> tmp;
+    tmp.swap (Ax);
+  }
+
+  umfpack_di_defaults (&control[0]);
+}
+
+
+
+void
+SparseDirectUMFPACK::
+initialize (const SparsityPattern &)
+{}
+
+
+
+void
+SparseDirectUMFPACK::
+factorize (const SparseMatrix<double> &matrix)
+{
+  Assert (matrix.m() == matrix.n(), ExcMatrixNotSquare())
+  
+  clear ();
+
+  const unsigned int N = matrix.m();
+
+                                   // copy over the data from the matrix to
+                                   // the data structures UMFPACK wants. note
+                                   // two things: first, UMFPACK wants
+                                   // compressed column storage whereas we
+                                   // always do compressed row storage; we
+                                   // work around this by, rather than
+                                   // shuffling things around, copy over the
+                                   // data we have, but then call the
+                                   // umfpack_di_solve function with the
+                                   // UMFPACK_At argument, meaning that we
+                                   // want to solve for the transpose system
+                                   //
+                                   // second: the data we have in the sparse
+                                   // matrices is "almost" right already;
+                                   // UMFPACK wants the entries in each row
+                                   // (i.e. really: column) to be sorted in
+                                   // ascending order. we almost have that,
+                                   // except that we usually store the
+                                   // diagonal first in each row to allow for
+                                   // some optimizations. thus, we have to
+                                   // resort things a little bit, but only
+                                   // within each row
+                                   //
+                                   // final note: if the matrix has entries in
+                                   // the sparsity pattern that are actually
+                                   // occupied by entries that have a zero
+                                   // numerical value, then we keep them
+                                   // anyway. people are supposed to provide
+                                   // accurate sparsity patterns.
+  Ap.resize (N+1);
+  Ai.resize (matrix.get_sparsity_pattern().n_nonzero_elements());
+  Ax.resize (matrix.get_sparsity_pattern().n_nonzero_elements());
+
+                                   // first fill row lengths array
+  Ap[0] = 0;
+  for (unsigned int row=1; row<=N; ++row)
+    Ap[row] = Ap[row-1] + matrix.get_sparsity_pattern().row_length(row-1);
+  Assert (static_cast<unsigned int>(Ap.back()) == Ai.size(),
+          ExcInternalError());
+  
+                                   // then copy over matrix elements
+  {
+    unsigned int index = 0;
+    for (SparseMatrix<double>::const_iterator p=matrix.begin();
+         p!=matrix.end(); ++p, ++index)
+      {
+        Ai[index] = p->column();
+        Ax[index] = p->value();
+      }
+    Assert (index == Ai.size(), ExcInternalError());
+  }
+
+                                   // finally do the copying around of entries
+                                   // so that the diagonal entry is in the
+                                   // right place. note that this is easy to
+                                   // detect: since all entries apart from the
+                                   // diagonal entry are sorted, we know that
+                                   // the diagonal entry is in the wrong place
+                                   // if and only if its column index is
+                                   // larger than the column index of the
+                                   // second entry in a row
+                                   //
+                                   // ignore rows with only one or no entry
+  {
+    for (unsigned int row=0; row<N; ++row)
+      {
+                                         // we may have to move some elements
+                                         // that are left of the diagonal but
+                                         // presently after the diagonal entry
+                                         // to the left, whereas the diagonal
+                                         // entry has to move to the right. we
+                                         // could first figure out where to
+                                         // move everything to, but for
+                                         // simplicity we just make a series
+                                         // of swaps instead (this is kind of
+                                         // a single run of bubble-sort, which
+                                         // gives us the desired result since
+                                         // the array is already "almost"
+                                         // sorted)
+                                         //
+                                         // in the first loop, the condition
+                                         // in the while-header also checks
+                                         // that the row has at least two
+                                         // entries and that the diagonal
+                                         // entry is really in the wrong place
+        int cursor = Ap[row];
+        while ((cursor < Ap[row+1]-1) &&
+               (Ai[cursor] > Ai[cursor+1]))
+          {
+            std::swap (Ai[cursor], Ai[cursor+1]);
+            std::swap (Ax[cursor], Ax[cursor+1]);
+            ++cursor;
+          }
+      }
+  }
+  
+  
+        
+  int status;
+
+  status = umfpack_di_symbolic (N, N,
+                                &Ap[0], &Ai[0], &Ax[0],
+                                &symbolic_decomposition,
+                                &control[0], 0);
+  AssertThrow (status == UMFPACK_OK, ExcUMFPACKError());
+  
+  status = umfpack_di_numeric (&Ap[0], &Ai[0], &Ax[0],
+                               symbolic_decomposition,
+                               &numeric_decomposition,
+                               &control[0], 0);
+  AssertThrow (status == UMFPACK_OK, ExcUMFPACKError());
+
+  umfpack_di_free_symbolic (&symbolic_decomposition) ;
+}
+
+
+
+void
+SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
+{
+                                   // make sure that some kind of factorize()
+                                   // call has happened before
+  Assert (Ap.size() != 0, ExcNotInitialized());
+  Assert (Ai.size() != 0, ExcNotInitialized());
+  Assert (Ai.size() == Ax.size(), ExcNotInitialized());
+  
+  Vector<double> rhs (rhs_and_solution.size());
+  rhs = rhs_and_solution;
+  
+                                   // solve the system. note that since
+                                   // UMFPACK wants compressed column storage
+                                   // instead of the compressed row storage
+                                   // format we use in deal.II's
+                                   // SparsityPattern classes, we solve for
+                                   // UMFPACK's A^T instead
+  const int status
+    = umfpack_di_solve (UMFPACK_At,
+                        &Ap[0], &Ai[0], &Ax[0],
+                        rhs_and_solution.begin(), rhs.begin(),
+                        numeric_decomposition,
+                        &control[0], 0);
+  AssertThrow (status == UMFPACK_OK, ExcUMFPACKError());
+}
+
+
+
+void
+SparseDirectUMFPACK::solve (const SparseMatrix<double> &matrix,
+                            Vector<double>             &rhs_and_solution)
+{
+  factorize (matrix);
+  solve (rhs_and_solution);
+}
+
+
+
+
+
+SparseDirectSuperLU::SparseDirectSuperLU ()
+                :
+                data (0)
+{}
+
+
+
+struct SparseDirectSuperLU::Data
+{
+    SuperMatrix A, X, L, U;
+    std::vector<int> perm_r;
+    std::vector<int> perm_c;
+    std::vector<double> solution;
+    std::vector<double> R,C;
+    std::vector<int> etree;
+    char equed[1];
+    void *work;
+    int lwork;
+
+    Data (const unsigned int N);
+    ~Data ();
+};
+
+
+SparseDirectSuperLU::Data::Data (const unsigned int N)
+                :
+                perm_r (N),
+                perm_c (N),
+                solution (N),
+                R (N),
+                C (N),
+                etree (N),
+                work (0),
+                lwork (0)
+{}
+
+
+SparseDirectSuperLU::Data::~Data ()
+{  
+  Destroy_SuperMatrix_Store(&A);
+  Destroy_SuperMatrix_Store(&X);
+  Destroy_SuperNode_Matrix(&L);
+  Destroy_CompCol_Matrix(&U);
+}
+
+
+
+SparseDirectSuperLU::~SparseDirectSuperLU ()
+{
+  clear ();
+}
+
+
+void
+SparseDirectSuperLU::clear ()
+{
+  if (data != 0)
+    delete data;
+  data = 0;
+}
+
+
+
+void
+SparseDirectSuperLU::
+initialize (const SparsityPattern &)
+{}
+
+
+
+void
+SparseDirectSuperLU::
+factorize (const SparseMatrix<double> &matrix)
+{
+  Assert (matrix.m() == matrix.n(), ExcMatrixNotSquare());
+
+                                   // delete old objects if there are any
+  clear ();
+
+  const unsigned int N = matrix.m();
+
+                                   // copy over the data from the matrix to
+                                   // the data structures SuperLU wants. note
+                                   // two things: first, SuperLU wants
+                                   // compressed column storage whereas we
+                                   // always do compressed row storage; we
+                                   // work around this by, rather than
+                                   // shuffling things around, copy over the
+                                   // data we have, but then call the
+                                   // umfpack_di_solve function with the
+                                   // SuperLU_At argument, meaning that we
+                                   // want to solve for the transpose system
+                                   //
+                                   // second: the data we have in the sparse
+                                   // matrices is "almost" right already;
+                                   // SuperLU wants the entries in each row
+                                   // (i.e. really: column) to be sorted in
+                                   // ascending order. we almost have that,
+                                   // except that we usually store the
+                                   // diagonal first in each row to allow for
+                                   // some optimizations. thus, we have to
+                                   // resort things a little bit, but only
+                                   // within each row
+                                   //
+                                   // final note: if the matrix has entries in
+                                   // the sparsity pattern that are actually
+                                   // occupied by entries that have a zero
+                                   // numerical value, then we keep them
+                                   // anyway. people are supposed to provide
+                                   // accurate sparsity patterns.
+  std::vector<int> Ap (N+1);
+  std::vector<int> Ai (matrix.get_sparsity_pattern().n_nonzero_elements());
+  std::vector<double> Ax (matrix.get_sparsity_pattern().n_nonzero_elements());
+
+                                   // first fill row lengths array
+  Ap[0] = 0;
+  for (unsigned int row=1; row<=N; ++row)
+    Ap[row] = Ap[row-1] + matrix.get_sparsity_pattern().row_length(row-1);
+  Assert (static_cast<unsigned int>(Ap.back()) == Ai.size(),
+          ExcInternalError());
+  
+                                   // then copy over matrix elements
+  {
+    unsigned int index = 0;
+    for (SparseMatrix<double>::const_iterator p=matrix.begin();
+         p!=matrix.end(); ++p, ++index)
+      {
+        Ai[index] = p->column();
+        Ax[index] = p->value();
+      }
+    Assert (index == Ai.size(), ExcInternalError());
+  }
+
+                                   // finally do the copying around of entries
+                                   // so that the diagonal entry is in the
+                                   // right place. note that this is easy to
+                                   // detect: since all entries apart from the
+                                   // diagonal entry are sorted, we know that
+                                   // the diagonal entry is in the wrong place
+                                   // if and only if its column index is
+                                   // larger than the column index of the
+                                   // second entry in a row
+                                   //
+                                   // ignore rows with only one or no entry
+  {
+    for (unsigned int row=0; row<N; ++row)
+      {
+                                         // we may have to move some elements
+                                         // that are left of the diagonal but
+                                         // presently after the diagonal entry
+                                         // to the left, whereas the diagonal
+                                         // entry has to move to the right. we
+                                         // could first figure out where to
+                                         // move everything to, but for
+                                         // simplicity we just make a series
+                                         // of swaps instead (this is kind of
+                                         // a single run of bubble-sort, which
+                                         // gives us the desired result since
+                                         // the array is already "almost"
+                                         // sorted)
+                                         //
+                                         // in the first loop, the condition
+                                         // in the while-header also checks
+                                         // that the row has at least two
+                                         // entries and that the diagonal
+                                         // entry is really in the wrong place
+        int cursor = Ap[row];
+        while ((cursor < Ap[row+1]-1) &&
+               (Ai[cursor] > Ai[cursor+1]))
+          {
+            std::swap (Ai[cursor], Ai[cursor+1]);
+            std::swap (Ax[cursor], Ax[cursor+1]);
+            ++cursor;
+          }
+      }
+  }
+
+
+                                   // now factorize the matrix. we need a
+                                   // dummy rhs vector as well as
+                                   // an object to hold the data we need
+  data = new Data (N);
+  
+  std::vector<double> dummy_rhs (N);  
+  SuperMatrix B;
+  
+  dCreate_CompRow_Matrix(&data->A, N, N, Ax.size(),
+                         &Ax[0], &Ai[0], &Ap[0], SLU_NC, SLU_D, SLU_GE);
+  
+  dCreate_Dense_Matrix(&B, N, 1, &dummy_rhs[0], N,
+                       SLU_DN, SLU_D, SLU_GE);
+  dCreate_Dense_Matrix(&data->X, N, 1, &data->solution[0], N,
+                       SLU_DN, SLU_D, SLU_GE);
+
+                                   // set options. note that just as with
+                                   // umfpack, we solve the transpose system,
+                                   // since we give compressed row storage and
+                                   // superlu wants compressed column storage
+  superlu_options_t options;
+  set_default_options(&options);
+  options.Trans = TRANS;
+
+                                   // this seems to be crucial. without we get
+                                   // atrocious performance
+  options.ColPerm = MMD_AT_PLUS_A;
+  options.SymmetricMode = YES;
+  
+                                   // indicate that we don't actually want to
+                                   // solve anything, just to factorize
+  B.ncol = 0;
+
+                                   // lots of unused output arguments of dgssvx
+  int info;
+  double rpg, rcond;
+  double ferr[1];
+  double berr[1];
+  mem_usage_t mem_usage;
+
+  SuperLUStat_t stat;
+  StatInit(&stat);
+  
+                                   // do the factorization
+  dgssvx(&options, &data->A, &data->perm_c[0], &data->perm_r[0],
+         &data->etree[0], data->equed, &data->R[0], &data->C[0],
+         &data->L, &data->U, data->work, data->lwork, &B,
+         &data->X, &rpg, &rcond, ferr, berr,
+         &mem_usage, &stat, &info);
+  AssertThrow (info == 0, ExcSuperLUError());
+
+                                   // delete temp vector again
+  Destroy_SuperMatrix_Store (&B);
+  StatFree(&stat);
+}
+
+
+
+void
+SparseDirectSuperLU::solve (Vector<double> &rhs_and_solution) const
+{
+  const unsigned int N = rhs_and_solution.size();
+
+                                   // create rhs vector
+  SuperMatrix B;
+  dCreate_Dense_Matrix(&B, N, 1, rhs_and_solution.begin(), N,
+                       SLU_DN, SLU_D, SLU_GE);
+
+                                   // set options. note that just as with
+                                   // umfpack, we solve the transpose system,
+                                   // since we give compressed row storage and
+                                   // superlu wants compressed column storage
+  superlu_options_t options;
+  set_default_options(&options);
+  options.Trans = TRANS;
+
+                                   // this seems to be crucial. without we get
+                                   // atrocious performance
+  options.ColPerm = MMD_AT_PLUS_A;
+  options.SymmetricMode = YES;
+
+                                   // indicate that the matrix has already
+                                   // been factorized
+  options.Fact = FACTORED;
+
+                                   // lots of unused output arguments of dgssvx
+  int info;
+  double rpg, rcond;
+  double ferr[1];
+  double berr[1];
+  mem_usage_t mem_usage;
+
+  SuperLUStat_t stat;
+  StatInit(&stat);
+
+                                   // do the solve
+  dgssvx(&options, &data->A, &data->perm_c[0], &data->perm_r[0],
+         &data->etree[0], data->equed, &data->R[0], &data->C[0],
+         &data->L, &data->U, data->work, data->lwork,
+         &B, &data->X, &rpg, &rcond, ferr, berr,
+         &mem_usage, &stat, &info);
+  AssertThrow (info == 0, ExcSuperLUError());
+
+                                   // copy result
+  std::copy ((double*) ((DNformat*) data->X.Store)->nzval,
+             (double*) ((DNformat*) data->X.Store)->nzval + N,
+             rhs_and_solution.begin());
+
+                                   // delete temp vectors
+  Destroy_SuperMatrix_Store(&B);
+  StatFree(&stat);
+}
+
+
+
+void
+SparseDirectSuperLU::solve (const SparseMatrix<double> &matrix,
+                            Vector<double>             &rhs_and_solution)
+{
+  factorize (matrix);
+  solve (rhs_and_solution);
+}
+
+
+
 // explicit instantiations
 template
 void

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.