Violator
08-03-2005, 14:58
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;
}
}
}
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;
}
}
}