Adding necessary changes for staggered scheme

src/MGroup.Solvers/AlgebraicModel/GlobalAlgebraicModel.cs

@@ -294,7 +294,7 @@ internal GlobalMatrix<TMatrix> CheckCompatibleMatrix(IGlobalMatrix matrix)
  // Casting inside here is usually safe since all global matrices should be created by this object
  if (matrix is GlobalMatrix<TMatrix> globalMatrix)
  {
- if (matrix.CheckForCompatibility == false || globalMatrix.Format == this.Format)
+ if (globalMatrix.CheckForCompatibility == false ||  globalMatrix.Format == this.Format)
  {
  return globalMatrix;
  }
@@ -310,7 +310,7 @@ internal GlobalVector CheckCompatibleVector(IGlobalVector vector)
  // Casting inside here is usually safe since all global vectors should be created by the this object
  if (vector is GlobalVector globalVector)
  {
- if (vector.CheckForCompatibility == false || globalVector.Format == this.Format)
+ if (globalVector.CheckForCompatibility == false || globalVector.Format == this.Format)
  {
  return globalVector;
  }

src/MGroup.Solvers/Direct/DenseMatrixNonSymmetricSolver.cs

@@ -0,0 +1,130 @@
+using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using MGroup.LinearAlgebra.Matrices;
+using MGroup.LinearAlgebra.Vectors;
+using MGroup.MSolve.Discretization;
+using MGroup.MSolve.Discretization.Entities;
+using MGroup.MSolve.Solution.LinearSystem;
+using MGroup.Solvers.AlgebraicModel;
+using MGroup.Solvers.Assemblers;
+using MGroup.Solvers.DofOrdering;
+using MGroup.Solvers.DofOrdering.Reordering;
+using MGroup.Solvers.LinearSystem;
+
+namespace MGroup.Solvers.Direct
+{
+ /// <summary>
+ /// Direct solver for models with only 1 subdomain. Uses Cholesky factorization on symmetric positive definite matrices
+ /// stored in full format. Its purpose is mainly for testing, since it is inefficient for large linear systems resulting 
+ /// from FEM . Cholesky or LU for Generic matrices
+ /// Authors: Serafeim Bakalakos
+ /// </summary>
+ public class DenseMatrixNonSymmetricSolver : SingleSubdomainSolverBase<Matrix>
+ {
+ private readonly bool isMatrixPositiveDefinite; //TODO: actually there should be 3 states: posDef, symmIndef, unsymm
+
+ private bool factorizeInPlace = true;
+ private bool mustInvert = true;
+ private Matrix inverse;
+
+ private DenseMatrixNonSymmetricSolver(GlobalAlgebraicModel<Matrix> model, bool isMatrixPositiveDefinite)
+ : base(model, "DenseMatrixNonSymmetricSolver")
+ {
+ this.isMatrixPositiveDefinite = isMatrixPositiveDefinite;
+ }
+
+ public override void HandleMatrixWillBeSet()
+ {
+ mustInvert = true;
+ inverse = null;
+ }
+
+ public override void Initialize() { }
+
+ public override void PreventFromOverwrittingSystemMatrices() => factorizeInPlace = false;
+
+ /// <summary>
+ /// Solves the linear system with back-forward substitution. If the matrix has been modified, it will be refactorized.
+ /// </summary>
+ public override void Solve()
+ {
+ var watch = new Stopwatch();
+ if (LinearSystem.Solution.SingleVector == null)
+ {
+ LinearSystem.Solution.SingleVector = Vector.CreateZero(LinearSystem.Matrix.SingleMatrix.NumRows);
+ }
+ else LinearSystem.Solution.Clear(); // no need to waste computational time on this in a direct solver
+
+ // Factorization
+ if (mustInvert)
+ {
+ watch.Start();
+ var matrix = LinearSystem.Matrix.SingleMatrix;
+ if (isMatrixPositiveDefinite)
+ {
+ inverse = matrix.FactorCholesky(factorizeInPlace).Invert(true);
+ }
+ else
+ {
+ inverse = matrix.FactorLU(factorizeInPlace).Invert(true);
+ }
+ watch.Stop();
+ Logger.LogTaskDuration("Matrix factorization", watch.ElapsedMilliseconds);
+ watch.Reset();
+ mustInvert = false;
+ }
+
+ // Substitutions
+ watch.Start();
+ inverse.MultiplyIntoResult(LinearSystem.RhsVector.SingleVector, LinearSystem.Solution.SingleVector);
+ watch.Stop();
+ Logger.LogTaskDuration("Back/forward substitutions", watch.ElapsedMilliseconds);
+ Logger.IncrementAnalysisStep();
+ }
+
+ protected override Matrix InverseSystemMatrixTimesOtherMatrix(IMatrixView otherMatrix)
+ {
+ var watch = new Stopwatch();
+ if (mustInvert)
+ {
+ watch.Start();
+ var matrix = LinearSystem.Matrix.SingleMatrix;
+ if (isMatrixPositiveDefinite)
+ {
+ inverse = matrix.FactorCholesky(factorizeInPlace).Invert(true);
+ }
+ else
+ {
+ inverse = matrix.FactorLU(factorizeInPlace).Invert(true);
+ }
+ watch.Stop();
+ Logger.LogTaskDuration("Matrix factorization", watch.ElapsedMilliseconds);
+ watch.Reset();
+ mustInvert = false;
+ }
+ watch.Start();
+ Matrix result = inverse.MultiplyRight(otherMatrix);
+ watch.Stop();
+ Logger.LogTaskDuration("Back/forward substitutions", watch.ElapsedMilliseconds);
+ Logger.IncrementAnalysisStep();
+ return result;
+ }
+
+ public class Factory
+ {
+ public Factory() { }
+
+ public IDofOrderer DofOrderer { get; set; }
+ = new DofOrderer(new NodeMajorDofOrderingStrategy(), new NullReordering());
+
+ public bool IsMatrixPositiveDefinite { get; set; } = false;
+
+ public DenseMatrixNonSymmetricSolver BuildSolver(GlobalAlgebraicModel<Matrix> model)
+ => new DenseMatrixNonSymmetricSolver(model, IsMatrixPositiveDefinite);
+
+ public GlobalAlgebraicModel<Matrix> BuildAlgebraicModel(IModel model)
+ => new GlobalAlgebraicModel<Matrix>(model, DofOrderer, new DenseMatrixAssembler());
+ }
+ }
+}

src/MGroup.Solvers/Direct/LUSolver.cs

@@ -0,0 +1,179 @@
+using System;
+using System.Diagnostics;
+using MGroup.LinearAlgebra.Matrices;
+using MGroup.LinearAlgebra.Triangulation;
+using MGroup.LinearAlgebra.Vectors;
+using MGroup.MSolve.Discretization;
+using MGroup.MSolve.Discretization.Entities;
+using MGroup.MSolve.Solution.LinearSystem;
+using MGroup.Solvers.AlgebraicModel;
+using MGroup.Solvers.Assemblers;
+using MGroup.Solvers.DofOrdering;
+using MGroup.Solvers.DofOrdering.Reordering;
+using MGroup.Solvers.LinearSystem;
+
+namespace MGroup.Solvers.Direct
+{
+ /// <summary>
+ /// Direct LU solver
+ /// Authors: Papas Orestis, Christodoulou Theofilos
+ /// </summary>
+ public class LUSolver : SingleSubdomainSolverBase<CscMatrix>, IDisposable
+ {
+ private const bool useSuperNodalFactorization = true; // For faster back/forward substitutions.
+ private readonly double factorizationPivotTolerance;
+
+ private bool mustFactorize = true;
+ private LUCSparseNet factorization;
+
+ private LUSolver(GlobalAlgebraicModel<CscMatrix> model, double factorizationPivotTolerance) 
+ : base(model, "LUSolver")
+ {
+ this.factorizationPivotTolerance = factorizationPivotTolerance;
+ }
+
+ ~LUSolver()
+ {
+ ReleaseResources();
+ }
+
+ public void Dispose()
+ {
+ ReleaseResources();
+ GC.SuppressFinalize(this);
+ }
+
+ public override void HandleMatrixWillBeSet()
+ {
+ mustFactorize = true;
+ if (factorization != null)
+ {
+ //factorization.Dispose(); //TODO: Maybe should open this again
+ factorization = null;
+ }
+ //TODO: make sure the native memory allocated has been cleared. We need all the available memory we can get.
+ }
+
+ public override void Initialize() { }
+
+ public override void PreventFromOverwrittingSystemMatrices()
+ {
+ // The factorization is done over different memory.
+ }
+
+ /// <summary>
+ /// Solves the linear system with back-forward substitution. If the matrix has been modified, it will be refactorized.
+ /// </summary>
+ public override void Solve()
+ {
+ var watch = new Stopwatch();
+ CscMatrix matrix = LinearSystem.Matrix.SingleMatrix;
+ int systemSize = matrix.NumRows;
+ if (LinearSystem.Solution.SingleVector == null)
+ {
+ LinearSystem.Solution.SingleVector = Vector.CreateZero(systemSize);
+ }
+ else LinearSystem.Solution.Clear();// no need to waste computational time on this in a direct solver
+
+ // Factorization
+ if (mustFactorize)
+ {
+ watch.Start();
+ factorization = LUCSparseNet.Factorize(matrix);
+ watch.Stop();
+ Logger.LogTaskDuration("Matrix factorization", watch.ElapsedMilliseconds);
+ watch.Reset();
+ mustFactorize = false;
+ }
+
+ // Substitutions
+ watch.Start();
+ factorization.SolveLinearSystem(LinearSystem.RhsVector.SingleVector, LinearSystem.Solution.SingleVector);
+ watch.Stop();
+ Logger.LogTaskDuration("Back/forward substitutions", watch.ElapsedMilliseconds);
+ Logger.IncrementAnalysisStep();
+ }
+
+ protected override Matrix InverseSystemMatrixTimesOtherMatrix(IMatrixView otherMatrix)
+ {
+ throw new NotImplementedException();
+ //var watch = new Stopwatch();
+
+ //// Factorization
+ //CscMatrix matrix = LinearSystem.Matrix.SingleMatrix;
+ //int systemSize = matrix.NumRows;
+ //if (mustFactorize)
+ //{
+ // watch.Start();
+ // factorization = LUCSparseNet.Factorize(matrix);
+ // watch.Stop();
+ // Logger.LogTaskDuration("Matrix factorization", watch.ElapsedMilliseconds);
+ // watch.Reset();
+ // mustFactorize = false;
+ //}
+
+ //// Substitutions
+ //watch.Start();
+ //Matrix solutionVectors;
+ //if (otherMatrix is Matrix otherDense)
+ //{
+ // return factorization.SolveLinearSystems(otherDense.);
+ //}
+ //else
+ //{
+ // try
+ // {
+ // // If there is enough memory, copy the RHS matrix to a dense one, to speed up computations. 
+ // //TODO: must be benchmarked, if it is actually more efficient than solving column by column.
+ // Matrix rhsVectors = otherMatrix.CopyToFullMatrix();
+ // solutionVectors = factorization.SolveLinearSystems(rhsVectors);
+ // }
+ // catch (InsufficientMemoryException) //TODO: what about OutOfMemoryException?
+ // {
+ // // Solution vectors
+ // int numRhs = otherMatrix.NumColumns;
+ // solutionVectors = Matrix.CreateZero(systemSize, numRhs);
+ // var solutionVector = Vector.CreateZero(systemSize);
+
+ // // Solve each linear system separately, to avoid copying the RHS matrix to a dense one.
+ // for (int j = 0; j < numRhs; ++j)
+ // {
+ // if (j != 0) solutionVector.Clear();
+ // Vector rhsVector = otherMatrix.GetColumn(j);
+ // factorization.SolveLinearSystem(rhsVector, solutionVector);
+ // solutionVectors.SetSubcolumn(j, solutionVector);
+ // }
+ // }
+ //}
+ //watch.Stop();
+ //Logger.LogTaskDuration("Back/forward substitutions", watch.ElapsedMilliseconds);
+ //Logger.IncrementAnalysisStep();
+ //return solutionVectors;
+ }
+
+ private void ReleaseResources()
+ {
+ if (factorization != null)
+ {
+ //factorization.Dispose(); //TODO: Maybe we should open this again
+ factorization = null;
+ }
+ }
+
+ public class Factory
+ {
+ public Factory() { }
+
+ public IDofOrderer DofOrderer { get; set; }
+ = new DofOrderer(new NodeMajorDofOrderingStrategy(), AmdReordering.CreateWithSuiteSparseAmd());
+
+ public double FactorizationPivotTolerance { get; set; } = 1E-15;
+
+ public LUSolver BuildSolver(GlobalAlgebraicModel<CscMatrix> model)
+ => new LUSolver(model, FactorizationPivotTolerance);
+
+ public GlobalAlgebraicModel<CscMatrix> BuildAlgebraicModel(IModel model)
+ => new GlobalAlgebraicModel<CscMatrix>(model, DofOrderer, new CscMatrixAssembler(false));
+ }
+ }
+}

src/MGroup.Solvers/LinearSystem/GlobalMatrix.cs

@@ -31,8 +31,7 @@ public GlobalMatrix(Guid format, Func<IGlobalVector, GlobalVector> checkCompatib
  internal Guid Format { get; }
 
  public TMatrix SingleMatrix { get; set; }
-
- public bool CheckForCompatibility { get; set; }
+ public bool CheckForCompatibility { get; set; } = true;
 
  public void AxpyIntoThis(IGlobalMatrix otherMatrix, double otherCoefficient)
  {

src/MGroup.Solvers/LinearSystem/GlobalVector.cs

@@ -20,6 +20,8 @@ public GlobalVector(Guid format, Func<IGlobalVector, GlobalVector> checkCompatib
  this.checkCompatibleVector = checkCompatibleVector;
  }
 
+ public bool CheckForCompatibility { get; set; } = true;
+
  internal Guid Format { get; }
 
  public int Length => SingleVector.Length;
@@ -28,8 +30,6 @@ public GlobalVector(Guid format, Func<IGlobalVector, GlobalVector> checkCompatib
  //TODO: I would rather this was internal, but it is needed by the test classes
  public Vector SingleVector { get; set; }
 
- public bool CheckForCompatibility { get; set; }
-
  public void AxpyIntoThis(IGlobalVector otherVector, double otherCoefficient)
  {
  GlobalVector globalOtherVector = checkCompatibleVector(otherVector);

src/MGroup.Solvers/MGroup.Solvers.csproj

@@ -18,4 +18,15 @@
  <PackageReference Include="MGroup.MSolve.Core" Version="0.2.0-unstable.48" />
  <PackageReference Include="Microsoft.SourceLink.GitHub" Version="1.1.1" PrivateAssets="All" />
  </ItemGroup>
+
+ <ItemGroup>
+ <PackageReference Update="Microsoft.CodeAnalysis.FxCopAnalyzers" Version="3.3.2">
+ <PrivateAssets>all</PrivateAssets>
+ <IncludeAssets>runtime; build; native; contentfiles; analyzers</IncludeAssets>
+ </PackageReference>
+ <PackageReference Update="StyleCop.Analyzers" Version="1.2.0-beta.435">
+ <PrivateAssets>all</PrivateAssets>
+ <IncludeAssets>runtime; build; native; contentfiles; analyzers</IncludeAssets>
+ </PackageReference>
+ </ItemGroup>
 </Project>