]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a detached mode for the MA27 solver.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 Sep 2002 15:45:32 +0000 (15:45 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 Sep 2002 15:45:32 +0000 (15:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@6509 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/contrib/hsl/Makefile
deal.II/contrib/hsl/source/detached_ma27.cc [new file with mode: 0644]
deal.II/lac/include/lac/sparse_direct.h
deal.II/lac/source/sparse_direct.cc

index febd87264ac03a00c8563b92b5f8ffe06ef53cc1..195ccfd79727444a4bfe024352efb9b78bc8f4b9 100644 (file)
@@ -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 (file)
index 0000000..1e072c8
--- /dev/null
@@ -0,0 +1,168 @@
+#include <cstdio>
+#include <unistd.h>
+#include <hsl/hsl.h>
+#include <vector>
+
+
+
+template <typename T>
+void put (const T *t, const size_t N, const char *info)
+{
+  write (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)
+{
+  unsigned int count = 0;
+  while (count < sizeof(T)*N)
+    {
+      int ret = read (0,
+                      reinterpret_cast<char *> (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<unsigned int> IRN, ICN, IW, IKEEP, IW1;
+  std::vector<double> 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<unsigned int> tmp1, tmp2, tmp3;
+                ICN.swap (tmp1);
+                IRN.swap (tmp2);
+                IKEEP.swap (tmp3);
+              };
+            
+                                             // finally return IFLAG
+            put (&IFLAG, 1, "IFLAG");
+
+            break;
+          };
+          
+
+          case '3':
+          {
+            std::vector<double> W(MAXFRT);
+            std::vector<double> 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<unsigned int>(action) << ")!!"
+                      << std::endl;
+            abort();
+          };
+        };
+    };
+};
+
+                
index 49a882b215dfc09d7dcac2c4396f5a9c36f06a66..d995aeb506c7c6d7a8a0ee392374d655293fb6c8 100644 (file)
  *
  * @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
 {
@@ -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<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);
 };
 
 
@@ -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<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);
 };
 
 
index c9213159a6a1f57a0ee7db7e3e43f85624d2ecc7..4f1d72339f3e6abbec7122c90d6fdb68fac249cc 100644 (file)
@@ -18,6 +18,8 @@
 #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,
@@ -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 <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,
@@ -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<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
@@ -287,9 +443,9 @@ SparseDirectMA27::initialize (const SparsityPattern &sp)
                                   // 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;
@@ -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<double> &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<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
@@ -479,10 +645,9 @@ SparseDirectMA27::solve (Vector<double> &rhs_and_solution) const
   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);
 };
 
 
@@ -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<double> &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<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;
@@ -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<double> &m)
 {
   Assert (initialize_called == false, ExcInitializeAlreadyCalled());
 
+                                   // some initialization stuff
+  call_ma47id (CNTL, ICNTL);
+  if (suppress_output)
+    ICNTL[0] = 0;
+
+                                   // then start with work
   matrix = &m;
   const SparsityPattern &sparsity_pattern = matrix->get_sparsity_pattern();
   
-  const unsigned int n_rows = sparsity_pattern.n_rows();
-  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<double> &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<double> &m)
   
                                   // 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);
   
@@ -706,15 +1047,13 @@ SparseDirectMA47::factorize (const SparseMatrix<double> &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<double> &rhs_and_solution)
   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]);
 };
 
 
@@ -883,3 +1221,76 @@ SparseDirectMA47::fill_A (const SparseMatrix<double> &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<double> W(*n_rows);
+  HSL::MA47::ma47cd_(n_rows, A, LA,
+                    IW, LIW, &W[0],
+                    rhs_and_solution, IW1, ICNTL);  
+};

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.