MKL gesv memory usage

[Gianluca_G_]

Testing LAPACK_zgesv, 

I saw that it allocates a large amount of memory, in particular this can crash my rogram when for example a matrix 5000x5000 is passed.
Is there any alternative solution or option in order to limitate this overuse of memory?

I'm using MKL 11.3.01 in C# if it can be helpfull.

Thanks

Gianluca

[Gennady_F_Intel]

yes, in your case > 40 GB of RAM would require for such cases. May be using only 1RHS may help? but the effect will be very negligible... 

[Gianluca_G_]

What do you mean 1RHS?

Thanks

Gianluca

[Gianluca_G_]

I think I'm already using (1 rhs) matrix.

[Gianluca_G_]

is it possible to avoid the iterative post refinement in order to reduce memory usage? 

[mecej4]

A double precision complex matrix of size 5000 X 5000 occupies just 400 MB, if allocated as a straight two-dimensional array. A single RHS vector of size 5000 occupies 80 kB. If you call any of the siblings of GESV other than DSGESV or ZCGESV, there is no iterative refinement done.

There is no routine with the name LAPACK_zgesv. There is ZGESV, and there is Lapacke_zgesv. Which one did you call? How did you allocate your matrices?

Thus, there are a number of discrepancies in the posts of this thread, and they need to be resolved.

[Gianluca_G_]

sorry, I use LAPACKE_zgesv().

below the example, it works but it use a lot of memory. 

[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
        internal static extern int LAPACKE_zgesv(
            int matrix_layout,
            int n,
            int nrhs,
            [In, Out] Complex[] input_matrix,
            int lda,
            [In, Out] int[] ipiv,
            [In, Out] Complex[] b,
            int ldb
        );

        public static Complex[,] MKLResolve(Complex[,] MSL, Complex[,] VTN)
        {
            Complex[] A = MKLMatrix(MSL);
            Complex[] b = MKLVector(VTN);

            int n = MSL.GetLength(0);
            int nrhs = 1;
            int lda = n;
            int ldb = 1;
            int info = 1;

            int[] ipiv = new int;

            info = MKLWrapper.LAPACKE_zgesv(LAPACK_ROW_MAJOR, n, nrhs, A, lda, ipiv, b, ldb);
            if (info == 0)
            {
                Complex[,] res = MatrixGeneric.Create<Complex>(n, true);
                for (int i = 0; i < b.Length; i++)
                {
                    res[i + 1, 1] = b;
                }
                return res;
            }
            else
            {
                return null;
            }
        }

Thks

[mecej4]

Line-15 of #7 makes a copy of the matrix, thereby doubling the memory footprint of the program. Such copying is fine when the matrices are small or when the program is being debugged, but is it really necessary to do this with a 5000 X 5000 matrix?

If you are saving a copy simply to be able to obtain solutions with other RHS vectors, instead of calling GESV you can call GETRF to do the factorization once and then call GETRS as many times as needed to obtain solutions.

[Gianluca_G_]

The copy is because your function accept a vector in input! the matrix is rappresented as vector.

So, I have a matrix Complex[,] and convert it in a Vector Complex[], otherwise I can't use the function, is there other possibilities, remeber that I'm working in C#.

[mecej4]

Gianluca G. wrote:

The copy is because your function accept a vector in input! the matrix is represented as vector.

Please note that Intel employees usually signal their status by having "(Intel)" as part of their name. I am a forum member with the same status as you. It is not clear which MKL function you refer to as "your function". The C-interface function, the Fortran interface function, or a C# interface function. If the last, I don't know anything about it, whether it came from Intel or somewhere else.

So, I have a matrix Complex[,] and convert it in a Vector Complex[], otherwise I can't use the function, is there other possibilities, remember that I'm working in C#.

If you look in the MKL-C reference manual, you will see that the type of the A and B matrix arguments is lapack_complex_double *. Therefore, what is passed is just the base address of the array. Whether that array is one or two-dimensional is immaterial. In fact, the LDA and LDB arguments are provided just to map double subscripts to offsets relative to the base of the array. If you now use a number of wrapper routines to conform to the stricter type rules of C# (as compared to C or Fortran), you may well box yourself into a corner where you are forced to make a copy of a huge matrix, costing you both memory and CPU time. I am only an occasional user of C#, and I do not know anything about the wrappers that you use, so I cannot provide specific advice. I have a hunch that it may be possible to typecast a reference to a two-dimensional array into a reference to a one-dimensional array, just as you can in C. Were it not so, C# would not be suitable for the kind of work that you have undertaken.

[Gianluca_G_]

sorry for the misunderstanding, I thought you were Intel employee. 

The question of representing a matrix as vector, was because when I tested the function I got the right result in that way, maybe I did something wrong, I will try without making a copy.

Thank you very much

[mecej4]

Gianluca G. wrote:

sorry for the misunderstanding, I thought you were Intel employee. 

No problem; new users often don't realize that this forum is visible to anyone, and that anyone who is interested may obtain an account and then post new messages and replies.

The question of representing a matrix as vector, was because when I tested the function I got the right result in that way, maybe I did something wrong, I will try without making a copy.

As I said above, I am only an occasional user of C#. However, since MKL is mostly Fortran 77, plus smaller parts in Fortran 95 and C, I felt that it should be possible to pass a contiguous array from C# just as it is possible in C and Fortran. A bit of experimentation shows that it works, so you should be able to get rid of the duplicate array and the work of copying back and forth. 

As a test case, I took the dgeev.cs example contained in the second zip file at https://software.intel.com/en-us/articles/using-intel-mkl-in-your-c-program , and changed the type of the array A from double [] to double [,]. Correspondingly, I changed the declared argument type of the function argument A in the prototype declarations. The modified program worked fine. I think that you could try the same way of keeping the compiler happy by telling it "white lies". In other words, Lapack routines only need the address of the base of the array, since they obtain the needed additional information about the shape and size of the array from other arguments (ldA, n).

Beware that this may be a risky and unapproved practice, and a C# expert may be able to provide a better/safer way.

[Gianluca_G_]

Hi mecej4,

I confirm that you were right, I tried again my tests with matrixs as input and this  works. So I was wrong with my previous ones tests. 

Now the memory usage is half!

just one more question. In a post before, you mentioned that there are two functions: ZGESV, and Lapacke_zgesv which one is better to use, in terms of memory?

Thank you very much

Gianluca

[mecej4]

Gianluca G. wrote:

In a post before, you mentioned that there are two functions: ZGESV, and Lapacke_zgesv which one is better to use, in terms of memory?

Lapacke_zgesv is just a wrapper around ZGESV. The actual work-horse is ZGESV (and other routines that it may call). In general, going through wrapper routines will entail a slightly higher memory consumption. In most cases, this overhead is negligible. On the other hand, there can be a significant CPU time overhead cost if a wrapper routine has to reorganize and restore a two-dimensional matrix from row-major storage to column-major storage.

My recommendation is that you use whichever interface is convenient and easy to use, and not worry about slight changes in memory or CPU time consumption.

For the history of LapackE please see https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/300188. ;For more information on LapackE and to obtain the source code, see http://netlib.org/lapack/lapacke.html .

[Gianluca_G_]

Thank you