Allow to run MA27 with float matrices as well.

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

diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h
index 32c63739ec1625af0fee4995ebb49079ffbd99cc..be8f691f7c185b87593759dd3c2bd8cb1ea1d570 100644 (file)
--- a/deal.II/lac/include/lac/sparse_direct.h
+++ b/deal.II/lac/include/lac/sparse_direct.h
@@ -266,7 +266,8 @@ class SparseDirectMA27 : public Subscriptor
                                      * called at the beginning of
                                      * this function.
                                      */
-    void factorize (const SparseMatrix<double> &matrix);
+    template <typename number>
+    void factorize (const SparseMatrix<number> &matrix);
 
                                     /**
                                      * Solve for a certain right hand
@@ -293,7 +294,8 @@ class SparseDirectMA27 : public Subscriptor
                                      * called, since we have no
                                      * access to the actual matrix.
                                      */
-    void solve (Vector<double> &rhs_and_solution) const;
+    template <typename number>
+    void solve (Vector<number> &rhs_and_solution) const;
 
                                     /**
                                      * Call the three functions above
@@ -306,7 +308,8 @@ class SparseDirectMA27 : public Subscriptor
                                      * in place of the right hand
                                      * side vector.
                                      */
-    void solve (const SparseMatrix<double> &matrix,
+    template <typename number>
+    void solve (const SparseMatrix<number> &matrix,
                Vector<double>             &rhs_and_solution);
 
                                     /**
@@ -510,7 +513,8 @@ class SparseDirectMA27 : public Subscriptor
                                      * symmetric part of the given
                                      * matrix.
                                      */
-    void fill_A (const SparseMatrix<double> &matrix);
+    template <typename number>    
+    void fill_A (const SparseMatrix<number> &matrix);
 
                                     /**
                                      * Call the respective function

diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc
index df2e1b77619f699fbc4ff79da784d44e046b0d6a..f8f67301ef98b0560620ea2cad64e84d920e5691 100644 (file)
--- a/deal.II/lac/source/sparse_direct.cc
+++ b/deal.II/lac/source/sparse_direct.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -702,8 +702,9 @@ SparseDirectMA27::initialize (const SparsityPattern &sp)
 
 
 
+template <typename number>
 void
-SparseDirectMA27::factorize (const SparseMatrix<double> &matrix)
+SparseDirectMA27::factorize (const SparseMatrix<number> &matrix)
 {
                                   // if necessary, initialize process
   if (initialize_called == false)
@@ -850,6 +851,7 @@ SparseDirectMA27::factorize (const SparseMatrix<double> &matrix)
 
 
 
+template <>
 void
 SparseDirectMA27::solve (Vector<double> &rhs_and_solution) const
 {
@@ -863,8 +865,31 @@ SparseDirectMA27::solve (Vector<double> &rhs_and_solution) const
 
 
 
+template <>
 void
-SparseDirectMA27::solve (const SparseMatrix<double> &matrix,
+SparseDirectMA27::solve (Vector<float> &rhs_and_solution) const
+{
+  Assert (factorize_called == true, ExcFactorizeNotCalled());
+
+                                   // first have to convert data type to
+                                   // doubles
+  Vector<double> tmp (rhs_and_solution.size());
+  tmp = rhs_and_solution;
+  
+  const unsigned int n_rows = rhs_and_solution.size();
+  call_ma27cd (&n_rows, &A[0], &LA,
+               &IW[0], &LIW, &MAXFRT,
+               &tmp(0), &IW1[0], &NSTEPS);
+
+                                   // then copy result back
+  rhs_and_solution = tmp;
+}
+
+
+
+template <typename number>
+void
+SparseDirectMA27::solve (const SparseMatrix<number> &matrix,
                         Vector<double>             &rhs_and_solution)
 {
   initialize (matrix.get_sparsity_pattern());
@@ -899,8 +924,9 @@ SparseDirectMA27::get_synchronisation_lock () const
 
 
 
+template <typename number>
 void
-SparseDirectMA27::fill_A (const SparseMatrix<double> &matrix)
+SparseDirectMA27::fill_A (const SparseMatrix<number> &matrix)
 {
   Assert (n_nonzero_elements <= A.size(), ExcInternalError());
 
@@ -1499,3 +1525,25 @@ call_ma47cd (const unsigned int *n_rows,           //scalar
                     IW, LIW, &W[0],
                     rhs_and_solution, IW1, ICNTL);  
 }
+
+
+
+// explicit instantiations
+template
+void
+SparseDirectMA27::factorize (const SparseMatrix<double> &matrix);
+
+template
+void
+SparseDirectMA27::factorize (const SparseMatrix<float> &matrix);
+
+template
+void
+SparseDirectMA27::solve (const SparseMatrix<double> &matrix,
+                        Vector<double>             &rhs_and_solution);
+
+template
+void
+SparseDirectMA27::solve (const SparseMatrix<float>  &matrix,
+                        Vector<double>             &rhs_and_solution);
+