From 567972692e1e70f1c38f58c66845dd1380b3d909 Mon Sep 17 00:00:00 2001 From: maier Date: Sun, 31 Mar 2013 20:38:00 +0000 Subject: [PATCH] Remove support for the obsolete MA27 and MA47 direct solvers of the HSL collection git-svn-id: https://svn.dealii.org/trunk@29122 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 + deal.II/include/deal.II/lac/sparse_direct.h | 1206 +------------- deal.II/source/lac/sparse_direct.cc | 1648 +------------------ 3 files changed, 156 insertions(+), 2705 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 3bb7acff64..a5cfafa161 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -25,6 +25,13 @@ inconvenience this causes.
    +
  1. Removed: The interfaces to the obsolete direct solvers MA27 and MA47 from +the Harwell Subroutine Library. Support for the HSL routines were not ported to +the new build system +
    +(Matthias Maier, 2013/04/01) + +
  2. Changed: The TimeDependent::end_sweep function with an argument indicating the number of threads has been removed. Use the corresponding function without an argument. Since the argument had a default value, few users will have used diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h index 5a5d9d673e..195ca99190 100644 --- a/deal.II/include/deal.II/lac/sparse_direct.h +++ b/deal.II/include/deal.II/lac/sparse_direct.h @@ -30,986 +30,6 @@ 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 Harwell Subroutine - * Library 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. - * - * - *

    Interface and Method

    - * - * For the meaning of the three functions initialize(), factorize(), - * and solve(), as well as for the method used in MA27, please see the - * documentation - * 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. - * - * - *

    Parameters to the constructor

    - * - * 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 LIW1, LIW2, and LA, - * where LIW1 denotes the size of the IW array in the call to - * MA27A, while LIW2 is the array size in the call to - * MA27B. 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 at least 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 LA, but also for LIW1). 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 LIW is only - * changed in the second call if the recommended value times - * LIW_factor_2 is larger than the array size already is from the - * call to MA27A; otherwise, LIW_factor_2 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 MA27C function writes out an error message - * if the value of LA 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. - * - * - *

    Note on parallelization

    - * - *

    Synchronisation

    - * - * 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. - * - *

    Detached mode

    - * - * 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 detached_ma27, that is started in its place - * through the execv 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. - * - * - *
    Internals of the detached mode
    - * - * The program that actually runs the detached solver is called - * detached_ma27, 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 ps 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 - void factorize (const SparseMatrix &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 - void solve (Vector &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 - void solve (const SparseMatrix &matrix, - Vector &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 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 row_numbers; - std::vector column_numbers; - - /** - * Array to hold the matrix - * elements, and later the - * elements of the factors. - */ - std::vector A; - - /** - * Length of the A 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 IW; - std::vector IKEEP; - std::vector 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 A array from the - * symmetric part of the given - * matrix. - */ - template - void fill_A (const SparseMatrix &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 Harwell - * Subroutine Library 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. - * - * - *

    Interface and Method

    - * - * For the meaning of the three functions initialize(), factorize(), - * and solve(), as well as for the method used in MA47, please see the - * documentation - * 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. - * - * - *

    Parameters to the constructor

    - * - * 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. - * - * - *

    Note on parallelization

    - * - * 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 - * MA47ID 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 &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 &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 &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 &matrix, - Vector &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,SparseDirectMA47> matrix; - - /** - * Number of nonzero elements in - * the sparsity pattern on and - * above the diagonal. - */ - unsigned int n_nonzero_elements; - - /** - * Control values set by MA47ID. - */ - 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 row_numbers; - std::vector column_numbers; - - /** - * Array to hold the matrix - * elements, and later the - * elements of the factors. - */ - std::vector A; - - /** - * Length of the A 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 IW; - std::vector KEEP; - std::vector IW1; - - /** - * Mutex for synchronising access - * to this class. - */ - static Threads::Mutex synchronisation_lock; - - /** - * Fill the A array from the - * symmetric part of the given - * matrix. - */ - void fill_A (const SparseMatrix &matrix); - - /** - * Call the ma47id function - * with the given args. - */ - void call_ma47id (double *CNTL, - unsigned int *ICNTL); - - /** - * Call the ma47ad 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 ma47bd 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 ma47bd function - * with the given args. - */ - void call_ma47cd (const unsigned int *n_rows, - const double *A, - const unsigned int *LA, - const unsigned int *IW, - const unsigned int *LIW, - double *rhs_and_solution, - unsigned int *IW1, - const unsigned int *ICNTL); -}; - - /** @@ -1051,8 +71,7 @@ private: * * There are instantiations of this class for SparseMatrix, * SparseMatrix, SparseMatrixEZ, SparseMatrixEZ, - * BlockSparseMatrix, and - * BlockSparseMatrix. + * BlockSparseMatrix, and BlockSparseMatrix. * * @ingroup Solvers Preconditioners * @@ -1062,19 +81,16 @@ class SparseDirectUMFPACK : public Subscriptor { public: /** - * Dummy class needed for the - * usual initalization interface - * of preconditioners. + * Dummy class needed for the usual initalization interface of + * preconditioners. */ class AdditionalData {}; /** - * Constructor. See the - * documentation of this class - * for the meaning of the - * parameters to this function. + * Constructor. See the documentation of this class for the meaning of + * the parameters to this function. */ SparseDirectUMFPACK (); @@ -1084,136 +100,94 @@ public: ~SparseDirectUMFPACK (); /** - * This function does nothing. It is only - * here to provide an interface that is - * consistent with that of the HSL MA27 - * and MA47 solver classes. + * This function does nothing. It is only here to provide a consistent + * interface. */ void initialize (const SparsityPattern &sparsity_pattern); /** - * Factorize the matrix. This function - * may be called multiple times for - * different matrices, after the object - * of this class has been initialized for - * a certain sparsity pattern. You may - * therefore save some computing time if - * you want to invert several matrices - * with the same sparsity - * pattern. However, note that the bulk - * of the computing time is actually - * spent in the factorization, so this - * functionality may not always be of - * large benefit. + * Factorize the matrix. This function may be called multiple times for + * different matrices, after the object of this class has been + * initialized for a certain sparsity pattern. You may therefore save + * some computing time if you want to invert several matrices with the + * same sparsity pattern. However, note that the bulk of the computing + * time is actually spent in the factorization, so this functionality may + * not always be of large benefit. * - * In contrast to the other direct solver - * classes, the initialisation method - * does nothing. Therefore initialise - * is not automatically called by this - * method, when - * the initialization step has - * not been performed yet. + * In contrast to the other direct solver classes, the initialisation + * method does nothing. Therefore initialise is not automatically called + * by this method, when the initialization step has not been performed + * yet. * - * This function copies the contents of - * the matrix into its own storage; the - * matrix can therefore be deleted after - * this operation, even if subsequent - * solves are required. + * This function copies the contents of the matrix into its own storage; + * the matrix can therefore be deleted after this operation, even if + * subsequent solves are required. */ template void factorize (const Matrix &matrix); /** - * Initialize memory and call - * SparseDirectUMFPACK::factorize. + * Initialize memory and call SparseDirectUMFPACK::factorize. */ template 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 Ax - * = b. 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 Ax = b. 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 &, const Vector &) const; /** - * Not implemented but necessary - * for compiling. + * Not implemented but necessary for compiling. */ void Tvmult (Vector &, const Vector &) 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 &, const Vector &) const; /** - * Not implemented but necessary - * for compiling. + * Not implemented but necessary for compiling. */ void Tvmult_add (Vector &, const Vector &) 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 &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 void solve (const Matrix &matrix, Vector &rhs_and_solution); /** - * One of the UMFPack routines - * threw an error. The error code - * is included in the output and - * can be looked up in the - * UMFPack user manual. The name - * of the routine is included for - * reference. + * One of the UMFPack routines threw an error. The error code is included + * in the output and can be looked up in the UMFPack user manual. The + * name of the routine is included for reference. */ DeclException2 (ExcUMFPACKError, char *, int, << "UMFPACK routine " << arg1 @@ -1223,33 +197,24 @@ public: private: /** - * The UMFPACK routines allocate objects - * in which they store information about - * symbolic and numeric values of the - * decomposition. The actual data type of - * these objects is opaque, and only - * passed around as void pointers. + * The UMFPACK routines allocate objects in which they store information + * about symbolic and numeric values of the decomposition. The actual + * data type of these objects is opaque, and only passed around as void + * pointers. */ void *symbolic_decomposition; void *numeric_decomposition; /** - * Free all memory that hasn't been freed - * yet. + * Free all memory that hasn't been freed yet. */ void clear (); - /** - * Make sure that the arrays Ai - * and Ap are sorted in each - * row. UMFPACK wants it this - * way. We need to have three - * versions of this function, one - * for the usual SparseMatrix, one - * for the SparseMatrixEZ, and - * one for the BlockSparseMatrix - * classes + * Make sure that the arrays Ai and Ap are sorted in each row. UMFPACK + * wants it this way. We need to have three versions of this function, + * one for the usual SparseMatrix, one for the SparseMatrixEZ, and one + * for the BlockSparseMatrix classes */ template void sort_arrays (const SparseMatrixEZ &); @@ -1261,32 +226,29 @@ private: void sort_arrays (const BlockSparseMatrix &); /** - * The arrays in which we store the data - * for the solver. + * The arrays in which we store the data for the solver. */ std::vector Ap; std::vector Ai; std::vector Ax; /** - * Control and info arrays for the solver - * routines. + * Control and info arrays for the solver routines. */ std::vector control; }; /** - * This class provides an interface to the parallel sparse direct - * solver MUMPS. 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 + * MUMPS. 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 MUMPS exists on your - * system and was detected during configuration of - * deal.II. + * @note This class is useable if and only if a working installation of MUMPS exists on your system and was + * detected during configuration of deal.II. * *

    Instantiations

    * @@ -1312,23 +274,20 @@ private: unsigned int nz; /** - * This function initializes a MUMPS instance - * and hands over the system's matrix - * matrix. + * This function initializes a MUMPS instance and hands over the system's + * matrix matrix. */ template 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 &vector); /** - * Flags storing whether the function - * initialize () has already been - * called. + * Flags storing whether the function initialize () has already + * been called. */ bool initialize_called; @@ -1350,35 +309,30 @@ public: DeclException0 (ExcInitializeAlreadyCalled); /** - * This function initializes a MUMPS instance - * and hands over the system's matrix - * matrix and right-hand side - * vector to the solver. + * This function initializes a MUMPS instance and hands over the system's + * matrix matrix and right-hand side vector to the + * solver. */ template void initialize (const Matrix &matrix, const Vector &vector); /** - * This function initializes a MUMPS instance - * and computes the factorization of the - * system's matrix matrix. + * This function initializes a MUMPS instance and computes the + * factorization of the system's matrix matrix. */ template void initialize (const Matrix &matrix); /** - * A function in which the linear system is - * solved and the solution vector is copied - * into the given vector. + * A function in which the linear system is solved and the solution + * vector is copied into the given vector. */ void solve (Vector &vector); /** - * A function in which the inverse of the - * matrix is applied to the input vector - * src and the solution is - * written into the output vector + * A function in which the inverse of the matrix is applied to the input + * vector src and the solution is written into the output vector * dst. */ void vmult (Vector &dst, const Vector &src); diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc index 2b10e88c2d..1328eb20ce 100644 --- a/deal.II/source/lac/sparse_direct.cc +++ b/deal.II/source/lac/sparse_direct.cc @@ -28,1455 +28,11 @@ 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 -#endif - // include UMFPACK file. #ifdef DEAL_II_WITH_UMFPACK # include #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 - 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 - 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 (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 - 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 (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(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; rowbegin(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; rowbegin(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((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(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 -void -SparseDirectMA27::factorize (const SparseMatrix &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(NRLNEC * LA_factor), - static_cast(n_nonzero_elements)); - A.resize (LA); - fill_A (matrix); - - // if necessary extend IW - if (LIW < NIRNEC * LIW_factor_2) - { - LIW = static_cast(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 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(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(LA * LA_increase_factor); - if (true) - { - std::vector 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 &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 &rhs_and_solution) const -{ - Assert (factorize_called == true, ExcFactorizeNotCalled()); - - // first have to convert data type to - // doubles - Vector 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 -void -SparseDirectMA27::solve (const SparseMatrix &matrix, - Vector &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 -void -SparseDirectMA27::fill_A (const SparseMatrix &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::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 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 &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::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::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((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(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 &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(INFO[5] * LA_factor), - static_cast(n_nonzero_elements)); - A.resize (LA); - fill_A (m); - - // if necessary extend IW - if (LIW < INFO[6] * LIW_factor_2) - { - LIW = static_cast(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(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(LA * LA_increase_factor); - if (true) - { - std::vector 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 &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 &matrix, - Vector &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 &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::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 W(*n_rows); - HSL::MA47::ma47cd_(n_rows, A, LA, - IW, LIW, &W[0], - rhs_and_solution, IW1, ICNTL); -} - SparseDirectUMFPACK::~SparseDirectUMFPACK () @@ -1507,8 +63,7 @@ SparseDirectUMFPACK::SparseDirectUMFPACK () void SparseDirectUMFPACK::clear () { - // delete objects that haven't been deleted - // yet + // delete objects that haven't been deleted yet if (symbolic_decomposition != 0) { umfpack_dl_free_symbolic (&symbolic_decomposition); @@ -1546,38 +101,26 @@ void SparseDirectUMFPACK:: sort_arrays (const SparseMatrix &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 Ai[cursor+1])) @@ -1617,37 +160,27 @@ void SparseDirectUMFPACK:: sort_arrays (const BlockSparseMatrix &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 Ai[element+1])) @@ -1673,35 +206,24 @@ factorize (const Matrix &matrix) const unsigned int N = matrix.m(); - // copy over the data from the matrix to - // the data structures UMFPACK wants. note - // two things: first, UMFPACK wants - // compressed column storage whereas we - // always do compressed row storage; we - // work around this by, rather than - // shuffling things around, copy over the - // data we have, but then call the - // umfpack_dl_solve function with the - // UMFPACK_At argument, meaning that we - // want to solve for the transpose system + // copy over the data from the matrix to the data structures UMFPACK + // wants. note two things: first, UMFPACK wants compressed column storage + // whereas we always do compressed row storage; we work around this by, + // rather than shuffling things around, copy over the data we have, but + // then call the umfpack_dl_solve function with the UMFPACK_At argument, + // meaning that we want to solve for the transpose system // - // second: the data we have in the sparse - // matrices is "almost" right already; - // UMFPACK wants the entries in each row - // (i.e. really: column) to be sorted in - // ascending order. we almost have that, - // except that we usually store the - // diagonal first in each row to allow for - // some optimizations. thus, we have to - // resort things a little bit, but only + // second: the data we have in the sparse matrices is "almost" right + // already; UMFPACK wants the entries in each row (i.e. really: column) + // to be sorted in ascending order. we almost have that, except that we + // usually store the diagonal first in each row to allow for some + // optimizations. thus, we have to resort things a little bit, but only // within each row // - // final note: if the matrix has entries in - // the sparsity pattern that are actually - // occupied by entries that have a zero - // numerical value, then we keep them - // anyway. people are supposed to provide - // accurate sparsity patterns. + // final note: if the matrix has entries in the sparsity pattern that are + // actually occupied by entries that have a zero numerical value, then we + // keep them anyway. people are supposed to provide accurate sparsity + // patterns. Ap.resize (N+1); Ai.resize (matrix.n_nonzero_elements()); Ax.resize (matrix.n_nonzero_elements()); @@ -1713,30 +235,23 @@ factorize (const Matrix &matrix) Assert (static_cast(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 row_pointers = Ap; - // loop over the elements of the matrix row by row, as suggested - // in the documentation of the sparse matrix iterator class + // loop over the elements of the matrix row by row, as suggested in the + // documentation of the sparse matrix iterator class for (unsigned int row = 0; row < matrix.m(); ++row) { for (typename Matrix::const_iterator p=matrix.begin(row); p!=matrix.end(row); ++p) { - // write entry into the first - // free one for this row + // write entry into the first free one for this row Ai[row_pointers[row]] = p->column(); Ax[row_pointers[row]] = p->value(); @@ -1745,17 +260,14 @@ factorize (const Matrix &matrix) } } - // at the end, we should have - // written all rows completely + // at the end, we should have written all rows completely for (unsigned int i=0; i &rhs_and_solution) const { - // make sure that some kind of factorize() - // call has happened before + // make sure that some kind of factorize() call has happened before Assert (Ap.size() != 0, ExcNotInitialized()); Assert (Ai.size() != 0, ExcNotInitialized()); Assert (Ai.size() == Ax.size(), ExcNotInitialized()); @@ -1790,12 +301,9 @@ SparseDirectUMFPACK::solve (Vector &rhs_and_solution) const Vector rhs (rhs_and_solution.size()); rhs = rhs_and_solution; - // solve the system. note that since - // UMFPACK wants compressed column storage - // instead of the compressed row storage - // format we use in deal.II's - // SparsityPattern classes, we solve for - // UMFPACK's A^T instead + // solve the system. note that since UMFPACK wants compressed column + // storage instead of the compressed row storage format we use in + // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead const int status = umfpack_dl_solve (UMFPACK_At, &Ap[0], &Ai[0], &Ax[0], @@ -1855,7 +363,6 @@ SparseDirectUMFPACK::solve (const Matrix &, AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); } - #endif @@ -1904,7 +411,10 @@ SparseDirectUMFPACK::Tvmult_add ( Assert(false, ExcNotImplemented()); } + + #ifdef DEAL_II_WITH_MUMPS + SparseDirectMUMPS::SparseDirectMUMPS () : initialize_called (false) @@ -1953,17 +463,16 @@ void SparseDirectMUMPS::initialize_matrix (const Matrix &matrix) // representation of the matrix a = new double[nz]; - // matrix indices pointing to the row and - // column dimensions respectively of the - // matrix representation above (a): ie. a[k] - // is the matrix element (irn[k], jcn[k]) + // matrix indices pointing to the row and column dimensions + // respectively of the matrix representation above (a): ie. a[k] is + // the matrix element (irn[k], jcn[k]) irn = new int[nz]; jcn = new int[nz]; unsigned int index = 0; - // loop over the elements of the matrix row by row, as suggested - // in the documentation of the sparse matrix iterator class + // loop over the elements of the matrix row by row, as suggested in + // the documentation of the sparse matrix iterator class for (unsigned int row = 0; row < matrix.m(); ++row) { for (typename Matrix::const_iterator ptr = matrix.begin (row); @@ -2037,12 +546,10 @@ void SparseDirectMUMPS::initialize (const Matrix &matrix) void SparseDirectMUMPS::solve (Vector &vector) { - // Check that the solver has been initialized - // by the routine above: + // Check that the solver has been initialized by the routine above: Assert (initialize_called == true, ExcNotInitialized()); - // and that the matrix has at least one - // nonzero element: + // and that the matrix has at least one nonzero element: Assert (nz != 0, ExcNotInitialized()); // Start solver @@ -2054,12 +561,10 @@ void SparseDirectMUMPS::solve (Vector &vector) void SparseDirectMUMPS::vmult (Vector &dst, const Vector &src) { - // Check that the solver has been initialized - // by the routine above: + // Check that the solver has been initialized by the routine above: Assert (initialize_called == true, ExcNotInitialized()); - // and that the matrix has at least one - // nonzero element: + // and that the matrix has at least one nonzero element: Assert (nz != 0, ExcNotInitialized()); // Hand over right-hand side @@ -2082,21 +587,6 @@ void SparseDirectMUMPS::vmult (Vector &dst, #endif // DEAL_II_WITH_MUMPS -// explicit instantiations for SparseMatrixMA27 -template -void SparseDirectMA27::factorize (const SparseMatrix &matrix); - -template -void SparseDirectMA27::factorize (const SparseMatrix &matrix); - -template -void SparseDirectMA27::solve (const SparseMatrix &matrix, - Vector &rhs_and_solution); - -template -void SparseDirectMA27::solve (const SparseMatrix &matrix, - Vector &rhs_and_solution); - // explicit instantiations for SparseMatrixUMFPACK #define InstantiateUMFPACK(MATRIX) \ -- 2.39.5