Ax=b Jacobi 100000x100000 sparse matrix

[hvmclrhu]

Hi I have a big problem. When we compile serial.c with gcc, I get this error
program is generating the sparse matrix
Segmentation fault

I think ı have to use malloc() but I don't know how to use and add this function in my program. Could you help me please? Thank you for your help...


[code=c]#include <stdio.h>
#include <math.h>

#define MAX_DIM 100000

typedef float MATRIX_T[MAX_DIM][MAX_DIM];

int Jacobi(MATRIX_T A, float x[], float b[], int n,
float tol, int max_iter);
void Read_matrix(cha r* prompt, MATRIX_T A, int n);
void Read_vector(cha r* prompt, float x[], int n);
void Print_vector(ch ar* title, float x[], int n);
void Print_matrix(ch ar* title, MATRIX_T A, int n);

main(int argc, char* argv) {
MATRIX_T A;
float x[MAX_DIM];
float b[MAX_DIM];
int n;
float tol;
int max_iter;
int converged;
n=100000;
tol=0.0000001;
max_iter=10000;
Read_matrix("pr ogram is generating the sparse matrix", A, n);
Read_vector("pr ogram is generating the sparse vector", b, n);

converged = Jacobi(A, x, b, n, tol, max_iter);

if (converged)
Print_vector("T he solution is", x, n);
else
printf("Failed to converge in %d iterations\n", max_iter);
} /* main */


/*************** *************** *************** *************** *********/
/* Return 1 if iteration converged, 0 otherwise */
/* MATRIX_T is just a 2-dimensional array */
int Jacobi(
MATRIX_T A /* in */,
float x[] /* out */,
float b[] /* in */,
int n /* in */,
float tol /* in */,
int max_iter /* in */) {
int i, j;
int iter_num;
float x_old[MAX_DIM];

float Distance(float x[], float y[], int n);

/* Initialize x */
for (i = 0; i < n; i++)
x[i] = b[i];

iter_num = 0;
do {
iter_num++;

for (i = 0; i < n; i++)
x_old[i] = x[i];

for (i = 0; i < n; i++){
x[i] = b[i];
for (j = 0; j < i; j++)
x[i] = x[i] - A[i][j]*x_old[j];
for (j = i+1; j < n; j++)
x[i] = x[i] - A[i][j]*x_old[j];
x[i] = x[i]/A[i][i];
}
} while ((iter_num < max_iter) &&
(Distance(x,x_o ld,n) >= tol));

if (Distance(x,x_o ld,n) < tol)
return 1;
else
return 0;
} /* Jacobi */


/*************** *************** *************** *************** *********/
float Distance(float x[], float y[], int n) {
int i;
float sum = 0.0;

for (i = 0; i < n; i++) {
sum = sum + (x[i] - y[i])*(x[i] - y[i]);
}
return sqrt(sum);
} /* Distance */


/*************** *************** *************** *************** *********/
void Read_matrix(
char* prompt /* in */,
MATRIX_T A /* out */,
int n /* in */) {
int i, j;
printf("%s\n", prompt);
for (i = 0; i < n; i++)
for (j = 0; j < n; j++){
if( rand()/2147483647.0 > 0.7 )
A[i][j] = rand() /2147483647.0;
else A[i][j] = 0;
}

} /* Read_matrix */


/*************** *************** *************** *************** *********/
void Read_vector(
char* prompt /* in */,
float x[] /* out */,
int n /* in */) {
int i;

printf("%s\n", prompt);
for (i = 0; i < n; i++){
if( rand()/2147483647.0 > 0.7 )
x[i] = rand() /2147483647.0;
else x[i] = 0;
}
} /* Read_vector */


/*************** *************** *************** *************** *********/
void Print_vector(
char* title /* in */,
float x[] /* in */,
int n /* in */) {
int i;

printf("%s\n", title);
for (i = 0; i < n; i++)
printf("%4.1f ", x[i]);
printf("\n");
} /* Print_vector */

/*************** *************** *************** *************** *********/
void Print_matrix(
char* title /* in */,
MATRIX_T A /* in */,
int n /* in */) {
int i,j;

printf("%s\n", title);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
printf("%f", A[i][j]);
printf("\n");
}
} /* Print_matrix */[/code]

[JosAH]

As far as I can tell there's nothing sparse about that 100,000x100,000 matrix,
i.e. every element is there and it's taking up too much space no matter what
you want to do about it as it is now.

kind regards,

Jos

[AdrianH]

Er, your stack has exploded. ;) You created an array that is about 4.6 TiB in size. That is over 4768 GiB!!! Good luck with that. :D

You defiantly need to implement a sparse array type. And you probably want to not pass it by value but by indirection (pass a pointer to the type).


Adrian

[weaknessforcats]

100000 x 100000 is probably too big for a stack variable. You will need to allocate this on the heap using malloc(). Check your documentation on how to use this function.

[AdrianH]

Oh, my bad. My calculation is wrong. Your actually allocating a little over 37 GiB. Not as bad as TiB ;).


Adrian

[Savage]

Oh, my bad. My calculation is wrong. Your actually allocating a little over 37 GiB. Not as bad as TiB ;).


Adrian

But still,very very bad.

Savage

[emaghero]

Jacobi Iteration with that size of Matrix could take about a week to converge especially since you're using sparse matrices. Get yourself a better solver or you will be waiting a long long time for a solution.

Another thing you could also do, if you insist using archaic numerical methods, is to start with a small example that you know the solution of, say 2*2, or 3*3, and then work your way up to 100000*100000.

If you know the sturcture of your matrices, tri-diagonal, band-diagonal etc you SHOULD find a more suitable solver than Jacobi Iteration.