Segmentation fault on encountering expression GSL_SET_COMPLEX

[Mehul Dixit]

I am writing a code in C++ to solve a problem in electrodynamics. The code is to diagonalize a matrix M with complex elements. I am using the GSL library functions for representing complex numbers and the diagonalization. The matrix M is defined in the code by the statement

gsl_complex M[dimension][dimension] //dimension is an int defined and initialized earlier in the code

Before diagonalization of the matrix M, I initialize each element of it to 0 using the statements
for (int p = 0; p <= (dimension - 1); p++){
for (int q = 0; q <= (dimension - 1); q++){
GSL_SET_COMPLEX(&M[p][q], 0, 0);}}

The program compiles without error but on execution gives a segmentation fault. Using GDB the segmentation fault is traced to the line GSL_SET_COMPLEX(&M[p][q],0,0); which sets the real and imaginary parts of the p,q-th element of M, M[p][q], to 0.

Wondering if anyone here has any clues on why the program is exiting with a segmentation fault. I seem to be using the correct syntax for setting the matrix elements (http://www.gnu.org/s/gsl/manual/html_node/Representation-of-complex-numbers.html)

Thanks

[weaknessforcats]

I think your inner loop needs to run from 0 to p-1 since a pq array has p elements that are q qrrays. Therefore each p is really a q array numbered from 0 to q-1 and not dimension-1.

[Mehul Dixit]

Hi there. Thanks for your reply. Here M is matrix of size (dimension)x(dimension) and is defined as so in the statement gsl_complex M[dimension][dimension]. The pq-th element of M in c++ is referenced as M[p][q] (see section on Multidimensional arrays in http://www.cplusplus.com/doc/tutorial/arrays/ for e.g.). Doesn't the for loop then as written access each element of the matrix and set it to 0 (think of p as the row index and q as the column index) ?

[weaknessforcats]

You are correct.

In fact I ran your code an there is no segmentation fault when I run the program.

Are you sure you have created the array?

Here's what I used:

struct Array
{
   int data;
};
#define dimension 10
 
void GSL_SET_COMPLEX(Array* arg, int real, int imag)
{
      arg->data = real;
}
 
int main()
{
Array M[dimension][dimension];
for (int p = 0; p <= (dimension - 1); p++){
 for (int q = 0; q <= (dimension - 1); q++){
 GSL_SET_COMPLEX(&M[p][q], 0, 0);}}
 
for (int p = 0; p <= (dimension - 1); p++){
 for (int q = 0; q <= (dimension - 1); q++){
 cout << M[p][q].data << endl;}}
 
 
}

[donbock]

As I understand it, GSL_SET_COMPLEX is a function-like macro. If it were a function, then its prototype would be:

void GSL_SET_COMPLEX(gsl_complex *p, double realvalue, double imaginaryvalue);

weaknessforcats: you passed the address of an Array rather than of a gsl_complex. Why did that work?
Oops. I see, you used a local definition of GSL_SET_COMPLEX.

M is declared as a two-dimensional array of gsl_complex; or more pedantically, as an array of one-dimensional arrays of gsl_complex. I don't use 2-dimensional arrays enough to remember: is &M[2][3] the address of one gsl_complex or is it the address of an array of gsl_complex? (I would have assumed this works the way you expect, but when something goes wrong you have to examine your assumptions.)
There's a great article in the Insights tab called Arrays Revealed, but I can't seem to call it up.

I suggest you verify you are passing proper arguments to GSL_SET_COMPLEX by adding a function prototype for it (but under a different name, perhaps MY_GSL_SET_COMPLEX) and then calling that new name instead. You will get compiler errors if any of the arguments are of the wrong type.

[donbock]

Are you sure dimension is greater than 0?

Take a look at gsl/gsl_complex.h and gsl/gsl_complex_math.h to find out how gsl_complex and GSL_SET_COMPLEX are defined.

[weaknessforcats]

@donbock: It worked because I wrote my own GSL_SET_COMPLEX function. I was trying to show that the loops did not cause the segmentation fault. That meant the fault was in the creation of the array or the processing of it.

[Mehul Dixit]

Thanks for the replies.

@weaknessforcats:
> Are you sure you have created the array?
Doesn't the declaration statement gsl_complex M[dimension][dimension]; create the array ?

@donblock:
>is &M[2][3] the address of one gsl_complex or is it the address of an array of gsl_complex?
I believe it is the address of one gsl_complex
>Are you sure dimension is greater than 0?
Yes. I printed out the value of dimension before the declaration of M[dimension][dimension] and the for loop to check whether dimension is being initialized correctly and indeed it is (in my program the value of dimension is 876)
> I suggest you verify you are passing proper arguments to GSL_SET_COMPLEX
I believe I am passing the proper arguments. I will try out your suggestion of defining a function prototype and post the results.

Thanks

[weaknessforcats]

If your dimension is 876 x 876 you will need 767376 elements. Probably this is too large to fit on the stack. I suggest you change your code to create the array on the heap using the new operator:

gsl_complex (* M)[dimension]= new gas_complex[dimension][dimension];
 

[Mehul Dixit]

Yes, the dimension indeed seems too large. I reduced the dimension to 243x243 and the code executes without segmentation fault.

I do not quite know what you mean by create the array on the heap ?

[weaknessforcats]

When you call a function, the compiler generates code to allocate memory for the variables the function needs. This memory is called a stack frame. If the called function itself calls a function another stack frame is created. These stack frames are managed using a data strcture called a heap.

So in programmerspeak,a variable is "on the stack" when the compiler creates it.

There may be some restrictions on the maximum size of a stack frame.

You may also create variables by allocating your own memory at run time. When you do this, the memory allocation is also maintained in a data structure called a heap. Hence, when you allocate your own memory, it is said to be "on the heap".

There are no size restrictions for allocating on the heap.

You use the C++ "new" operator to allocate your own memory.

Just remember that any memory you allocate you must personally release when you no longer need it. Use the C++ "delete" operator to do that.

The code sinppet on my last poat creates your array "on the heap" because it uses the "new" operator.

Lastly, stack and heap are compiler design concepts and not part of the C++ language. All C++ requires is that you and the compiler must be able to create variables.