Project cleanup: the Matrix class

[wo80]

As mentioned in https://github.com/BriefFiniteElementNet/BriefFiniteElement.Net/issues/36#issuecomment-673732028, here are some issues I see regarding the Matrix class:

1. There are approximately 70-80 public methods in the Matrix class and 30-40 of them are never used anywhere in the whole solution. That's about 50%. Since BFE.NET is not a linear algebra package (meaning the Matrix class is rather an internal implementation detail), all of this code should be removed (so it doesn't have to be maintained and tested).

2. Implementing the IEnumerable<double> interface doesn't make too much sense. It's only used in DkqHelper (line 309) for LINQ, but [matrix].Max(i=>Math.Abs(i)) could easily be written as [matrix].CoreArray.Max(i=>Math.Abs(i)).

3. There are a few methods dealing with matrix inversion (Determinant, Inverse, Inverse2, Solve). This could all be handled with LU factorization. Why are there two Invert methods? Missing documentation!

4. Code duplication:

   • line 457: Matrix Multiply(Matrix m1, Matrix m2) vs void Multiply(Matrix m1, Matrix m2, Matrix result)
   • line 1596: Matrix operator +(Matrix mat1, Matrix mat2) vs void Plus(Matrix mat1, Matrix mat2, Matrix result)

5. Avoid memory allocations or at least provide an alternative. For example line 977, Matrix ExtractRow(int i) should have a twin void ExtractRow(int i, Matrix target). Using a double array would be even better, why is it a Matrix? Missing documentation!

6. Line 612: public static double[] Multiply(Matrix m, double[] vec). Why is this a static method and not a class method?

7. There's a mysterious #region nonzero Pattern region at the end of the file. It is used in SparseMatrixMultiplyValidation, but why is it part of the Matrix class?

8. Line 684: Dictionary<int, int> dists = new Dictionary<int, int>()? Not used anywhere!

Many of the points above could be labeled messy code and are not related to the matrix/linear algebra subject (which I chose for obvious reasons). I haven't looked at any other classes, but the problem I see is that a user or potential contributor cannot guess the intent of the code, so it will be hard to clean it up (for any other person than you @epsi1on ).

Regarding the Matrix class: I'd suggest to replace it with the DenseMatrix class of CSparse. Missing functionality could be added using extension methods or by moving the code to CSparse. I'm also thinking about adding some dense factorizations to CSparse, which would address point 3 from above list.

[epsi1on]

Hi, Yes i had plan to exclude builtin Matrix and use CSparse.DenseMatrix. Also Currently an array pooling system is implemented to reuse Matrix.Values arrays. I have plan to use BriefFiniteElementNet.Common.DisposableDenseMatrix class which is inherited from CSparse.DenseMatrix and have a reference to an ArrayPool<double> to return the DenseMatrix.Values array to pool when disposed (see code here). The usage would be something like this:

using(var matrix = AollocateDense(3,3))//allocate a 3x3 matrix from pool
{
       //do stuff
}

or simply matrix.Dispose() will return array to pool. I have plan to replace the BriefFiniteElementNet.Common.Matrix with newer class of BriefFiniteElementNet.Common.DisposableDenseMatrix. Do you have any suggestion about it? I have tested and found out having a Pooling system have noticeable affect on speed, also reduce pressure on GC,

Thanks

[wo80]

You might want to consider using the System.Buffers package and not write your own array pool.

Whenever possible, I prefer allocating memory upfront and then hand it down to methods that need it. You also want to avoid allocating arrays in loops. There are a couple of places in BriefFiniteElementNet, where this could easily be applied, so I'm not sure, if a pool make sense for those cases (pooling is no free lunch, Memory Pools in .NET).

I have tested and found out having a Pooling system have noticeable affect on speed

Do have some code comparing the two cases solving the same model with and without pool?

[epsi1on]

Do have some code comparing the two cases solving the same model with and without pool?

today i did tested that again. This time i did not see any noticeable difference in using memory pool. let me test it tommorow once more with .net profiler. For a usual model to solve, There can be ~1M matrix allocations (size from 2x2 to 12x12) with current architecture.

[wo80]

For a usual model to solve, There can be ~1M matrix allocations (size from 2x2 to 12x12) with current architecture.

You're talking about assembling the element matrices into the global system matrix. If you choose to use an array pool, you will have to do the benchmarks, because it's not clear if using a pool (with appropriate locking) is faster than just allocating new memory for each element.

A different approach: In MatrixAssemblerUtil.AssembleFullStiffnessMatrix you have

var maxNodePerElement = model.Elements.Any() ? model.Elements.Select(i => i.Nodes.Length).Max() : 1;

You could do the same to get the maximum memory needed for the local matrices and hand that array over to elm.GetGlobalStifnessMatrix() and use it for the returned matrix. Then, of course, you will have to clear that memory for each element, so again you have to do the benchmarks and check if this is actually faster.

[wo80]

I did some memory profiling of Benchmark1.cs and it's a perfect example of what I mean with allocating memory upfront. In FrameElement2Node.MultSubMatrix(...) you allocate two small double arrays. Compare the following analysis results: In the second run, both arrays are passed as function parameters.

BTW: since FrameElement2Node is marked obsolete, you should update the benchmark code.

[epsi1on]

You're talking about assembling the element matrices into the global system matrix. If you choose to use an array pool, you will have to do the benchmarks, because it's not clear if using a pool (with appropriate locking) is faster than just allocating new memory for each element.

Not exactly, before assembling full stiffness matrix, stiffness matrix of each individual element should be calculated separately.

For a model with 3400 node (~20k DoF) and ~10k BarElements, there is about 1M matrix allocation I did calculate the 1M allocation with increasing a counter on Matrix class constructor. this is the code to create the model:

var model = StructureGenerator.Generate3DBarElementGrid(15, 15, 15);

model.Trace.Listeners.Add(new ConsoleTraceListener());

model.Solve_MPC();

This could be our reference for doing benchmarks.

I think we should calculate the number of allocation using a counter and size of each allocation time for above model. As allocations for small arrays more expensive than borrowing from a pool.

Also FrameElement2Node is deprecated and BarElement is the successor.

[wo80]

You're talking about assembling the element matrices into the global system matrix. If you choose to use an array pool, you will have to do the benchmarks, because it's not clear if using a pool (with appropriate locking) is faster than just allocating new memory for each element.

Not exactly, before assembling full stiffness matrix, stiffness matrix of each individual element should be calculated separately.

That's what I meant. I haven't looked at all the element implementations, so please correct me if you think this won't work, but for me a preferable way to reduce allocations is to provide memory via function parameters. Take for example the MatrixAssemblerUtil.AssembleFullStiffnessMatrix() method I already mentioned above. Currently, in your Element class you have

public abstract Matrix GetGlobalStifnessMatrix();

That should become

public abstract Matrix GetGlobalStiffnessMatrix(double[] work);
public abstract int GetMaxWorkSize(); // The caller has to figure out how much memory is needed!

The work array could be re-used for each element in the model. This would save all the allocations or renting and returning memory from a pool for each element. The same could be done for mass and damping matrix etc.

I think we should calculate the number of allocation using a counter and size of each allocation time for above model.

The above pictures were taken using the memory profiler included in Visual Studio Community edition. Do you see any problem using this to track allocations? As for the timing: a Stopwatch could be started at the beginning of StaticLinearAnalysisResult.AddAnalysisResult() and reported before the call to the linear solver.

---

Important

Changing the Matrix class will affect code all across the project. At the moment only a small part of this code is covered by unit tests, so there's a risk that something might break, without me noticing. How should we deal with this?

A possible solution: provide simple, but non-trivial integration tests for all element types.

Since there also is a request for examples, this could be combined. Let me give an example: in the CodeProject examples, instead of just printing the results

Console.WriteLine("Support reaction of n0: {0}", model.Nodes[0].GetSupportReaction());
Console.WriteLine("Support reaction of n1: {0}", model.Nodes[1].GetSupportReaction());

do

CheckResult("Support reaction of n0", model.Nodes[0].GetSupportReaction(), expectedSupportReaction0);
CheckResult("Support reaction of n1", model.Nodes[1].GetSupportReaction(), expectedSupportReaction1);

with a method CheckResult(string message, object actual, object expected) that shows an error or a success message. Same for nodal displacement, internal force etc. if applicable.

[epsi1on]

You're talking about assembling the element matrices into the global system matrix. If you choose to use an array pool, you will have to do the benchmarks, because it's not clear if using a pool (with appropriate locking) is faster than just allocating new memory for each element.

Not exactly, before assembling full stiffness matrix, stiffness matrix of each individual element should be calculated separately.

That's what I meant. I haven't looked at all the element implementations, so please correct me if you think this won't work, but for me a preferable way to reduce allocations is to provide memory via function parameters. Take for example the MatrixAssemblerUtil.AssembleFullStiffnessMatrix() method I already mentioned above. Currently, in your Element class you have

public abstract Matrix GetGlobalStifnessMatrix();

That should become

public abstract Matrix GetGlobalStiffnessMatrix(double[] work);
public abstract int GetMaxWorkSize(); // The caller has to figure out how much memory is needed!

The work array could be re-used for each element in the model. This would save all the allocations or renting and returning memory from a pool for each element. The same could be done for mass and damping matrix etc.

I think we should calculate the number of allocation using a counter and size of each allocation time for above model.

The above pictures were taken using the memory profiler included in Visual Studio Community edition. Do you see any problem using this to track allocations? As for the timing: a Stopwatch could be started at the beginning of StaticLinearAnalysisResult.AddAnalysisResult() and reported before the call to the linear solver.

Important

Changing the Matrix class will affect code all across the project. At the moment only a small part of this code is covered by unit tests, so there's a risk that something might break, without me noticing. How should we deal with this?

A possible solution: provide simple, but non-trivial integration tests for all element types.

Since there also is a request for examples, this could be combined. Let me give an example: in the CodeProject examples, instead of just printing the results

Console.WriteLine("Support reaction of n0: {0}", model.Nodes[0].GetSupportReaction());
Console.WriteLine("Support reaction of n1: {0}", model.Nodes[1].GetSupportReaction());

do

CheckResult("Support reaction of n0", model.Nodes[0].GetSupportReaction(), expectedSupportReaction0);
CheckResult("Support reaction of n1", model.Nodes[1].GetSupportReaction(), expectedSupportReaction1);

with a method CheckResult(string message, object actual, object expected) that shows an error or a success message. Same for nodal displacement, internal force etc. if applicable.

the double[] work array would not work (at least simply) because each element should pass it to its helpers (for example BarElement have three helper, then each helper will make 3 matrices (named D and B and J) to calculate K matrix at one gauss sampling point, then this samplings should be added together and returned to element from helper and then next helper and so on. so the worker array approach is not (at least) simply possible to implement.

About integral tests instead of unit tests, the whole BriefFiniteElement.Validation project is about integral tests, but they are not labeled with test attributes and project type is not test. they should run manually from a console project

[wo80]

About integral tests instead of unit tests, the whole BriefFiniteElement.Validation project is about integral tests, but they are not labeled with test attributes and project type is not test. they should run manually from a console project

It would be great to have some entry point to run all the integration tests that are passing.

[epsi1on]

I did checked again, and i've seen that external applications are used for validation. i'm afraid appveyor do not have these external applications, thus build will not pass and appveyor's red flag will rise up...

[epsi1on]

This would work as a simple small integration test which is good for checking Matrix operations:

using BriefFiniteElementNet.Elements;
using BriefFiniteElementNet.Materials;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Net.Configuration;
using System.Text;
using System.Threading.Tasks;

namespace BriefFiniteElementNet.CodeProjectExamples
{
    public class SimpleCantilever
    {
        public static void Run()
        {
            var model = new Model();

            var l = 5;

            var n1 = new Node(0, 0, 0);
            var n2 = new Node(0, 0, l);

            var axialLoad = 1000;
            var horizontalLoad = 1000;

            var f = new Force(horizontalLoad, 0, axialLoad, 0, 0, 0);

            /**/
            var h = 0.1;
            var w = 0.05;

            var a = h * w;
            var iy = h * h * h * w / 12;
            var iz = w * w * w * h / 12;
            var j = iy + iz;
            var e = 210e9;
            var nu = 0.3;

            var g = e / (2 * 1 + nu);
            /**/

            var sec = new Sections.UniformParametric1DSection(a, iy, iz, j);
            var mat = UniformIsotropicMaterial.CreateFromYoungShear(e, g);

            var belm = new BarElement(n1, n2) { Material = mat, Section = sec, Behavior = BarElementBehaviours.FullFrame };
            model.Elements.Add(belm);
            model.Nodes.Add(n1, n2);

            n1.Constraints = Constraints.Fixed;

            n2.Loads.Add(new NodalLoad(f));

            model.Solve_MPC();

            var d = model.Nodes[1].GetNodalDisplacement();

            var expectedDx = (horizontalLoad * l * l * l) / (3 * e * iy);
            var expectedRy = (horizontalLoad * l * l) / (2 * e * iy);
            var expectedDz = axialLoad * l / (e * a);

            var epsilon = 1e-6;

            if (Math.Abs(d.DX - expectedDx) > epsilon) throw new Exception("wrong value");
            if (Math.Abs(d.RY - expectedRy) > epsilon) throw new Exception("wrong value");
            if (Math.Abs(d.DZ - expectedDz) > epsilon) throw new Exception("wrong value");

        }
    }
}

if there is any problem with matrix or integration operations, then one of three lines at the end of method will throw exception with message wrong value. Is it applicable?

[epsi1on]

Are you working on this? I think I understand what you mean by sending each matrix down into each method. but it needs extensive changes. do you think it worth?

[wo80]

Are you working on this?

I started playing around with a local copy (outside of source control). I'm still not feeling comfortable with having no proper testing methods. Maybe you can write up some more for the remaining element types. As mentioned, integration tests can be seen as examples, so this would also address the request for examples in #36 (I'd suggest you add detailed comments to make clear what those tests/examples are doing). If you don't think this is too much of an issue, I will start working on this and do a pull request. Since this will affect code all across the project, the PR will be kind of blocking (at least if we want to avoid merge conflicts).

I think I understand what you mean by sending each matrix down into each method. but it needs extensive changes. do you think it worth?

I'm not sure. Lets stick with the pool for now and discuss this in a separate issue.

[epsi1on]

How is everything going on? Is there any thing I can engage?

[wo80]

How is everything going on?

The changes I did outside of source control are working, unit tests are passing and your bar element example from above completes without exception.

Is there any thing I can engage?

As I said, having more elements covered by examples/integration tests would be very welcome. The easiest way would be to go through the BriefFiniteElementNet.Validation code and wrap up all tests that are supposed to succeed into a single calling class, that could then be called in the BriefFiniteElementNet.TestConsole. Since I don't know which of the available tests are supposed to work, that's up to you.

If you don't think that more testing is needed at this point, I can start a pull request.

If there are any commits pending from your side, you should merge them now to reduce conflicts.

[epsi1on]

How is everything going on?

The changes I did outside of source control are working, unit tests are passing and your bar element example from above completes without exception.

Is there any thing I can engage?

As I said, having more elements covered by examples/integration tests would be very welcome. The easiest way would be to go through the BriefFiniteElementNet.Validation code and wrap up all tests that are supposed to succeed into a single calling class, that could then be called in the BriefFiniteElementNet.TestConsole. Since I don't know which of the available tests are supposed to work, that's up to you.

If you don't think that more testing is needed at this point, I can start a pull request.

If there are any commits pending from your side, you should merge them now to reduce conflicts.

Hi, Actually there are no changes on my side. if you send the pull request i can run validations on my side and tell you if there are any problems. So please make the pull request. Thanks ...