-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 ();
-};
-
/*@}*/
}
#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
-
-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
projects / dealii-svn.git / commitdiff