]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Compute inverses at beginning.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 17 Jan 2000 20:21:25 +0000 (20:21 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 17 Jan 2000 20:21:25 +0000 (20:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@2241 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_vanka.h
deal.II/lac/include/lac/sparse_vanka.templates.h

index 1191b6ea0c0eadf20fd1dbd3d7498c9d8561d6d7..718f47fa2a4aeb151bf63562dd73c578db06830c 100644 (file)
@@ -98,7 +98,7 @@
  * whether this might pose some problems in the inversion of the local matrices.
  * Maybe someone would like to check this.
  *
- * @author Guido Kanschat, documentation mostly by Wolfgang Bangerth; 1999
+ * @author Guido Kanschat, documentation and extensions by Wolfgang Bangerth; 1999, 2000
  */
 template<typename number>
 class SparseVanka
@@ -114,16 +114,50 @@ class SparseVanka
                                      * longer than the Vanka
                                      * object. The same is true for
                                      * the matrix.
+                                     *
+                                     * The matrix #M# which is passed
+                                     * here may or may not be the
+                                     * same matrix for which this
+                                     * object shall act as
+                                     * preconditioner. In particular,
+                                     * it is conceivable that the
+                                     * preconditioner is build up for
+                                     * one matrix once, but is used
+                                     * for subsequent steps in a
+                                     * nonlinear process as well,
+                                     * where the matrix changes in
+                                     * each step slightly.
+                                     *
+                                     * If #conserve_mem# is #false#,
+                                     * then the inverses of the local
+                                     * systems are computed now, if
+                                     * the flag is #true#, then they
+                                     * are computed every time the
+                                     * preconditioner is
+                                     * applied. This saves some
+                                     * memory, but makes
+                                     * preconditioning very
+                                     * slow. Note also, that if the
+                                     * flag is #false#, the the
+                                     * contents of the matrix #M# at
+                                     * the time of calling this
+                                     * constructor are used, while if
+                                     * the flag is #true#, then the
+                                     * values in #M# at the time of
+                                     * preconditioning are used. This
+                                     * may lead to different results,
+                                     * obviously, of #M# changes.
                                      */
     SparseVanka(const SparseMatrix<number> &M,
-               const vector<bool>         &selected);
+               const vector<bool>         &selected,
+               const bool                  conserve_memory = false);
     
                                     /**
                                      * Destructor.
                                      * Delete all allocated matrices.
                                      */
     ~SparseVanka();
-    
+
                                     /**
                                      * Do the preconditioning.
                                      * This function takes the residual
@@ -134,18 +168,6 @@ class SparseVanka
     void operator() (Vector<number2>       &dst,
                     const Vector<number2> &src) const;
 
-                                    /**
-                                     * Minimize memory consumption.
-                                     * Activating this option reduces
-                                     * memory needs of the Vanka object
-                                     * to nearly zero. You pay for this
-                                     * by a high increase of computing
-                                     * time, since all local matrices
-                                     * are built up and inverted every
-                                     * time the Vanka operator is applied.
-                                     */
-    void conserve_memory();
-
                                     /**
                                      * Exception
                                      */
@@ -175,6 +197,12 @@ class SparseVanka
                                      * that are tagged in #selected#.
                                      */
     mutable vector<SmartPointer<FullMatrix<float> > > inverses;
+
+                                    /**
+                                     * Compute the inverses of all
+                                     * selected diagonal elements.
+                                     */
+    void compute_inverses ();
 };
 
 
index 3591447074f9f7e61e409f492ce431e5f8bb8b50..8b5ad7ca044209182b1d1621f4fc81fd60352b0a 100644 (file)
 
 template<typename number>
 SparseVanka<number>::SparseVanka(const SparseMatrix<number> &M,
-                                const vector<bool>         &selected)
+                                const vector<bool>         &selected,
+                                const bool                  conserve_mem)
                :
                matrix(&M),
                selected(selected),
-               conserve_mem(false),
+               conserve_mem(conserve_mem),
                inverses(M.m(),0)
 {
   Assert (M.m() == M.n(),
          ExcMatrixNotSquare ());
+
+  if (conserve_mem == false)
+    compute_inverses ();
 }
 
 
@@ -28,9 +32,9 @@ template<typename number>
 SparseVanka<number>::~SparseVanka()
 {
   vector<SmartPointer<FullMatrix<float> > >::iterator i;
-  for(i=inverses.begin();i!=inverses.end();++i)
+  for(i=inverses.begin(); i!=inverses.end(); ++i)
     {
-      FullMatrix<float>p = *i;
+      FullMatrix<float> *p = *i;
       *i = 0;
       if (p != 0) delete p;
     }
@@ -38,6 +42,95 @@ SparseVanka<number>::~SparseVanka()
 
 
 
+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 SparseMatrixStruct &structure
+    = matrix->get_sparsity_pattern();
+
+                                  // 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
+                                        //   sparsematrixstruct 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#.
+       map<unsigned int, unsigned int> local_index;
+       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);
+           
+                                            // copy rhs
+//         b(i) = src(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 not, then still use
+                                                // this dof to modify the rhs
+                                                //
+                                                // note that if so, we already
+                                                // have copied the entry above
+               if (js == local_index.end())
+;//b(i) -= matrix->raw_entry(irow,j) * dst(col);
+               else
+                                                    // if so, then build the
+                                                    // matrix out of it
+                 (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j);
+             };
+         };
+       
+                                        // Compute new values
+       inverses[row]->gauss_jordan();
+      };
+};
+
+
+
 template<typename number>
 template<typename number2>
 void
@@ -68,24 +161,8 @@ SparseVanka<number>::operator ()(Vector<number2>       &dst,
   for (unsigned int row=0; row< matrix->m() ; ++row)
     if (selected[row] == true)
       {
-                                        // shall we reconstruct the system
-                                        // of equations for this DoF?
-       bool build_matrix = true;
        const unsigned int row_length = structure.row_length(row);
        
-       if (!conserve_mem)
-         {
-           if (inverses[row] == 0)
-                                              // inverse not yet built
-             {
-               FullMatrix<float> *p = new FullMatrix<float>;
-               inverses[row] = p;
-             } else {
-                                                // inverse already built
-               build_matrix = false;
-             }
-         }
-       
                                         // if we don't store the
                                         // inverse matrices, then alias
                                         // the entry in the global
@@ -95,10 +172,7 @@ SparseVanka<number>::operator ()(Vector<number2>       &dst,
          {
            inverses[row] = &local_matrix;
            inverses[row]->reinit (row_length, row_length);
-         }
-       else
-         if (build_matrix)
-           inverses[row]->reinit (row_length, row_length);
+         };
        
        b.reinit (row_length);
        x.reinit (row_length);
@@ -163,14 +237,16 @@ SparseVanka<number>::operator ()(Vector<number2>       &dst,
                else
                                                     // if so, then build the
                                                     // matrix out of it
-                 if (build_matrix)
+                 if (conserve_mem == true)
                    (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j);
              };
          };
        
                                         // Compute new values
-       if (build_matrix)
+       if (conserve_mem == true)
          inverses[row]->gauss_jordan();
+
+                                        // apply preconditioner
        inverses[row]->vmult(x,b);
        
                                         // Distribute new values

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.