]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove support for the obsolete MA27 and MA47 direct solvers of the HSL collection
authormaier <maier@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 31 Mar 2013 20:38:00 +0000 (20:38 +0000)
committermaier <maier@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 31 Mar 2013 20:38:00 +0000 (20:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@29122 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/sparse_direct.h
deal.II/source/lac/sparse_direct.cc

index 3bb7acff64ea3afa4c7b8d2cffe274eccaa3b1e3..a5cfafa16115f82c1378dc0f22ccec457a1e0e86 100644 (file)
@@ -25,6 +25,13 @@ inconvenience this causes.
 
 <ol>
 
+<li> Removed: The interfaces to the obsolete direct solvers MA27 and MA47 from
+the Harwell Subroutine Library. Support for the HSL routines were not ported to
+the new build system
+<br>
+(Matthias Maier, 2013/04/01)
+
+
 <li> Changed: The TimeDependent::end_sweep function with an argument indicating
 the number of threads has been removed. Use the corresponding function without
 an argument. Since the argument had a default value, few users will have used
index 5a5d9d673e3e4d93f6867c8f655b63bf28d03ee5..195ca99190f4a6288b43c80acacd575987bd790e 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-/**
- * This class provides an interface to the sparse direct solver MA27
- * from the Harwell Subroutine Library. MA27 is a direct solver
- * specialized for sparse symmetric indefinite systems of linear
- * equations and uses a modified form of Gauss elimination. It is
- * included in the <a
- * href="http://www.cse.clrc.ac.uk/Activity/HSL">Harwell Subroutine
- * Library</a> and is written in Fortran. The present class only
- * transforms the data stored in SparseMatrix objects into the
- * form which is required by the functions resembling MA27, calls
- * these Fortran functions, and interprets some of the returned values
- * indicating error codes, etc. It also manages allocation of the
- * right amount of temporary storage required by these functions.
- *
- * Note that this class only works if configuration of the deal.II library has
- * detected the presence of this solver. Please read the README file on what
- * the configure script is looking for and how to provide it.
- *
- *
- * <h3>Interface and Method</h3>
- *
- * For the meaning of the three functions initialize(), factorize(),
- * and solve(), as well as for the method used in MA27, please see the
- * <a href="http://www.cse.clrc.ac.uk/Activity/HSL">documentation</a>
- * of these functions. In practice, you will most often call the
- * second solve() function, which solves the linear system for a
- * given right hand side, but one can as well call the three functions
- * separately if, for example, one would like to solve the same matrix
- * for several right hand side vectors; the MA27 solver can do this
- * efficiently, as it computes a decomposition of the matrix, so that
- * subsequent solves only amount to a forward-backward substitution
- * which is significantly less costly than the decomposition process.
- *
- *
- * <h3>Parameters to the constructor</h3>
- *
- * The constructor of this class takes several arguments. The meaning
- * is the following: the MA27 functions require the user to allocate
- * and pass a certain amount of memory for temporary variables or for
- * data to be passed to subsequent functions. The sizes of these
- * arrays are denoted by the variables <tt>LIW1</tt>, <tt>LIW2</tt>, and <tt>LA</tt>,
- * where <tt>LIW1</tt> denotes the size of the <tt>IW</tt> array in the call to
- * <tt>MA27A</tt>, while <tt>LIW2</tt> is the array size in the call to
- * <tt>MA27B</tt>. The documentation of the MA27 functions gives ways to
- * obtain estimates for their values, e.g. by evaluating values
- * returned by functions called before. However, the documentation
- * only states that the values have to be <b>at least</b> as large as
- * the estimates, a hint that is not very useful oftentimes (in my
- * humble opinion, the lack of dynamic memory allocation mechanism is
- * a good reason not to program in Fortran 77 :-).
- *
- * In our experience, it is often necessary to go beyond the proposed
- * values (most often for <tt>LA</tt>, but also for <tt>LIW1</tt>). The first
- * three parameters of the constructor denote by which factor the
- * initial estimates shall be increased. The default values are 1.2
- * (the documentation recommends this value, 1, and 1.5, values which
- * have often worked for us. Note that the value of <tt>LIW</tt> is only
- * changed in the second call if the recommended value times
- * <tt>LIW_factor_2</tt> is larger than the array size already is from the
- * call to <tt>MA27A</tt>; otherwise, <tt>LIW_factor_2</tt> is ignored.
- *
- * If the values thus constructed fail to work, we try to restart the
- * called function with larger values until the calls succeed. The
- * second triple of values passed to the constructor denotes by which
- * factor we shall increase the array sizes. If the increment factors
- * are less than or equal to one, then we only try to call the
- * respective calls to the functions once and abort by throwing an
- * error. Note that the <tt>MA27C</tt> function writes out an error message
- * if the value of <tt>LA</tt> is too small and gives an indication to
- * which size it should be increased. However, most often the
- * indicated value is far too small and can not be relied upon.
- *
- *
- * <h3>Note on parallelization</h3>
- *
- * <h4>Synchronisation</h4>
- *
- * Due to the use of global variables through COMMON blocks, the calls
- * to the sparse direct solver routines are not multithreading-safe,
- * i.e. at each time there may only be one call to these functions
- * active. You have to synchronise your calls to the functions
- * provided by this class using mutexes (see the Threads
- * namespace for such classes) to avoid multiple active calls at the
- * same time if you use multithreading. Since you may use this class
- * in different parts of your program, and may not want to use a
- * global variable for locking, this class has a lock as static member
- * variable, which may be accessed using the
- * get_synchronisation_lock() function. Note however, that this class
- * does not perform the synchronisation for you within its member
- * functions. The reason is that you will usually want to synchronise
- * over the calls to initialize() and factorize(), since there should
- * probably not be a call to one of these function with another matrix
- * between the calls for one matrix. (The author does not really know
- * whether this is true, but it is probably safe to assume that.)
- * Since such cross-function synchronisation can only be performed
- * from outside, it is left to the user of this class to do so.
- *
- * <h4>Detached mode</h4>
- *
- * As an alternative, you can call the function set_detached_mode()
- * right after calling the constructor. This lets the program fork, so
- * that we now have two programs that communicate via pipes.  The
- * forked copy of the program then actually replaces itself by a
- * program called <tt>detached_ma27</tt>, that is started in its place
- * through the <tt>execv</tt> system call. Now every time you call one of
- * the functions of this class, it relays the data to the other
- * program and lets it execute the respective function. The results
- * are then transferred back. Since the MA27 functions are only called
- * in the detached program, they will now no longer interfere with the
- * respective calls to other functions with different data, so no
- * synchronisation is necessary any more.
- *
- * The advantage of this approach is that as many instances of this
- * class may be active at any time as you want. This is handy, if your
- * programs spens a significant amount of time in them, and you are
- * using many threads, for example in a machine with 4 or more
- * processors. The disadvantage, of course, is that the data has to
- * copied to and from the detached program, which might make things
- * slower (though, as we use block writes, this should not be so much
- * of a factor).
- *
- * Since no more synchronisation is necessary, the
- * get_synchronisation_lock() returns a reference to a member
- * variable when the detached mode is set. Thus, you need not change
- * your program: you can still acquire and release the lock as before,
- * it will only have no effect now, since different objects of this
- * class no longer share the lock, i.e. you will get it always without
- * waiting. On the other hand, it will prevent that you call functions
- * of this object multiply in parallel at the same time, which is what
- * you probably wanted.
- *
- *
- * <h5>Internals of the detached mode</h5>
- *
- * The program that actually runs the detached solver is called
- * <tt>detached_ma27</tt>, and will show up under this name in the process
- * list. It communicates with the main program through a pipe.
- *
- * Since the solver and the main program are two separated processes,
- * the solver program will not be notified if the main program dies,
- * for example because it is aborted with Control-C, because an
- * exception is raised and not caught, or some other reason. It will
- * just not get any new jobs, but will happily wait until the end of
- * times. For this reason, the detached solver has a second thread
- * running in parallel that simply checks in regular intervals whether
- * the main program is still alive, using the <tt>ps</tt> program. If this
- * is no longer the case, the detached solver exits as well.
- *
- * Since the intervals between two such checks are a couple of second,
- * it may happen that the detached solver survives the main program by
- * some time. Presently, the check interval is once every 20
- * seconds. After that time, the detached solver should have noticed
- * the main programs demise.
- *
- * @ingroup Solvers Preconditioners
- *
- * @author Wolfgang Bangerth, 2000, 2001, 2002
- */
-class SparseDirectMA27 : public Subscriptor
-{
-public:
-  /**
-   * Constructor. See the
-   * documentation of this class
-   * for the meaning of the
-   * parameters to this function.
-   */
-  SparseDirectMA27 (const double LIW_factor_1          = 1.2,
-                    const double LIW_factor_2          = 1,
-                    const double LA_factor             = 1.5,
-                    const double LIW_increase_factor_1 = 1.2,
-                    const double LIW_increase_factor_2 = 1.2,
-                    const double LA_increase_factor    = 1.2,
-                    const bool   suppress_output       = true);
-
-  /**
-   * Destructor.
-   */
-  ~SparseDirectMA27 ();
-
-  /**
-   * Set the detached mode (see the
-   * general class documentation
-   * for a description of what this
-   * is).
-   *
-   * This function must not be
-   * called after initialize()
-   * (or the two-argument solve()
-   * function has been called. If
-   * it is to be called, then only
-   * right after construction of
-   * the object, and before first
-   * use.
-   */
-  void set_detached_mode ();
-
-  /**
-   * Return whether the detached
-   * mode is set.
-   */
-  bool detached_mode_set () const;
-
-  /**
-   * Initialize some data
-   * structures. This function
-   * computes symbolically some
-   * information based on the
-   * sparsity pattern, but does not
-   * actually use the values of the
-   * matrix, so only the sparsity
-   * pattern has to be passed as
-   * argument.
-   */
-  void initialize (const SparsityPattern &sparsity_pattern);
-
-  /**
-   * Actually 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.
-   */
-  template <typename number>
-  void factorize (const SparseMatrix<number> &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.
-   */
-  template <typename number>
-  void solve (Vector<number> &rhs_and_solution) const;
-
-  /**
-   * Call the three functions
-   * initialize, factorize and solve
-   * 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.
-   */
-  template <typename number>
-  void solve (const SparseMatrix<number> &matrix,
-              Vector<double>             &rhs_and_solution);
-
-  /**
-   * Return an estimate of the
-   * memory used by this class.
-   */
-  std::size_t memory_consumption () const;
-
-  /**
-   * Get a reference to the
-   * synchronisation lock which can
-   * be used for this class. See
-   * the general description of
-   * this class for more
-   * information.
-   */
-  Threads::Mutex &get_synchronisation_lock () const;
-
-  /** @addtogroup Exceptions
-   * @{ */
-
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcMA27AFailed,
-                  int,
-                  << "The function MA27A failed with an exit code of " << arg1);
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcMA27BFailed,
-                  int,
-                  << "The function MA27B failed with an exit code of " << arg1);
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcMA27CFailed,
-                  int,
-                  << "The function MA27C failed with an exit code of " << arg1);
-  /**
-   * Exception
-   */
-  DeclException0 (ExcInitializeAlreadyCalled);
-
-  /**
-   * Exception
-   */
-  DeclException0 (ExcFactorizeNotCalled);
-
-  /**
-   * Exception
-   */
-  DeclException0 (ExcDifferentSparsityPatterns);
-
-  /**
-   * Exception
-   */
-  DeclException2 (ExcReadError,
-                  int, int,
-                  << "Error while reading in detached mode. Return value "
-                  << "for 'read' was " << arg1
-                  << ", errno has value " << arg2);
-  /**
-   * Exception
-   */
-  DeclException0 (ExcMatrixNotSymmetric);
-  //@}
-private:
-  /**
-   * Declare a local type which
-   * will store the data necessary
-   * to communicate with a detached
-   * solver. To avoid adding
-   * various system include files,
-   * the actual declaration of this
-   * class is in the implementation
-   * file.
-   */
-  struct DetachedModeData;
-
-  /**
-   * Store in the constructor
-   * whether the MA27 routines
-   * shall deliver output to stdout
-   * or not.
-   */
-  const bool suppress_output;
-
-  /**
-   * Store whether
-   * set_detached_mode() has been
-   * called.
-   */
-  bool detached_mode;
-
-  /**
-   * Pointer to a structure that
-   * will hold the data necessary
-   * to uphold communication with a
-   * detached solver.
-   */
-  DetachedModeData *detached_mode_data;
-
-  /**
-   * Store the three values passed
-   * to the cinstructor. See the
-   * documentation of this class
-   * for the meaning of these
-   * variables.
-   */
-  const double LIW_factor_1;
-  const double LIW_factor_2;
-  const double LA_factor;
-
-  /**
-   * Increase factors in case a
-   * call to a function fails.
-   */
-  const double LIW_increase_factor_1;
-  const double LIW_increase_factor_2;
-  const double LA_increase_factor;
-
-  /**
-   * Flags storing whether the
-   * first two functions have
-   * already been called.
-   */
-  bool initialize_called;
-  bool factorize_called;
-
-  /**
-   * Store a pointer to the
-   * sparsity pattern, to make sure
-   * that we use the same thing for
-   * all calls.
-   */
-  SmartPointer<const SparsityPattern,SparseDirectMA27> sparsity_pattern;
-
-  /**
-   * Number of nonzero elements in
-   * the sparsity pattern on and
-   * above the diagonal.
-   */
-  unsigned int n_nonzero_elements;
-
-  /**
-   * Arrays holding row and column
-   * indices.
-   */
-  std::vector<unsigned int> row_numbers;
-  std::vector<unsigned int> column_numbers;
-
-  /**
-   * Array to hold the matrix
-   * elements, and later the
-   * elements of the factors.
-   */
-  std::vector<double>       A;
-
-  /**
-   * Length of the <tt>A</tt> array.
-   */
-  unsigned int         LA;
-
-  /**
-   * Scratch arrays and variables
-   * used by the MA27 functions. We
-   * keep to the names introduced
-   * in the documentation of these
-   * functions, in all uppercase
-   * letters as is usual in
-   * Fortran.
-   */
-  unsigned int LIW;
-  std::vector<unsigned int> IW;
-  std::vector<unsigned int> IKEEP;
-  std::vector<unsigned int> IW1;
-
-  unsigned int NSTEPS;
-  unsigned int MAXFRT;
-
-  /**
-   * Two values that live inside a
-   * COMMON block of the Fortran
-   * code and are mirrored at these
-   * locations. They are used to
-   * transport information about
-   * the required length of arrays
-   * from the Fortran functions to
-   * the outside world.
-   */
-  unsigned int NRLNEC;
-  unsigned int NIRNEC;
-
-  /**
-   * Flag indicating the level of
-   * output desired and returning
-   * error values if error occured.
-   */
-  int IFLAG;
-
-  /**
-   * Mutexes for synchronising access
-   * to this class.
-   */
-  static Threads::Mutex static_synchronisation_lock;
-  mutable Threads::Mutex non_static_synchronisation_lock;
-
-  /**
-   * Fill the <tt>A</tt> array from the
-   * symmetric part of the given
-   * matrix.
-   */
-  template <typename number>
-  void fill_A (const SparseMatrix<number> &matrix);
-
-  /**
-   * Call the respective function
-   * with the given args, either
-   * locally or remote.
-   */
-  void call_ma27ad (const unsigned int *N,
-                    const unsigned int *NZ,
-                    const unsigned int *IRN,
-                    const unsigned int *ICN,
-                    unsigned int       *IW,
-                    const unsigned int *LIW,
-                    unsigned int       *IKEEP,
-                    unsigned int       *IW1,
-                    unsigned int       *NSTEPS,
-                    int                *IFLAG);
-
-  /**
-   * Call the respective function
-   * with the given args, either
-   * locally or remote.
-   */
-  void call_ma27bd (const unsigned int *N,
-                    const unsigned int *NZ,
-                    const unsigned int *IRN,
-                    const unsigned int *ICN,
-                    double             *A,
-                    const unsigned int *LA,
-                    unsigned int       *IW,
-                    const unsigned int *LIW,
-                    const unsigned int *IKEEP,
-                    const unsigned int *NSTEPS,
-                    unsigned int       *MAXFRT,
-                    unsigned int       *IW1,
-                    int                *IFLAG);
-
-  /**
-   * Call the respective function
-   * with the given args, either
-   * locally or remote.
-   */
-  void call_ma27cd (const unsigned int *N,
-                    const double       *A,
-                    const unsigned int *LA,
-                    const unsigned int *IW,
-                    const unsigned int *LIW,
-                    const unsigned int *MAXFRT,
-                    double             *RHS,
-                    const unsigned int *IW1,
-                    const unsigned int *NSTEPS) const;
-
-  /**
-   * Call the respective function
-   * with the given args, either
-   * locally or remote.
-   */
-  void call_ma27x1 (unsigned int *NRLNEC);
-
-  /**
-   * Call the respective function
-   * with the given args, either
-   * locally or remote.
-   */
-  void call_ma27x2 (unsigned int *NIRNEC);
-
-  /**
-   * Call the respective function
-   * with the given args, either
-   * locally or remote.
-   */
-  void call_ma27x3 (const unsigned int *LP);
-};
-
-
-
-
-
-
-/**
- * This class provides an interface to the sparse direct solver MA47
- * from the Harwell Subroutine Library. MA47 is a direct solver
- * specialized for sparse symmetric indefinite systems of linear
- * equations and uses a frontal elimination method. It is included in
- * the <a href="http://www.cse.clrc.ac.uk/Activity/HSL">Harwell
- * Subroutine Library</a> and is written in Fortran. The present class
- * only transforms the data stored in SparseMatrix objects into
- * the form which is required by the functions resembling MA47, calls
- * these Fortran functions, and interprets some of the returned values
- * indicating error codes, etc. It also manages allocation of the
- * right amount of temporary storage required by these functions.
- *
- * Note that this class only works if configuration of the deal.II library has
- * detected the presence of this solver. Please read the README file on what
- * the configure script is looking for and how to provide it.
- *
- *
- * <h3>Interface and Method</h3>
- *
- * For the meaning of the three functions initialize(), factorize(),
- * and solve(), as well as for the method used in MA47, please see the
- * <a href="http://www.cse.clrc.ac.uk/Activity/HSL">documentation</a>
- * of these functions. In practice, one will most often call the
- * second solve() function, which solves the linear system for a given
- * right hand side, but one can as well call the three functions
- * separately if, for example, one would like to solve the same matrix
- * for several right hand side vectors; the MA47 solver can do this
- * efficiently, as it computes a decomposition of the matrix, so that
- * subsequent solves only amount to a forward-backward substitution
- * which is significantly less costly than the decomposition process.
- *
- *
- * <h3>Parameters to the constructor</h3>
- *
- * The constructor of this class takes several arguments. Their
- * meaning is equivalent to those of the constructor of the
- * SparseDirectMA27 class; see there for more information.
- *
- *
- * <h3>Note on parallelization</h3>
- *
- * Due to the use of global variables through COMMON blocks, the calls
- * to the sparse direct solver routines is not multithreading-capable,
- * i.e. at each time there may only be one call to these functions
- * active. You have to synchronise your calls to the functions
- * provided by this class using mutexes (see the Threads
- * namespace for such classes) to avoid multiple active calls at the
- * same time if you use multithreading. Since you may use this class
- * in different parts of your program, and may not want to use a
- * global variable for locking, this class has a lock as static member
- * variable, which may be accessed using the
- * get_synchronisation_lock() function. Note however, that this class
- * does not perform the synchronisation for you within its member
- * functions. The reason is that you will usually want to synchronise
- * over the calls to initialize() and factorize(), since there should
- * probably not be a call to one of these function with another matrix
- * between the calls for one matrix. (The author does not really know
- * whether this is true, but it is probably safe to assume that.)
- * Since such cross-function synchronisation can only be performed
- * from outside, it is left to the user of this class to do so.
- *
- * A detached mode as for MA27 has not yet been implemented for this
- * class.
- *
- *
- * @ingroup Solvers Preconditioners
- *
- * @author Wolfgang Bangerth, 2000, 2001
- */
-class SparseDirectMA47 : public Subscriptor
-{
-public:
-  /**
-   * Constructor. See the
-   * documentation of this class
-   * for the meaning of the
-   * parameters to this function.
-   *
-   * This function already calls
-   * the initialization function
-   * <tt>MA47ID</tt> to set up some
-   * values.
-   */
-  SparseDirectMA47 (const double LIW_factor_1 = 1.4,
-                    const double LIW_factor_2 = 1,
-                    const double LA_factor    = 3,
-                    const double LIW_increase_factor_1 = 1.2,
-                    const double LIW_increase_factor_2 = 1.2,
-                    const double LA_increase_factor    = 1.2,
-                    const bool   suppress_output       = true);
-
-  /**
-   * Initialize some data
-   * structures. This function
-   * computes symbolically some
-   * information based on the
-   * sparsity pattern, but does not
-   * actually use the values of the
-   * matrix, so only the sparsity
-   * pattern has to be passed as
-   * argument.
-   *
-   * Since the MA47 solver requires
-   * us to omit zero-entries in the
-   * matrix (even if they are in
-   * the sparsity pattern), we have
-   * to actually use the matrix
-   * here, as opposed to the MA27
-   * solver that only required the
-   * sparsity pattern.
-   */
-  void initialize (const SparseMatrix<double> &matrix);
-
-  /**
-   * Actually factorize the
-   * matrix. Unlike for the MA27
-   * solver, this function may not
-   * be called multiple times for
-   * different matrices, since we
-   * have eliminated entries from
-   * the sparsity pattern where
-   * matrix entries happen to be
-   * zero. Since this is likely to
-   * change between matrices
-   * although they have the same
-   * sparsity pattern.
-   *
-   * If the initialization step has
-   * not been performed yet, then
-   * the initialize() function is
-   * called at the beginning of
-   * this function.
-   */
-  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);
-
-  /**
-   * Call the three functions
-   * initialize, factorize and
-   * solve
-   * 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);
-
-  /**
-   * Return an estimate of the
-   * memory used by this class.
-   */
-  std::size_t memory_consumption () const;
-
-  /**
-   * Get a reference to the
-   * synchronisation lock which can
-   * be used for this class. See
-   * the general description of
-   * this class for more
-   * information.
-   */
-  Threads::Mutex &get_synchronisation_lock () const;
-
-  /** @addtogroup Exceptions
-   * @{ */
-
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcMA47AFailed,
-                  int,
-                  << "The function MA47A failed with an exit code of " << arg1);
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcMA47BFailed,
-                  int,
-                  << "The function MA47B failed with an exit code of " << arg1);
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcMA47CFailed,
-                  int,
-                  << "The function MA47C failed with an exit code of " << arg1);
-  /**
-   * Exception
-   */
-  DeclException0 (ExcInitializeAlreadyCalled);
-
-  /**
-   * Exception
-   */
-  DeclException0 (ExcFactorizeNotCalled);
-
-  /**
-   * Exception
-   */
-  DeclException0 (ExcCantFactorizeAgain);
-
-  /**
-   * Exception
-   */
-  DeclException0 (ExcDifferentMatrices);
-  /**
-   * Exception
-   */
-  DeclException0 (ExcMatrixNotSymmetric);
-  //@}
-private:
-  /**
-   * Store in the constructor
-   * whether the MA47 routines
-   * shall deliver output to stdout
-   * or not.
-   */
-  const bool suppress_output;
-
-  /**
-   * Store the three values passed
-   * to the cinstructor. See the
-   * documentation of this class
-   * for the meaning of these
-   * variables.
-   */
-  const double LIW_factor_1;
-  const double LIW_factor_2;
-  const double LA_factor;
-
-  /**
-   * Increase factors in case a
-   * call to a function fails.
-   */
-  const double LIW_increase_factor_1;
-  const double LIW_increase_factor_2;
-  const double LA_increase_factor;
-
-  /**
-   * Flags storing whether the
-   * first two functions have
-   * already been called.
-   */
-  bool initialize_called;
-  bool factorize_called;
-
-  /**
-   * Store a pointer to the matrix,
-   * to make sure that we use the
-   * same thing for all calls.
-   */
-  SmartPointer<const SparseMatrix<double>,SparseDirectMA47> matrix;
-
-  /**
-   * Number of nonzero elements in
-   * the sparsity pattern on and
-   * above the diagonal.
-   */
-  unsigned int n_nonzero_elements;
-
-  /**
-   * Control values set by <tt>MA47ID</tt>.
-   */
-  double       CNTL[2];
-  unsigned int ICNTL[7];
-
-  /**
-   * Info field filled by the MA47
-   * functions and (partially) used
-   * for subsequent MA47 calls.
-   */
-  int          INFO[24];
-
-  /**
-   * Arrays holding row and column
-   * indices.
-   */
-  std::vector<unsigned int> row_numbers;
-  std::vector<unsigned int> column_numbers;
-
-  /**
-   * Array to hold the matrix
-   * elements, and later the
-   * elements of the factors.
-   */
-  std::vector<double>       A;
-
-  /**
-   * Length of the <tt>A</tt> array.
-   */
-  unsigned int         LA;
-
-  /**
-   * Scratch arrays and variables
-   * used by the MA47 functions. We
-   * keep to the names introduced
-   * in the documentation of these
-   * functions, in all uppercase
-   * letters as is usual in
-   * Fortran.
-   */
-  unsigned int LIW;
-  std::vector<unsigned int> IW;
-  std::vector<unsigned int> KEEP;
-  std::vector<unsigned int> IW1;
-
-  /**
-   * Mutex for synchronising access
-   * to this class.
-   */
-  static Threads::Mutex synchronisation_lock;
-
-  /**
-   * Fill the <tt>A</tt> array from the
-   * symmetric part of the given
-   * matrix.
-   */
-  void fill_A (const SparseMatrix<double> &matrix);
-
-  /**
-   * Call the <tt>ma47id</tt> function
-   * with the given args.
-   */
-  void call_ma47id (double       *CNTL,
-                    unsigned int *ICNTL);
-
-  /**
-   * Call the <tt>ma47ad</tt> function
-   * with the given args.
-   */
-  void call_ma47ad (const unsigned int *n_rows,
-                    const unsigned int *n_nonzero_elements,
-                    unsigned int       *row_numbers,
-                    unsigned int       *column_numbers,
-                    unsigned int       *IW,
-                    const unsigned int *LIW,
-                    unsigned int       *KEEP,
-                    const unsigned int *ICNTL,
-                    int                *INFO);
-
-  /**
-   * Call the <tt>ma47bd</tt> function
-   * with the given args.
-   */
-  void call_ma47bd (const unsigned int *n_rows,
-                    const unsigned int *n_nonzero_elements,
-                    const unsigned int *column_numbers,
-                    double             *A,
-                    const unsigned int *LA,
-                    unsigned int       *IW,
-                    const unsigned int *LIW,
-                    const unsigned int *KEEP,
-                    const double       *CNTL,
-                    const unsigned int *ICNTL,
-                    unsigned int       *IW1,
-                    int                *INFO);
-
-  /**
-   * Call the <tt>ma47bd</tt> function
-   * with the given args.
-   */
-  void call_ma47cd (const unsigned int *n_rows,
-                    const double       *A,
-                    const unsigned int *LA,
-                    const unsigned int *IW,
-                    const unsigned int *LIW,
-                    double             *rhs_and_solution,
-                    unsigned int       *IW1,
-                    const unsigned int *ICNTL);
-};
-
-
 
 
 /**
@@ -1051,8 +71,7 @@ private:
  *
  * There are instantiations of this class for SparseMatrix<double>,
  * SparseMatrix<float>, SparseMatrixEZ<float>, SparseMatrixEZ<double>,
- * BlockSparseMatrix<double>, and
- * BlockSparseMatrix<float>.
+ * BlockSparseMatrix<double>, and BlockSparseMatrix<float>.
  *
  * @ingroup Solvers Preconditioners
  *
@@ -1062,19 +81,16 @@ class SparseDirectUMFPACK : public Subscriptor
 {
 public:
   /**
-   * Dummy class needed for the
-   * usual initalization interface
-   * of preconditioners.
+   * Dummy class needed for the usual initalization interface of
+   * preconditioners.
    */
   class AdditionalData
 {};
 
 
   /**
-   * Constructor. See the
-   * documentation of this class
-   * for the meaning of the
-   * parameters to this function.
+   * Constructor. See the documentation of this class for the meaning of
+   * the parameters to this function.
    */
   SparseDirectUMFPACK ();
 
@@ -1084,136 +100,94 @@ public:
   ~SparseDirectUMFPACK ();
 
   /**
-   * 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.
+   * This function does nothing. It is only here to provide a consistent
+   * interface.
    */
   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.
+   * 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.
    *
-   * In contrast to the other direct solver
-   * classes, the initialisation method
-   * does nothing. Therefore initialise
-   * is not automatically called by this
-   * method, when
-   * the initialization step has
-   * not been performed yet.
+   * In contrast to the other direct solver classes, the initialisation
+   * method does nothing. Therefore initialise is not automatically called
+   * by this method, when the initialization step has not been performed
+   * yet.
    *
-   * 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.
+   * 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.
    */
   template <class Matrix>
   void factorize (const Matrix &matrix);
 
   /**
-   * Initialize memory and call
-   * SparseDirectUMFPACK::factorize.
+   * Initialize memory and call SparseDirectUMFPACK::factorize.
    */
   template <class Matrix>
   void initialize(const Matrix &matrix,
                   const AdditionalData additional_data = AdditionalData());
 
   /**
-   * Preconditioner interface
-   * function. Usually, given the source
-   * vector, this method returns an
-   * approximated solution of <i>Ax
-   * = b</i>. As this class provides a
-   * wrapper to a direct solver, here
-   * it is actually the exact solution
-   * (exact within the range of numerical
-   * accuracy of course).
+   * Preconditioner interface function. Usually, given the source vector,
+   * this method returns an approximated solution of <i>Ax = b</i>. As this
+   * class provides a wrapper to a direct solver, here it is actually the
+   * exact solution (exact within the range of numerical accuracy of
+   * course).
    */
   void vmult (Vector<double> &, const Vector<double> &) const;
 
   /**
-   * Not implemented but necessary
-   * for compiling.
+   * Not implemented but necessary for compiling.
    */
   void Tvmult (Vector<double> &, const Vector<double> &) const;
 
   /**
-   * Same as vmult(), but adding to
-   * the previous solution. Not
-   * implemented yet.
+   * Same as vmult(), but adding to the previous solution. Not implemented
+   * yet.
    */
   void vmult_add (Vector<double> &, const Vector<double> &) const;
 
   /**
-   * Not implemented but necessary
-   * for compiling.
+   * Not implemented but necessary for compiling.
    */
   void Tvmult_add (Vector<double> &, const Vector<double> &) const;
 
   /**
-   * 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.
+   * 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.
+   * 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.
+   * 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 two functions
-   * factorize and solve
-   * in that order, i.e. perform
-   * the whole solution process for
-   * the given right hand side
-   * vector.
+   * Call the two functions factorize and solve 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.
+   * The solution will be returned in place of the right hand side vector.
    */
   template <class Matrix>
   void solve (const Matrix   &matrix,
               Vector<double> &rhs_and_solution);
 
   /**
-   * One of the UMFPack routines
-   * threw an error. The error code
-   * is included in the output and
-   * can be looked up in the
-   * UMFPack user manual. The name
-   * of the routine is included for
-   * reference.
+   * One of the UMFPack routines threw an error. The error code is included
+   * in the output and can be looked up in the UMFPack user manual. The
+   * name of the routine is included for reference.
    */
   DeclException2 (ExcUMFPACKError, char *, int,
                   << "UMFPACK routine " << arg1
@@ -1223,33 +197,24 @@ public:
 
 private:
   /**
-   * The UMFPACK routines allocate objects
-   * in which they store information about
-   * symbolic and numeric values of the
-   * decomposition. The actual data type of
-   * these objects is opaque, and only
-   * passed around as void pointers.
+   * The UMFPACK routines allocate objects in which they store information
+   * about symbolic and numeric values of the decomposition. The actual
+   * data type of these objects is opaque, and only passed around as void
+   * pointers.
    */
   void *symbolic_decomposition;
   void *numeric_decomposition;
 
   /**
-   * Free all memory that hasn't been freed
-   * yet.
+   * Free all memory that hasn't been freed yet.
    */
   void clear ();
 
-
   /**
-   * Make sure that the arrays Ai
-   * and Ap are sorted in each
-   * row. UMFPACK wants it this
-   * way. We need to have three
-   * versions of this function, one
-   * for the usual SparseMatrix, one
-   * for the SparseMatrixEZ, and
-   * one for the BlockSparseMatrix
-   * classes
+   * Make sure that the arrays Ai and Ap are sorted in each row. UMFPACK
+   * wants it this way. We need to have three versions of this function,
+   * one for the usual SparseMatrix, one for the SparseMatrixEZ, and one
+   * for the BlockSparseMatrix classes
    */
   template <typename number>
   void sort_arrays (const SparseMatrixEZ<number> &);
@@ -1261,32 +226,29 @@ private:
   void sort_arrays (const BlockSparseMatrix<number> &);
 
   /**
-   * The arrays in which we store the data
-   * for the solver.
+   * The arrays in which we store the data for the solver.
    */
   std::vector<long int> Ap;
   std::vector<long int> Ai;
   std::vector<double> Ax;
 
   /**
-   * Control and info arrays for the solver
-   * routines.
+   * Control and info arrays for the solver routines.
    */
   std::vector<double> control;
 };
 
 
 /**
- * This class provides an interface to the parallel sparse direct
- * solver <a href="http://mumps.enseeiht.fr">MUMPS</a>. MUMPS is
- * direct method based on a multifrontal approach, which performs a
- * direct LU factorization. The matrix coming in may have either
- * symmetric or nonsymmetric sparsity pattern.
+ * This class provides an interface to the parallel sparse direct solver
+ * <a href="http://mumps.enseeiht.fr">MUMPS</a>. MUMPS is direct method
+ * based on a multifrontal approach, which performs a direct LU
+ * factorization. The matrix coming in may have either symmetric or
+ * nonsymmetric sparsity pattern.
  *
- * @note This class is useable if and only if a working installation
- * of <a href="http://mumps.enseeiht.fr">MUMPS</a> exists on your
- * system and was detected during configuration of
- * <code>deal.II</code>.
+ * @note This class is useable if and only if a working installation of <a
+ * href="http://mumps.enseeiht.fr">MUMPS</a> exists on your system and was
+ * detected during configuration of <code>deal.II</code>.
  *
  * <h4>Instantiations</h4>
  *
@@ -1312,23 +274,20 @@ private:
   unsigned int  nz;
 
   /**
-   * This function initializes a MUMPS instance
-   * and hands over the system's matrix
-   * <tt>matrix</tt>.
+   * This function initializes a MUMPS instance and hands over the system's
+   * matrix <tt>matrix</tt>.
    */
   template<class Matrix>
   void initialize_matrix (const Matrix &matrix);
 
   /**
-   * Copy the computed solution into the
-   * solution vector.
+   * Copy the computed solution into the solution vector.
    */
   void copy_solution (Vector<double> &vector);
 
   /**
-   * Flags storing whether the function
-   * <tt>initialize ()</tt> has already been
-   * called.
+   * Flags storing whether the function <tt>initialize ()</tt> has already
+   * been called.
    */
   bool initialize_called;
 
@@ -1350,35 +309,30 @@ public:
   DeclException0 (ExcInitializeAlreadyCalled);
 
   /**
-   * This function initializes a MUMPS instance
-   * and hands over the system's matrix
-   * <tt>matrix</tt> and right-hand side
-   * <tt>vector</tt> to the solver.
+   * This function initializes a MUMPS instance and hands over the system's
+   * matrix <tt>matrix</tt> and right-hand side <tt>vector</tt> to the
+   * solver.
    */
   template <class Matrix>
   void initialize (const Matrix &matrix,
                    const Vector<double>       &vector);
 
   /**
-   * This function initializes a MUMPS instance
-   * and computes the factorization of the
-   * system's matrix <tt>matrix</tt>.
+   * This function initializes a MUMPS instance and computes the
+   * factorization of the system's matrix <tt>matrix</tt>.
    */
   template <class Matrix>
   void initialize (const Matrix &matrix);
 
   /**
-   * A function in which the linear system is
-   * solved and the solution vector is copied
-   * into the given <tt>vector</tt>.
+   * A function in which the linear system is solved and the solution
+   * vector is copied into the given <tt>vector</tt>.
    */
   void solve (Vector<double> &vector);
 
   /**
-   * A function in which the inverse of the
-   * matrix is applied to the input vector
-   * <tt>src</tt> and the solution is
-   * written into the output vector
+   * A function in which the inverse of the matrix is applied to the input
+   * vector <tt>src</tt> and the solution is written into the output vector
    * <tt>dst</tt>.
    */
   void vmult (Vector<double> &dst, const Vector<double> &src);
index 2b10e88c2daeffd2157d5c065cbebdda0c58420e..1328eb20ce4c10c261a6091d019c79aac88dd98b 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 
-// if we know that at least one of the HSL functions are there,
-// include the respective include file. Otherwise save some CPU cycles
-// in the compiler
-#if defined(HAVE_HSL_MA27) || defined(HAVE_HSL_MA47)
-#  include <hsl/hsl.h>
-#endif
-
 // include UMFPACK file.
 #ifdef DEAL_II_WITH_UMFPACK
 #  include <umfpack.h>
 #endif
 
-// if the HSL functions are not there, define them empty and throw an
-// exception
-#ifndef HAVE_HSL_MA27
-namespace HSL
-{
-  namespace MA27
-  {
-    extern "C"
-    void ma27ad_ (const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  unsigned int *,
-                  const unsigned int *,
-                  unsigned int *,
-                  unsigned int *,
-                  unsigned int *,
-                  int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C"
-    void ma27bd_ (const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  double *,
-                  const unsigned int *,
-                  unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  unsigned int *,
-                  unsigned int *,
-                  int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C"
-    void ma27cd_ (const unsigned int *,
-                  const double *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  double *,
-                  const unsigned int *,
-                  double *,
-                  const unsigned int *,
-                  const unsigned int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C" void ma27x1_ (unsigned int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C" void ma27x2_ (unsigned int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C" void ma27x3_ (const unsigned int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-  }
-}
-#endif  // ifndef HAVE_HSL_MA27
-
-
-#ifndef HAVE_HSL_MA47
-namespace HSL
-{
-  namespace MA47
-  {
-    extern "C"
-    void ma47id_ (double *,
-                  unsigned int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C"
-    void ma47ad_ (const unsigned int *,
-                  const unsigned int *,
-                  unsigned int *,
-                  unsigned int *,
-                  unsigned int *,
-                  const unsigned int *,
-                  unsigned int *,
-                  const unsigned int *,
-                  double *,
-                  int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C"
-    void ma47bd_ (const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  double *,
-                  const unsigned int *,
-                  unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  const double *,
-                  const unsigned int *,
-                  unsigned int *,
-                  double *,
-                  int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-
-
-    extern "C"
-    void ma47cd_ (const unsigned int *,
-                  const double *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  const unsigned int *,
-                  double *,
-                  double *,
-                  unsigned int *,
-                  const unsigned int *)
-    {
-      AssertThrow (false,
-                   ExcMessage("You can only use the HSL functions after putting "
-                              "the respective files in the right place, "
-                              "re-configuring the library and re-building it!"));
-    }
-  }
-}
-#endif   // ifndef HAVE_HSL_MA47
-
-
-
-
-namespace
-{
-  /**
-   * Output an error message and terminate the program.
-   */
-  void die (const std::string &text,
-            const pid_t child)
-  {
-    std::cerr << "+++++ detached_ma27 driver(" << child << "): " << text
-              << std::endl;
-    std::abort ();
-  }
-
-
-  /**
-   * Output an error message and terminate the program. Write two error
-   * codes.
-   */
-  template <typename T1, typename T2>
-  void die (const std::string &text,
-            const T1 t1,
-            const T2 t2,
-            const pid_t child)
-  {
-    std::cerr << "+++++ detached_ma27 driver(" << child << "): " << text
-              << " code1=" << t1 << ", code2=" << t2
-              << std::endl;
-    std::abort ();
-  }
-}
-
-
-
-/* -------------------------- MA27 ---------------------------- */
-
-Threads::Mutex SparseDirectMA27::static_synchronisation_lock;
-
-
-struct SparseDirectMA27::DetachedModeData
-{
-  /**
-   * Mutex to assure that only one
-   * thread is currently talking
-   * through the pipe.
-   */
-  Threads::Mutex mutex;
-
-  /**
-   * File handles for the pipe
-   * between server (computing
-   * process) and client (display
-   * process).
-   */
-  int server_client_pipe[2];
-  int client_server_pipe[2];
-
-  /**
-   * PID of the forked child
-   * process.
-   */
-  pid_t child_pid;
-
-  /**
-   * Put a message from the server
-   * to the client program. Obey
-   * all the rules the operating
-   * system sets, and create a log
-   * entry for this communication
-   */
-  template <typename T>
-  void put (const T *t,
-            const std::size_t N,
-            const char * /*debug_info*/) const
-  {
-    unsigned int count = 0;
-    while (count < sizeof(T)*N)
-      {
-        // repeat writing until
-        // syscall is not
-        // interrupted
-        int ret = -1;
-#ifndef DEAL_II_MSVC
-        do
-          ret = write (server_client_pipe[1],
-                       reinterpret_cast<const char *> (t) + count,
-                       sizeof(T) * N - count);
-        while ((ret<0) && (errno==EINTR));
-#else
-        Assert (false,
-                ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-        if (ret < 0)
-          die ("error on client side in 'put'", ret, errno, child_pid);
-
-        count += ret;
-      };
-
-    std::fflush (NULL);
-  }
-
-
-  /**
-   * Get a message from the client
-   * program. Obey all the rules
-   * the operating system sets, and
-   * create a log entry for this
-   * communication
-   */
-  template <typename T>
-  void get (T *t,
-            const std::size_t N,
-            const char * /*debug_info*/) const
-  {
-    unsigned int count = 0;
-    while (count < sizeof(T)*N)
-      {
-        int ret = -1;
-#ifndef DEAL_II_MSVC
-        do
-          ret = write (server_client_pipe[1],
-                       reinterpret_cast<const char *> (t) + count,
-                       sizeof(T) * N - count);
-        while ((ret<0) && (errno==EINTR));
-#else
-        Assert (false,
-                ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-
-        if (ret < 0)
-          die ("error on client side in 'get'", ret, errno, child_pid);
-
-        count += ret;
-      }
-  }
-};
-
-
-
-SparseDirectMA27::SparseDirectMA27 (const double LIW_factor_1,
-                                    const double LIW_factor_2,
-                                    const double LA_factor,
-                                    const double LIW_increase_factor_1,
-                                    const double LIW_increase_factor_2,
-                                    const double LA_increase_factor,
-                                    const bool   suppress_output)
-  :
-  suppress_output (suppress_output),
-  detached_mode (false),
-  detached_mode_data (0),
-  LIW_factor_1 (LIW_factor_1),
-  LIW_factor_2 (LIW_factor_2),
-  LA_factor (LA_factor),
-  LIW_increase_factor_1 (LIW_increase_factor_1),
-  LIW_increase_factor_2 (LIW_increase_factor_2),
-  LA_increase_factor (LA_increase_factor),
-  initialize_called (false),
-  factorize_called (false),
-  sparsity_pattern (0, typeid(*this).name())
-{}
-
-
-
-SparseDirectMA27::~SparseDirectMA27()
-{
-  if (detached_mode)
-    if (detached_mode_data != 0)
-      {
-        // close down client
-        Threads::Mutex::ScopedLock lock (detached_mode_data->mutex);
-        // Assign the result of write
-        // and reset the variable to
-        // avoid compiler warnings
-#ifndef DEAL_II_MSVC
-//TODO:[WB] Shouldn't t be used to trace errors?
-        ssize_t t = write (detached_mode_data->server_client_pipe[1], "7", 1);
-        (void)t;
-#else
-        Assert (false,
-                ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-        // then also delete data
-        delete detached_mode_data;
-        detached_mode_data = 0;
-      }
-}
-
-
-
-void
-SparseDirectMA27::set_detached_mode ()
-{
-  Assert (initialize_called == false,
-          ExcInitializeAlreadyCalled());
-  detached_mode = true;
-}
-
-
-
-bool
-SparseDirectMA27::detached_mode_set () const
-{
-  return detached_mode;
-}
-
-
-
-void
-SparseDirectMA27::initialize (const SparsityPattern &sp)
-{
-  Assert (initialize_called == false,
-          ExcInitializeAlreadyCalled());
-
-
-  // first thing is: if detached mode
-  // is requested, then we need to
-  // spawn an instance of the
-  // detached solver and open
-  // communication channels with it
-  if (detached_mode_set())
-    {
-      Assert (detached_mode_data == 0, ExcInternalError());
-      detached_mode_data = new DetachedModeData();
-
-      // create pipes to which we can
-      // write and from which the
-      // slave process will read its
-      // stdin
-
-      // Assign the return value to a
-      // variable to avoid compiler
-      // warnings
-#ifndef DEAL_II_MSVC
-//TODO:[WB] Use t to trace errors?
-      int t = pipe(detached_mode_data->server_client_pipe);
-      (void)t;
-#else
-      Assert (false,
-              ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-      // fflush(NULL) is said to be a
-      // good idea before fork()
-      std::fflush(NULL);
-
-      // now fork and create child
-      // process
-#ifndef DEAL_II_MSVC
-      // BG comment out until pipes are implemented in MSVC
-      detached_mode_data->child_pid = fork();
-#else
-      Assert (false,
-              ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-
-      if (detached_mode_data->child_pid == 0)
-        // child process starts here
-        {
-          // copy read end of input
-          // pipe to stdin, and
-          // likewise with write end
-          // of pipe to stdout
-#ifndef DEAL_II_MSVC
-          dup2(detached_mode_data->server_client_pipe[0], 0);
-          close(detached_mode_data->server_client_pipe[0]);
-
-          dup2(detached_mode_data->client_server_pipe[1], 1);
-          close(detached_mode_data->client_server_pipe[1]);
-
-          // then dispose of this
-          // copy of the program, and
-          // run the detached solver
-          // slave instead
-          /*
-           * TODO: Does this invocation work with just the name of the
-           * executable? Maier, 2012
-           */
-          const char *const program_name = "detached_ma27";
-          const char *const child_argv[] = { program_name, NULL };
-          execv(program_name, const_cast<char *const *>(child_argv));
-
-
-          // usually execv does not
-          // return. if it does, then an
-          // error happened and we report it
-          // herewith:
-          AssertThrow (false,
-                       ExcMessage ("execv returned, which it is not supposed to do!"));
-          std::exit(1);
-
-#else
-          Assert (false,
-                  ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-        };
-      // parent process continues
-      // here.  first thing is to
-      // send the process id of the
-      // present process. this is
-      // used to make sure that the
-      // client can end itself when
-      // it finds that the master
-      // process was somehow
-      // terminated without sending
-      // him this information
-#ifndef DEAL_II_MSVC
-      const pid_t parent_pid = getpid();
-      detached_mode_data->put (&parent_pid, 1, "parent_pid");
-#else
-      Assert (false,
-              ExcMessage ("Detached mode isn't currently implemented on Windows"));
-#endif
-    };
-
-
-  // suppress error output if
-  // requested
-  if (suppress_output)
-    {
-      const unsigned int LP = 0;
-      call_ma27x3 (&LP);
-    };
-
-  sparsity_pattern = &sp;
-
-  const unsigned int
-  n_rows           = sparsity_pattern->n_rows();
-
-  // first count number of nonzero elements in the upper right part. the
-  // matrix is symmetric, so this suffices
-  n_nonzero_elements = 0;
-  for (unsigned int row=0; row<n_rows; ++row)
-    for (SparsityPattern::iterator col = sparsity_pattern->begin(row);
-         col < sparsity_pattern->end(row); ++col)
-      if (row <= col->column())
-        ++n_nonzero_elements;
-
-
-  // fill the row numbers and column numbers arrays from the sparsity
-  // pattern. note that we have Fortran convention, i.e. indices need to be
-  // 1-base, as opposed to C's 0-based convention!
-  row_numbers.resize (n_nonzero_elements);
-  column_numbers.resize (n_nonzero_elements);
-
-  unsigned int global_index = 0;
-  for (unsigned int row=0; row<n_rows; ++row)
-    for (SparsityPattern::iterator col = sparsity_pattern->begin(row);
-         col < sparsity_pattern->end(row); ++col)
-      // note that the matrix must be
-      // symmetric, so only treat the
-      // upper right part
-      if (row <= col->column())
-        {
-          Assert (global_index < n_nonzero_elements, ExcInternalError());
-
-          row_numbers[global_index] = row+1;
-          column_numbers[global_index] = col->column()+1;
-          ++global_index;
-        };
-  Assert (global_index == n_nonzero_elements, ExcInternalError());
-
-  // initialize scratch arrays and
-  // variables
-  LIW = static_cast<unsigned int>((2*n_nonzero_elements + 3*n_rows + 1) *
-                                  LIW_factor_1);
-  IW.resize    (detached_mode_set() ? 0 : LIW);
-  IKEEP.resize (detached_mode_set() ? 0 : 3*n_rows);
-  IW1.resize   (detached_mode_set() ? 0 : 2*n_rows);
-
-  // no output please
-  IFLAG = 0;
-
-  // loop until memory requirements
-  // are satisfied or we are not
-  // allowed to allocate more memory
-  // no more
-  bool call_succeeded = true;
-  do
-    {
-      call_ma27ad (&n_rows, &n_nonzero_elements,
-                   &row_numbers[0], &column_numbers[0],
-                   &IW[0], &LIW, &IKEEP[0],
-                   &IW1[0], &NSTEPS, &IFLAG);
-      call_succeeded = (IFLAG==0);
-
-      // if enough memory or no
-      // increase allowed: exit loop
-      if (call_succeeded || (LIW_increase_factor_1 <= 1))
-        break;
-
-      // otherwise: increase LIW and retry
-      LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_1);
-      IW.resize (LIW);
-    }
-  while (true);
-
-  // if we were not allowed to
-  // allocate more memory, then throw
-  // an exception
-  AssertThrow (call_succeeded, ExcMA27AFailed(IFLAG));
-
-  // catch returned values from the
-  // COMMON block. we need these
-  // values in order to set array
-  // sizes in the next function
-  call_ma27x1 (&NRLNEC);
-  call_ma27x2 (&NIRNEC);
-
-  // note that we have already been
-  // in this function
-  initialize_called = true;
-}
-
-
-
-template <typename number>
-void
-SparseDirectMA27::factorize (const SparseMatrix<number> &matrix)
-{
-  // if necessary, initialize process
-  if (initialize_called == false)
-    initialize (matrix.get_sparsity_pattern());
-
-  // make sure the sparsity patterns
-  // are the same
-  Assert (sparsity_pattern == &matrix.get_sparsity_pattern(),
-          ExcDifferentSparsityPatterns());
-
-
-  // set LA and fill the A array of
-  // values
-  LA = std::max (static_cast<int>(NRLNEC * LA_factor),
-                 static_cast<int>(n_nonzero_elements));
-  A.resize (LA);
-  fill_A (matrix);
-
-  // if necessary extend IW
-  if (LIW < NIRNEC * LIW_factor_2)
-    {
-      LIW = static_cast<unsigned int>(NIRNEC * LIW_factor_2);
-      IW.resize (LIW);
-    };
-
-  const unsigned int n_rows = matrix.get_sparsity_pattern().n_rows();
-
-  // loop until memory requirements
-  // are satisfied or we are not
-  // allowed to allocate more memory
-  // no more
-  bool call_succeeded = true;
-  do
-    {
-      call_ma27bd (&n_rows, &n_nonzero_elements,
-                   &row_numbers[0], &column_numbers[0],
-                   &A[0], &LA,
-                   &IW[0], &LIW, &IKEEP[0], &NSTEPS, &MAXFRT,
-                   &IW1[0], &IFLAG);
-      call_succeeded = (IFLAG==0);
-
-      // if enough memory or no
-      // increase allowed: exit
-      // loop. delete data that is no
-      // more used
-      if (call_succeeded)
-        {
-          std::vector<unsigned int> tmp1, tmp2, tmp3;
-          row_numbers.swap (tmp1);
-          column_numbers.swap (tmp2);
-          IKEEP.swap (tmp3);
-
-          break;
-        };
-
-
-      // otherwise: increase LIW or
-      // LA if that is allowed and
-      // retry
-      switch (IFLAG)
-        {
-        case -3:
-        {
-          if (LIW_increase_factor_2 <= 1)
-            goto exit_loop;
-
-          LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_2);
-          IW.resize (LIW);
-          break;
-        };
-
-        case -4:
-        {
-          if (LA_increase_factor <= 1)
-            goto exit_loop;
-          // increase A. note that
-          // since the function has
-          // already part of the
-          // array @p{A}, we have
-          // to re-fill it with the
-          // original values. minor
-          // clue: since the old
-          // entries are no more
-          // needed, we can discard
-          // them; we use this to
-          // first release all
-          // memory (through the
-          // call to @p{swap} and
-          // the subsequent call to
-          // the destructor of the
-          // @p{tmp} object) and
-          // only then re-allocate
-          // it. If we called
-          // @p{resize} directly,
-          // this would first
-          // allocate more memory,
-          // then copy the old
-          // contents, and only
-          // then release the old
-          // memory, but keeping
-          // both memory regions at
-          // the same time could
-          // sometimes be more than
-          // we can do, leading to
-          // an exception on the
-          // allocation.
-          std::cout << "<*>" << std::flush;
-
-          LA  = static_cast<unsigned int>(LA * LA_increase_factor);
-          if (true)
-            {
-              std::vector<double> tmp;
-              A.swap (tmp);
-            };
-
-          A.resize (LA);
-          fill_A (matrix);
-
-          break;
-        };
-
-        // ups, other return
-        // value, don't know
-        // what to do here
-        default:
-          AssertThrow (false, ExcMA27BFailed(IFLAG));
-        };
-      continue;
-
-exit_loop:
-      break;
-    }
-  while (true);
-
-  AssertThrow (call_succeeded, ExcMA27BFailed(IFLAG));
-
-  // note that we have been here
-  // already and release the sparsity
-  // pattern object, since we won't
-  // need it any more
-  factorize_called = true;
-  sparsity_pattern = 0;
-}
-
-
-
-template <>
-void
-SparseDirectMA27::solve (Vector<double> &rhs_and_solution) const
-{
-  Assert (factorize_called == true, ExcFactorizeNotCalled());
-
-  const unsigned int n_rows = rhs_and_solution.size();
-  call_ma27cd (&n_rows, &A[0], &LA,
-               &IW[0], &LIW, &MAXFRT,
-               &rhs_and_solution(0), &IW1[0], &NSTEPS);
-}
-
-
-
-template <>
-void
-SparseDirectMA27::solve (Vector<float> &rhs_and_solution) const
-{
-  Assert (factorize_called == true, ExcFactorizeNotCalled());
-
-  // first have to convert data type to
-  // doubles
-  Vector<double> tmp (rhs_and_solution.size());
-  tmp = rhs_and_solution;
-
-  const unsigned int n_rows = rhs_and_solution.size();
-  call_ma27cd (&n_rows, &A[0], &LA,
-               &IW[0], &LIW, &MAXFRT,
-               &tmp(0), &IW1[0], &NSTEPS);
-
-  // then copy result back
-  rhs_and_solution = tmp;
-}
-
-
-
-template <typename number>
-void
-SparseDirectMA27::solve (const SparseMatrix<number> &matrix,
-                         Vector<double>             &rhs_and_solution)
-{
-  initialize (matrix.get_sparsity_pattern());
-  factorize (matrix);
-  solve (rhs_and_solution);
-}
-
-
-
-std::size_t
-SparseDirectMA27::memory_consumption () const
-{
-  return (sizeof(*this) +
-          MemoryConsumption::memory_consumption (row_numbers) +
-          MemoryConsumption::memory_consumption (column_numbers) +
-          MemoryConsumption::memory_consumption (A) +
-          MemoryConsumption::memory_consumption (IW) +
-          MemoryConsumption::memory_consumption (IKEEP) +
-          MemoryConsumption::memory_consumption (IW1));
-}
-
-
-
-Threads::Mutex &
-SparseDirectMA27::get_synchronisation_lock () const
-{
-  if (detached_mode)
-    return non_static_synchronisation_lock;
-  else
-    return static_synchronisation_lock;
-}
-
-
-
-template <typename number>
-void
-SparseDirectMA27::fill_A (const SparseMatrix<number> &matrix)
-{
-  Assert (n_nonzero_elements <= A.size(), ExcInternalError());
-
-  const SparsityPattern &sparsity_pattern = matrix.get_sparsity_pattern ();
-
-  const unsigned int n_rows = sparsity_pattern.n_rows();
-
-  unsigned int global_index = 0;
-  for (unsigned int row=0; row<n_rows; ++row)
-    for (typename SparseMatrix<number>::const_iterator col=matrix.begin(row);
-         col < matrix.end(row); ++col)
-      // note that the matrix must be
-      // symmetric, so only treat the
-      // upper right part
-      if (row <= col->column())
-        {
-          Assert (global_index < n_nonzero_elements, ExcInternalError());
-
-          A[global_index] = col->value();
-          ++global_index;
-
-          // make sure that the symmetric
-          // entry exists and has the same
-          // value, unless this one is zero
-          Assert ((col->value() == 0)
-                  ||
-                  (std::fabs(col->value() - matrix(col->column(),row))
-                   <= 1e-15 * std::fabs (col->value())),
-                  ExcMatrixNotSymmetric());
-        }
-      else
-        // lower left part. just check
-        // symmetry
-        Assert ((col->value() == 0)
-                ||
-                (std::fabs(col->value() - matrix(col->column(),row))
-                 <= 1e-15 * std::fabs (col->value())),
-                ExcMatrixNotSymmetric());
-
-  Assert (global_index == n_nonzero_elements, ExcInternalError());
-}
-
-
-
-
-void SparseDirectMA27::call_ma27ad (const unsigned int *N,
-                                    const unsigned int *NZ,
-                                    const unsigned int *IRN,
-                                    const unsigned int *ICN,
-                                    unsigned int       *IW,
-                                    const unsigned int *LIW,
-                                    unsigned int       *IKEEP,
-                                    unsigned int       *IW1,
-                                    unsigned int       *NSTEPS,
-                                    int                *IFLAG)
-{
-  if (detached_mode_set() == false)
-    HSL::MA27::ma27ad_ (N, NZ, IRN, ICN, IW, LIW,
-                        IKEEP, IW1, NSTEPS, IFLAG);
-  else
-    {
-      Threads::Mutex::ScopedLock lock (detached_mode_data->mutex);
-      // first write the data we have
-      // to push over, i.e. first
-      // function index, then array
-      // sizes, then arrays
-      detached_mode_data->put ("1", 1, "ACTION 1");
-
-      detached_mode_data->put (N,   1, "N");
-      detached_mode_data->put (NZ,  1, "NZ");
-      detached_mode_data->put (IRN, *NZ, "IRN");
-      detached_mode_data->put (ICN, *NZ, "ICN");
-      detached_mode_data->put (LIW, 1, "LIW");
-      detached_mode_data->put (IFLAG, 1, "IFLAG");
-
-      // all other fields are kept at
-      // the client. array should not
-      // be in used on this side
-      Assert (this->IKEEP.size() == 0, ExcInternalError());
-      Assert (this->IW1.size() == 0, ExcInternalError());
-
-      // next get back what we need
-      // to know
-      detached_mode_data->get (IFLAG, 1, "IFLAG");
-    };
-}
-
-
-
-void SparseDirectMA27::call_ma27bd (const unsigned int *N,
-                                    const unsigned int *NZ,
-                                    const unsigned int *IRN,
-                                    const unsigned int *ICN,
-                                    double             *A,
-                                    const unsigned int *LA,
-                                    unsigned int       *IW,
-                                    const unsigned int *LIW,
-                                    const unsigned int *IKEEP,
-                                    const unsigned int *NSTEPS,
-                                    unsigned int       *MAXFRT,
-                                    unsigned int       *IW1,
-                                    int                *IFLAG)
-{
-  if (detached_mode_set() == false)
-    HSL::MA27::ma27bd_ (N, NZ, IRN, ICN, A, LA, IW, LIW,
-                        IKEEP, NSTEPS, MAXFRT, IW1, IFLAG);
-  else
-    {
-      // basically, everything is
-      // already over the line,
-      // except for A and LA
-      Threads::Mutex::ScopedLock lock (detached_mode_data->mutex);
-      detached_mode_data->put ("2", 1, "ACTION 2");
-
-      detached_mode_data->put (LA, 1, "LA");
-      detached_mode_data->put (A,  *LA, "A");
-
-      // next get back what we need
-      // to know
-      detached_mode_data->get (IFLAG, 1, "IFLAG");
-    };
-}
-
-
-
-void SparseDirectMA27::call_ma27cd (const unsigned int *N,
-                                    const double       *A,
-                                    const unsigned int *LA,
-                                    const unsigned int *IW,
-                                    const unsigned int *LIW,
-                                    const unsigned int *MAXFRT,
-                                    double             *RHS,
-                                    const unsigned int *IW1,
-                                    const unsigned int *NSTEPS) const
-{
-  if (detached_mode_set() == false)
-    {
-      std::vector<double> W(*MAXFRT);
-      HSL::MA27::ma27cd_ (N, A, LA, IW, LIW, &W[0], MAXFRT, RHS, IW1, NSTEPS);
-    }
-  else
-    {
-      detached_mode_data->put ("3", 1, "ACTION 3");
-
-      // we only have to push and get
-      // the rhs vector
-      detached_mode_data->put (RHS, *N, "RHS");
-      detached_mode_data->get (RHS, *N, "RHS");
-    };
-}
-
-
-
-void SparseDirectMA27::call_ma27x1 (unsigned int *NRLNEC)
-{
-  if (detached_mode_set() == false)
-    HSL::MA27::ma27x1_ (NRLNEC);
-  else
-    {
-      Threads::Mutex::ScopedLock lock (detached_mode_data->mutex);
-      // ma27x1 only reads data, so
-      // don't send anything except
-      // for the id
-      detached_mode_data->put ("4", 1, "ACTION 4");
-      detached_mode_data->get (NRLNEC, 1, "NRLNEC");
-    };
-}
-
-
-
-void SparseDirectMA27::call_ma27x2 (unsigned int *NIRNEC)
-{
-  if (detached_mode_set() == false)
-    HSL::MA27::ma27x2_ (NIRNEC);
-  else
-    {
-      Threads::Mutex::ScopedLock lock (detached_mode_data->mutex);
-      // ma27x2 only reads data, so
-      // don't send anything except
-      // for the id
-      detached_mode_data->put ("5", 1, "ACTION 5");
-      detached_mode_data->get (NIRNEC, 1, "NIRNEC");
-    };
-}
-
-
-
-void SparseDirectMA27::call_ma27x3 (const unsigned int *LP)
-{
-  if (detached_mode_set() == false)
-    HSL::MA27::ma27x3_ (LP);
-  else
-    {
-      Threads::Mutex::ScopedLock lock (detached_mode_data->mutex);
-      // ma27x2 only reads data, so
-      // don't send anything except
-      // for the id
-      detached_mode_data->put ("6", 1, "ACTION 6");
-      detached_mode_data->put (LP, 1, "LP");
-    };
-}
-
-
-
-
-
-/* -------------------------- MA47 ---------------------------- */
-
-Threads::Mutex SparseDirectMA47::synchronisation_lock;
-
-
-SparseDirectMA47::SparseDirectMA47 (const double LIW_factor_1,
-                                    const double LIW_factor_2,
-                                    const double LA_factor,
-                                    const double LIW_increase_factor_1,
-                                    const double LIW_increase_factor_2,
-                                    const double LA_increase_factor,
-                                    const bool   suppress_output)
-  :
-  suppress_output (suppress_output),
-  LIW_factor_1 (LIW_factor_1),
-  LIW_factor_2 (LIW_factor_2),
-  LA_factor (LA_factor),
-  LIW_increase_factor_1 (LIW_increase_factor_1),
-  LIW_increase_factor_2 (LIW_increase_factor_2),
-  LA_increase_factor (LA_increase_factor),
-  initialize_called (false),
-  factorize_called (false),
-  matrix (0, typeid(*this).name())
-{}
-
-
-
-void
-SparseDirectMA47::initialize (const SparseMatrix<double> &m)
-{
-  Assert (initialize_called == false,
-          ExcInitializeAlreadyCalled());
-
-  // some initialization stuff
-  call_ma47id (CNTL, ICNTL);
-  if (suppress_output)
-    ICNTL[0] = 0;
-
-  // then start with work
-  matrix = &m;
-  const SparsityPattern &sparsity_pattern = matrix->get_sparsity_pattern();
-
-  const unsigned int
-  n_rows           = sparsity_pattern.n_rows();
-
-  // first count number of nonzero
-  // elements in the upper right
-  // part. the matrix is symmetric,
-  // so this suffices
-  n_nonzero_elements = 0;
-  for (unsigned int row=0; row<n_rows; ++row)
-    for (SparseMatrix<double>::const_iterator col = m.begin(row);
-         col < m.end(row); ++col)
-      // skip zero elements, as required by the docs of MA47
-      if (row <= col->column() && col->value() != 0)
-        ++n_nonzero_elements;
-
-
-  // fill the row numbers and column
-  // numbers arrays from the sparsity
-  // pattern. note that we have
-  // Fortran convention, i.e. indices
-  // need to be 1-base, as opposed to
-  // C's 0-based convention!
-  row_numbers.resize (n_nonzero_elements);
-  column_numbers.resize (n_nonzero_elements);
-
-  unsigned int global_index = 0;
-  for (unsigned int row=0; row<n_rows; ++row)
-    for (SparseMatrix<double>::const_iterator col = m.begin(row);
-         col < m.end(row); ++col)
-      // note that the matrix must be
-      // symmetric, so only treat the
-      // upper right part
-      if ((row <= col->column()) && (col->value() != 0))
-        {
-          Assert (global_index < n_nonzero_elements, ExcInternalError());
-
-          row_numbers[global_index] = row+1;
-          column_numbers[global_index] = col->column()+1;
-          ++global_index;
-        };
-  Assert (global_index == n_nonzero_elements, ExcInternalError());
-
-  // initialize scratch arrays and
-  // variables
-  LIW = static_cast<unsigned int>((2*n_nonzero_elements + 5*n_rows + 4) *
-                                  LIW_factor_1);
-  IW.resize (LIW);
-  KEEP.resize (n_nonzero_elements + 5*n_rows + 2);
-
-  // declare output info fields
-  bool call_succeeded;
-  do
-    {
-      call_ma47ad(&n_rows, &n_nonzero_elements,
-                  &row_numbers[0], &column_numbers[0],
-                  &IW[0], &LIW, &KEEP[0],
-                  &ICNTL[0], &INFO[0]);
-      call_succeeded = (INFO[0] == 0);
-
-      // if enough memory or no
-      // increase allowed: exit loop
-      if (call_succeeded || (LIW_increase_factor_1 <= 1))
-        break;
-
-      // otherwise: increase LIW and retry
-      LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_1);
-      IW.resize (LIW);
-    }
-  while (true);
-
-  AssertThrow (call_succeeded, ExcMA47AFailed(INFO[0]));
-
-  // note that we have already been
-  // in this function
-  initialize_called = true;
-}
-
-
-
-void
-SparseDirectMA47::factorize (const SparseMatrix<double> &m)
-{
-  Assert (factorize_called == false,
-          ExcCantFactorizeAgain());
-
-  // if necessary, initialize process
-  if (initialize_called == false)
-    initialize (m);
-
-  // make sure the matrices
-  // are the same
-  Assert (matrix == &m, ExcDifferentMatrices());
-
-
-  // set LA and fill the A array of
-  // values
-  LA = std::max (static_cast<int>(INFO[5] * LA_factor),
-                 static_cast<int>(n_nonzero_elements));
-  A.resize (LA);
-  fill_A (m);
-
-  // if necessary extend IW
-  if (LIW < INFO[6] * LIW_factor_2)
-    {
-      LIW = static_cast<unsigned int>(INFO[6] * LIW_factor_2);
-      IW.resize (LIW);
-    };
-
-  const unsigned int n_rows = m.get_sparsity_pattern().n_rows();
-  IW1.resize (2*n_rows+2);
-
-  // output info flags
-  bool call_succeeded;
-  do
-    {
-      call_ma47bd (&n_rows, &n_nonzero_elements, &column_numbers[0],
-                   &A[0], &LA,
-                   &IW[0], &LIW, &KEEP[0], &CNTL[0], &ICNTL[0],
-                   &IW1[0], &INFO[0]);
-      call_succeeded = (INFO[0] == 0);
-
-      // if enough memory or no
-      // increase allowed: exit loop
-      if (call_succeeded)
-        break;
-
-      // otherwise: increase LIW or
-      // LA if that is allowed and
-      // retry
-      switch (INFO[0])
-        {
-        case -3:
-        {
-          if (LIW_increase_factor_2 <= 1)
-            goto exit_loop;
-
-          LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_2);
-          IW.resize (LIW);
-          break;
-        };
-
-        case -4:
-        {
-          if (LA_increase_factor <= 1)
-            goto exit_loop;
-          // increase A. note that
-          // since the function has
-          // already part of the
-          // array @p{A}, we have
-          // to re-fill it with the
-          // original values. minor
-          // clue: since the old
-          // entries are no more
-          // needed, we can discard
-          // them; we use this to
-          // first release all
-          // memory (through the
-          // call to @p{swap} and
-          // the subsequent call to
-          // the destructor of the
-          // @p{tmp} object) and
-          // only then re-allocate
-          // it. If we called
-          // @p{resize} directly,
-          // this would first
-          // allocate more memory,
-          // then copy the old
-          // contents, and only
-          // then release the old
-          // memory, but keeping
-          // both memory regions at
-          // the same time could
-          // sometimes be more than
-          // we can do, leading to
-          // an exception on the
-          // allocation.
-          std::cout << "<*>" << std::flush;
-
-          LA  = static_cast<unsigned int>(LA * LA_increase_factor);
-          if (true)
-            {
-              std::vector<double> tmp;
-              A.swap (tmp);
-            };
-
-          A.resize (LA);
-          fill_A (m);
-
-          break;
-        };
-
-        // ups, other return
-        // value, don't know
-        // what to do here
-        default:
-          AssertThrow (false, ExcMA47BFailed(INFO[0]));
-        };
-      continue;
-
-exit_loop:
-      break;
-    }
-  while (true);
-
-  AssertThrow (call_succeeded, ExcMA47BFailed(INFO[0]));
-
-  // note that we have been here
-  // already
-  factorize_called = true;
-}
-
-
-
-void
-SparseDirectMA47::solve (Vector<double> &rhs_and_solution)
-{
-  Assert (factorize_called == true, ExcFactorizeNotCalled());
-
-  const unsigned int n_rows = rhs_and_solution.size();
-  call_ma47cd (&n_rows, &A[0], &LA,
-               &IW[0], &LIW,
-               &rhs_and_solution(0), &IW1[0], &ICNTL[0]);
-}
-
-
-
-void
-SparseDirectMA47::solve (const SparseMatrix<double> &matrix,
-                         Vector<double>             &rhs_and_solution)
-{
-  initialize (matrix);
-  factorize (matrix);
-  solve (rhs_and_solution);
-}
-
-
-
-std::size_t
-SparseDirectMA47::memory_consumption () const
-{
-  return (sizeof(*this) +
-          MemoryConsumption::memory_consumption (row_numbers) +
-          MemoryConsumption::memory_consumption (column_numbers) +
-          MemoryConsumption::memory_consumption (A) +
-          MemoryConsumption::memory_consumption (IW) +
-          MemoryConsumption::memory_consumption (KEEP) +
-          MemoryConsumption::memory_consumption (IW1));
-}
-
-
-
-Threads::Mutex &
-SparseDirectMA47::get_synchronisation_lock () const
-{
-  return synchronisation_lock;
-}
-
-
-
-void
-SparseDirectMA47::fill_A (const SparseMatrix<double> &matrix)
-{
-  Assert (n_nonzero_elements <= A.size(), ExcInternalError());
-
-  const SparsityPattern &sparsity_pattern = matrix.get_sparsity_pattern ();
-
-  const unsigned int n_rows = sparsity_pattern.n_rows();
-
-  unsigned int global_index = 0;
-  for (unsigned int row=0; row<n_rows; ++row)
-    for (SparseMatrix<double>::const_iterator col=matrix.begin(row);
-         col < matrix.end(row); ++col)
-      // note that the matrix must be
-      // symmetric, so only treat the
-      // upper right part
-      if ((row <= col->column()) && (col->value() != 0))
-        {
-          Assert (global_index < n_nonzero_elements, ExcInternalError());
-
-          A[global_index] = col->value();
-          ++global_index;
-
-          // make sure that the symmetric
-          // entry exists and has the same
-          // value, unless this one is zero
-          Assert ((col->value() == 0)
-                  ||
-                  (col->value() == matrix(col->column(),row)),
-                  ExcMatrixNotSymmetric());
-        }
-      else
-        // lower left part. just check
-        // symmetry
-        Assert ((col->value() == 0)
-                ||
-                (col->value() == matrix(col->column(),row)),
-                ExcMatrixNotSymmetric());
-
-  Assert (global_index == n_nonzero_elements, ExcInternalError());
-}
-
-
-
-void
-SparseDirectMA47::call_ma47id (double       *CNTL,   // length 2
-                               unsigned int *ICNTL)  // length 7
-{
-  HSL::MA47::ma47id_ (CNTL, ICNTL);
-}
-
-
-
-void
-SparseDirectMA47::
-call_ma47ad (const unsigned int *n_rows,             //scalar
-             const unsigned int *n_nonzero_elements, //scalar
-             unsigned int       *row_numbers,        //length n_nonzero
-             unsigned int       *column_numbers,     //length n_nonzero
-             unsigned int       *IW,                 //length LIW
-             const unsigned int *LIW,                //scalar
-             unsigned int       *KEEP,               //n_nonzero+5*n_rows+2
-             const unsigned int *ICNTL,              //length 7
-             int                *INFO)               //length 24
-{
-  double RINFO[4];
-  HSL::MA47::ma47ad_(n_rows, n_nonzero_elements,
-                     row_numbers, column_numbers,
-                     IW, LIW, KEEP,
-                     ICNTL, &RINFO[0], INFO);
-}
-
-
-
-void
-SparseDirectMA47::
-call_ma47bd (const unsigned int *n_rows,             //scalar
-             const unsigned int *n_nonzero_elements, //scalar
-             const unsigned int *column_numbers,     //length n_nonzero
-             double             *A,                  //length LA
-             const unsigned int *LA,                 //scalar
-             unsigned int       *IW,                 //length LIW
-             const unsigned int *LIW,                //scalar
-             const unsigned int *KEEP,               //n_nonzero+5*n_rows+2
-             const double       *CNTL,               //length 2
-             const unsigned int *ICNTL,              //length 7
-             unsigned int       *IW1,                //2*n_rows+2
-             int                *INFO)               //length 24
-{
-  double RINFO[4];
-  HSL::MA47::ma47bd_(n_rows, n_nonzero_elements, column_numbers,
-                     A, LA,
-                     IW, LIW, KEEP, CNTL, ICNTL,
-                     IW1, &RINFO[0], INFO);
-}
-
-
-
-void
-SparseDirectMA47::
-call_ma47cd (const unsigned int *n_rows,           //scalar
-             const double       *A,                //length LA
-             const unsigned int *LA,               //scalar
-             const unsigned int *IW,               //length LIW
-             const unsigned int *LIW,              //scalar
-             double             *rhs_and_solution, //length n_rows
-             unsigned int       *IW1,              //length 2*n_rows+2
-             const unsigned int *ICNTL)            //length 7
-{
-  std::vector<double> W(*n_rows);
-  HSL::MA47::ma47cd_(n_rows, A, LA,
-                     IW, LIW, &W[0],
-                     rhs_and_solution, IW1, ICNTL);
-}
-
 
 
 SparseDirectUMFPACK::~SparseDirectUMFPACK ()
@@ -1507,8 +63,7 @@ SparseDirectUMFPACK::SparseDirectUMFPACK ()
 void
 SparseDirectUMFPACK::clear ()
 {
-  // delete objects that haven't been deleted
-  // yet
+  // delete objects that haven't been deleted yet
   if (symbolic_decomposition != 0)
     {
       umfpack_dl_free_symbolic (&symbolic_decomposition);
@@ -1546,38 +101,26 @@ void
 SparseDirectUMFPACK::
 sort_arrays (const SparseMatrix<number> &matrix)
 {
-  // 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
+  // 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<matrix.m(); ++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)
+      // 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
+      // 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
       long int cursor = Ap[row];
       while ((cursor < Ap[row+1]-1) &&
              (Ai[cursor] > Ai[cursor+1]))
@@ -1617,37 +160,27 @@ void
 SparseDirectUMFPACK::
 sort_arrays (const BlockSparseMatrix<number> &matrix)
 {
-  // the case for block matrices is a
-  // bit more difficult, since all we
-  // know is that *within each
-  // block*, the diagonal of that
-  // block may come first. however,
-  // that means that there may be as
-  // many entries per row in the
-  // wrong place as there are block
-  // columns. we can do the same
-  // thing as above, but we have to
-  // do it multiple times
+  // the case for block matrices is a bit more difficult, since all we know
+  // is that *within each block*, the diagonal of that block may come
+  // first. however, that means that there may be as many entries per row
+  // in the wrong place as there are block columns. we can do the same
+  // thing as above, but we have to do it multiple times
   for (unsigned int row=0; row<matrix.m(); ++row)
     {
       long int cursor = Ap[row];
       for (unsigned int block=0; block<matrix.n_block_cols(); ++block)
         {
-
-          // find the next
-          // out-of-order element
+          // find the next out-of-order element
           while ((cursor < Ap[row+1]-1) &&
                  (Ai[cursor] < Ai[cursor+1]))
             ++cursor;
 
-          // if there is none, then
-          // just go on
+          // if there is none, then just go on
           if (cursor == Ap[row+1]-1)
             break;
 
-          // otherwise swap this entry
-          // with successive ones as
-          // long as necessary
+          // otherwise swap this entry with successive ones as long as
+          // necessary
           long int element = cursor;
           while ((element < Ap[row+1]-1) &&
                  (Ai[element] > Ai[element+1]))
@@ -1673,35 +206,24 @@ factorize (const Matrix &matrix)
 
   const unsigned int N = matrix.m();
 
-  // copy over the data from the matrix to
-  // the data structures UMFPACK wants. note
-  // two things: first, UMFPACK 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_dl_solve function with the
-  // UMFPACK_At argument, meaning that we
-  // want to solve for the transpose system
+  // copy over the data from the matrix to the data structures UMFPACK
+  // wants. note two things: first, UMFPACK 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_dl_solve function with the UMFPACK_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;
-  // UMFPACK 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
+  // second: the data we have in the sparse matrices is "almost" right
+  // already; UMFPACK 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.
+  // 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.
   Ap.resize (N+1);
   Ai.resize (matrix.n_nonzero_elements());
   Ax.resize (matrix.n_nonzero_elements());
@@ -1713,30 +235,23 @@ factorize (const Matrix &matrix)
   Assert (static_cast<unsigned int>(Ap.back()) == Ai.size(),
           ExcInternalError());
 
-  // then copy over matrix
-  // elements. note that for sparse
-  // matrices, iterators are sorted
-  // so that they traverse each row
-  // from start to end before moving
-  // on to the next row. however,
-  // this isn't true for block
-  // matrices, so we have to do a bit
-  // of book keeping
+  // then copy over matrix elements. note that for sparse matrices,
+  // iterators are sorted so that they traverse each row from start to end
+  // before moving on to the next row. however, this isn't true for block
+  // matrices, so we have to do a bit of book keeping
   {
-    // have an array that for each
-    // row points to the first entry
-    // not yet written to
+    // have an array that for each row points to the first entry not yet
+    // written to
     std::vector<long int> row_pointers = Ap;
 
-    // loop over the elements of the matrix row by row, as suggested
-    // in the documentation of the sparse matrix iterator class
+    // loop over the elements of the matrix row by row, as suggested in the
+    // documentation of the sparse matrix iterator class
     for (unsigned int row = 0; row < matrix.m(); ++row)
       {
         for (typename Matrix::const_iterator p=matrix.begin(row);
             p!=matrix.end(row); ++p)
           {
-            // write entry into the first
-            // free one for this row
+            // write entry into the first free one for this row
             Ai[row_pointers[row]] = p->column();
             Ax[row_pointers[row]] = p->value();
 
@@ -1745,17 +260,14 @@ factorize (const Matrix &matrix)
           }
       }
 
-    // at the end, we should have
-    // written all rows completely
+    // at the end, we should have written all rows completely
     for (unsigned int i=0; i<Ap.size()-1; ++i)
       Assert (row_pointers[i] == Ap[i+1], ExcInternalError());
   }
 
-  // make sure that the elements in
-  // each row are sorted. we have to
-  // be more careful for block sparse
-  // matrices, so ship this task out
-  // to a different function
+  // make sure that the elements in each row are sorted. we have to be more
+  // careful for block sparse matrices, so ship this task out to a
+  // different function
   sort_arrays (matrix);
 
   int status;
@@ -1781,8 +293,7 @@ factorize (const Matrix &matrix)
 void
 SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
 {
-  // make sure that some kind of factorize()
-  // call has happened before
+  // make sure that some kind of factorize() call has happened before
   Assert (Ap.size() != 0, ExcNotInitialized());
   Assert (Ai.size() != 0, ExcNotInitialized());
   Assert (Ai.size() == Ax.size(), ExcNotInitialized());
@@ -1790,12 +301,9 @@ SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
   Vector<double> rhs (rhs_and_solution.size());
   rhs = rhs_and_solution;
 
-  // solve the system. note that since
-  // UMFPACK wants compressed column storage
-  // instead of the compressed row storage
-  // format we use in deal.II's
-  // SparsityPattern classes, we solve for
-  // UMFPACK's A^T instead
+  // solve the system. note that since UMFPACK wants compressed column
+  // storage instead of the compressed row storage format we use in
+  // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
   const int status
     = umfpack_dl_solve (UMFPACK_At,
                         &Ap[0], &Ai[0], &Ax[0],
@@ -1855,7 +363,6 @@ SparseDirectUMFPACK::solve (const Matrix &,
   AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
 }
 
-
 #endif
 
 
@@ -1904,7 +411,10 @@ SparseDirectUMFPACK::Tvmult_add (
   Assert(false, ExcNotImplemented());
 }
 
+
+
 #ifdef DEAL_II_WITH_MUMPS
+
 SparseDirectMUMPS::SparseDirectMUMPS ()
   :
   initialize_called (false)
@@ -1953,17 +463,16 @@ void SparseDirectMUMPS::initialize_matrix (const Matrix &matrix)
       // representation of the matrix
       a   = new double[nz];
 
-      // matrix indices pointing to the row and
-      // column dimensions respectively of the
-      // matrix representation above (a): ie. a[k]
-      // is the matrix element (irn[k], jcn[k])
+      // matrix indices pointing to the row and column dimensions
+      // respectively of the matrix representation above (a): ie. a[k] is
+      // the matrix element (irn[k], jcn[k])
       irn = new int[nz];
       jcn = new int[nz];
 
       unsigned int index = 0;
 
-      // loop over the elements of the matrix row by row, as suggested
-      // in the documentation of the sparse matrix iterator class
+      // loop over the elements of the matrix row by row, as suggested in
+      // the documentation of the sparse matrix iterator class
       for (unsigned int row = 0; row < matrix.m(); ++row)
         {
           for (typename Matrix::const_iterator ptr = matrix.begin (row);
@@ -2037,12 +546,10 @@ void SparseDirectMUMPS::initialize (const Matrix &matrix)
 
 void SparseDirectMUMPS::solve (Vector<double> &vector)
 {
-  // Check that the solver has been initialized
-  // by the routine above:
+  // Check that the solver has been initialized by the routine above:
   Assert (initialize_called == true, ExcNotInitialized());
 
-  // and that the matrix has at least one
-  // nonzero element:
+  // and that the matrix has at least one nonzero element:
   Assert (nz != 0, ExcNotInitialized());
 
   // Start solver
@@ -2054,12 +561,10 @@ void SparseDirectMUMPS::solve (Vector<double> &vector)
 void SparseDirectMUMPS::vmult (Vector<double>       &dst,
                                const Vector<double> &src)
 {
-  // Check that the solver has been initialized
-  // by the routine above:
+  // Check that the solver has been initialized by the routine above:
   Assert (initialize_called == true, ExcNotInitialized());
 
-  // and that the matrix has at least one
-  // nonzero element:
+  // and that the matrix has at least one nonzero element:
   Assert (nz != 0, ExcNotInitialized());
 
   // Hand over right-hand side
@@ -2082,21 +587,6 @@ void SparseDirectMUMPS::vmult (Vector<double>       &dst,
 
 #endif // DEAL_II_WITH_MUMPS
 
-// explicit instantiations for SparseMatrixMA27
-template
-void SparseDirectMA27::factorize (const SparseMatrix<double> &matrix);
-
-template
-void SparseDirectMA27::factorize (const SparseMatrix<float> &matrix);
-
-template
-void SparseDirectMA27::solve (const SparseMatrix<double> &matrix,
-                              Vector<double>             &rhs_and_solution);
-
-template
-void SparseDirectMA27::solve (const SparseMatrix<float>  &matrix,
-                              Vector<double>             &rhs_and_solution);
-
 
 // explicit instantiations for SparseMatrixUMFPACK
 #define InstantiateUMFPACK(MATRIX) \

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.