]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove superlu, but it is at least now in the archive.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 20 Aug 2004 20:47:19 +0000 (20:47 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 20 Aug 2004 20:47:19 +0000 (20:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@9567 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 1e56ead71c5e3031354d2f28d4393638b6615580..a3281797bc2d8d1c39c232e2586c99b392c8d3c2 100644 (file)
@@ -1143,134 +1143,6 @@ class SparseDirectUMFPACK : public Subscriptor
 
 
 
-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 cfe1eae35024147a3a3497e08d2b898b9c38e651..2081d1a4f077e660e545d9213e269994ec088a4d 100644 (file)
@@ -51,7 +51,6 @@ extern "C" {
 }
 #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
@@ -1814,310 +1813,6 @@ SparseDirectUMFPACK::solve (const SparseMatrix<double> &matrix,
 
 
 
-
-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.