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.
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.
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.
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.
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 everytime 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 transfered 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.
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.
SparseDirectMA27::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 | |||
) |
Constructor. See the documentation of this class for the meaning of the parameters to this function.
SparseDirectMA27::~SparseDirectMA27 | ( | ) |
Destructor.
void SparseDirectMA27::set_detached_mode | ( | ) |
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.
bool SparseDirectMA27::detached_mode_set | ( | ) | const |
Return whether the detached mode is set.
void SparseDirectMA27::initialize | ( | const SparsityPattern & | sparsity_pattern | ) |
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 SparseDirectMA27::factorize | ( | const SparseMatrix< number > & | matrix | ) | [inline] |
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.
void SparseDirectMA27::solve | ( | Vector< number > & | rhs_and_solution | ) | const [inline] |
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 SparseDirectMA27::solve | ( | const SparseMatrix< number > & | matrix, | |
Vector< double > & | rhs_and_solution | |||
) | [inline] |
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.
Return an estimate of the memory used by this class.
Threads::ThreadMutex& SparseDirectMA27::get_synchronisation_lock | ( | ) | 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.
void SparseDirectMA27::fill_A | ( | const SparseMatrix< number > & | matrix | ) | [inline, private] |
Fill the A
array from the symmetric part of the given matrix.
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 | |||
) | [private] |
Call the respective function with the given args, either locally or remote.
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 | |||
) | [private] |
Call the respective function with the given args, either locally or remote.
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 [private] |
Call the respective function with the given args, either locally or remote.
Call the respective function with the given args, either locally or remote.
Call the respective function with the given args, either locally or remote.
Call the respective function with the given args, either locally or remote.
const bool SparseDirectMA27::suppress_output [private] |
Store in the constructor whether the MA27 routines shall deliver output to stdout or not.
bool SparseDirectMA27::detached_mode [private] |
Store whether set_detached_mode() has been called.
DetachedModeData* SparseDirectMA27::detached_mode_data [private] |
Pointer to a structure that will hold the data necessary to uphold communication with a detached solver.
const double SparseDirectMA27::LIW_factor_1 [private] |
Store the three values passed to the cinstructor. See the documentation of this class for the meaning of these variables.
const double SparseDirectMA27::LIW_factor_2 [private] |
const double SparseDirectMA27::LA_factor [private] |
const double SparseDirectMA27::LIW_increase_factor_1 [private] |
Increase factors in case a call to a function fails.
const double SparseDirectMA27::LIW_increase_factor_2 [private] |
const double SparseDirectMA27::LA_increase_factor [private] |
bool SparseDirectMA27::initialize_called [private] |
Flags storing whether the first two functions have already been called.
bool SparseDirectMA27::factorize_called [private] |
SmartPointer<const SparsityPattern> SparseDirectMA27::sparsity_pattern [private] |
Store a pointer to the sparsity pattern, to make sure that we use the same thing for all calls.
unsigned int SparseDirectMA27::n_nonzero_elements [private] |
Number of nonzero elements in the sparsity pattern on and above the diagonal.
std::vector<unsigned int> SparseDirectMA27::row_numbers [private] |
Arrays holding row and column indices.
std::vector<unsigned int> SparseDirectMA27::column_numbers [private] |
std::vector<double> SparseDirectMA27::A [private] |
Array to hold the matrix elements, and later the elements of the factors.
unsigned int SparseDirectMA27::LA [private] |
Length of the A
array.
unsigned int SparseDirectMA27::LIW [private] |
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.
std::vector<unsigned int> SparseDirectMA27::IW [private] |
std::vector<unsigned int> SparseDirectMA27::IKEEP [private] |
std::vector<unsigned int> SparseDirectMA27::IW1 [private] |
unsigned int SparseDirectMA27::NSTEPS [private] |
unsigned int SparseDirectMA27::MAXFRT [private] |
unsigned int SparseDirectMA27::NRLNEC [private] |
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 SparseDirectMA27::NIRNEC [private] |
int SparseDirectMA27::IFLAG [private] |
Flag indicating the level of output desired and returning error values if error occured.
Threads::ThreadMutex SparseDirectMA27::static_synchronisation_lock [static, private] |
Mutexes for synchronising access to this class.
Threads::ThreadMutex SparseDirectMA27::non_static_synchronisation_lock [mutable, private] |