]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
improvements in copy_from, conjugate_add
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 11 Oct 2002 12:13:12 +0000 (12:13 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 11 Oct 2002 12:13:12 +0000 (12:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@6631 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix_ez.h

index c892790bd8466deae96100751c13d3dce0e1c50d..340b79ccb6b67e10252389a7cc09b4a3f58801a8 100644 (file)
@@ -644,16 +644,16 @@ class SparseMatrixEZ : public Subscriptor
                                      * Obviously, the matrix needs to
                                      * be square for this operation.
                                      */
-    template <typename somenumber>
-    somenumber matrix_norm_square (const Vector<somenumber> &v) const;
+//    template <typename somenumber>
+//    somenumber matrix_norm_square (const Vector<somenumber> &v) const;
 
                                     /**
                                      * Compute the matrix scalar
                                      * product $\left(u,Mv\right)$.
                                      */
-    template <typename somenumber>
-    somenumber matrix_scalar_product (const Vector<somenumber> &u,
-                                     const Vector<somenumber> &v) const;
+//    template <typename somenumber>
+//    somenumber matrix_scalar_product (const Vector<somenumber> &u,
+//                                   const Vector<somenumber> &v) const;
 
                                     /**
                                      * Frobenius-norm of the matrix.
@@ -731,6 +731,15 @@ class SparseMatrixEZ : public Subscriptor
                            const Vector<somenumber> &src,
                            const number              om = 1.) const;
 
+                                    /**
+                                     * Add the matrix @p{A}
+                                     * conjugated by @p{B}, that is,
+                                     * $B A B^T$ to this object.
+                                     */
+    template <class MATRIXA, class MATRIXB>
+    void conjugate_add (const MATRIXA& A,
+                       const MATRIXB& B);
+    
                                     /**
                                      * STL-like iterator with the
                                      * first entry.
@@ -774,7 +783,7 @@ class SparseMatrixEZ : public Subscriptor
                                      * nonzero entry of the matrix
                                      * per line.
                                      */
-    void print (std::ostream &out) const;
+//    void print (std::ostream &out) const;
 
                                     /**
                                      * Print the matrix in the usual
@@ -815,12 +824,12 @@ class SparseMatrixEZ : public Subscriptor
                                      * may produce @em{large} amounts of
                                      * output if applied to a large matrix!
                                      */
-    void print_formatted (std::ostream       &out,
-                         const unsigned int  precision   = 3,
-                         const bool          scientific  = true,
-                         const unsigned int  width       = 0,
-                         const char         *zero_string = " ",
-                         const double        denominator = 1.) const;
+//      void print_formatted (std::ostream       &out,
+//                       const unsigned int  precision   = 3,
+//                       const bool          scientific  = true,
+//                       const unsigned int  width       = 0,
+//                       const char         *zero_string = " ",
+//                       const double        denominator = 1.) const;
 
                                     /**
                                      * Determine an estimate for the
@@ -914,24 +923,6 @@ class SparseMatrixEZ : public Subscriptor
                                         const unsigned int        end_row,
                                         somenumber               *partial_sum) const;
 
-                                    /**
-                                     * Version of @p{residual} which
-                                     * only performs its actions on
-                                     * the region defined by
-                                     * @p{[begin_row,end_row)} (these
-                                     * numbers are the components of
-                                     * @p{interval}). This function is
-                                     * called by @p{residual} in the
-                                     * case of enabled
-                                     * multithreading.
-                                     */
-    template <typename somenumber>
-    void threaded_residual (Vector<somenumber>       &dst,
-                           const Vector<somenumber> &u,
-                           const Vector<somenumber> &b,
-                           const std::pair<unsigned int,unsigned int> interval,
-                           somenumber               *partial_norm) const;
-
                                     /**
                                      * Number of columns. This is
                                      * used to check vector
@@ -1346,7 +1337,8 @@ SparseMatrixEZ<number>::copy_from (const MATRIX& M)
 
   while (start != final)
     {
-      set(start->row(), start->column(), start->value());
+      if (start->value() != 0.)
+       set(start->row(), start->column(), start->value());
       ++start;
     }
   return *this;
@@ -1361,16 +1353,52 @@ SparseMatrixEZ<number>::add_scaled (number factor,
 {
   Assert (M.m() == m(), ExcDimensionMismatch(M.m(), m()));
   Assert (M.n() == n(), ExcDimensionMismatch(M.n(), n()));
+
+  if (factor == 0.)
+    return;
   
   typename MATRIX::const_iterator start = M.begin();
   const typename MATRIX::const_iterator final = M.end();
 
   while (start != final)
     {
-      add(start->row(), start->column(), factor * start->value());
+      if (start->value() != 0.)
+       add(start->row(), start->column(), factor * start->value());
       ++start;
     }
 }
 
+
+template<typename number>
+template <class MATRIXA, class MATRIXB>
+void
+SparseMatrixEZ<number>::conjugate_add (const MATRIXA& A,
+                                      const MATRIXB& B)
+{
+// Compute the result
+// r_ij = \sum_kl b_ik b_jl a_kl
+  
+  typename MATRIXB::const_iterator b1 = B.begin();
+  typename MATRIXB::const_iterator b2 = B.begin();
+  const typename MATRIXB::const_iterator b_final = B.end();
+  while (b1 != b_final)
+    {
+      const unsigned int i = b->row();
+      const unsigned int k = b->column();
+      while (b2 != b_final)
+       {
+         const unsigned int j = b->row();
+         const unsigned int l = b->column();
+         
+         const typename MATRIXA::value_type a = A.el(k,l);
+         
+         if (a != 0.)
+           add (i, j, a * b1->value() * b2->value());
+         ++b2;
+       }
+      ++b1;
+    }
+}
+
 #endif
 /*----------------------------   sparse_matrix.h     ---------------------------*/

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.