00001
00002
00003
00004
00005
00006
00007
00008
00009
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