*
* @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}
* 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
{
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
* 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
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
* matrix.
*/
void fill_A (const SparseMatrix<double> &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);
};
* 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
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
* matrix.
*/
void fill_A (const SparseMatrix<double> &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);
};
#include <lac/vector.h>
#include <iostream>
+#include <sys/wait.h>
+#include <unistd.h>
// if we know that at least one of the HSL functions are there,
/* -------------------------- 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 <typename T>
+ void put (const T *t, const size_t N, const char *info) const
+ {
+ write (server_client_pipe[1],
+ reinterpret_cast<const char *> (t),
+ sizeof(T) * N);
+ fflush (NULL);
+ };
+
+ template <typename T>
+ 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<char *> (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,
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),
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;
};
{
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<char * const *>(child_argv));
+
+ // usually execv does not
+ // return. if it does, then an
+ // error happened and we report it
+ // herewith:
+ AssertThrow (false,
+ ExcMessage ("execv returned, which it is not supposed to do!"));
+ 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
// variables
LIW = static_cast<unsigned int>((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;
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
// 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
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<unsigned int> tmp1, tmp2, tmp3;
+ row_numbers.swap (tmp1);
+ column_numbers.swap (tmp2);
+ IKEEP.swap (tmp3);
+
+ break;
+ };
+
// otherwise: increase LIW or
// LA if that is allowed and
Assert (factorize_called == true, ExcFactorizeNotCalled());
const unsigned int n_rows = rhs_and_solution.size();
- std::vector<double> 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);
};
Threads::ThreadMutex &
SparseDirectMA27::get_synchronisation_lock () const
{
- return synchronisation_lock;
+ if (detached_mode)
+ return non_static_synchronisation_lock;
+ else
+ return static_synchronisation_lock;
};
+
+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<double> W(*MAXFRT);
+ HSL::MA27::ma27cd_ (N, A, LA, IW, LIW, &W[0], MAXFRT, RHS, IW1, NSTEPS);
+ }
+ else
+ {
+ detached_mode_data->put ("3", 1, "ACTION 3");
+
+ // we only have to push and get
+ // the rhs vector
+ detached_mode_data->put (RHS, *N, "RHS");
+ detached_mode_data->get (RHS, *N, "RHS");
+ };
+};
+
+
+
+void SparseDirectMA27::call_ma27x1 (unsigned int *NRLNEC)
+{
+ if (detached_mode_set() == false)
+ HSL::MA27::ma27x1_ (NRLNEC);
+ else
+ {
+ 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;
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),
initialize_called (false),
factorize_called (false),
matrix (0)
-{
- HSL::MA47::ma47id_ (CNTL, ICNTL);
-
- // suppress error output if
- // requested
- if (suppress_output)
- ICNTL[0] = 0;
-};
+{};
{
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
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
// set LA and fill the A array of
// values
- LA=static_cast<int>(INFO[5] * LA_factor);
+ LA = static_cast<int>(INFO[5] * LA_factor);
A.resize (LA);
fill_A (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
Assert (factorize_called == true, ExcFactorizeNotCalled());
const unsigned int n_rows = rhs_and_solution.size();
- std::vector<double> 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]);
};
};
Assert (global_index == n_nonzero_elements, ExcInternalError());
};
+
+
+
+void
+SparseDirectMA47::call_ma47id (double *CNTL, // length 2
+ unsigned int *ICNTL) // length 7
+{
+ HSL::MA47::ma47id_ (CNTL, ICNTL);
+};
+
+
+
+void
+SparseDirectMA47::
+call_ma47ad (const unsigned int *n_rows, //scalar
+ const unsigned int *n_nonzero_elements, //scalar
+ unsigned int *row_numbers, //length n_nonzero
+ unsigned int *column_numbers, //length n_nonzero
+ unsigned int *IW, //length LIW
+ const unsigned int *LIW, //scalar
+ unsigned int *KEEP, //n_nonzero+5*n_rows+2
+ const unsigned int *ICNTL, //length 7
+ int *INFO) //length 24
+{
+ double RINFO[4];
+ HSL::MA47::ma47ad_(n_rows, n_nonzero_elements,
+ row_numbers, column_numbers,
+ IW, LIW, KEEP,
+ ICNTL, &RINFO[0], INFO);
+};
+
+
+
+void
+SparseDirectMA47::
+call_ma47bd (const unsigned int *n_rows, //scalar
+ const unsigned int *n_nonzero_elements, //scalar
+ const unsigned int *column_numbers, //length n_nonzero
+ double *A, //length LA
+ const unsigned int *LA, //scalar
+ unsigned int *IW, //length LIW
+ const unsigned int *LIW, //scalar
+ const unsigned int *KEEP, //n_nonzero+5*n_rows+2
+ const double *CNTL, //length 2
+ const unsigned int *ICNTL, //length 7
+ unsigned int *IW1, //2*n_rows+2
+ int *INFO) //length 24
+{
+ double RINFO[4];
+ HSL::MA47::ma47bd_(n_rows, n_nonzero_elements, column_numbers,
+ A, LA,
+ IW, LIW, KEEP, CNTL, ICNTL,
+ IW1, &RINFO[0], INFO);
+};
+
+
+
+void
+SparseDirectMA47::
+call_ma47cd (const unsigned int *n_rows, //scalar
+ const double *A, //length LA
+ const unsigned int *LA, //scalar
+ const unsigned int *IW, //length LIW
+ const unsigned int *LIW, //scalar
+ double *rhs_and_solution, //length n_rows
+ unsigned int *IW1, //length 2*n_rows+2
+ const unsigned int *ICNTL) //length 7
+{
+ std::vector<double> W(*n_rows);
+ HSL::MA47::ma47cd_(n_rows, A, LA,
+ IW, LIW, &W[0],
+ rhs_and_solution, IW1, ICNTL);
+};