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);
-};
-
-
/**
*
* 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
*
{
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 ();
~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
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> &);
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>
*
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;
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);
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 ()
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);
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]))
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]))
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());
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();
}
}
- // 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;
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());
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],
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
Assert(false, ExcNotImplemented());
}
+
+
#ifdef DEAL_II_WITH_MUMPS
+
SparseDirectMUMPS::SparseDirectMUMPS ()
:
initialize_called (false)
// 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);
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
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
#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) \
projects / dealii-svn.git / commitdiff