From 1af9c734f7b20fbc9a5e07dda4e4b31521b2d489 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 20 Aug 2004 20:47:19 +0000 Subject: [PATCH] Remove superlu, but it is at least now in the archive. git-svn-id: https://svn.dealii.org/trunk@9567 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_direct.h | 128 ---------- deal.II/lac/source/sparse_direct.cc | 305 ------------------------ 2 files changed, 433 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h index 1e56ead71c..a3281797bc 100644 --- a/deal.II/lac/include/lac/sparse_direct.h +++ b/deal.II/lac/include/lac/sparse_direct.h @@ -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 &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 &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 &matrix, - Vector &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 (); -}; - /*@}*/ diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc index cfe1eae350..2081d1a4f0 100644 --- a/deal.II/lac/source/sparse_direct.cc +++ b/deal.II/lac/source/sparse_direct.cc @@ -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 &matrix, - -SparseDirectSuperLU::SparseDirectSuperLU () - : - data (0) -{} - - - -struct SparseDirectSuperLU::Data -{ - SuperMatrix A, X, L, U; - std::vector perm_r; - std::vector perm_c; - std::vector solution; - std::vector R,C; - std::vector 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 &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 Ap (N+1); - std::vector Ai (matrix.get_sparsity_pattern().n_nonzero_elements()); - std::vector 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(Ap.back()) == Ai.size(), - ExcInternalError()); - - // then copy over matrix elements - { - unsigned int index = 0; - for (SparseMatrix::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 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 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 &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 &matrix, - Vector &rhs_and_solution) -{ - factorize (matrix); - solve (rhs_and_solution); -} - - - // explicit instantiations template void -- 2.39.5