I am kind of new to C and I was doing a physics simulation in python which was taking forever to finish so I decided to switch to C , and I don t understand how to make a function which will return a 2D array (or 3D array)

#include <iostream>
#include <cmath>
// #include <complex> //

using namespace std;

double** psiinit(int L, int n, double alpha){
    double yj[400][400] = {};
    for (int j = 0; j < n; j  )
    {
        double xi[400] = {};
        for (int i = 0; i < n; i  )
        {
            xi[i] = exp(-(pow((i-(L/4)), 2)   (pow((j-(L/4)), 2)))/alpha) / (sqrt(2)*3.14159*alpha);
        };
        yj[j] = xi;
    };
    return yj;
}

int main(){
    int L = 10;
    int n = 400;
    int nt = 200*n;
    double alpha = 1;
    double m = 1;
    double hbar = 1;

    double x[n] = {};
    double y[n] = {};
    double t[nt] = {};

    double psi[nt][n][n] = {};
    psi[0] = psiinit(L, n, alpha);

    cout << psi <<endl;
    return 0;
}

I have look for answers but it doesn't seems to be for my kind of problems

Thanks

CodePudding user response：

If you're new to c you should read about the concepts of heap and stack, and about stack frames. There are a ton of good resources for that.

In short, when you declare a C-style array (such as yj), it is created in the stack frame of the function, and therefore there are no guarantees about it once you exit the frame, and your program invokes undefined behavior when it references that returned array.

There are 3 options:

1. Pass the array to the function as an output parameter (very C-style and not recommended).
2. Wrap the array in a class (like std::array already does for you), in which case it remains on the stack and is copied to the calling frame when returned, but then its size has to be known at compile time.
3. Allocate the array on the heap and return it, which seems to me to best suit your case. std::vector does that for you:

std::vector<std::vector<double>> psiinit(int L, int n, double alpha){
    std::vector<std::vector<double>> yj;
    for (int j = 0; j < n; j  )
    {
        std::vector<double> xi;
        for (int i = 0; i < n; i  )
        {
            const int value = exp(-(pow((i-(L/4)), 2)   (pow((j-(L/4)), 2)))/alpha) / (sqrt(2)*3.14159*alpha);
            xi.push_back(value);
        }

        yj.push_back(xi);
    }
    return yj;
}

If you're concerned with performance and all of your inner vectors are of a fixed size N, it might be better to use std::vector<std::array<double, N>>.

CodePudding user response：

Either make a wrapper as said above, or use a vector of vectors.

#include <vector>
#include <iostream>

auto get_2d_array()
{
    // use std::vector since it will allocate (the large amount of) data on the heap
    // construct a vector of 400 vectors with 400 doubles each
    std::vector<std::vector<double>> arr(400, std::vector<double>(400));
    arr[100][100] = 3.14;
    return arr;
}

int main()
{
    auto arr = get_2d_array();
    std::cout << arr[100][100];
}

CodePudding user response：

Your understanding of arrays, pointers and return values is incomplete. I cannot write you a whole tutorial on the topic but I recommend you read up on this.

In the mean time, I recommend you use std::vector instead of C-style arrays and treat your multidimensional arrays as 1D vectors with proper indexing, e.g. cell = vector[row * cols col]

Something like this:

#include <cmath>
// using std::exp, M_PI, M_SQRT2
#include <vector>

std::vector<double> psiinit(int L, int n, double alpha) {
    std::vector<double> yj(n * n);
    double div = M_SQRT2 * M_PI * alpha;
    for (int j = 0; j < n; j  )
    {
        double jval = j - L/4;
        jval = jval * jval;
        for (int i = 0; i < n; i  )
        {
            double ival = i - L/4;
            ival = ival * ival;
            yj[j * n   i] = std::exp(-(ival   jval) / alpha) / div;
        }
    }
    return yj;
}

Addendum: There are also specialized libraries to support matrices better and faster. For example Eigen https://eigen.tuxfamily.org/dox/GettingStarted.html

CodePudding user response：

heap allocating and returning that pointer will also work... instead of

double yj[400][400] = {}; 

do,

double** yj; 
yj = new double*[400]; 
yj[i] = new double[400];

then just,

return yj;