* resize resets the matrix

* added an Interface to create a SuperMatrix
  (only available if HAVE_SUPERLU is defined... autoconf test will follow)

[[Imported from SVN: r764]]

fem/feop/spmatrix.cc

@@ -13,6 +13,7 @@ namespace Dune
   /*****************************/
   /*  Constructor(s)           */
   /*****************************/
+
   template <class T>
   SparseRowMatrix<T>::SparseRowMatrix()
   {
@@ -27,10 +28,6 @@ namespace Dune
   SparseRowMatrix<T>::~SparseRowMatrix()
   {}
 
-  /***********************************/
-  /*  Construct from storage vectors */
-  /***********************************/
-
   template <class T>
   SparseRowMatrix<T>::SparseRowMatrix(int rows, int cols, int nz)
   {
@@ -58,7 +55,7 @@ namespace Dune
   /* Iterator creation          */
   /******************************/
 
-  // Return iterator refering to first nonzero element in row
+  //! Return iterator refering to first nonzero element in row
   template <class T>
   typename SparseRowMatrix<T>::ColumnIterator SparseRowMatrix<T>::rbegin (int row)
   {
@@ -74,7 +71,7 @@ namespace Dune
     return it;
   }
 
-  // Return iterator refering to one past the last nonzero element of row
+  //! Return iterator refering to one past the last nonzero element of row
   template <class T>
   typename SparseRowMatrix<T>::ColumnIterator SparseRowMatrix<T>::rend (int row)
   {
@@ -86,9 +83,7 @@ namespace Dune
     return it;
   }
 
-  /***************************/
-  /* Assignment operator...  (deep copy) */
-  /***************************/
+  //! Assignment operator...  (deep copy)
 
   template <class T>
   SparseRowMatrix<T>& SparseRowMatrix<T>::operator=(const SparseRowMatrix<T>& other)
@@ -106,9 +101,7 @@ namespace Dune
   }
 
 
-  /***************************/
-  /* resize()                */
-  /***************************/
+  //! resize the Matrix and reset the content
 
   template <class T>
   void SparseRowMatrix<T>::resize(int M, int N, int nz)
@@ -121,8 +114,13 @@ namespace Dune
     // resize data fields
     values_.resize(dim_[0]*nz_);
     col_.resize(dim_[0]*nz_);
+
+    values_.set(0.0);
+    col_.set(-1);
   }
 
+  //! resize the Matrix and reset the content
+
   template <class T>
   void SparseRowMatrix<T>::resize(int M, int N)
   {
@@ -133,6 +131,9 @@ namespace Dune
     // resize data fields
     values_.resize(dim_[0]*nz_);
     col_.resize(dim_[0]*nz_);
+
+    values_.set(0.0);
+    col_.set(-1);
   }
 
   template <class T>
@@ -240,6 +241,7 @@ namespace Dune
   /***************************************/
   /*  Matrix-MV_Vector multiplication    */
   /***************************************/
+
   template <class T> template <class VECtype>
   void SparseRowMatrix<T>::mult(VECtype &ret, const VECtype &x) const
   {
@@ -486,4 +488,41 @@ namespace Dune
         set(col,col,1.0);
   }
 
+  template <class T>
+  void SparseRowMatrix<T>::createSuperMatrix(SuperMatrix & A)
+  {
+    // non zero values before this line
+    nzval_.resize(dim_[0]+1);
+    for (int i=0; i<nzval_.size(); i++)
+      nzval_[i] = 3*i;
+    // fill missing entries
+    for (int row = 0; row < dim_[0]; row++)
+    {
+      int whichCol;
+      while((whichCol = colIndex(row,-2)) != -1)
+      {
+        int col;
+        for (col = 0; col < dim_[1]; col++)
+        {
+          bool newcol = true;
+          for (int i=0; i<nz_; i++)
+            if (col_[row*nz_ +i] == col)
+              newcol = false;
+          if (newcol) break;
+        }
+        col_[row*nz_ + whichCol] = col;
+      }
+    }
+    // create matrix
+    dCreate_CompCol_Matrix(&A, dim_[0], dim_[1], values_.size(),
+                           values_.raw(), col_.raw(), nzval_.raw(),
+                           SLU_NR, SLU_D, SLU_GE);
+  }
+
+  template <class T>
+  void SparseRowMatrix<T>::destroySuperMatrix(SuperMatrix & A)
+  {
+    nzval_.resize(0);
+  }
+
 } // end namespace Dune

fem/feop/spmatrix.hh

@@ -3,6 +3,10 @@
 #ifndef __DUNE_SPMATRIX_HH
 #define __DUNE_SPMATRIX_HH
 
+#ifdef HAVE_SUPERLU
+#include <dsp_defs.h>
+#endif
+
 #include <dune/common/simplevector.hh>
 
 namespace Dune
@@ -215,6 +219,15 @@ namespace Dune
 
     //! Makes a given column a unit column
     void unitCol(int col);
+
+#ifdef HAVE_SUPERLU
+  private:
+    Array<int> nzval_;
+  public:
+    void createSuperMatrix(SuperMatrix & A);
+    void destroySuperMatrix(SuperMatrix & A);
+#endif
+
   private:
 
     //! Always contains zero.  It's only here so the index operator