//---------------------------------------------------------------------------
// $Id$
//
--// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
++// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
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:
- /**
- * Declare type for container size.
- */
- typedef types::global_dof_index size_type;
-
- /**
- * 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.
- */
- size_type n_nonzero_elements;
-
- /**
- * Arrays holding row and column
- * indices.
- */
- std::vector<size_type> row_numbers;
- std::vector<size_type> 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.
- */
- size_type 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.
- */
- size_type LIW;
- std::vector<size_type> IW;
- std::vector<size_type> IKEEP;
- std::vector<size_type> IW1;
-
- size_type NSTEPS;
- size_type 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.
- */
- size_type NRLNEC;
- size_type 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 size_type *N,
- const size_type *NZ,
- const size_type *IRN,
- const size_type *ICN,
- size_type *IW,
- const size_type *LIW,
- size_type *IKEEP,
- size_type *IW1,
- size_type *NSTEPS,
- int *IFLAG);
-
- /**
- * Call the respective function
- * with the given args, either
- * locally or remote.
- */
- void call_ma27bd (const size_type *N,
- const size_type *NZ,
- const size_type *IRN,
- const size_type *ICN,
- double *A,
- const size_type *LA,
- size_type *IW,
- const size_type *LIW,
- const size_type *IKEEP,
- const size_type *NSTEPS,
- size_type *MAXFRT,
- size_type *IW1,
- int *IFLAG);
-
- /**
- * Call the respective function
- * with the given args, either
- * locally or remote.
- */
- void call_ma27cd (const size_type *N,
- const double *A,
- const size_type *LA,
- const size_type *IW,
- const size_type *LIW,
- const size_type *MAXFRT,
- double *RHS,
- const size_type *IW1,
- const size_type *NSTEPS) const;
-
- /**
- * Call the respective function
- * with the given args, either
- * locally or remote.
- */
- void call_ma27x1 (size_type *NRLNEC);
-
- /**
- * Call the respective function
- * with the given args, either
- * locally or remote.
- */
- void call_ma27x2 (size_type *NIRNEC);
-
- /**
- * Call the respective function
- * with the given args, either
- * locally or remote.
- */
- void call_ma27x3 (const size_type *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:
- /**
- * Declare type for container size.
- */
- typedef types::global_dof_index size_type;
-
- /**
- * 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.
- */
- size_type 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<size_type> row_numbers;
- std::vector<size_type> 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.
- */
- size_type 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.
- */
- size_type LIW;
- std::vector<size_type> IW;
- std::vector<size_type> KEEP;
- std::vector<size_type> 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 size_type *n_rows,
- const size_type *n_nonzero_elements,
- size_type *row_numbers,
- size_type *column_numbers,
- size_type *IW,
- const size_type *LIW,
- size_type *KEEP,
- const unsigned int *ICNTL,
- int *INFO);
-
- /**
- * Call the <tt>ma47bd</tt> function
- * with the given args.
- */
- void call_ma47bd (const size_type *n_rows,
- const size_type *n_nonzero_elements,
- const size_type *column_numbers,
- double *A,
- const size_type *LA,
- size_type *IW,
- const size_type *LIW,
- const size_type *KEEP,
- const double *CNTL,
- const unsigned int *ICNTL,
- size_type *IW1,
- int *INFO);
-
- /**
- * Call the <tt>ma47bd</tt> function
- * with the given args.
- */
- void call_ma47cd (const size_type *n_rows,
- const double *A,
- const size_type *LA,
- const size_type *IW,
- const size_type *LIW,
- double *rhs_and_solution,
- size_type *IW1,
- const unsigned int *ICNTL);
- };
-
-
--
--
/**
* This class provides an interface to the sparse direct solver
* UMFPACK (see <a
class SparseDirectUMFPACK : public Subscriptor
{
public:
+ /**
+ * Declare type for container size.
+ */
+ typedef types::global_dof_index size_type;
+
/**
- * Dummy class needed for the
- * usual initalization interface
- * of preconditioners.
+ * Dummy class needed for the usual initalization interface of
+ * preconditioners.
*/
class AdditionalData
{};
DMUMPS_STRUC_C id;
#endif // DEAL_II_WITH_MUMPS
- double *a;
- double *rhs;
- int *irn;
- int *jcn;
- unsigned int n;
- unsigned int nz;
+ double *a;
+ double *rhs;
+ int *irn;
+ int *jcn;
+ types::global_dof_index n;
+ types::global_dof_index 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);
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)
+ for (size_type 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]))
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)
+ for (size_type row=0; row<matrix.m(); ++row)
{
long int cursor = Ap[row];
- for (unsigned int block=0; block<matrix.n_block_cols(); ++block)
+ for (size_type 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;
clear ();
- const unsigned int N = matrix.m();
+ const size_type 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());
// first fill row lengths array
Ap[0] = 0;
- for (unsigned int row=1; row<=N; ++row)
+ for (size_type row=1; row<=N; ++row)
Ap[row] = Ap[row-1] + matrix.get_row_length(row-1);
- Assert (static_cast<unsigned int>(Ap.back()) == Ai.size(),
+ Assert (static_cast<size_type>(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 (size_type row = 0; row < matrix.m(); ++row)
{
for (typename Matrix::const_iterator p=matrix.begin(row);
p!=matrix.end(row); ++p)
}
}
- // 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)
+ for (size_type i=0; i<Ap.size()-1; ++i)
Assert (row_pointers[i] == Ap[i+1], ExcInternalError());
}
irn = new int[nz];
jcn = new int[nz];
- unsigned int index = 0;
+ size_type 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 (size_type row = 0; row < matrix.m(); ++row)
{
for (typename Matrix::const_iterator ptr = matrix.begin (row);
ptr != matrix.end (row); ++ptr)