changed: preinit sparsity pattern if we are running with multiple threads

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1403 e10b68d5-8a6e-419e-a041-bce267b0401d

src/LinAlg/SAM.C

@@ -602,6 +602,30 @@ bool SAM::getElmEqns (IntVec& meen, int iel, int nedof) const
 }
 
 
+int SAM::getNoElmEqns (int iel) const
+{
+  if (iel < 1 || iel > nel)
+  {
+    std::cerr <<"SAM::getElmEqns: Element "<< iel <<" is out of range [1,"
+	      << nel <<"]"<< std::endl;
+    return false;
+  }
+
+  int ip = mpmnpc[iel-1];
+  int nenod = mpmnpc[iel] - ip;
+
+  int result=0;
+#ifndef USE_F77SAM 
+  // TODO F77SAM?
+  for (int i = 0; i < nenod; i++, ip++) {
+    int node = mmnpc[ip-1];
+    result += madof[node]-madof[node-1];
+  }
+#endif
+  return result;
+}
+
+
 bool SAM::getNodeEqns (IntVec& mnen, int inod) const
 {
   if (inod < 1 || inod > nnod)

src/LinAlg/SAM.h

@@ -168,6 +168,11 @@ public:
   //! (used for internal consistency checking, unless zero)
   virtual bool getElmEqns(IntVec& meen, int iel, int nedof = 0) const;
 
+  //! \brief Finds the number equations for an element.
+  //! \param[in] iel Identifier for the element to get the equation numbers for
+  //! \return The number of equations for this element
+  virtual int getNoElmEqns(int iel) const;
+
   //! \brief Finds the matrix of equation numbers for a node.
   //! \param[out] mnen Matrix of node equation numbers
   //! \param[in] inod Identifier for the node to get the equation numbers for

src/LinAlg/SparseMatrix.C

@@ -23,6 +23,10 @@
 #endif
 #include <algorithm>
 
+#ifdef USE_OPENMP
+#include <omp.h>
+#endif
+
 #if defined(HAS_SUPERLU_MT)
 #define sluop_t superlumt_options_t
 #elif defined(HAS_SUPERLU)
@@ -557,6 +561,22 @@ static void assemSparse (const RealArray& V, SparseMatrix& SM, size_t col,
 void SparseMatrix::initAssembly (const SAM& sam)
 {
   this->resize(sam.neq,sam.neq);
+#ifdef USE_OPENMP
+  // dummy assembly loop to avoid matrix resizes during assembly
+  if (omp_get_max_threads() > 1) {
+    for (int i=0;i<sam.getNoElms();++i) {
+      int siz = sam.getNoElmEqns(i+1);
+      Matrix Ek;
+      Ek.resize(siz,siz);
+      assemble(Ek,sam,i+1);
+    }
+  }
+  switch (solver) {
+    case SUPERLU: optimiseSLU(); break;
+    case S_A_M_G: optimiseSAMG(); break;
+    default: break;
+  }
+#endif
 }