This project has moved and is read-only. For the latest updates, please go here.

Memory doubled when using LUFactor

Jun 17, 2014 at 1:20 PM

I'm solving a large complex dense matrix (10k*10k) by using "void LUFactor(Complex[] data, int order, Int32[] ipiv)".

If I understand correctly, LUFactor should be pointed to zgetf2 or zgetrf Fortran function which will do a in-place LU-decomposition. So there should be no extra copies produced.

But the memory usage becomes doubled (from 1.6G to more than 3G) when LUFactor once be called. Can someone kindly tell me the reason? It would be great if I can figure out a solution to avoid the extra memory usage.

Thanks alot in advance
Jun 17, 2014 at 1:54 PM

It does use zgetrf.
Are you calling the LUFactor routine directly from the MKL native provider?

Jun 17, 2014 at 2:13 PM
Hi Christoph,

Here is my small test code:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics.Algorithms.LinearAlgebra.Mkl;
using System.Numerics;

namespace TestMathNet
    class Program
        static void Main(string[] args)
            int n1 = 10000;
            int n2 = n1 * n1;
            Complex[] A = new Complex[n2];
            Complex[] b = new Complex[n1];

            for (int i = 0; i < n1; i++)
                for (int j = 0; j < n1; j++)
                    int ij = j * n1 + i; //column major
                    A[ij] = new Complex(i + 1, j + 1);
                b[i] = new Complex(i + 1, i + 2);

            int[] ipiv = new int[n1];

            MklLinearAlgebraProvider MKL = new MklLinearAlgebraProvider();
            MKL.LUFactor(A, n1, ipiv);
The memory usage will be doubled when the line "MKL.LUFactor(...);" is stepped over. But if I call zgetrf directly in Fortran the memory will not increase.

Thanks for reply
Yi Luo
Jun 17, 2014 at 3:09 PM
Edited Jun 17, 2014 at 3:13 PM
Thanks for the repro.

I can confirm:
  • The process private memory is doubled while LUFactor is running, but freed right before it returns
  • The additional memory does not seem to be managed memory.
  • This only happens with Complex and Complex32, but not with double or float.
It may be caused by the marshaling of the Complex/Complex32 structure.
This may be tweakable by adjusting the structure layout and marshaling (maybe manual pinning and pointer/unsafe).

Update: would be great if we could avoid any copying at marshaling!
Jun 17, 2014 at 3:22 PM
Hi Christoph,

Thanks a lot for the reply.

I'm not so good at low-level or native programming. Could you explain a bit more about how to avoid this problem (maybe some lines of code)?

Thanks again
Yi Luo
Jun 17, 2014 at 9:36 PM
Just to clarify, this would need to be fixed within Math.NET Numerics itself and the native provider.

Some relevant links: