0

I am solving sparse system of linear equations with Math.Net Numerics in C#. I'm trying to export this code to dll (using Robert Giesecke's Unmanaged Exports) and run it from Delphi code.

Here is my C# code:

namespace MathNetSolverLib
{
    public class ParticleInterpolator
    {
        public ParticleInterpolator(int eqnumber)
        {
            n = eqnumber;
            b = Vector.Build.DenseOfArray(new double[eqnumber]);
            result = Vector.Build.DenseOfArray(new double[eqnumber]);
            input_points = new List<Tuple<int, int, double>>();
        }

        public void solve()
        {
            using (var writer = new StreamWriter("C:\\users\\s\\desktop\\logb.txt"))
            {
                for (int i = 0; i < n; ++i)
                {
                    writer.WriteLine(b[i]);
                }
            }
            using (var writer = new StreamWriter("C:\\users\\s\\desktop\\loginput.txt"))
            {
                for (int i = 0; i < input_points.Count(); ++i)
                {
                    writer.WriteLine(input_points[i]);
                }
            }

            TFQMR solver = new TFQMR();

            MathNet.Numerics.LinearAlgebra.Solvers.Iterator<double> iter = 
                new MathNet.Numerics.LinearAlgebra.Solvers.Iterator<double>();

            SparseMatrix sd = SparseMatrix.OfIndexed(n, n, input_points);
            solver.Solve(sd, b, result, iter, new DiagonalPreconditioner());



            using (var writer = new StreamWriter("C:\\users\\s\\desktop\\logres.txt"))
            {
                for (int i = 0; i < n; ++i)
                {
                    writer.WriteLine(result[i]);
                }
            }

        }
        public int n;
        public List<Tuple<int, int, double>> input_points;
        public Vector<double> result;
        public Vector<double> b;
    }

    public class Class1
    {
        [DllExport]
        public static ParticleInterpolator create(int eqnumber)
        {           
            return new ParticleInterpolator(eqnumber);
        }

        [DllExport]
        public static void addValue(int row, int col, double value, ParticleInterpolator ins)
        {
            ins.input_points.Add(new Tuple<int, int, double>(row, col, value));
        }

        [DllExport]
        public static void setRhs(int row, double value, ParticleInterpolator ins)
        {
            ins.b[row] = value;
        }

        [DllExport]
        public static void solve(ParticleInterpolator ins)
        {      
            ins.solve();
        }

        [DllExport]
        public static double getSolutionByIndex(int index, ParticleInterpolator ins)
        {
            return ins.result[index];
        }
    }
}

I have imported this to C# dll like this:

 [DllImport("MathNetSolverLib.dll", EntryPoint = "create")]
        public static extern int create(int eqnumber);

        [DllImport("MathNetSolverLib.dll", EntryPoint = "addValue")]
        public static extern void addValue(int row, int col, double value, int ins);

        [DllImport("MathNetSolverLib.dll", EntryPoint = "setRhs")]
        public static extern void setRhs(int row, double value, int ins);

        [DllImport("MathNetSolverLib.dll", EntryPoint = "solve")]
        public static extern void solve(int ins);

        [DllImport("MathNetSolverLib.dll", EntryPoint = "getSolutionByIndex")]
        public static extern double getSolutionByIndex(int index, int ins);

and it works fine.

But when I am trying to import this in Delphi:

function create(n: integer): pointer; stdcall; external MathNetSolverLib;
procedure addValue(row, col: Integer; value:double; p:pointer); stdcall; external MathNetSolverLib;
procedure setRhs(row: Integer; value: double; p:pointer); stdcall; external MathNetSolverLib;
procedure solve(p: pointer); stdcall; external MathNetSolverLib;
function getSolutionByIndex(index: Integer; p:pointer): double; stdcall; external MathNetSolverLib;

the TFQMR solver returns wrong result. It seems that algorithm inside TFQMR has different behaviour.

It is very strange because all input data (variables A, b, n and input_points) are remain the same when calling dll from C# or Delphi (I have checked it).

Does anyone have an idea how to explain this?

TKireev
  • 43
  • 7
  • It's quite possible that the difference is down to floating point control state differences. And then there's the fact that you pass a pointer to the function from Delphi, but ignore it on the C# side. Is this your real code? Personally this seems like an inelegant way to add a sparse solver to your Delphi program. I would recommend one of Tim Davis's sparse solvers, written in C. Easy enough to wrap up in a DLL or even compile to .obj and link directly. – David Heffernan Jun 09 '15 at 11:52
  • Anyway, I think that if you want detailed help you would do well to add the missing detail. We've no way to know whether or not A, b and n really are the same. You've hidden away the code that sets them. – David Heffernan Jun 09 '15 at 12:05
  • You'll need to add the code in an edit to the question. – David Heffernan Jun 09 '15 at 12:09
  • I can't see either the input or the output. I still think that floating point control status could be relevant. Your design is broken because there's no reason for the GC to keep the `ParticleInterpolator` instance alive. So far as it knows it is collectable. I'd still expect the performance to be pretty shaky with all those managed to unmanaged transitions. I strongly recommend that you don't try to implement it this way. I'd keep .net out of the picture completely. But if you had to use .net then you'd pass all the inputs in a single call. – David Heffernan Jun 09 '15 at 12:19
  • Also, a pointer is `IntPtr` in C#, not `int`. – David Heffernan Jun 09 '15 at 12:20
  • Furthermore, swapping pointers around with managed code, afaik, is also dangerous - pointers are not fixed to specific memory in .NET and your object can move around if the runtime memory manager decides to put it somewhere else -- even if it doesn't get GC'd you may end up with a bunk reference. To anchor a pointer you have to write the C# code in an `unsafe` section and use the `fixed` statement to force the variable into a fixed pointer. https://msdn.microsoft.com/en-US/library/f58wzh21%28v=vs.80%29.aspx – J... Jun 09 '15 at 12:21
  • Thank you a lot! Anyway, because the instance of ParticleInterpolator is collectable, even if this code works well, I shouldn't use .net. It is not necessary for me. – TKireev Jun 09 '15 at 12:23
  • 1
    The sparse solver that I use is one of Tim Davis's direct solvers. As used by Matlab. I commend his solvers highly. – David Heffernan Jun 09 '15 at 12:25
  • I agree with @DavidHeffernan that mixing managed code into a Delphi application should probably be avoided given that there are ample native alternatives. If you want to pursue this project, however, the key to take away from this is that you need to keep managed objects inside the managed code - you should not be passing object references back and forth. – J... Jun 09 '15 at 12:33
  • Don't you know where can I get Tim Davis's direct sparse solver? – TKireev Jun 09 '15 at 12:40
  • 1
    See : https://github.com/PetterS/SuiteSparse/tree/master/CSparse Also : https://books.google.ca/books?hl=en&lr=&id=aTLYrafw3vUC&oi=fnd&pg=PR1&dq=Tim+Davis%27s+direct+solvers&ots=QIPsoNvYSr&sig=RNYElSvLnrGSYfB3oUFu6poY6uY#v=onepage&q=Tim%20Davis%27s%20direct%20solvers&f=false – J... Jun 09 '15 at 13:09
  • @J... Thanks. There's also UMFPack. which I believe is the solver that MATLAB uses. – David Heffernan Jun 09 '15 at 13:29
  • @TKireev A web search would have found the answers. That's presumably how J... did it. Don't forget to use search engines! You may find that there are other solvers you prefer. You've picked an iterative one in the question. Perhaps for good reason. In which case you would not use one of Tim's direct solvers. – David Heffernan Jun 09 '15 at 13:30
  • i'd expect C3 to use SSE rather than x87 which perhaps would lessen dependance upon x87 state register. I wonder what wold happen if to try to make Delphi Win64 app - perhaps it can contrast the problems with pointers/datatypes mismatch that otherwise go unnoticed and cause hidden damage – Arioch 'The Jun 10 '15 at 00:13
  • @Arioch SSE has floating point control state with all the same options as x87 – David Heffernan Jun 10 '15 at 06:26

0 Answers0