From c093abeb807f8505e95ba61551e4a6dc0964f365 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 25 Sep 2002 15:45:32 +0000 Subject: [PATCH] Implement a detached mode for the MA27 solver. git-svn-id: https://svn.dealii.org/trunk@6509 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/contrib/hsl/Makefile | 15 +- deal.II/contrib/hsl/source/detached_ma27.cc | 168 ++++++ deal.II/lac/include/lac/sparse_direct.h | 240 ++++++++- deal.II/lac/source/sparse_direct.cc | 539 +++++++++++++++++--- 4 files changed, 888 insertions(+), 74 deletions(-) create mode 100644 deal.II/contrib/hsl/source/detached_ma27.cc diff --git a/deal.II/contrib/hsl/Makefile b/deal.II/contrib/hsl/Makefile index febd87264a..195ccfd797 100644 --- a/deal.II/contrib/hsl/Makefile +++ b/deal.II/contrib/hsl/Makefile @@ -27,24 +27,31 @@ $(LIBDIR)/contrib/hsl/%.o : # rules ifeq ($(enable-shared),yes) - lib: $(LIBDIR)/libhsl.so + lib: $(LIBDIR)/libhsl.so $(LIBDIR)/bin/detached_ma27 else - lib: $(LIBDIR)/libhsl.a + lib: $(LIBDIR)/libhsl.a $(LIBDIR)/bin/detached_ma27 endif $(LIBDIR)/libhsl.a: $(forward-declarations) $(o-files) - @echo =====hsl==============optimized==$(MT)== Linking library: $@ + @echo =====hsl==============optimized==$(MT)== Linking library: $(@F) @ar ruv $@ $(o-files) # link the object files to a shared library. since this is not C++, # simply use the C compiler, rather than $(SHLIBLD) $(LIBDIR)/libhsl.so: $(forward-declarations) $(o-files) - @echo =====hsl==============optimized==$(MT)== Linking library: $@ + @echo =====hsl==============optimized==$(MT)== Linking library: $(@F) @$(CC) $(LDFLAGS) -shared -o $@ $(o-files) +# rule to make the program that runs an MA27 solver detached from the +# main program, and communicates through a pipe +$(LIBDIR)/bin/detached_ma27: source/detached_ma27.cc include/hsl/hsl.h + @echo =====hsl=========================$(MT)== Making $(@F) + @$(CXX) $(CXXFLAGS.g) $< -o $@ $(lib-contrib-hsl) $(F77LIBS) + + clean: -rm -f *~ lib/lib* lib/o/*.o lib/Makefile.dep diff --git a/deal.II/contrib/hsl/source/detached_ma27.cc b/deal.II/contrib/hsl/source/detached_ma27.cc new file mode 100644 index 0000000000..1e072c83fc --- /dev/null +++ b/deal.II/contrib/hsl/source/detached_ma27.cc @@ -0,0 +1,168 @@ +#include +#include +#include +#include + + + +template +void put (const T *t, const size_t N, const char *info) +{ + write (1, + reinterpret_cast (t), + sizeof(T) * N); + fflush (NULL); +}; + + +template +void get (T *t, const size_t N, const char *info) +{ + unsigned int count = 0; + while (count < sizeof(T)*N) + { + int ret = read (0, + reinterpret_cast (t) + count, + sizeof(T) * N - count); + if (ret < 0) + { + abort (); + std::cerr << "------ error " << ret << " on client side!" + << std::endl; + } + else + count += ret; + }; +}; + + + +int main () +{ + unsigned int N, NZ, NSTEPS, LA, MAXFRT, LIW; + int IFLAG; + std::vector IRN, ICN, IW, IKEEP, IW1; + std::vector A; + + while (true) + { + char action; + get(&action, 1, "ACTION"); + + switch (action) + { + case '1': + { + get (&N, 1, "N"); + get (&NZ, 1, "NZ"); + + IRN.resize (NZ); + ICN.resize (NZ); + + get (&IRN[0], NZ, "IRN"); + get (&ICN[0], NZ, "ICN"); + get (&LIW, 1, "LIW"); + get (&IFLAG, 1, "IFLAG"); + + IW.resize (LIW); + IKEEP.resize (3*N); + IW1.resize (2*N); + + // next call function + HSL::MA27::ma27ad_ (&N, &NZ, &IRN[0], &ICN[0], &IW[0], &LIW, + &IKEEP[0], &IW1[0], &NSTEPS, &IFLAG); + + // then return IFLAG + put (&IFLAG, 1, "IFLAG"); + + sleep (1); + break; + }; + + case '2': + { + get (&LA, 1, "LA"); + A.resize (LA); + get (&A[0], LA, "A"); + + HSL::MA27::ma27bd_ (&N, &NZ, &IRN[0], &ICN[0], &A[0], &LA, + &IW[0], &LIW, &IKEEP[0], + &NSTEPS, &MAXFRT, &IW1[0], &IFLAG); + + // if IFLAG==0, then the + // call succeeded and we + // won't need ICN, IRN, + // and IKEEP any more + if (IFLAG==0) + { + std::vector tmp1, tmp2, tmp3; + ICN.swap (tmp1); + IRN.swap (tmp2); + IKEEP.swap (tmp3); + }; + + // finally return IFLAG + put (&IFLAG, 1, "IFLAG"); + + break; + }; + + + case '3': + { + std::vector W(MAXFRT); + std::vector rhs(N); + get (&rhs[0], N, "RHS"); + + HSL::MA27::ma27cd_ (&N, &A[0], &LA, &IW[0], + &LIW, &W[0], &MAXFRT, &rhs[0], + &IW1[0], &NSTEPS); + + put (&rhs[0], N, "RHS"); + + break; + }; + + + case '4': + { + unsigned int NRLNEC; + HSL::MA27::ma27x1_ (&NRLNEC); + put (&NRLNEC, 1, "NRLNEC"); + break; + }; + + case '5': + { + unsigned int NIRNEC; + HSL::MA27::ma27x2_ (&NIRNEC); + put (&NIRNEC, 1, "NIRNEC"); + break; + }; + + case '6': + { + unsigned int LP; + get (&LP, 1, "LP"); + HSL::MA27::ma27x3_ (&LP); + break; + }; + + case '7': + { + exit (0); + break; + }; + + default: + { + std::cerr << "Invalid action key ('" << action + << "'=" << static_cast(action) << ")!!" + << std::endl; + abort(); + }; + }; + }; +}; + + diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h index 49a882b215..d995aeb506 100644 --- a/deal.II/lac/include/lac/sparse_direct.h +++ b/deal.II/lac/include/lac/sparse_direct.h @@ -99,8 +99,10 @@ * * @sect3{Note on parallelization} * + * @sect4{Synchronisation} + * * Due to the use of global variables through COMMON blocks, the calls - * to the sparse direct solver routines is not multithreading-capable, + * 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 @ref{Threads} @@ -120,7 +122,36 @@ * be performed from outside, it is left to the user of this class to * do so. * - * @author Wolfgang Bangerth, 2000, 2001 + * @sect4{Detached mode} + * + * As an alternative, you can call the function @p{set_detached_mode} + * right after calling the constructor. This lets the program fork, so + * that we now have two programs that communicate via pipes. 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 + * @p{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. + * + * @author Wolfgang Bangerth, 2000, 2001, 2002 */ class SparseDirectMA27 : public Subscriptor { @@ -139,6 +170,34 @@ class SparseDirectMA27 : public Subscriptor 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 @p{initialize} + * (or the two-argument @p{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 @@ -265,8 +324,43 @@ class SparseDirectMA27 : public Subscriptor * Exception */ DeclException0 (ExcDifferentSparsityPatterns); - + 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 + * @p{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 @@ -366,10 +460,11 @@ class SparseDirectMA27 : public Subscriptor int IFLAG; /** - * Mutex for synchronising access + * Mutexes for synchronising access * to this class. */ - static Threads::ThreadMutex synchronisation_lock; + static Threads::ThreadMutex static_synchronisation_lock; + mutable Threads::ThreadMutex non_static_synchronisation_lock; /** * Fill the @p{A} array from the @@ -377,6 +472,77 @@ class SparseDirectMA27 : public Subscriptor * matrix. */ 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); }; @@ -445,6 +611,9 @@ class SparseDirectMA27 : public Subscriptor * 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. + * * @author Wolfgang Bangerth, 2000, 2001 */ class SparseDirectMA47 : public Subscriptor @@ -606,7 +775,15 @@ class SparseDirectMA47 : public Subscriptor DeclException0 (ExcDifferentMatrices); 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 @@ -705,6 +882,57 @@ class SparseDirectMA47 : public Subscriptor * matrix. */ void fill_A (const SparseMatrix &matrix); + + /** + * Call the @p{ma47id} function + * with the given args. + */ + void call_ma47id (double *CNTL, + unsigned int *ICNTL); + + /** + * Call the @p{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 @p{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 *IW, + int *INFO); + + /** + * Call the @p{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); }; diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc index c9213159a6..4f1d72339f 100644 --- a/deal.II/lac/source/sparse_direct.cc +++ b/deal.II/lac/source/sparse_direct.cc @@ -18,6 +18,8 @@ #include #include +#include +#include // if we know that at least one of the HSL functions are there, @@ -197,8 +199,64 @@ namespace HSL /* -------------------------- MA27 ---------------------------- */ +Threads::ThreadMutex SparseDirectMA27::static_synchronisation_lock; + + +struct SparseDirectMA27::DetachedModeData +{ + /** + * Mutex to assure that only one + * thread is currently talking + * through the pipe. + */ + Threads::ThreadMutex 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; + + template + void put (const T *t, const size_t N, const char *info) const + { + write (server_client_pipe[1], + reinterpret_cast (t), + sizeof(T) * N); + fflush (NULL); + }; + + template + void get (T *t, const size_t N, const char *info) const + { + unsigned int count = 0; + while (count < sizeof(T)*N) + { + int ret = read (client_server_pipe[0], + reinterpret_cast (t) + count, + sizeof(T) * N - count); + if (ret < 0) + { + abort (); + std::cerr << "****** error " << ret << " on server side!" + << std::endl; + } + else + count += ret; + }; + }; + +}; -Threads::ThreadMutex SparseDirectMA27::synchronisation_lock; SparseDirectMA27::SparseDirectMA27 (const double LIW_factor_1, @@ -207,7 +265,11 @@ SparseDirectMA27::SparseDirectMA27 (const double LIW_factor_1, const double LIW_increase_factor_1, const double LIW_increase_factor_2, const double LA_increase_factor, - const bool suppress_output) : + 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), @@ -217,14 +279,42 @@ SparseDirectMA27::SparseDirectMA27 (const double LIW_factor_1, initialize_called (false), factorize_called (false), sparsity_pattern (0) +{}; + + + +SparseDirectMA27::~SparseDirectMA27() { - // suppress error output if - // requested - if (suppress_output) - { - const unsigned int LP = 0; - HSL::MA27::ma27x3_ (&LP); - }; + if (detached_mode) + if (detached_mode_data != 0) + { + + // first close down client + detached_mode_data->mutex.acquire (); + write (detached_mode_data->server_client_pipe[1], "7", 1); + detached_mode_data->mutex.acquire (); + + // 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; }; @@ -234,14 +324,80 @@ 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 + pipe(detached_mode_data->server_client_pipe); + pipe(detached_mode_data->client_server_pipe); + // fflush(NULL) is said to be a + // good idea before fork() + fflush(NULL); + + // now fork and create child + // process + detached_mode_data->child_pid = fork(); + 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 + 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 + const char * const program_name = DEAL_II_PATH"/lib/bin/detached_ma27"; + const char * const child_argv[] = { program_name, 0 }; + 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!")); + exit(1); + }; + // parent process continues here. + // close unneeded end of pipe + }; + + + // 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(); - const unsigned int - *rowstart_indices = sparsity_pattern->get_rowstart_indices(); - const unsigned int - *col_nums = sparsity_pattern->get_column_numbers(); + n_rows = sparsity_pattern->n_rows(); + const unsigned int * const + rowstart_indices = sparsity_pattern->get_rowstart_indices(); + const unsigned int * const + col_nums = sparsity_pattern->get_column_numbers(); // first count number of nonzero // elements in the upper right @@ -287,9 +443,9 @@ SparseDirectMA27::initialize (const SparsityPattern &sp) // variables LIW = static_cast((2*n_nonzero_elements + 3*n_rows + 1) * LIW_factor_1); - IW.resize (LIW); - IKEEP.resize (3*n_rows); - IW1.resize (2*n_rows); + 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; @@ -301,10 +457,10 @@ SparseDirectMA27::initialize (const SparsityPattern &sp) bool call_succeeded = true; do { - HSL::MA27::ma27ad_(&n_rows, &n_nonzero_elements, - &row_numbers[0], &column_numbers[0], - &IW[0], &LIW, &IKEEP[0], - &IW1[0], &NSTEPS, &IFLAG); + 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 @@ -328,8 +484,8 @@ SparseDirectMA27::initialize (const SparsityPattern &sp) // COMMON block. we need these // values in order to set array // sizes in the next function - HSL::MA27::ma27x1_(&NRLNEC); - HSL::MA27::ma27x2_(&NIRNEC); + call_ma27x1 (&NRLNEC); + call_ma27x2 (&NIRNEC); // note that we have already been // in this function @@ -373,17 +529,27 @@ SparseDirectMA27::factorize (const SparseMatrix &matrix) bool call_succeeded = true; do { - HSL::MA27::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_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 + // increase allowed: exit + // loop. delete data that is no + // more used if (call_succeeded) - break; + { + 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 @@ -479,10 +645,9 @@ SparseDirectMA27::solve (Vector &rhs_and_solution) const Assert (factorize_called == true, ExcFactorizeNotCalled()); const unsigned int n_rows = rhs_and_solution.size(); - std::vector W(MAXFRT); - HSL::MA27::ma27cd_(&n_rows, &A[0], &LA, - &IW[0], &LIW, &W[0], &MAXFRT, - &rhs_and_solution(0), &IW1[0], &NSTEPS); + call_ma27cd (&n_rows, &A[0], &LA, + &IW[0], &LIW, &MAXFRT, + &rhs_and_solution(0), &IW1[0], &NSTEPS); }; @@ -515,7 +680,10 @@ SparseDirectMA27::memory_consumption () const Threads::ThreadMutex & SparseDirectMA27::get_synchronisation_lock () const { - return synchronisation_lock; + if (detached_mode) + return non_static_synchronisation_lock; + else + return static_synchronisation_lock; }; @@ -550,6 +718,176 @@ SparseDirectMA27::fill_A (const SparseMatrix &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) +{ + if (detached_mode_set() == false) + HSL::MA27::ma27ad_ (N, NZ, IRN, ICN, IW, LIW, + IKEEP, IW1, NSTEPS, IFLAG); + else + { + detached_mode_data->mutex.acquire (); + // 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 (IW == 0, ExcInternalError()); + Assert (IKEEP == 0, ExcInternalError()); + Assert (IW == 0, ExcInternalError()); + + // next get back what we need + // to know + detached_mode_data->get (IFLAG, 1, "IFLAG"); + + detached_mode_data->mutex.release (); + }; +}; + + + +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 + detached_mode_data->mutex.acquire (); + 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"); + + detached_mode_data->mutex.release (); + }; +}; + + + +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 + { + detached_mode_data->mutex.acquire (); + // 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"); + detached_mode_data->mutex.release (); + }; +}; + + + +void SparseDirectMA27::call_ma27x2 (unsigned int *NIRNEC) +{ + if (detached_mode_set() == false) + HSL::MA27::ma27x2_ (NIRNEC); + else + { + detached_mode_data->mutex.acquire (); + // 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"); + detached_mode_data->mutex.release (); + }; +}; + + + +void SparseDirectMA27::call_ma27x3 (const unsigned int *LP) +{ + if (detached_mode_set() == false) + HSL::MA27::ma27x3_ (LP); + else + { + detached_mode_data->mutex.acquire (); + // 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"); + detached_mode_data->mutex.release (); + }; +}; + + + + + /* -------------------------- MA47 ---------------------------- */ Threads::ThreadMutex SparseDirectMA47::synchronisation_lock; @@ -561,7 +899,9 @@ SparseDirectMA47::SparseDirectMA47 (const double LIW_factor_1, const double LIW_increase_factor_1, const double LIW_increase_factor_2, const double LA_increase_factor, - const bool suppress_output) : + const bool suppress_output) + : + suppress_output (suppress_output), LIW_factor_1 (LIW_factor_1), LIW_factor_2 (LIW_factor_2), LA_factor (LA_factor), @@ -571,14 +911,7 @@ SparseDirectMA47::SparseDirectMA47 (const double LIW_factor_1, initialize_called (false), factorize_called (false), matrix (0) -{ - HSL::MA47::ma47id_ (CNTL, ICNTL); - - // suppress error output if - // requested - if (suppress_output) - ICNTL[0] = 0; -}; +{}; @@ -587,12 +920,21 @@ 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(); - const unsigned int *rowstart_indices = sparsity_pattern.get_rowstart_indices(); - const unsigned int *col_nums = sparsity_pattern.get_column_numbers(); + const unsigned int + n_rows = sparsity_pattern.n_rows(); + const unsigned int * const + rowstart_indices = sparsity_pattern.get_rowstart_indices(); + const unsigned int * const + col_nums = sparsity_pattern.get_column_numbers(); // first count number of nonzero // elements in the upper right @@ -644,14 +986,13 @@ SparseDirectMA47::initialize (const SparseMatrix &m) KEEP.resize (n_nonzero_elements + 5*n_rows + 2); // declare output info fields - double RINFO[4]; bool call_succeeded; do { - HSL::MA47::ma47ad_(&n_rows, &n_nonzero_elements, - &row_numbers[0], &column_numbers[0], - &IW[0], &LIW, &KEEP[0], - &ICNTL[0], &RINFO[0], &INFO[0]); + 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 @@ -691,7 +1032,7 @@ SparseDirectMA47::factorize (const SparseMatrix &m) // set LA and fill the A array of // values - LA=static_cast(INFO[5] * LA_factor); + LA = static_cast(INFO[5] * LA_factor); A.resize (LA); fill_A (m); @@ -706,15 +1047,13 @@ SparseDirectMA47::factorize (const SparseMatrix &m) IW1.resize (2*n_rows+2); // output info flags - double RINFO[4]; - bool call_succeeded; do - { - HSL::MA47::ma47bd_(&n_rows, &n_nonzero_elements, &column_numbers[0], - &A[0], &LA, - &IW[0], &LIW, &KEEP[0], &CNTL[0], &ICNTL[0], - &IW1[0], &RINFO[0], &INFO[0]); + { + 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 @@ -816,10 +1155,9 @@ SparseDirectMA47::solve (Vector &rhs_and_solution) Assert (factorize_called == true, ExcFactorizeNotCalled()); const unsigned int n_rows = rhs_and_solution.size(); - std::vector W(n_rows); - HSL::MA47::ma47cd_(&n_rows, &A[0], &LA, - &IW[0], &LIW, &W[0], - &rhs_and_solution(0), &IW1[0], &ICNTL[0]); + call_ma47cd (&n_rows, &A[0], &LA, + &IW[0], &LIW, + &rhs_and_solution(0), &IW1[0], &ICNTL[0]); }; @@ -883,3 +1221,76 @@ SparseDirectMA47::fill_A (const SparseMatrix &matrix) }; 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); +}; -- 2.39.5