What is the best structure to store matrices for Finite Differences Methods?

[makubex85]

Hi,
I implemented the so-called Finite Differences Methods which consist on computing a value of a matrix recursevely. The issue I am having is a stack overflow since I use 4 dimension table of table with something like 50x200x200x200 values so it's over the capacity of the machine. Once my matrices created, there is a recursive procedure where each node in the time coordinate nT is related to the one before nT-1 but where the spatial coordinates as decremented by a value that relates any coordinate smaller than the one in nT. (see the code below). My issue is how optimize the code in order not to surpass the machine capacity in terms of memory. The code roughly looks like this loops over the nodes of the dimensions and a recursion which links a point to its predecessors.
v_new = creatArray(NT+1, NGMIBB+1, NAV+1 ,NRP+1, NRUBB+1, NTaxFreeBase+1);
v_old = creatArray(NT+1, NGMIBB+1, NAV+1, NRP+1, NRUBB+1, NTaxFreeBase+1);
w_opt = creatArray(NT+1, NGMIBB+1, NAV+1, NRP+1, NRUBB+1, NTaxFreeBase+1);
tau = creatArray(NT+1, NGMIBB+1, NAV+1, NRP+1, NRUBB+1, NTaxFreeBase+1);




for(nT = NT-1; nT >= 0; nT--)
{
for(nGMIBB = 0; nGMIBB<=NGMIBB; nGMIBB++)
{
for(nAV = nGMIBB+1; nAV<=NAV; nAV++)
{
for(nRP=0; nRP<=NRP; nRP++)
{
for(nRUBB=0; nRUBB<=NRUBB; nRUBB++)
{
for(nTaxFreeBase=0; nTaxFreeBase<=NTaxFreeBase; nTaxFreeBase++)
{
v_old[nT][nGMIBB][nAV][nRP][nRUBB][nTaxFreeBase] = v_old[nT][nAV-10][nAV-80][-30][nRUBB-10][nTaxFreeBase-10];
}
}
}
}
}
}

[helios]

Assuming you're using doubles, 50*200^3*8 = 2.98 GiB. It's not an insane amount of memory. Unless you have specific reasons that require you to use a 32-bit process or you have other memory constraints, I see little point in applying special techniques to optimize memory usage. As long as the process is 64-bit, any modern system should let you allocate an array of that size regardless of how much physical memory the system has. The application will simply not perform as well because the OS will be forced to swap pages out to disk.

[ne555]

> Finite Differences Methods
you have an sparse matrix, make use of that.


> each node in the time coordinate nT is related to the one before nT-1
so you need to store only two time instants, you may dump the rest.

[makubex85]

I actually encountered a stack overflow exception when allocating my first variable (i.e creatArray(NT+1, NGMIBB+1, NAV+1 ,NRP+1, NRUBB+1, NTaxFreeBase+1) with NT = 30, NGMIBB = NAV = NRP = NRUBB = 200) and might need up to 50*300^5 values but not all are needed for the calculations that is why.

Yes I need to store two instants for my v_old but the thing is, my recursion is a cost optimization problem. It's more like finding X such that

v_old[nT][nGMIBB][nAV][nRP][nRUBB][nTaxFreeBase] = min_{X \in [0,min(nAV,nGMIBB,nRP,nRUBB]}v_old[nT][nAV-X][nGMIBB-X][nRP-X][nRUBB-X][nTaxFreeBase-X];

and I need to stock the values of X for all instants and space state variables.

I will check the sparse matrix. Thank you.

[helios]

Why would you get a stack overflow? Obviously matrices of this size are allocated on the heap, since most stacks are not larger than 1 MiB.

[makubex85]

My CreateArray function is defined as the following:

double ******* creatArray( int l, int m, int n, int o, int p, int q, int r){
	double *******A = new double ******[l];
	for(int i=0; i<l; i++){
		A[i] = new double *****[m];
		for(int j=0; j<m; j++){
			A[i][j] = new double ****[n];
			for(int k=0; k<n; k++){
				A[i][j][k] = new double ***[o];
				for(int ii = 0; ii<o; ii++)
				{
					A[i][j][k][ii] = new double **[p];
					for(int jj = 0; jj<p; jj++)
					{
						A[i][j][k][ii][jj] = new double *[q];
						for(int kk=0; kk<q; kk++){
							A[i][j][k][ii][jj][kk] = new double [r];
						}
					}
				}
			}
		}
	}
	return A;
}

The stack overflow came when I chose each dimension to be large enough in visual studio. I increased the memory but kept having this problem. When I tried to see it on Xcode creating the first matrix already used 12 out of 16 GB of memory. Can it come from this?

[helios]

That's not a stack overflow, it's just an out of memory exception.

Allocating the matrix as an array of pointers to arrays of pointers to arrays of pointers ... to doubles increases the memory requirements because a sizeable portion (dependent on the relative sizes of the dimensions) is used up by pointers.
Try allocating the matrix as a single array of size l * m * n * o * p * q * r. This won't help for the 50*300^5 case, though.

[makubex85]

Thanx. Ok. I will give on monday the number of the exception visual studio gives me.
The truth is I was trying to avoid this solution since I have complicated indices change formulas on my recursion. Since I am used to math softwares like matlab I was hoping to use a more intuitive structure but I will see if it's that complicated :)

[48T7M4Gy]

There are a couple of ways of being more intuitive with C++ on this problem - eg classes/vectors and purpose designed exception handling - which all provide more comfort than basic arrays. There is nothing wrong with arrays though and the linear approach definitely makes the array better tuned to what happens in memory but has a downside in that it places more complexity with the programmer and hence far less intuitiveness.

Unless you use specialised matrix, linear (or otherwise) C++ libraries, and there are plenty around, you won't get anything like the intuitive approach in mathlab, cad, sage or maxima.

[JLBorges]

ne555 > you have a sparse matrix, make use of that.

+1

Start with something simple (along these lines, perhaps):

#include <iostream>
#include <map>
#include <array>

template< std::size_t NDIMENSIONS > using indices = std::array< std::size_t, NDIMENSIONS > ;

template< std::size_t NDIMENSIONS, typename T > using sparse_array = std::map< indices<NDIMENSIONS>, T > ;

template< std::size_t NDIMENSIONS, typename T >
T get( const sparse_array< NDIMENSIONS, T >& arr, const indices<NDIMENSIONS>& pos )
{
    const auto iter = arr.find(pos) ;
    return iter == arr.end() ? T{} : iter->second ;
}

int main()
{
    constexpr std::size_t NDIMS = 4 ;
    sparse_array< NDIMS, double > matrix ;

    matrix[ { { 35, 123, 7, 196 } } ] = 7.89 ;
    matrix[ { { 40, 8, 78, 104 } } ] = 1.2345 ;
    matrix[ { { 0, 0, 0, 0 } } ] = -5.682 ;
    matrix[ { { 99, 999, 999, 999 } } ] += 3.34 ; // 0.0 + 3.34

    std::cout << get( matrix, { { 40, 8, 78, 104 } } ) << '\n' // 1.2345
              << get( matrix, { { 0, 0, 0, 0 } } ) << '\n' // -5.682
              << get( matrix, { { 99, 999, 999, 999 } } ) << '\n' // 3.34
              << get( matrix, { { 15, 100, 62, 4 } } ) << '\n' ; // 0 (not present)
}

http://coliru.stacked-crooked.com/a/a15c3709eebc592c

[48T7M4Gy]

Whether or not the modelling involves sparse matrices is beside the point especially given the earlier estimates on memory capacity required of about 2GB. Finite difference methods may or may not indicate sparse matrices and to take advantage of them (ie memory and maybe speed) if they occur only adds complexity and again loss of intuitiveness which is what the OP is asking about.

It might be simpler to create 4 arrays separately and manipulate them as the differential calculation proceeds rather than doing what appears to be all at once in a complex loop structure. it's not even clear from the above what the recursion is, possibly a misnomer maybe?

OO approach is ideal for this. Write a matrix class and manipulate the matrix objects at will. There are plenty of other finite difference/like modelling exercises done this way without a huge deal of bother.

[cire]

Whether or not the modelling involves sparse matrices is beside the point especially given the earlier estimates on memory capacity required of about 2GB.

Seems cogent to me, given the OP is running out of memory (or at least running out of large enough chunks.)

it's not even clear from the above what the recursion is, possibly a misnomer maybe?

It's pretty clear from what's posted that the manipulation is not the current problem (which is what the OP used the term recursion to qualify,) so the recursion seems to be "beside the point."

[48T7M4Gy]

Thanks, Cire, that's certainly something to consider in the final mix.