sorry, forgot attachment

[Thorsten Bagdonat <t dot bagdonat at tu-bs dot de>]

I have attached a C++ code which I'm currently working on.
If it is compiled using -finline-functions it won't work
correctly. Otherwise it does fine. 

I figured out, that this has to have something to do with
the class 'vfield', especially the 
double* vfield::operator () (int x,int y,int z) {
  return data+x*ydim*zdim*3+y*zdim*3+z*3;
}
function. Whenever this is processed using -finline-functions
the program behaves strangely.

I have to say that I'm doing C programming for quite a while, but C++
only for some weeks, so maybe, there is something simple I overlooked. 

Is this a bug?

        Thorsten Bagdonat
        Inst. f. Theoretical Physics
        TU Braunschweig
        D-38102 Braunschweig

        Tel.: +49-(0)531-391-5187
        email: t.bagdonat@tu-bs.de
#include <stdio.h>
//#include <matpack.h>
#include <time.h>
#include <math.h>
#include <stdlib.h>

long  total_mem_use = 0;
clock_t  old_time;
FILE*  prot=stderr;

void start_timer() {
  old_time = clock();
}

void eval_timer(FILE* outfile) {
  clock_t time_used;
  
  time_used=clock()-old_time;
  fprintf (outfile," %4.5f sec.\n",(double)(1.0*time_used/CLOCKS_PER_SEC));
}

void errtxt (char *fname, char *txt) {
  fprintf (stderr,"%s: %s\n",fname,txt);
  exit (1);
}

class vfield {
private:
  double *data;   // 3-dim array of 3-comp. vectors 
  int    xdim,ydim,zdim;
public:

  void    init (int, int, int);
  void    done (void);
  double* operator () (int,int,int);
  void    show (FILE*);
};

void  vfield::show(FILE* outfile) {
  int x,y,z,c;
  double *v;

  for (z=0; z<zdim; z++) {
    for (y=0; y<ydim; y++) {
      for (x=0; x<xdim; x++) {
	v = data+x*ydim*zdim*3+y*zdim*3+z*3;
	fprintf (outfile,"(%3.2f,%3.2f,%3.2f) ",v[0],v[1],v[2]);
      }
      fprintf (outfile,"\n");
    } 
    fprintf (outfile,"\n");
  }
}

double* vfield::operator () (int x,int y,int z) {
  return data+x*ydim*zdim*3+y*zdim*3+z*3;
}

void vfield::init (int Xdim, int Ydim, int Zdim) {
  char *fname="vfield::Init";

  xdim=Xdim;
  ydim=Ydim;
  zdim=Zdim;
  if ((data=(double*)malloc(sizeof(double)*xdim*ydim*zdim*3)) == NULL) 
    errtxt (fname,"not enough memory");
  fprintf (prot,"%s: allocated %d bytes of memory\n",fname,
	   sizeof(double)*xdim*ydim*zdim*3);
  total_mem_use += sizeof(double)*xdim*ydim*zdim*3;
}

void vfield::done () {
  char *fname="vfield:free";
  
  free (data);
}

// vector routines
void vec_sub(double *res, double *v1, double *v2) {
  int i;

  for (i=0; i<3; i++) 
    res[i]=v1[i]-v2[i];
}

double vec_scalar(double *v1, double *v2) {
  int i;
  double erg=0;
  
  for (i=0; i<3; i++)
    erg+=v1[i]*v2[i];
  return erg;
}

void vec_cross(double *res, double *v1, double *v2) {
  res[0] = v1[1]*v2[2]-v1[2]*v2[1];
  res[1] = v1[2]*v2[0]-v1[0]*v2[2];
  res[2] = v1[0]*v2[1]-v1[1]*v2[0];
}

void vec_copy(double *from, double *to) {
  int i;
  for (i=0; i<3; i++)
    to[i]=from[i];
}

// grid data
vfield gr;           // even integer grid points
vfield cob[3];         // co-variant basis-vectors
int    g_ngp[3]={4,4,4};  // number of even integer grid points
double qmin[3]={5,0,0};  // min,max values for curvilinear grid coordinates
double qmax[3]={10,M_PI,M_PI};

// magnetic field
vfield A,B;

void init_vfields() {
  char *fname="init_vfields()";
  int  i;

  fprintf (prot,"%s:initializing vector field 'gr' (%dx%dx%d).\n",
	   fname,g_ngp[0],g_ngp[1],g_ngp[2]);
  gr.init (g_ngp[0],g_ngp[1],g_ngp[2]);

  for (i=0; i<3; i++) {
    fprintf (prot,"%s:initializing vector field 'cob[%d]' (%dx%dx%d).\n",
	     fname,i,g_ngp[0],g_ngp[1],g_ngp[2]);
    cob[i].init (g_ngp[0],g_ngp[1],g_ngp[2]);
  }

  fprintf (prot,"%s:initializing vector fields 'A','B' (%dx%dx%d).\n",
	   fname,g_ngp[0],g_ngp[1],g_ngp[2]);
  B.init (g_ngp[0],g_ngp[1],g_ngp[2]);
  A.init (g_ngp[0],g_ngp[1],g_ngp[2]); 

  fprintf (prot,"%s:total mem usage: %d bytes.\n",fname,total_mem_use);
}

void free_vfields() {
  char *fname="free_vfields()";
  int  i;

  fprintf (prot,"%s:free memory of vector field 'gr'.\n",fname);
  gr.done;

  for (i=0; i<3; i++) {
    fprintf (prot,"%s:free memory of vector field 'cob[%d]'.\n",fname,i);
    cob[i].done;
  }

  fprintf (prot,"%s:free memory of vector field 'A','B'.\n",fname);
  B.done;
  A.done;
}

// These functions declare the coordinate transformation
// The drx(), dry(), drz() are the partial differentiations
// qisel selects the variable after which the differentiation is to be made
double rx(double *q) {
  return (q[0]*cos(q[1])*sin(q[2]));
}

double ry(double *q) {
  return (q[0]*sin(q[1])*sin(q[2]));
}

double rz(double *q) {
  return (q[0]*cos(q[2]));
}

//  double rx(double *qi) {
//    return (qi[0]);
//  }

//  double ry(double *qi) {
//    return (qi[1]);
//  }

//  double rz(double *qi) {
//    return (qi[2]);
//  }

void generate_grid() {
  char   *fname="generate_grid()";
  int    qn[3],i,j;
  double qi[3],qileft[3],qiright[3],grleft[3],grright[3];

  fprintf (prot,"%s:generating even integer grid points.",fname);
  start_timer();
  for (qn[0]=0; qn[0]<g_ngp[0]; qn[0]++)
    for (qn[1]=0; qn[1]<g_ngp[1]; qn[1]++)
      for (qn[2]=0; qn[2]<g_ngp[2]; qn[2]++) {
	for (i=0; i<3; i++) 
	  qi[i] = (qmax[i]-qmin[i])/(g_ngp[i]-1)*qn[i]+qmin[i];
	gr(qn[0],qn[1],qn[2])[0] = rx(qi);
	gr(qn[0],qn[1],qn[2])[1] = ry(qi);
	gr(qn[0],qn[1],qn[2])[2] = rz(qi);
      }
  eval_timer(prot);

  fprintf (prot,"%s:calculating covariante basis vectors.",fname);
  start_timer();

  for (qn[0]=0; qn[0]<g_ngp[0]; qn[0]++)
    for (qn[1]=0; qn[1]<g_ngp[1]; qn[1]++)
      for (qn[2]=0; qn[2]<g_ngp[2]; qn[2]++) {
	for (i=0; i<3; i++) {
	  // i indicates the actual direction to be differentiated
	  // we need the cartesian coordinates of the halg integer
	  // grid points in direction i
	  for (j=0; j<3; j++) {
	    if (i==j) {
	      // for the diff. direction the half integer qn have to be chosen
	      qileft[j] = (qmax[j]-qmin[j])/(g_ngp[j]-1)*(qn[j]-0.5)+qmin[j];
	      qiright[j] = (qmax[j]-qmin[j])/(g_ngp[j]-1)*(qn[j]+0.5)+qmin[j];
	    } else {
	      // for the other directions the full integer qn are chosen
	      qileft[j] = qiright[j] = 
		(qmax[j]-qmin[j])/(g_ngp[j]-1)*(qn[j])+qmin[j];
	    }
	  }
	  grright[0] = rx(qiright);
	  grright[1] = ry(qiright);  // to cartesian coordinates
	  grright[2] = rz(qiright);
	  grleft[0] = rx(qileft);
	  grleft[1] = ry(qileft);
	  grleft[2] = rz(qileft);
	  
	  vec_sub(cob[i](qn[0],qn[1],qn[2]),grright,grleft);
	  // the basis vector in direction i is the difference of the two
	  // half integer grid points.
	}
      }
  eval_timer(prot);
}

double Bx(double *x) {
  return (x[1]);
}
double By(double *x) {
  return (-x[0]);
}
double Bz(double *x) {
  return (0.0);
}

void generate_Bfield() {
  // generates a field in cartesian components for each grid point
  // from analytic functions
  int    qn[3];
  double *pos;

  for (qn[0]=0; qn[0]<g_ngp[0]; qn[0]++)
    for (qn[1]=0; qn[1]<g_ngp[1]; qn[1]++)
      for (qn[2]=0; qn[2]<g_ngp[2]; qn[2]++) {
	pos = gr(qn[0],qn[1],qn[2]);  // cartesian position. indexed 1..3
	B(qn[0],qn[1],qn[2])[0] = Bx(pos);
	B(qn[0],qn[1],qn[2])[1] = By(pos);
	B(qn[0],qn[1],qn[2])[2] = Bz(pos);
      }
}

void get_co_basis (int ijk[3],double* basev[3]) {
  // calc the three covariant basis vectors basev[i][comp] at grid point 
  // ijk[comp]
  int i;

  for (i=0; i<3; i++)
    basev[i] = cob[i](ijk[0],ijk[1],ijk[2]);
}

void cart2cov (vfield a) {
  // assuming vector field 'a' given in cartesian components, 'a' is
  // converted to covariant components with respect to the coordinates
  // defined by the covariante basis vectors cob[]
  double* b[3];   // three covariant basis vectors
  double scratch[3];
  int    qn[3],i;

  for (qn[0]=0; qn[0]<g_ngp[0]; qn[0]++)
    for (qn[1]=0; qn[1]<g_ngp[1]; qn[1]++)
      for (qn[2]=0; qn[2]<g_ngp[2]; qn[2]++) {
	get_co_basis (qn,b);
	// covariant components simply by scalar multiplication
	for (i=0; i<3; i++)
	  scratch[i] = vec_scalar (a(qn[0],qn[1],qn[2]),b[i]);
	// put covariant components back in vector field
	vec_copy (scratch,a(qn[0],qn[1],qn[2]));
      }
}

void contra2cart (vfield a) {
  // converts vector field a from contravariant coordinates to cartesian
  // coordinates
  double* b[3];   // three covariant basis vectors
  double scratch[3];
  int    qn[3],i,j;

  for (qn[0]=0; qn[0]<g_ngp[0]; qn[0]++)
    for (qn[1]=0; qn[1]<g_ngp[1]; qn[1]++)
      for (qn[2]=0; qn[2]<g_ngp[2]; qn[2]++) {
	get_co_basis (qn,b);
	// remember contravariant components
	vec_copy (a(qn[0],qn[1],qn[2]),scratch);
	// A_cart,i = A^1 b_1i + A^2 b_2i + A^3 b_3i
	for (i=0; i<3; i++)
	  a(qn[0],qn[1],qn[2])[i] = 
	    scratch[0]*b[0][i]+scratch[1]*b[1][i]+scratch[2]*b[2][i];
      }
}

double cell_volume (int ijk[3]) {
  // calculates b1*(b2xb3);
  double* b[3];
  double b2xb3[3];

  get_co_basis (ijk,b);
  vec_cross (b2xb3,b[1],b[2]);
  return vec_scalar(b[0],b2xb3);
}

void  rot(vfield a, vfield b) {
  // calculates a=curl b . b must be given in covariant 
  // components. The result a is given in contravariant components.
  // The calculation is done only for the grid points 1..dim-1 for
  // each dimension. The other values have to be applied using some boundary
  // condition.
  double scratch[3];
  double S;
  int    qn[3],i;

  for (qn[0]=1; qn[0]<g_ngp[0]-1; qn[0]++)
    for (qn[1]=1; qn[1]<g_ngp[1]-1; qn[1]++)
      for (qn[2]=1; qn[2]<g_ngp[2]-1; qn[2]++) { 
	a(qn[0],qn[1],qn[2])[0] =
	  0.5*(b(qn[0],qn[1]+1,qn[2])[2]-b(qn[0],qn[1]-1,qn[2])[2]-
	       b(qn[0],qn[1],qn[2]+1)[1]+b(qn[0],qn[1],qn[2]-1)[1]);
	a(qn[0],qn[1],qn[2])[1] = 
	  0.5*(b(qn[0],qn[1],qn[2]+1)[0]-b(qn[0],qn[1],qn[2]-1)[0]-
	       b(qn[0]+1,qn[1],qn[2])[2]+b(qn[0]-1,qn[1],qn[2])[2]);
	a(qn[0],qn[1],qn[2])[2] = 
	  0.5*(b(qn[0]+1,qn[1],qn[2])[1]-b(qn[0]-1,qn[1],qn[2])[1]-
	       b(qn[0],qn[1]+1,qn[2])[0]+b(qn[0],qn[1]-1,qn[2])[0]);
	S = cell_volume(qn);
	for (i=0; i<3; i++)  a(qn[0],qn[1],qn[2])[i] /= S;
      }
}

main() {
  prot=stderr;
  init_vfields();

  generate_grid();
  // gr.show(stdout);
  // grh.show(stdout);
  cob[0].show(stdout);
  cob[1].show(stdout);
  cob[2].show(stdout);

  generate_Bfield();
  B.show(stdout);
  fprintf (prot,"cart2cov(B):");
  start_timer();
  cart2cov(B);
  eval_timer(prot);
  B.show(stdout);
  fprintf (prot,"rot(A,B):");
  start_timer();
  rot(A,B);
  eval_timer(prot);
  A.show(stdout);
  fprintf (prot,"contra2cart (A):");
  start_timer();
  contra2cart (A);
  eval_timer(prot);
  A.show(stdout);

  free_vfields();
}