Changed: Let the second argument of SystemMatrix::initAssembly()

be a char to enable forced preassembly in addition to delayed
sparsity pattern lock. Forced preassembly is activated by
num_threads_SLU equal to -2 (maybe a bit hacklish, but...).
This is needed for the single-thread or non/optimized builds
where the preassembly step normally is skipped.

Author: kmokstad

src/ASM/AlgEqSystem.C

@@ -34,7 +34,8 @@ bool AlgEqSystem::init (LinAlg::MatrixType mtype, const LinSolParams* spar,
                         bool withReactions, int num_threads_SLU)
 {
   // Using the sign of the num_threads_SLU argument to flag this (convenience)
-  bool dontLockSparsityPattern = num_threads_SLU < 0;
+  char paFlag = num_threads_SLU == -2 ? 'f' : (num_threads_SLU < 0 ? 'd' : 0);
+  if (num_threads_SLU < 0) num_threads_SLU = 1;
 
   size_t i;
   for (i = nmat; i < A.size(); i++)
@@ -55,11 +56,11 @@ bool AlgEqSystem::init (LinAlg::MatrixType mtype, const LinSolParams* spar,
       if (spar)
         A[i]._A = SystemMatrix::create(adm,mtype,*spar);
       else
-        A[i]._A = SystemMatrix::create(adm,mtype,abs(num_threads_SLU));
+        A[i]._A = SystemMatrix::create(adm,mtype,num_threads_SLU);
       if (!A[i]._A) return false;
     }
 
-    A[i]._A->initAssembly(sam,dontLockSparsityPattern);
+    A[i]._A->initAssembly(sam,paFlag);
     A[i]._b = nullptr;
   }
 

src/LinAlg/DenseMatrix.C

@@ -74,7 +74,7 @@ DenseMatrix::DenseMatrix (const Matrix& A, bool s) : myMat(A)
 }
 
 
-void DenseMatrix::initAssembly (const SAM& sam, bool)
+void DenseMatrix::initAssembly (const SAM& sam, char)
 {
   myMat.resize(sam.neq,sam.neq,true);
 }

src/LinAlg/DenseMatrix.h

@@ -71,7 +71,7 @@ class DenseMatrix : public SystemMatrix
 
   //! \brief Initializes the element assembly process.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  virtual void initAssembly(const SAM& sam, bool);
+  virtual void initAssembly(const SAM& sam, char);
 
   //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
   virtual void init();

src/LinAlg/DiagMatrix.C

@@ -24,7 +24,7 @@ DiagMatrix::DiagMatrix (const RealArray& data, size_t nrows)
 }
 
 
-void DiagMatrix::initAssembly (const SAM& sam, bool)
+void DiagMatrix::initAssembly (const SAM& sam, char)
 {
   myMat.resize(sam.neq,true);
 }

src/LinAlg/DiagMatrix.h

@@ -58,9 +58,8 @@ class DiagMatrix : public SystemMatrix
                     const char* label);
 
   //! \brief Initializes the element assembly process.
-  //! \details Must be called once before the element assembly loop.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  virtual void initAssembly(const SAM& sam, bool);
+  virtual void initAssembly(const SAM& sam, char);
 
   //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
   virtual void init() { myMat.fill(Real(0)); }

src/LinAlg/ISTLMatrix.C

@@ -117,18 +117,18 @@ ISTLMatrix::~ISTLMatrix ()
 }
 
 
-void ISTLMatrix::initAssembly (const SAM& sam, bool delayLocking)
+void ISTLMatrix::initAssembly (const SAM& sam, char)
 {
-  SparseMatrix::initAssembly(sam, delayLocking);
-  SparseMatrix::preAssemble(sam, delayLocking);
+  this->resize(sam.getNoEquations());
+  this->preAssemble(sam,false);
 
-  std::vector<std::set<int>> dofc;
+  std::vector<IntSet> dofc;
   sam.getDofCouplings(dofc);
 
   // Set correct number of rows and columns for matrix.
   size_t sum = 0;
-  for (const auto& it : dofc)
-    sum += it.size();
+  for (const IntSet& dofs : dofc)
+    sum += dofs.size();
 
   iA.setSize(rows(), cols(), sum);
   iA.setBuildMode(ISTL::Mat::random);
@@ -138,8 +138,8 @@ void ISTLMatrix::initAssembly (const SAM& sam, bool delayLocking)
   iA.endrowsizes();
 
   for (size_t i = 0; i < dofc.size(); ++i)
-    for (const auto& it : dofc[i])
-      iA.addindex(i, it-1);
+    for (int dof : dofc[i])
+      iA.addindex(i,dof-1);
 
   iA.endindices();
 
@@ -159,7 +159,7 @@ bool ISTLMatrix::endAssembly()
 
 void ISTLMatrix::init ()
 {
-  SparseMatrix::init();
+  this->SparseMatrix::init();
 
   // Set all matrix elements to zero
   iA = 0;

src/LinAlg/ISTLMatrix.h

@@ -102,12 +102,9 @@ class ISTLMatrix : public SparseMatrix
   SystemMatrix* copy() const override { return new ISTLMatrix(*this); }
 
   //! \brief Initializes the element assembly process.
-  //! \details Must be called once before the element assembly loop.
-  //! The PETSc data structures are initialized and the all symbolic operations
   //! that are needed before the actual assembly can start are performed.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  //! \param[in] delayLocking If \e true, do not lock the sparsity pattern yet
-  void initAssembly(const SAM& sam, bool delayLocking) override;
+  void initAssembly(const SAM& sam, char) override;
 
   //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
   void init() override;

src/LinAlg/PETScMatrix.C

@@ -169,14 +169,11 @@ PETScMatrix::~PETScMatrix ()
 }
 
 
-void PETScMatrix::initAssembly (const SAM& sam, bool delayLocking)
+void PETScMatrix::initAssembly (const SAM& sam, char)
 {
-  if (adm.dd.isPartitioned())
-    this->resize(sam.neq,sam.neq);
-  else {
-    this->SparseMatrix::initAssembly(sam, delayLocking);
-    this->SparseMatrix::preAssemble(sam, delayLocking);
-  }
+  this->resize(sam.neq,sam.neq);
+  if (!adm.dd.isPartitioned())
+    this->preAssemble(sam,false);
 
   // Get number of local equations in linear system
   PetscInt neq = adm.dd.getMaxEq() - adm.dd.getMinEq() + 1;
@@ -220,7 +217,7 @@ void PETScMatrix::initAssembly (const SAM& sam, bool delayLocking)
     isvec.resize(adm.dd.getNoBlocks());
     // index sets
     for (size_t i = 0; i < isvec.size(); ++i) {
-      std::vector<int> blockEq;
+      IntVec blockEq;
       blockEq.reserve(adm.dd.getMaxEq(i+1)-adm.dd.getMinEq(i+1)+1);
       for (int leq : adm.dd.getBlockEqs(i)) {
         int eq = adm.dd.getGlobalEq(leq);

src/LinAlg/PETScMatrix.h

@@ -100,12 +100,12 @@ class PETScMatrix : public SparseMatrix
   LinAlg::MatrixType getType() const override { return LinAlg::PETSC; }
 
   //! \brief Initializes the element assembly process.
+  //! \param[in] sam Auxiliary data describing the FE model topology, etc.
+  //!
   //! \details Must be called once before the element assembly loop.
   //! The PETSc data structures are initialized and the all symbolic operations
   //! that are needed before the actual assembly can start are performed.
-  //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  //! \param[in] delayLocking If \e true, do not lock the sparsity pattern yet
-  void initAssembly(const SAM& sam, bool delayLocking) override;
+  void initAssembly(const SAM& sam, char) override;
 
   //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
   void init() override;

src/LinAlg/SPRMatrix.C

@@ -196,7 +196,7 @@ size_t SPRMatrix::dim (int idim) const
   that are need before the actual assembly can start are performed.
 */
 
-void SPRMatrix::initAssembly (const SAM& sam, bool)
+void SPRMatrix::initAssembly (const SAM& sam, char)
 {
   memset(mpar,0,NS*sizeof(int));
   msica = msifa = mtrees = mvarnc = nullptr;

src/LinAlg/SPRMatrix.h

@@ -48,7 +48,7 @@ class SPRMatrix : public SystemMatrix
 
   //! \brief Initializes the element assembly process.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  virtual void initAssembly(const SAM& sam, bool);
+  virtual void initAssembly(const SAM& sam, char);
 
   //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
   virtual void init();

src/LinAlg/SparseMatrix.C

@@ -145,10 +145,8 @@ SparseMatrix::SparseMatrix (SparseSolver eqSolver, int nt)
   nrow = ncol = 0;
   solver = eqSolver;
   numThreads = nt;
-#ifdef HAS_UMFPACK
-  umfSymbolic = nullptr;
-#endif
   slu = nullptr;
+  umfSymbolic = nullptr;
 }
 
 
@@ -161,9 +159,7 @@ SparseMatrix::SparseMatrix (size_t m, size_t n)
   solver = NONE;
   numThreads = 0;
   slu = nullptr;
-#ifdef HAS_UMFPACK
   umfSymbolic = nullptr;
-#endif
 }
 
 
@@ -177,9 +173,7 @@ SparseMatrix::SparseMatrix (const SparseMatrix& B)
   solver = B.solver;
   numThreads = B.numThreads;
   slu = nullptr; // The SuperLU data (if any) is not copied
-#ifdef HAS_UMFPACK
-  umfSymbolic = nullptr;
-#endif
+  umfSymbolic = nullptr; // The UMFPACK data (if any) is not copied
 }
 
 
@@ -210,6 +204,7 @@ bool SparseMatrix::lockPattern (bool doLock)
 void SparseMatrix::resize (size_t r, size_t c, bool forceEditable)
 {
   factored = false;
+  if (c == 0) c = r;
   if (r == nrow && c == ncol && !forceEditable)
   {
     // Clear the matrix content but retain its sparsity pattern
@@ -227,15 +222,14 @@ void SparseMatrix::resize (size_t r, size_t c, bool forceEditable)
   A.clear();
 
   nrow = r;
-  ncol = c > 0 ? c : r;
+  ncol = c;
 
   delete slu;
   slu = nullptr;
 #ifdef HAS_UMFPACK
-  if (umfSymbolic) {
+  if (umfSymbolic)
     umfpack_di_free_symbolic(&umfSymbolic);
-    umfSymbolic = nullptr;
-  }
+  umfSymbolic = nullptr;
 #endif
 }
 
@@ -477,8 +471,8 @@ bool SparseMatrix::augment (const SystemMatrix& B, size_t r0, size_t c0)
 
   for (const ValueMap::value_type& val : Bptr->elem)
   {
-    elem[std::make_pair(r0+val.first.first,c0+val.first.second)] += val.second;
-    elem[std::make_pair(c0+val.first.second,r0+val.first.first)] += val.second;
+    elem[{ r0+val.first.first,  c0+val.first.second }] += val.second;
+    elem[{ c0+val.first.second, r0+val.first.first  }] += val.second;
   }
 
   return true;
@@ -736,12 +730,14 @@ static void assemSparse (const RealArray& V, SparseMatrix& SM, size_t col,
 }
 
 
-void SparseMatrix::initAssembly (const SAM& sam, bool delayLocking)
+void SparseMatrix::initAssembly (const SAM& sam, char preAssembly)
 {
   this->resize(sam.neq,sam.neq);
+  if (preAssembly == 'f')
+    this->preAssemble(sam,false);
 #ifdef USE_OPENMP
-  if (omp_get_max_threads() > 1)
-    this->preAssemble(sam,delayLocking);
+  else if (omp_get_max_threads() > 1)
+    this->preAssemble(sam,preAssembly=='d');
 #endif
 }
 
@@ -758,7 +754,7 @@ void SparseMatrix::preAssemble (const SAM& sam, bool delayLocking)
 
   // If we are not locking the sparsity pattern yet, the index pair map over
   // the non-zero matrix elements needs to be initialized before the assembly.
-  // This is used when SAM::getDofCouplings does not return all connectivities
+  // This is used when SAM::getDofCouplings() does not return all connectivities
   if (delayLocking) // that will exist in the final matrix.
     for (size_t i = 0; i < dofc.size(); i++)
       for (int j : dofc[i])

src/LinAlg/SparseMatrix.h

@@ -105,8 +105,12 @@ class SparseMatrix : public SystemMatrix
 
   //! \brief Initializes the element assembly process.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  //! \param[in] delayLocking If \e true, do not lock the sparsity pattern yet
-  virtual void initAssembly(const SAM& sam, bool delayLocking);
+  //! \param[in] preAssembly Flag for doing preassembly
+  //! - f : Force preassembly always, and lock the sparsity pattern
+  //! - d : Do preassembly only if multi-threaded assembly,
+  //!       but delay locking the sparsity pattern
+  //! - 0 : Do preassembly if multi-threading, and lock parsity pattern
+  virtual void initAssembly(const SAM& sam, char preAssembly);
 
   //! \brief Initializes the element sparsity pattern based on node connections.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
@@ -275,10 +279,7 @@ class SparseMatrix : public SystemMatrix
   SparseSolver solver; //!< Which equation solver to use
   SuperLUdata*    slu; //!< Matrix data for the SuperLU equation solver
   int      numThreads; //!< Number of threads to use for the SuperLU_MT solver
-
-#ifdef HAS_UMFPACK
-  void* umfSymbolic; //!< Symbolically factored matrix for UMFPACK
-#endif
+  void*   umfSymbolic; //!< Symbolically factored matrix for UMFPACK
 
 protected:
   bool factored; //!< Set to \e true when the matrix is factorized

src/LinAlg/SystemMatrix.h

@@ -230,10 +230,9 @@ class SystemMatrix
 
   //! \brief Initializes the element assembly process.
   //! \param[in] sam Auxiliary data describing the FE model topology, etc.
-  //! \param[in] delayLocking If \e true, do not lock the sparsity pattern yet
   //!
   //! \details This method must be called once before the element assembly loop.
-  virtual void initAssembly(const SAM& sam, bool delayLocking = false) = 0;
+  virtual void initAssembly(const SAM& sam, char = 0) = 0;
 
   //! \brief Initializes the matrix to zero assuming it is properly dimensioned.
   virtual void init() = 0;

src/SIM/SIMinput.C

@@ -368,6 +368,12 @@ bool SIMinput::parseGeometryTag (const tinyxml2::XMLElement* elem)
       IFEM::cout <<"\tActivation function for P"<< patch
                  <<": "<< funcdef << std::endl;
       pch->setElementActivator(utl::parseIntFunc(funcdef,type));
+#ifndef USE_OPENMP
+      // Force pre-assembly also for Debug and Nopt builds, such that
+      // the final sparsity pattern is established up front when not
+      // all elements are active in the first step/iteration
+      opt.num_threads_SLU = -2;
+#endif
     }
   }