Single out the function that computes the inverses in a first step to parallelize...

git-svn-id: https://svn.dealii.org/trunk@2295 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/lac/include/lac/sparse_vanka.h b/deal.II/lac/include/lac/sparse_vanka.h
index 718f47fa2a4aeb151bf63562dd73c578db06830c..c33ae98435d1c2f6414804289ef5086438e76aa3 100644 (file)
--- a/deal.II/lac/include/lac/sparse_vanka.h
+++ b/deal.II/lac/include/lac/sparse_vanka.h
@@ -11,6 +11,7 @@
 #include <lac/forward_declarations.h>
 
 #include <vector>
+#include <map>
 
 
 /**
@@ -203,6 +204,23 @@ class SparseVanka
                                      * selected diagonal elements.
                                      */
     void compute_inverses ();
+
+                                    /**
+                                     * Compute the inverse of the
+                                     * block located at position
+                                     * #row#. Since the map is used
+                                     * quite often, it is generated
+                                     * only once in the caller of
+                                     * this function and passed to
+                                     * this function which first
+                                     * clears it. Reusing the map
+                                     * makes the process
+                                     * significantly faster than in
+                                     * the case where this function
+                                     * re-creates it each time.
+                                     */
+    void compute_inverse (const unsigned int               row,
+                         map<unsigned int, unsigned int> &local_index);
 };
 
 

diff --git a/deal.II/lac/include/lac/sparse_vanka.templates.h b/deal.II/lac/include/lac/sparse_vanka.templates.h
index 8ed7e4e1ea9acb73498316c99615790ce32a34dc..71f97593ca00ded9b3ebfa2a31c6d7775320dcff 100644 (file)
--- a/deal.II/lac/include/lac/sparse_vanka.templates.h
+++ b/deal.II/lac/include/lac/sparse_vanka.templates.h
@@ -46,82 +46,91 @@ template <typename number>
 void
 SparseVanka<number>::compute_inverses () 
 {
-                                  // first define an alias to the sparsity
-                                  // pattern of the matrix, since this
-                                  // will be used quite often
-  const SparsityPattern &structure
-    = matrix->get_sparsity_pattern();
-
   map<unsigned int, unsigned int> local_index;
 
                                   // traverse all rows of the matrix
                                   // which are selected
   for (unsigned int row=0; row< matrix->m() ; ++row)
     if (selected[row] == true)
-      {
-       const unsigned int row_length = structure.row_length(row);
-       inverses[row] = new FullMatrix<float> (row_length,
-                                              row_length); 
-                                        // mapping between:
-                                        // 1 column number of all
-                                        //   entries in this row, and
-                                        // 2 the position within this
-                                        //   row (as stored in the
-                                        //   SparsityPattern object
-                                        //
-                                        // since we do not explicitely
-                                        // consider nonsysmmetric sparsity
-                                        // patterns, the first element
-                                        // of each entry simply denotes
-                                        // all degrees of freedom that
-                                        // couple with #row#.
-       local_index.clear ();
-       for (unsigned int i=0; i<row_length; ++i)
-         local_index.insert(pair<unsigned int, unsigned int>
-                            (structure.column_number(row, i), i));
-       
-                                        // Build local matrix and rhs
-       for (map<unsigned int, unsigned int>::const_iterator is=local_index.begin();
-            is!=local_index.end(); ++is)
-         {
-                                            // irow loops over all DoFs that
-                                            // couple with the present DoF
-           const unsigned int irow = is->first;
-                                            // index of DoF irow in the matrix
-                                            // row corresponding to DoF #row#.
-                                            // runs between 0 and row_length
-           const unsigned int i = is->second;
-                                            // number of DoFs coupling to
-                                            // irow (including irow itself)
-           const unsigned int irow_length = structure.row_length(irow);
-           
-                                            // for all the DoFs that irow
-                                            // couples with
-           for (unsigned int j=0; j<irow_length; ++j)
-             {
-                                                // col is the number of
-                                                // this dof
-               const unsigned int col = structure.column_number(irow, j);
-                                                // find out whether this DoF
-                                                // (that couples with #irow#,
-                                                // which itself couples with
-                                                // #row#) also couples with
-                                                // #row#.
-               const map<unsigned int, unsigned int>::const_iterator js
-                 = local_index.find(col);
+      compute_inverse (row, local_index);
+};
 
-               if (js != local_index.end())
-                 (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j);
-             };
-         };
-       
-                                        // Compute new values
-       inverses[row]->gauss_jordan();
-      };
+
+
+template <typename number>
+void
+SparseVanka<number>::compute_inverse (const unsigned int               row,
+                                     map<unsigned int, unsigned int> &local_index) 
+{
+                                  // first define an alias to the sparsity
+                                  // pattern of the matrix, since this
+                                  // will be used quite often
+  const SparsityPattern &structure
+    = matrix->get_sparsity_pattern();
+
+  const unsigned int row_length = structure.row_length(row);
+  inverses[row] = new FullMatrix<float> (row_length,
+                                        row_length); 
+                                  // mapping between:
+                                  // 1 column number of all
+                                  //   entries in this row, and
+                                  // 2 the position within this
+                                  //   row (as stored in the
+                                  //   SparsityPattern object
+                                  //
+                                  // since we do not explicitely
+                                  // consider nonsysmmetric sparsity
+                                  // patterns, the first element
+                                  // of each entry simply denotes
+                                  // all degrees of freedom that
+                                  // couple with #row#.
+  local_index.clear ();
+  for (unsigned int i=0; i<row_length; ++i)
+    local_index.insert(pair<unsigned int, unsigned int>
+                      (structure.column_number(row, i), i));
+  
+                                  // Build local matrix and rhs
+  for (map<unsigned int, unsigned int>::const_iterator is=local_index.begin();
+       is!=local_index.end(); ++is)
+    {
+                                      // irow loops over all DoFs that
+                                      // couple with the present DoF
+      const unsigned int irow = is->first;
+                                      // index of DoF irow in the matrix
+                                      // row corresponding to DoF #row#.
+                                      // runs between 0 and row_length
+      const unsigned int i = is->second;
+                                      // number of DoFs coupling to
+                                      // irow (including irow itself)
+      const unsigned int irow_length = structure.row_length(irow);
+      
+                                      // for all the DoFs that irow
+                                      // couples with
+      for (unsigned int j=0; j<irow_length; ++j)
+       {
+                                          // col is the number of
+                                          // this dof
+         const unsigned int col = structure.column_number(irow, j);
+                                          // find out whether this DoF
+                                          // (that couples with #irow#,
+                                          // which itself couples with
+                                          // #row#) also couples with
+                                          // #row#.
+         const map<unsigned int, unsigned int>::const_iterator js
+           = local_index.find(col);
+         
+         if (js != local_index.end())
+           (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j);
+       };
+    };
+  
+                                  // Compute new values
+  inverses[row]->gauss_jordan();
 };
 
 
 
+
 template<typename number>
 template<typename number2>
 void