include/deal.II/lac/sparse_direct.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: sparse_direct.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__sparse_direct_h
00013 #define __deal2__sparse_direct_h
00014 
00015 
00016 
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/base/exceptions.h>
00019 #include <deal.II/base/subscriptor.h>
00020 #include <deal.II/base/thread_management.h>
00021 #include <deal.II/lac/vector.h>
00022 #include <deal.II/lac/sparse_matrix.h>
00023 #include <deal.II/lac/sparse_matrix_ez.h>
00024 #include <deal.II/lac/block_sparse_matrix.h>
00025 
00026 #ifdef DEAL_II_USE_MUMPS
00027 #  include <deal.II/base/utilities.h>
00028 #  include <dmumps_c.h>
00029 #endif
00030 
00031 DEAL_II_NAMESPACE_OPEN
00032 
00191 class SparseDirectMA27 : public Subscriptor
00192 {
00193   public:
00200     SparseDirectMA27 (const double LIW_factor_1          = 1.2,
00201                       const double LIW_factor_2          = 1,
00202                       const double LA_factor             = 1.5,
00203                       const double LIW_increase_factor_1 = 1.2,
00204                       const double LIW_increase_factor_2 = 1.2,
00205                       const double LA_increase_factor    = 1.2,
00206                       const bool   suppress_output       = true);
00207 
00211     ~SparseDirectMA27 ();
00212 
00228     void set_detached_mode ();
00229 
00234     bool detached_mode_set () const;
00235 
00247     void initialize (const SparsityPattern &sparsity_pattern);
00248 
00273     template <typename number>
00274     void factorize (const SparseMatrix<number> &matrix);
00275 
00301     template <typename number>
00302     void solve (Vector<number> &rhs_and_solution) const;
00303 
00316     template <typename number>
00317     void solve (const SparseMatrix<number> &matrix,
00318                 Vector<double>             &rhs_and_solution);
00319 
00324     std::size_t memory_consumption () const;
00325 
00334     Threads::ThreadMutex & get_synchronisation_lock () const;
00335 
00342     DeclException1 (ExcMA27AFailed,
00343                     int,
00344                     << "The function MA27A failed with an exit code of " << arg1);
00348     DeclException1 (ExcMA27BFailed,
00349                     int,
00350                     << "The function MA27B failed with an exit code of " << arg1);
00354     DeclException1 (ExcMA27CFailed,
00355                     int,
00356                     << "The function MA27C failed with an exit code of " << arg1);
00360     DeclException0 (ExcInitializeAlreadyCalled);
00361 
00365     DeclException0 (ExcFactorizeNotCalled);
00366 
00370     DeclException0 (ExcDifferentSparsityPatterns);
00371 
00375     DeclException2 (ExcReadError,
00376                     int, int,
00377                     << "Error while reading in detached mode. Return value "
00378                     << "for 'read' was " << arg1
00379                     << ", errno has value " << arg2);
00383     DeclException0 (ExcMatrixNotSymmetric);
00385   private:
00396     struct DetachedModeData;
00397 
00404     const bool suppress_output;
00405 
00411     bool detached_mode;
00412 
00419     DetachedModeData *detached_mode_data;
00420 
00428     const double LIW_factor_1;
00429     const double LIW_factor_2;
00430     const double LA_factor;
00431 
00436     const double LIW_increase_factor_1;
00437     const double LIW_increase_factor_2;
00438     const double LA_increase_factor;
00439 
00445     bool initialize_called;
00446     bool factorize_called;
00447 
00454     SmartPointer<const SparsityPattern,SparseDirectMA27> sparsity_pattern;
00455 
00461     unsigned int n_nonzero_elements;
00462 
00467     std::vector<unsigned int> row_numbers;
00468     std::vector<unsigned int> column_numbers;
00469 
00475     std::vector<double>       A;
00476 
00480     unsigned int         LA;
00481 
00491     unsigned int LIW;
00492     std::vector<unsigned int> IW;
00493     std::vector<unsigned int> IKEEP;
00494     std::vector<unsigned int> IW1;
00495 
00496     unsigned int NSTEPS;
00497     unsigned int MAXFRT;
00498 
00509     unsigned int NRLNEC;
00510     unsigned int NIRNEC;
00511 
00517     int IFLAG;
00518 
00523     static Threads::ThreadMutex static_synchronisation_lock;
00524     mutable Threads::ThreadMutex non_static_synchronisation_lock;
00525 
00531     template <typename number>
00532     void fill_A (const SparseMatrix<number> &matrix);
00533 
00539     void call_ma27ad (const unsigned int *N,
00540                       const unsigned int *NZ,
00541                       const unsigned int *IRN,
00542                       const unsigned int *ICN,
00543                       unsigned int       *IW,
00544                       const unsigned int *LIW,
00545                       unsigned int       *IKEEP,
00546                       unsigned int       *IW1,
00547                       unsigned int       *NSTEPS,
00548                       int                *IFLAG);
00549 
00555     void call_ma27bd (const unsigned int *N,
00556                       const unsigned int *NZ,
00557                       const unsigned int *IRN,
00558                       const unsigned int *ICN,
00559                       double             *A,
00560                       const unsigned int *LA,
00561                       unsigned int       *IW,
00562                       const unsigned int *LIW,
00563                       const unsigned int *IKEEP,
00564                       const unsigned int *NSTEPS,
00565                       unsigned int       *MAXFRT,
00566                       unsigned int       *IW1,
00567                       int                *IFLAG);
00568 
00574     void call_ma27cd (const unsigned int *N,
00575                       const double       *A,
00576                       const unsigned int *LA,
00577                       const unsigned int *IW,
00578                       const unsigned int *LIW,
00579                       const unsigned int *MAXFRT,
00580                       double             *RHS,
00581                       const unsigned int *IW1,
00582                       const unsigned int *NSTEPS) const;
00583 
00589     void call_ma27x1 (unsigned int *NRLNEC);
00590 
00596     void call_ma27x2 (unsigned int *NIRNEC);
00597 
00603     void call_ma27x3 (const unsigned int *LP);
00604 };
00605 
00606 
00607 
00608 
00609 
00610 
00681 class SparseDirectMA47 : public Subscriptor
00682 {
00683   public:
00695     SparseDirectMA47 (const double LIW_factor_1 = 1.4,
00696                       const double LIW_factor_2 = 1,
00697                       const double LA_factor    = 3,
00698                       const double LIW_increase_factor_1 = 1.2,
00699                       const double LIW_increase_factor_2 = 1.2,
00700                       const double LA_increase_factor    = 1.2,
00701                       const bool   suppress_output       = true);
00702 
00723     void initialize (const SparseMatrix<double> &matrix);
00724 
00745     void factorize (const SparseMatrix<double> &matrix);
00746 
00772     void solve (Vector<double> &rhs_and_solution);
00773 
00787     void solve (const SparseMatrix<double> &matrix,
00788                 Vector<double>             &rhs_and_solution);
00789 
00794     std::size_t memory_consumption () const;
00795 
00804     Threads::ThreadMutex & get_synchronisation_lock () const;
00805 
00812     DeclException1 (ExcMA47AFailed,
00813                     int,
00814                     << "The function MA47A failed with an exit code of " << arg1);
00818     DeclException1 (ExcMA47BFailed,
00819                     int,
00820                     << "The function MA47B failed with an exit code of " << arg1);
00824     DeclException1 (ExcMA47CFailed,
00825                     int,
00826                     << "The function MA47C failed with an exit code of " << arg1);
00830     DeclException0 (ExcInitializeAlreadyCalled);
00831 
00835     DeclException0 (ExcFactorizeNotCalled);
00836 
00840     DeclException0 (ExcCantFactorizeAgain);
00841 
00845     DeclException0 (ExcDifferentMatrices);
00849     DeclException0 (ExcMatrixNotSymmetric);
00851   private:
00858     const bool suppress_output;
00859 
00867     const double LIW_factor_1;
00868     const double LIW_factor_2;
00869     const double LA_factor;
00870 
00875     const double LIW_increase_factor_1;
00876     const double LIW_increase_factor_2;
00877     const double LA_increase_factor;
00878 
00884     bool initialize_called;
00885     bool factorize_called;
00886 
00892     SmartPointer<const SparseMatrix<double>,SparseDirectMA47> matrix;
00893 
00899     unsigned int n_nonzero_elements;
00900 
00904     double       CNTL[2];
00905     unsigned int ICNTL[7];
00906 
00912     int          INFO[24];
00913 
00918     std::vector<unsigned int> row_numbers;
00919     std::vector<unsigned int> column_numbers;
00920 
00926     std::vector<double>       A;
00927 
00931     unsigned int         LA;
00932 
00942     unsigned int LIW;
00943     std::vector<unsigned int> IW;
00944     std::vector<unsigned int> KEEP;
00945     std::vector<unsigned int> IW1;
00946 
00951     static Threads::ThreadMutex synchronisation_lock;
00952 
00958     void fill_A (const SparseMatrix<double> &matrix);
00959 
00964     void call_ma47id (double       *CNTL,
00965                       unsigned int *ICNTL);
00966 
00971     void call_ma47ad (const unsigned int *n_rows,
00972                       const unsigned int *n_nonzero_elements,
00973                       unsigned int       *row_numbers,
00974                       unsigned int       *column_numbers,
00975                       unsigned int       *IW,
00976                       const unsigned int *LIW,
00977                       unsigned int       *KEEP,
00978                       const unsigned int *ICNTL,
00979                       int                *INFO);
00980 
00985     void call_ma47bd (const unsigned int *n_rows,
00986                       const unsigned int *n_nonzero_elements,
00987                       const unsigned int *column_numbers,
00988                       double             *A,
00989                       const unsigned int *LA,
00990                       unsigned int       *IW,
00991                       const unsigned int *LIW,
00992                       const unsigned int *KEEP,
00993                       const double       *CNTL,
00994                       const unsigned int *ICNTL,
00995                       unsigned int       *IW1,
00996                       int                *INFO);
00997 
01002     void call_ma47cd (const unsigned int *n_rows,
01003                       const double       *A,
01004                       const unsigned int *LA,
01005                       const unsigned int *IW,
01006                       const unsigned int *LIW,
01007                       double             *rhs_and_solution,
01008                       unsigned int       *IW1,
01009                       const unsigned int *ICNTL);
01010 };
01011 
01012 
01013 
01014 
01060 class SparseDirectUMFPACK : public Subscriptor
01061 {
01062   public:
01068     class AdditionalData
01069     {};
01070 
01071 
01078     SparseDirectUMFPACK ();
01079 
01083     ~SparseDirectUMFPACK ();
01084 
01091     void initialize (const SparsityPattern &sparsity_pattern);
01092 
01122     template <class Matrix>
01123     void factorize (const Matrix &matrix);
01124 
01129     template <class Matrix>
01130     void initialize(const Matrix &matrix,
01131                     const AdditionalData additional_data = AdditionalData());
01132 
01144     void vmult (Vector<double>&, const Vector<double>&) const;
01145 
01150     void Tvmult (Vector<double>&, const Vector<double>&) const;
01151 
01157     void vmult_add (Vector<double>&, const Vector<double>&) const;
01158 
01163     void Tvmult_add (Vector<double>&, const Vector<double>&) const;
01164 
01190     void solve (Vector<double> &rhs_and_solution) const;
01191 
01204     template <class Matrix>
01205     void solve (const Matrix   &matrix,
01206                 Vector<double> &rhs_and_solution);
01207 
01217     DeclException2 (ExcUMFPACKError, char*, int,
01218                     << "UMFPACK routine " << arg1
01219                     << " returned error status " << arg2
01220                     << ". See the file <contrib/umfpack/UMFPACK/Include/umfpack.h>"
01221                     << " for a description of error codes.");
01222 
01223   private:
01232     void *symbolic_decomposition;
01233     void *numeric_decomposition;
01234 
01239     void clear ();
01240 
01241 
01253     template <typename number>
01254     void sort_arrays (const SparseMatrixEZ<number> &);
01255 
01256     template <typename number>
01257     void sort_arrays (const SparseMatrix<number> &);
01258 
01259     template <typename number>
01260     void sort_arrays (const BlockSparseMatrix<number> &);
01261 
01266     std::vector<long int> Ap;
01267     std::vector<long int> Ai;
01268     std::vector<double> Ax;
01269 
01274     std::vector<double> control;
01275 };
01276 
01277 
01298 class SparseDirectMUMPS
01299 {
01300  private:
01301 
01302 #ifdef DEAL_II_USE_MUMPS
01303   DMUMPS_STRUC_C id;
01304 #endif // DEAL_II_USE_MUMPS
01305 
01306   double       *a;
01307   double       *rhs;
01308   unsigned int *irn;
01309   unsigned int *jcn;
01310   unsigned int  n;
01311   unsigned int  nz;
01312 
01318   bool initialize_called;
01319 
01320  public:
01321 
01325   SparseDirectMUMPS ();
01326 
01330   ~SparseDirectMUMPS ();
01331 
01335   DeclException0 (ExcInitializeAlreadyCalled);
01336 
01343   template <class Matrix>
01344     void initialize (const Matrix& matrix,
01345                      const Vector<double>      & vector);
01346 
01352   void solve (Vector<double>& vector);
01353 };
01354 
01355 DEAL_II_NAMESPACE_CLOSE
01356 
01357 #endif // __deal2__sparse_direct_h
01358 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3