C++ VEC函数代码示例

int
main(int argc, char **argv)
{
  const int n = 3;

  // init our vectors of size n
  struct vector *x = vector_create(n);
  struct vector *y = vector_create(n);
  struct vector *z = vector_create(n);
  struct vector *z_ref = vector_create(n);

  // initialize values for testing
  for (int i = 0; i < n; i++) {
    VEC(x, i) = i + 1.f;
    VEC(y, i) = i + 2.f;
    VEC(z_ref, i) = 2*i + 3.f;
  }

  // test the vector_add()
  vector_add(x, y, z);

  // and make sure it equals the reference result
  assert(vector_is_equal(z, z_ref));

  // clean up
  vector_destroy(x);
  vector_destroy(y);
  vector_destroy(z);
  vector_destroy(z_ref);

  return 0;
}

/*
 * write_data
 *
 * Takes a FILE*, a data vector, a grid vector 
 * a size and a mode.
 * mode 0: write real values
 * mode 1: write complex values
 * mode 2: write both real and comples values
 */
void write_data( FILE *fp, vector *v, vector *grid, 
                      const int n, int mode )
{
  int i;
  for( i = 0; i < n; i++ )
    fprintf( fp, "%g %g\n", creal(VEC( grid, i )), creal(VEC( v, i )) );
}

void vv_remove(void *xx, int size, int i, int n) {
    char *x = xx;
    assert(size == VEC(x)->size);

    memmove(x+i*size, x+(i+n)*size, size * (VEC(x)->nelt-i-n));
    VEC(x)->nelt -= n;
}

void *vv_insert(void *xx, int size, int i, void *e, int n V_DEBUGPARAMS) {
    char *x = xx;
    assert(!x || size == VEC(x)->size);

    x=vv_growk(x, size, n V_DEBUGARGS);
    memmove(x+(i+n)*size, x+i*size, size*(VEC(x)->nelt-n-i));
    memcpy(x+i*size, e, n*size);
    return x;
}

static CVECTOR *VECTOR_copy(CVECTOR *_object)
{
	CVECTOR *copy = VECTOR_create(SIZE(THIS), COMPLEX(THIS), FALSE);
	if (!COMPLEX(THIS))
		gsl_vector_memcpy(VEC(copy), VEC(THIS));
	else
		gsl_vector_complex_memcpy(CVEC(copy), CVEC(THIS));
	
	return copy;
}

void pdu_pdt( Vector *u, Vector *s, double dx )
{
    int i;
    double dx_3 = pow(dx,3);
    #pragma omp parallel for
    for( i = 0; i < u -> N; i++ )
    { 
        VEC(s,i) = -6 *VEC(u,i) * du_x(u,i,dx) -du_xxx(u,i,dx_3);
    }
}

/*
 * fft
 *
 * Takes an input vector, output vector and a number of points
 * Performs a Fast Fourier Transform on the input vector
 * and stores it in the output vector.
 */
void fft( vector *in, vector *out, int n )
{
  int k;
  if( n == 1 )
    VEC( out, 0 ) = VEC( in, 0 );
  else
  {
    complex double wn = cexp( -2*M_PI*I / n );
    vector *x, *y, *p, *q;
    create_vector( &x, n/2 );
    create_vector( &y, n/2 );
    create_vector( &p, n/2 );
    create_vector( &q, n/2 );
    for( k = 0; k < n/2; k++ )
    {
      VEC( x, k ) = VEC( in, 2*k );
      VEC( y, k ) = VEC( in, 2*k + 1 );
    }
    fft( x, p, n/2 );
    fft( y, q, n/2 );
    for( k = 0; k < n; k++ )
    {
      VEC( out, k ) = VEC( p, k % (n/2) ) + cpow(wn,k)*VEC( q, k % (n/2) );
    }
    destroy_vector( x );
    destroy_vector( y );
    destroy_vector( p );
    destroy_vector( q );
  }
}

static int _equal(CVECTOR *a, CVECTOR *b, bool invert)
{
	if (COMPLEX(a) || COMPLEX(b))
	{
		VECTOR_ensure_complex(a);
		VECTOR_ensure_complex(b);
		return gsl_vector_complex_equal(CVEC(a), CVEC(b));
	}
	else
		return gsl_vector_equal(VEC(a), VEC(b));
}

double
vector_dot(struct vector *x, struct vector *y)
{
    // can only dot vectors of equal length
    assert(x->n == y->n);

    double sum = 0.;
    for (int i = 0; i < x->n; i++) {
        sum += VEC(x, i) * VEC(y, i);
    }
    return sum;
}

static void
setup_test_vectors(struct vector *x, struct vector *y_ref)
{
  for (int i = 0; i < x->n; i++) {
    VEC(x, i) = i;
  }
  int len = y_ref->n;
  if (x->n < len) {
    len = x->n;
  }
  for (int i = 0; i < len; i++) {
    VEC(y_ref, i) = i * i;
  }
}

static CVECTOR *_sub(CVECTOR *a, CVECTOR *b, bool invert)
{
	CVECTOR *v = VECTOR_make(a);
	
	if (COMPLEX(v) || COMPLEX(b))
	{
		VECTOR_ensure_complex(v);
		VECTOR_ensure_complex(b);
		gsl_vector_complex_sub(CVEC(v), CVEC(b));
	}
	else
		gsl_vector_sub(VEC(v), VEC(b));
	
	return v;
}

VEC VEC::Vcross(VEC a)
  {
  double aa = v[1]*a.v[2] - v[2]*a.v[1];
  double b = v[2]*a.v[0] - v[0]*a.v[2];
  double c = v[0]*a.v[1] - v[1]*a.v[0];
  return VEC(aa,b,c);
  }

int main(int argc, char** argv) {
    int64 nside=4096;
    double rad_radians=-1;
    int rad_in_file=0;
    process_args(argc, argv, &nside, &rad_radians, &rad_in_file);

    struct healpix* hpix = hpix_new(nside);

    VEC(int64) pixlist = VEC_NEW(int64);

    double ra=0, dec=0;
    while (2 == fscanf(stdin,"%lf %lf", &ra, &dec)) {
        if (rad_in_file) {
            if (1 != fscanf(stdin,"%lf", &rad_radians)) {
                fprintf(stderr,"failed to read radius\n");
                exit(EXIT_FAILURE);
            }
            rad_radians *= D2R/3600.;
        }

        hpix_disc_intersect_radec(hpix, ra, dec, rad_radians, pixlist);

        printf("%.16g %.16g %lu", ra, dec, VEC_SIZE(pixlist);
        VEC_FOREACH(pix_ptr, pixlist) {
            printf(" %ld", *pix_ptr);
        }
        printf("\n");
    }

static char *_to_string(CVECTOR *_object, bool local)
{
	char *result = NULL;
	int i;
	int size = SIZE(THIS);
	char *str;
	int len;
	
	result = GB.AddChar(result, '[');
	
	for (i = 0; i < size; i++)
	{
		if (i)
			result = GB.AddChar(result, local ? ' ' : ',');
		
		if (!COMPLEX(THIS))
		{
			GB.NumberToString(local, gsl_vector_get(VEC(THIS), i), NULL, &str, &len);
			result = GB.AddString(result, str, len);
		}
		else
		{
			str = COMPLEX_to_string(gsl_vector_complex_get(CVEC(THIS), i), local);
			result = GB.AddString(result, str, GB.StringLength(str));
			GB.FreeString(&str);
		}
	}
	
	result = GB.AddChar(result, ']');
	
	return result;
}

/*
 * vec_fill_grid_mpi
 *
 * Takes a vector and a scalar and fills the vector
 * to be an evenly spaced grid from 0 ... scale.
 * Scaling with the rank from MPI.
 */
void vec_fill_grid_mpi( vector *v, val_type scale, int rank, int size )
{
  int i;
  for( i = 0; i < v->n; i++ )
    VEC( v, i ) = (double)((rank*size)+i) / (double)v->n;
    //VEC( v, i ) = ((double)i / (double)v->n) * (double)(rank+1);
}