[C] Metodo Jacobi in linguaggio C

[eele1]

Ciao a tutti! Per un progetto universitario devo convertire alcuni programmi scritti in MatLab nei corsi di Analisi Numerica in linguaggio C.
Momentaneamente sto lavorando al metodo di risoluzione dei sistemi lineari tramite Jacobi.
Ho scritto il codice in C ma non riesco a farlo funzionare. Tenete conto che sono alle primissime armi con il linguaggio C. Vi posto il codice, se qualcuno ha suggerimenti sono ben accetti!
Il compilatore comunque non da errori di sintassi!

//Algoritmo di Jacobi [Chiamata MatLab: [x,its]=jacobi(A,b,x_0,maxit,tol)]
/* L'algoritmo deve prendere in input una matrice quadrata A con elementi della diagonale non nulli,
un vettore (termine noto) b, una approssimazione della soluzione iniziale x_0, e i parametri d'arresto
dati da maxit= numero massimo di iterezioni consentite e tol=tolleranza errore.
L'algoritmo deve restituire il vettore soluzione x e il parametro its delle iterazioni compiute.*/

#include<stdio.h>
#include<math.h>
#define SIZE 50

double molt(size_t dim, double mat[dim][dim],double vet[dim]);
double norma(size_t dim,double vet[dim]);
double molt_mat(size_t dim,double mat1[dim][dim],double mat2[dim][dim]);
double somma_vet(size_t dim, double vec1[dim],double vec2[dim]);

int main(){

size_t dim;
int i,j,its,maxit;
double tol,norm;


printf("-----METODO JACOBI-----\n");
printf("Si inseriscano i dati nella maniera richiesta.\n");
printf("Si inserisca la dimensione della matrice A.\n");
scanf("%d", &dim);

double A[dim][dim];
double b[dim];
double x0[dim];
double x[dim];
double err[dim];

printf("Si inseriscano gli elementi della matrice A nell'ordine rischiesto.\n");
for (i=0; i<dim; i++){
        for(j=0;j<dim;j++){
    printf("Inserisci l'elemento di posto %d , %d : \n",i+1,j+1);
    scanf("%lf", &A[i][j]);
        }
}

//controllo sugli elementi di A che nn ci siano elementi della diag nulli
for (i=0;i<dim;i++){
    if(A[i][i]==0.00){
        printf("Errore: non e' possibile applicare il metodo di Jacobi.");
        return 0;
    }
}



printf("Si inseriscano gli elementi del termine noto b nell'ordine richiesto.\n");
for (i=0; i<dim; i++){
    printf("Inserisci l'elemento di posto %d: \n",i+1);
    scanf("%lf", &b[i]);
}

printf("Si inseriscano gli elementi del vettore iniziale x_0.\n");
for (i=0; i<dim; i++){
    printf("Inserisci l'elemento di posto %d: \n",i+1);
    scanf("%lf", &x0[i]);
}

printf("Si inserisca il numero di iterazioni massime.\n");
scanf("%d",&maxit);

printf("Si inserisca la tolleranza.\n");
scanf("%g", &tol);
    printf("toll=%lg",&tol);

//---------------------------------------------
//creo i pezzi per l'iterazione

//inizializzo x a x0
for(i=0;i<dim;i++){
    x[i]=x0[i];
}

//calcolo la norma iniziale
for(i=0;i<dim;i++){
    err[i]=b[i]-molt(dim,A,x);
}
norm=norma(dim,err);
printf("norma %g",norm);

//creo matrice diagonale

double D[dim][dim];
for(i=0;i<dim;i++){
    for(j=0;j<dim;j++){
            if(i==j){
                D[i][j]=A[i][i];
            }
            else{
                D[i][j]=0.00;
            }
        }
}

//Creo D^-1
double D_inv[dim][dim];
for (i=0;i<dim;i++){
    D_inv[i][i]=1/D[i][i];
}

//creo N=A-D
double N[dim][dim];
for (i=0;i<dim;i++){
    for(j=0;j<dim;j++){
        if(i==j){
            N[i][j]=0.00;
        }
        else{
            N[i][j]=A[i][j];
        }
    }
}
//creo T=D^-1*N
double T[dim][dim];
T[dim][dim]=molt_mat(dim,D_inv,N);

//creo S=D_inv*b
double S[dim];
S[dim]=molt(dim,D_inv,b);


//ITERAZIONE JACOBI
int iter=0;
while(iter<maxit && norm>tol){
double temp[dim];
temp[dim]=molt(dim,T,x);
x[dim]=somma_vet(dim,temp,S);

iter++;
for(i=0;i<dim;i++){
    err[i]=b[i]-molt(dim,A,x);
}
norm=norma(dim,err);

}


//presentazione risultati


printf("La soluzione calcolata e':\n");
for (i=0;i<dim;i++){
    printf("%f ",x[i]);
    printf("\n");
}
printf("\n");

printf("le iterazioni fatte sono %d \n",iter);
return 0;
}







//_________FUNZIONI A SUPPORTO___________

double molt(size_t dim, double mat[dim][dim], double vet[dim]){

int i,j;
double sum;
double RIS[dim];

for(i=0;i<dim;i++){
    sum=0.00;
    for(j=0;j<dim;j++){
        sum+=mat[i][j]*vet[j];
    }
    RIS[i]=sum;
}

return RIS[dim];
}

//____________

double norma(size_t dim, double vet[dim]){

int i;
double sum,norm;

sum=0.00;
for(i=0;i<dim;i++){
    sum+=pow(vet[i],2);
}

norm=fabs(sqrt(sum));
return norm;
}

//__________
double molt_mat(size_t dim,double mat1[dim][dim],double mat2[dim][dim]){

double RIS[dim][dim];
double sum;
int n,i,j,k;


for(i=0;i<dim;i++){
        for(j=0;j<dim;j++){
            sum=0.00;
            for(n=0;n<dim;n++){
                sum+=mat1[i][n]*mat2[n][j];
            }
            RIS[i][j]=sum;
        }
}
return RIS[dim][dim];
}

//___________
double somma_vet(size_t dim, double vec1[dim],double vec2[dim]){

double RIS[dim];
int i;


for(i=0;i<dim;i++){
    RIS[i]=vec1[i]+vec2[i];
    }

return RIS[dim];
}

[apatriarca]

Ho guardato solo velocemente, ma il valore di ritorno delle funzioni di supporto è al di fuori dell'array a cui appartiene (e le conseguenze possono quindi andare dalla lettura di un valore casuale ad un crash del programma). Che valore stai effettivamente cercando di restituire?

[eele1]

Ogni funzione esterna dovrebbe restituirmi il risultato dell'operazione. Quindi, ad esempio la matrice risultante dalla moltiplicazione tra due matrici.
Infatti per ora il programma gira e non da errori ma se faccio stampare le matrici create con le operazioni implementate ci sono dei valori casuali!

Di preciso cosa intendi per "al di fuori dell'array a cui appartengono?"

[apatriarca]

Stai restituendo un singolo valore, non una matrice e questo valore non appartiene alla matrice.

[eele1]

Sono d'accordo, il problema è lì... Però quando implemento le funzioni a supporto gli dico di ritornare il risultato nella dimensione matriciale o vettoriale dandogli un comando tipo

 return RIS[dim][dim]; 

e assegnando poi il risultato a una variabile oppurtunamente definita.

Dove sbaglio?
Grazie!

[apatriarca]

Sbagli che quella istruzione non restituisce una matrice, ma il valore della matrice agli indici (dim, dim) che non è valido. Non è possibile restituire una matrice in C. Puoi restituire un puntatore (ma devi allora allocare la matrice dinamicamente) o passare il risultato come argomento della funzione.

[vict85]

Per tradurre un codice da un linguaggio ad un altro devi conoscere molto bene il linguaggio di partenza e soprattutto devi conoscere ancora meglio quello di arrivo. Per scrivere un buon codice di questo tipo devi anche conoscere l'algoritmo.

Detto questo non tradurrei da Matlab, ma reimplementerei l'algoritmo. A meno di sapere molto bene cosa fa matlab al suo interno convertire un codice matlab in c produce un pessimo codice c (in termini di performance e/o di leggibilita).

[eele1]

Grazie Mille a tutti per i suggerimenti, sono riuscita a risolvere! Avevo sbagliato a passare i dati da un pezzo di codice all'altro...!

Sono riuscita a togliere gli errori e ad alleggerire molto il codice!

Alla prossima!