[AIUTO] Differenza C/C++

[Violator]

Sono un principiante, ho fatto questo algoritmo per il calcolo degli autovalori di una matrice 3x3.
Sotto windows (Microsoft Visual C++) viene compilato e funziona correttamente, mentre sotto linux (gcc) lo compila ma poi l'algoritmo si pianta.
E' possibile che gcc compili erroneamente il codice perchè questo è scritto in C++ e non in C?
Il codice è riportato di seguito, mi potete dire se è C o C++ e se ci sono errori vistosi.
Grazie mille a TUTTI.



/* --- Include file */
#include "nrutil.h"
#include "udf.h"
#include "mem.h"
#include <math.h>

/* --- Local definition */

#define DIM_MAT 4


/* --- Global variables */

float um[DIM_MAT][DIM_MAT], sm[DIM_MAT][DIM_MAT], om[DIM_MAT][DIM_MAT], a[DIM_MAT][DIM_MAT], z[DIM_MAT][DIM_MAT];

/*define function*/
void tred2(/*float **a,*/ int n, float d[], float e[]);
float pythag (float a, float b);
void tqli(float d[], float e[], int n/*, float **z*/);
void eigsrt ( float d[], int n );

/*
// ---[First Function]--- Reductions of a Symmetric Matrix to Tridiagonal Form (Householder reduction)
// Input:
// a[1..n][1..n] --> real symmetric matrix
// Output:
// d[1..n] --> returns the diagonal elements of the tridiagonal matrix
// e[1..n] --> returns the subdiagonal elements of the tridiagonal matrix
*/

DEFINE_ON_DEMAND(compute_streaks)

/* --- Main program */


{
float d[DIM_MAT], e[DIM_MAT];
int n=DIM_MAT;
int i=0;


/* declare domain pointer */

Domain *domain;

/* declare global thread pointer */

Thread *t;

/* declare global cell pointer */

cell_t c;

/* get domain */

domain = Get_Domain(1);

/* Print start of computing */

Message("\n\nComputing streaks...\n");

/* Loop over all cell threads in the domain */

thread_loop_c(t,domain)

{

/* Loop over all cells */

begin_c_loop(c,t)

{
i++;
um [1][1] = C_DUDX(c,t);
um [1][2] = C_DUDY(c,t);
um [1][3] = C_DUDZ(c,t);
um [2][1] = C_DVDX(c,t);
um [2][2] = C_DVDY(c,t);
um [2][3] = C_DVDZ(c,t);
um [3][1] = C_DWDX(c,t);
um [3][2] = C_DWDY(c,t);
um [3][3] = C_DWDZ(c,t);
/* Input
um [1][1] = 2.0;
um [1][2] = 0.0;
um [1][3] = 4.0;
um [2][1] = -2.0;
um [2][2] = 6.0;
um [2][3] = 2.0;
um [3][1] = 4.0;
um [3][2] = 8.0;
um [3][3] = 12.0;
*/

sm [1][1] = (um[1][1] + um[1][1])*0.5;
sm [1][2] = (um[2][1] + um[1][2])*0.5;
sm [1][3] = (um[3][1] + um[1][3])*0.5;
sm [2][1] = (um[1][2] + um[2][1])*0.5;
sm [2][2] = (um[2][2] + um[2][2])*0.5;
sm [2][3] = (um[3][2] + um[2][3])*0.5;
sm [3][1] = (um[1][3] + um[3][1])*0.5;
sm [3][2] = (um[2][3] + um[3][2])*0.5;
sm [3][3] = (um[3][3] + um[3][3])*0.5;

/*Message("\n sm"); */

om [1][1] = (um[1][1] - um[1][1])*0.5;
om [1][2] = (um[2][1] - um[1][2])*0.5;
om [1][3] = (um[3][1] - um[1][3])*0.5;
om [2][1] = (um[1][2] - um[2][1])*0.5;
om [2][2] = (um[2][2] - um[2][2])*0.5;
om [2][3] = (um[3][2] - um[2][3])*0.5;
om [3][1] = (um[1][3] - um[3][1])*0.5;
om [3][2] = (um[2][3] - um[3][2])*0.5;
om [3][3] = (um[3][3] - um[3][3])*0.5;

/*Message("\n om");*/

a [1][1] = (sm[1][1] * sm[1][1] + om[1][1] * om[1][1]) + (sm[1][2] * sm[2][1] + om[1][2] * om[2][1]) + (sm[1][3] * sm[3][1] + om[1][3] * om[3][1]);
a [1][2] = (sm[1][1] * sm[1][2] + om[1][1] * om[1][2]) + (sm[1][2] * sm[2][2] + om[1][2] * om[2][2]) + (sm[1][3] * sm[3][2] + om[1][3] * om[3][2]);
a [1][3] = (sm[1][1] * sm[1][3] + om[1][1] * om[1][3]) + (sm[1][2] * sm[2][3] + om[1][2] * om[2][3]) + (sm[1][3] * sm[3][3] + om[1][3] * om[3][3]);
a [2][1] = (sm[2][1] * sm[1][1] + om[2][1] * om[1][1]) + (sm[2][2] * sm[2][1] + om[2][2] * om[2][1]) + (sm[2][3] * sm[3][1] + om[2][3] * om[3][1]);
a [2][2] = (sm[2][1] * sm[1][2] + om[2][1] * om[1][2]) + (sm[2][2] * sm[2][2] + om[2][2] * om[2][2]) + (sm[2][3] * sm[3][2] + om[2][3] * om[3][2]);
a [2][3] = (sm[2][1] * sm[1][3] + om[2][1] * om[1][3]) + (sm[2][2] * sm[2][3] + om[2][2] * om[2][3]) + (sm[2][3] * sm[3][3] + om[2][3] * om[3][3]);
a [3][1] = (sm[3][1] * sm[1][1] + om[3][1] * om[1][1]) + (sm[3][2] * sm[2][1] + om[3][2] * om[2][1]) + (sm[3][3] * sm[3][1] + om[3][3] * om[3][1]);
a [3][2] = (sm[3][1] * sm[1][2] + om[3][1] * om[1][2]) + (sm[3][2] * sm[2][2] + om[3][2] * om[2][2]) + (sm[3][3] * sm[3][2] + om[3][3] * om[3][2]);
a [3][3] = (sm[3][1] * sm[1][3] + om[3][1] * om[1][3]) + (sm[3][2] * sm[2][3] + om[3][2] * om[2][3]) + (sm[3][3] * sm[3][3] + om[3][3] * om[3][3]);

/*Message("\n a");*/

/*First function*/
tred2( n-1, d, e );
/*tred2( a, n, d, e );*/

/*Message("\n First function");

Second function*/
tqli( d, e, n-1 );
/*tqli(d,e,n,z);*/

/*Message("\n Second function");*/

/*Third Function*/
eigsrt( d, n-1 );

/*Message("\n Third function");
*/


/*Output*/
/*if(i==1)
{
Message("\nGradient velocity tensor:");
Message("\n%3.0f %6.0f %6.0f", um[1][1], um[1][2], um[1][3]);
Message("\n%3.0f %6.0f %6.0f", um[2][1], um[2][2], um[2][3]);
Message("\n%3.0f %6.0f %6.0f\n", um[3][1], um[3][2], um[3][3]);

Message("\nEigenvalues:");
Message("\nd(1)=%.5lf", d[1]);
Message("\nd(2)=%.5lf", d[2]);
Message("\nd(3)=%.5lf\n", d[3]);
}
*/

C_UDMI(c,t,0)=d[1];
C_UDMI(c,t,1)=d[2];
C_UDMI(c,t,2)=d[3];

}

end_c_loop(c,t)

}
Message("Done.\n\n");
return( 0 );
}











void tred2(/*float **a,*/ int n, float d[], float e[])
{
int l,k,j,i;
float scale,hh,h,g,f;

for (i=n;i>=2;i--) {
l=i-1;
h=scale=0.0;
if (l > 1) {
for (k=1;k<=l;k++)
scale += (float)fabs(a[i][k]);
if (scale == 0.0)
e[i]=a[i][l];
else {
for (k=1;k<=l;k++) {
a[i][k] /= scale;
h += a[i][k]*a[i][k];
}
f=a[i][l];
g=(f >= 0.0 ? (float)-sqrt(h) : (float)sqrt(h));
e[i]=scale*g;
h -= f*g;
a[i][l]=f-g;
f=0.0;
for (j=1;j<=l;j++) {
/* Next statement can be omitted if eigenvectors not wanted */
a[j][i]=a[i][j]/h;
g=0.0;
for (k=1;k<=j;k++)
g += a[j][k]*a[i][k];
for (k=j+1;k<=l;k++)
g += a[k][j]*a[i][k];
e[j]=g/h;
f += e[j]*a[i][j];
}
hh=f/(h+h);
for (j=1;j<=l;j++) {
f=a[i][j];
e[j]=g=e[j]-hh*f;
for (k=1;k<=j;k++)
a[j][k] -= (f*e[k]+g*a[i][k]);
}
}
} else
e[i]=a[i][l];
d[i]=h;
}
/* Next statement can be omitted if eigenvectors not wanted */
d[1]=0.0;
e[1]=0.0;
/* Contents of this loop can be omitted if eigenvectors not wanted except for statement d[i]=a[i][i]; */
for (i=1;i<=n;i++) {
l=i-1;
if (d[i]) {
for (j=1;j<=l;j++) {
g=0.0;
for (k=1;k<=l;k++)
g += a[i][k]*a[k][j];
for (k=1;k<=l;k++)
a[k][j] -= g*a[k][i];
}
}
d[i]=a[i][i];
a[i][i]=1.0;
for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;
}
}

/* Temporary function !!! */

float pythag (float a, float b)
{
return(float)( sqrt(a*a + b*b) );
}

/*
// ---[Second Function]--- Eigenvalues and Eigenvectors of a Tridiagonal Matrix (QL Algorithm with Implicit Shifts)
// Input:
// d[1..n] --> diagonal elements of the tridiagonal matrix
// e[1..n] --> subdiagonal elements of the tridiagonal matrix
// Output:
// d[1..n] --> returns the eigenvalues (sorted into descending order) of matrix a
*/

void tqli(float d[], float e[], int n/*, float **z*/)
{
int m,l,iter,i,k;
float s,r,p,g,f,dd,c,b;

for (i=2;i<=n;i++) e[i-1]=e[i];
e[n]=0.0;
for (l=1;l<=n;l++) {
iter=0;
do {
for (m=l;m<=n-1;m++) {
dd=(float)fabs(d[m])+(float)fabs(d[m+1]);
if ((float)(fabs(e[m])+dd) == dd) break;
}
if (m != l) {
if (iter++ == 30) nrerror("Too many iterations in tqli", "", "");
g = (float)(( d[l+1] - d[l] ) / ( 2.0 * e[l] ));
r=pythag(g,1.0);
g=(float)(d[m]-d[l]+e[l]/(g+SIGN(r,g)));
s=c=1.0;
p=0.0;
for (i=m-1;i>=l;i--) {
f=s*e[i];
b=c*e[i];
e[i+1]=(r=pythag(f,g));
if (r == 0.0) {
d[i+1] -= p;
e[m]=0.0;
break;
}
s=f/r;
c=g/r;
g=d[i+1]-p;
r=(float)((d[i]-g)*s+2.0*c*b);
d[i+1]=g+(p=s*r);
g=c*r-b;
/* Next loop can be omitted if eigenvectors not wanted*/
for (k=1;k<=n;k++) {
f=z[k][i+1];
z[k][i+1]=s*z[k][i]+c*f;
z[k][i]=c*z[k][i]-s*f;
}
}
if (r == 0.0 && i >= l) continue;
d[l] -= p;
e[l]=g;
e[m]=0.0;
}
} while (m != l);
}
}


/* ---[Third Function]--- Reordering output
Input:
d[1..n] --> eigenvalues of matrix a
Output:
d[1..n] --> returns the eigenvalues (sorted into descending order) of matrix a
*/

void eigsrt ( float d[], int n )
{
int k,j,i;
float p;

for( i = 1; i < n; i++ )
{
p = d[ k = i ];

for( j = i + 1; j <= n; j++ )
if( d[ j ] >= p ) p = d[ k = j ];

if( k != i)
{
d[ k ] = d[ i ];
d[ i ] = p;
}
}
}

[cipi]

Originariamente inviato da Violator
[edit]...ho fatto questo algoritmo... [edit] mi potete dire se è C o C++

Infatti non penso che l'abbia fatto tu... comunque... programmi una routine da implementare in Fluent? Interessante....

Secondo me ci starebbe bene un bel int main() tipo:

 etc. etc. etc.
int main()
{
float d[DIM_MAT], e[DIM_MAT];
int n=DIM_MAT;
int i=0;
etc. etc. etc.

poi per il resto mi pare ok... (chiedo conferma a chi ne sa di più...) ricorda di mettere tutte le librerie (Numerical Recepies.... etc.) e di nominarlo miofile.cpp Poi usa gcc -lm miofile.cpp e se vuoi metti qualche ottimizzazione tipo -O3
Mandi!

[Violator]

Errata corrige: l'algoritmo non l'ho fatto io, peròà l'ho dovuto modificare in qua e in la...
Anche tu usi Fluent e carichi le udf al suo interno?

Faccio la modifica e ti so dire.

Grazie mille!