Commit 14559f53 authored by Takhir Fakhrutdinov's avatar Takhir Fakhrutdinov

Добавили папки gravity,memcached

parent f9673ee9
......@@ -2,10 +2,13 @@
!.gitignore
!makefile
!readme.md
!configure.sh
!Makefile
!include/
!fson/
!sofa_c/
!sofa/
!rdjson/
!memcached/
!gravity/
!fte/
#!/bin/sh
CDIR=`pwd`
for i in `find zve/* -type dir -maxdepth 0` ; do
cd $CDIR/$i
../buildinfo
done
!*.c
!*.cpp
!*.h
!*.cl
!*.f
!PZ90_2
!PZ90_02_2.rtf
!EGM96
!EGM2008
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
# /System/Library/Frameworks/OpenCL.framework/Libraries/openclc -x cl -cl-std=CL1.1 -cl-auto-vectorize-enable -emit-gcl gravity.cl
PLATFORM=$(shell uname -s |awk -F- '{print $$1}')
BIN = ../build/$(PLATFORM)/bin
LIBS = ../build/$(PLATFORM)/lib
BUILD = ../build/$(PLATFORM)/obj/gravity
SHARE = ../build/$(PLATFORM)/share
INCLUDE = ../build/$(PLATFORM)/include
MODDIR = ../build/$(PLATFORM)/obj/mod
TARGET=gravity
TARGETCL=cl_gravity
TARGLIB=$(LIBS)/libgravity.a
TARGF=test
TARGETICL=openclinfo
F=.f
C=.c
SLV=solvers
OBJ=$(BUILD)/main.o $(BUILD)/wtime.o
OBJ1=$(BUILD)/cl_gravity.o $(BUILD)/wtime.o $(BUILD)/cl_hlp.o
OBJLIB=$(BUILD)/gravity.o
OBJF=$(BUILD)/test.o
LD=ld
GCC=gcc
FC=gfortran
AR=ar
GCCFLAGS += -g -Wall -c -O -I/usr/local/include -I$(INCLUDE)
GCCFLAGS += $(FCEXTRA)
FCFLAGS += -Wall -ffree-line-length-none -c -O -J$(MODDIR) -cpp
FCFLAGS += $(FCEXTRA)
ifeq ($(PLATFORM),$(filter $(PLATFORM),FreeBSD))
LDFLAGS += -lexecinfo
endif
ifeq ($(PLATFORM),$(filter $(PLATFORM),Darwin Linux))
LDFLAGS += -lgfortran -lstdc++ -lgravity -lm -L$(LIBS) -L/opt/local/lib -L/opt/local/lib/gcc48
else
LDFLAGS += -lgfortran -lstdc++ -lgravity -lm -L$(LIBS) -L/usr/local/lib -L/usr/local/lib/gcc48 -rpath /usr/local/lib/gcc48
endif
ifeq ($(PLATFORM),Darwin)
GCC=clang
GCCFLAGS += -framework OpenCL -std=c99
LDFLAGS += -framework OpenCL
else
ifeq ($(PLATFORM),Linux)
LDFLAGS += -lOpenCL
endif
endif
LDFLAGS += $(FCEXTRA)
RM=rm
all: dirs $(TARGLIB)
ifeq ($(PLATFORM),FreeBSD)
test: $(TARGET) $(TARGF)
else
test: $(TARGET) $(TARGETCL) $(TARGF) $(TARGETICL)
endif
dirs:
mkdir -p $(BIN) $(LIBS) $(BUILD) $(SHARE) $(MODDIR) $(INCLUDE)
cp EGM2008 $(SHARE)
cp EGM96 $(SHARE)
cp PZ90_2 $(SHARE)
cp gravity.h $(INCLUDE)
cp ../memcached/memc_gate.h $(INCLUDE)
$(TARGLIB) : $(OBJLIB)
$(AR) rvs $@ $(OBJLIB)
$(TARGET) : $(OBJ) $(OBJLIB)
$(GCC) -O -o $@ $(OBJ) $(LDFLAGS) -lmemc_gate -lmemcached -lsofa_c -lsofa
$(TARGETCL) : $(OBJ1)
$(GCC) -O -o $@ $(OBJ1) $(LDFLAGS)
$(TARGETICL) : $(TARGETICL).o
$(GCC) -O -o $@ $< $(LDFLAGS)
$(TARGF) : $(OBJF)
$(FC) -O -o $@ $(OBJF) $(LDFLAGS) -lmemc_gate -lmemcached -lsofa_c -lsofa -lsolv
$(BUILD)/%.o : %.c
$(GCC) $(GCCFLAGS) -c $< -o $@
$(BUILD)/%.o : %.cpp
$(GCC) $(GCCFLAGS) -c $< -o $@
$(BUILD)/%.o : %.f
$(FC) $(FCFLAGS) -c $< -o $@
clean:
rm -rf $(BUILD) $(TARGET) $(TARGETICL) $(TARGETCL) $(TARGLIB) $(TARGF) $(SHARE)/EGM2008 $(SHARE)/EGM96 $(SHARE)/PZ90_2
#all: memc_gate.c
# gcc -Wall -I/usr/local/include -L/usr/local/lib -L/usr/home/viking/develop/solvers/build/FreeBSD/lib -lmemcached -lsofa_c -lm -lexecinfo -lstdc++ -lc -o memc_gate memc_gate.c
This diff is collapsed.
This diff is collapsed.
#include <stdio.h>
#include <stdlib.h>
#ifdef __APPLE__
#include <OpenCL/cl.h>
#else
#include <CL/cl.h>
#endif
#define VECTOR_SIZE 4
#define VECTORS 180
#define MAX(a,b) ((a>=b)?a:b)
#define MIN(a,b) ((a<=b)?a:b)
#define MAINDEV CL_DEVICE_TYPE_GPU
extern double wtime();
extern char *getKernelSource(char *);
extern int output_device_info(cl_device_id);
extern void error(char*, cl_int);
extern char* err_code (cl_int);
double fRand(double fMin, double fMax)
{
double f = (double)rand() / RAND_MAX;
return fMin + f * (fMax - fMin);
}
int main(void) {
int i;
const char *k_harmo = getKernelSource("gravity.cl");
// Allocate space for vectors A, B and C
double *A = (double*)malloc(sizeof(double)*VECTOR_SIZE*VECTORS);
double *B = (double*)malloc(sizeof(double)*VECTOR_SIZE*VECTORS);
for (i=0; i<VECTORS; i++) {
A[i*VECTOR_SIZE+0] = fRand(-35000,35000);
A[i*VECTOR_SIZE+1] = fRand(-35000,35000);
A[i*VECTOR_SIZE+2] = fRand(-35000,35000);
A[i*VECTOR_SIZE+3] = 0;
}
/*
A[0] = 30544.803906240781;
A[1] = 28708.707913342118;
A[2] = 4777.8592942243504;
*/
// for(i = 0; i < VECTOR_SIZE; i++) B[i] = 0;
// Get platform and device information
cl_platform_id * platforms = NULL;
cl_uint num_platforms;
cl_int clStatus = clGetPlatformIDs(0, NULL, &num_platforms);
if (clStatus != CL_SUCCESS || num_platforms <= 0) error("Failed to find a platforms.",clStatus);
platforms = (cl_platform_id *) malloc(sizeof(cl_platform_id)*num_platforms);
clStatus = clGetPlatformIDs(num_platforms, platforms, NULL);
if (clStatus != CL_SUCCESS || num_platforms <= 0) error("Failed to get the platforms.",clStatus);
//Get the devices list and choose the device you want to run on
cl_device_id *device_list = NULL;
cl_uint num_devices;
clStatus = clGetDeviceIDs( platforms[0], MAINDEV, 0, NULL, &num_devices);
device_list = (cl_device_id *) malloc(sizeof(cl_device_id)*num_devices);
clStatus = clGetDeviceIDs( platforms[0], MAINDEV, num_devices, device_list, NULL);
// for (i=0; i<num_devices; i++)
// output_device_info(device_list[i]);
cl_context context = clCreateContext( NULL, num_devices, device_list, NULL, NULL, &clStatus);
if (!context) error("Failed to create a compute context.", clStatus);
cl_command_queue command_queue = clCreateCommandQueue( context, device_list[0], 0, &clStatus);
if (!command_queue) error("Failed to create a command queue.", clStatus);
// Create a program from the kernel source
cl_program program = clCreateProgramWithSource(context, 1, (const char **)&k_harmo, 0, &clStatus);
if (clStatus != CL_SUCCESS) error("Error: Failed to create compute program.", clStatus);
// Build the program
// clStatus = clBuildProgram(program, 1, device_list, NULL, NULL, NULL);
clStatus = clBuildProgram(program, 1, device_list, "-DHARMONIX=64", NULL, NULL);
if (clStatus != CL_SUCCESS) {
size_t len;
char buffer[2048];
clGetProgramBuildInfo(program, NULL, CL_PROGRAM_BUILD_LOG, sizeof(buffer), buffer, &len);
printf("%s\n", buffer);
error("Error: Failed to build program executable.", clStatus);
}
// Create the OpenCL kernel
cl_kernel kernel = clCreateKernel(program, "harmo", &clStatus);
if (clStatus != CL_SUCCESS) error("Failed to translate kernel.",clStatus);
// Create memory buffers on the device for each vector
cl_mem A_clmem = clCreateBuffer(context, CL_MEM_READ_ONLY, VECTOR_SIZE * VECTORS * sizeof(double), NULL, &clStatus);
if (clStatus != CL_SUCCESS) error("Failed to create clBuffer.",clStatus);
cl_mem B_clmem = clCreateBuffer(context, CL_MEM_WRITE_ONLY, VECTOR_SIZE * VECTORS * sizeof(double), NULL, &clStatus);
if (clStatus != CL_SUCCESS) error("Failed to create clBuffer.",clStatus);
double rtime = wtime();
// Copy the Buffer A to the device
clStatus = clEnqueueWriteBuffer(command_queue, A_clmem, CL_TRUE, 0, VECTOR_SIZE * VECTORS * sizeof(double), A, 0, NULL, NULL);
if (clStatus != CL_SUCCESS) error("clEnqueueWriteBuffer failed.", clStatus);
// Set the arguments of the kernel
clStatus = clSetKernelArg(kernel, 0, sizeof(cl_mem), (void *)&A_clmem);
clStatus = clSetKernelArg(kernel, 1, sizeof(cl_mem), (void *)&B_clmem);
// Execute the OpenCL kernel on the list
size_t global_size = VECTORS; // Process the entire lists
size_t local_size = MAX(VECTORS,2048); // Process items at a time
clStatus = clEnqueueNDRangeKernel(command_queue, kernel, 1, NULL, &global_size, &local_size, 0, NULL, NULL);
// Read the cl memory C_clmem on device to the host variable C
clStatus = clEnqueueReadBuffer(command_queue, B_clmem, CL_TRUE, 0, VECTOR_SIZE * VECTORS * sizeof(double), B, 0, NULL, NULL);
// Clean up and wait for all the comands to complete.
clStatus = clFlush(command_queue);
clStatus = clFinish(command_queue);
rtime = wtime() - rtime;
// Display the result to the screen
for(i = 0; i < VECTORS; i++) printf("%+10.3lf %+10.3lf %+10.3lf %+.16e %+.16e %+.16e\n", A[i*VECTOR_SIZE],A[i*VECTOR_SIZE+1],A[i*VECTOR_SIZE+2],B[i*VECTOR_SIZE],B[i*VECTOR_SIZE+1],B[i*VECTOR_SIZE+2]);
printf("%lf\n", rtime);
// Finally release all OpenCL allocated objects and host buffers.
clStatus = clReleaseKernel(kernel);
clStatus = clReleaseProgram(program);
clStatus = clReleaseMemObject(A_clmem);
clStatus = clReleaseMemObject(B_clmem);
clStatus = clReleaseCommandQueue(command_queue);
clStatus = clReleaseContext(context);
free(A);
free(B);
free(platforms);
free(device_list);
return 0;
}
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define vImage_Utilities_h
#define vImage_CVUtilities_h
#ifdef __APPLE__
#include <OpenCL/opencl.h>
#else
#include <CL/cl.h>
#endif
//
// define VERBOSE if you want to print info about work groups sizes
#define VERBOSE 1
//------------------------------------------------------------------------------
//
// Name: err_code()
//
// Purpose: Function to output descriptions of errors for an input error code
//
//
// RETURN: echoes the input error code
//
// HISTORY: Written by Tim Mattson, June 2010
// This version automatically produced by genErrCode.py
// script written by Tom Deakin, August 2013
//
//------------------------------------------------------------------------------
char *err_code (cl_int err_in)
{
switch (err_in) {
case CL_SUCCESS :
return (char*)" CL_SUCCESS ";
case CL_DEVICE_NOT_FOUND :
return (char*)" CL_DEVICE_NOT_FOUND ";
case CL_DEVICE_NOT_AVAILABLE :
return (char*)" CL_DEVICE_NOT_AVAILABLE ";
case CL_COMPILER_NOT_AVAILABLE :
return (char*)" CL_COMPILER_NOT_AVAILABLE ";
case CL_MEM_OBJECT_ALLOCATION_FAILURE :
return (char*)" CL_MEM_OBJECT_ALLOCATION_FAILURE ";
case CL_OUT_OF_RESOURCES :
return (char*)" CL_OUT_OF_RESOURCES ";
case CL_OUT_OF_HOST_MEMORY :
return (char*)" CL_OUT_OF_HOST_MEMORY ";
case CL_PROFILING_INFO_NOT_AVAILABLE :
return (char*)" CL_PROFILING_INFO_NOT_AVAILABLE ";
case CL_MEM_COPY_OVERLAP :
return (char*)" CL_MEM_COPY_OVERLAP ";
case CL_IMAGE_FORMAT_MISMATCH :
return (char*)" CL_IMAGE_FORMAT_MISMATCH ";
case CL_IMAGE_FORMAT_NOT_SUPPORTED :
return (char*)" CL_IMAGE_FORMAT_NOT_SUPPORTED ";
case CL_BUILD_PROGRAM_FAILURE :
return (char*)" CL_BUILD_PROGRAM_FAILURE ";
case CL_MAP_FAILURE :
return (char*)" CL_MAP_FAILURE ";
case CL_MISALIGNED_SUB_BUFFER_OFFSET :
return (char*)" CL_MISALIGNED_SUB_BUFFER_OFFSET ";
case CL_EXEC_STATUS_ERROR_FOR_EVENTS_IN_WAIT_LIST :
return (char*)" CL_EXEC_STATUS_ERROR_FOR_EVENTS_IN_WAIT_LIST ";
case CL_INVALID_VALUE :
return (char*)" CL_INVALID_VALUE ";
case CL_INVALID_DEVICE_TYPE :
return (char*)" CL_INVALID_DEVICE_TYPE ";
case CL_INVALID_PLATFORM :
return (char*)" CL_INVALID_PLATFORM ";
case CL_INVALID_DEVICE :
return (char*)" CL_INVALID_DEVICE ";
case CL_INVALID_CONTEXT :
return (char*)" CL_INVALID_CONTEXT ";
case CL_INVALID_QUEUE_PROPERTIES :
return (char*)" CL_INVALID_QUEUE_PROPERTIES ";
case CL_INVALID_COMMAND_QUEUE :
return (char*)" CL_INVALID_COMMAND_QUEUE ";
case CL_INVALID_HOST_PTR :
return (char*)" CL_INVALID_HOST_PTR ";
case CL_INVALID_MEM_OBJECT :
return (char*)" CL_INVALID_MEM_OBJECT ";
case CL_INVALID_IMAGE_FORMAT_DESCRIPTOR :
return (char*)" CL_INVALID_IMAGE_FORMAT_DESCRIPTOR ";
case CL_INVALID_IMAGE_SIZE :
return (char*)" CL_INVALID_IMAGE_SIZE ";
case CL_INVALID_SAMPLER :
return (char*)" CL_INVALID_SAMPLER ";
case CL_INVALID_BINARY :
return (char*)" CL_INVALID_BINARY ";
case CL_INVALID_BUILD_OPTIONS :
return (char*)" CL_INVALID_BUILD_OPTIONS ";
case CL_INVALID_PROGRAM :
return (char*)" CL_INVALID_PROGRAM ";
case CL_INVALID_PROGRAM_EXECUTABLE :
return (char*)" CL_INVALID_PROGRAM_EXECUTABLE ";
case CL_INVALID_KERNEL_NAME :
return (char*)" CL_INVALID_KERNEL_NAME ";
case CL_INVALID_KERNEL_DEFINITION :
return (char*)" CL_INVALID_KERNEL_DEFINITION ";
case CL_INVALID_KERNEL :
return (char*)" CL_INVALID_KERNEL ";
case CL_INVALID_ARG_INDEX :
return (char*)" CL_INVALID_ARG_INDEX ";
case CL_INVALID_ARG_VALUE :
return (char*)" CL_INVALID_ARG_VALUE ";
case CL_INVALID_ARG_SIZE :
return (char*)" CL_INVALID_ARG_SIZE ";
case CL_INVALID_KERNEL_ARGS :
return (char*)" CL_INVALID_KERNEL_ARGS ";
case CL_INVALID_WORK_DIMENSION :
return (char*)" CL_INVALID_WORK_DIMENSION ";
case CL_INVALID_WORK_GROUP_SIZE :
return (char*)" CL_INVALID_WORK_GROUP_SIZE ";
case CL_INVALID_WORK_ITEM_SIZE :
return (char*)" CL_INVALID_WORK_ITEM_SIZE ";
case CL_INVALID_GLOBAL_OFFSET :
return (char*)" CL_INVALID_GLOBAL_OFFSET ";
case CL_INVALID_EVENT_WAIT_LIST :
return (char*)" CL_INVALID_EVENT_WAIT_LIST ";
case CL_INVALID_EVENT :
return (char*)" CL_INVALID_EVENT ";
case CL_INVALID_OPERATION :
return (char*)" CL_INVALID_OPERATION ";
case CL_INVALID_GL_OBJECT :
return (char*)" CL_INVALID_GL_OBJECT ";
case CL_INVALID_BUFFER_SIZE :
return (char*)" CL_INVALID_BUFFER_SIZE ";
case CL_INVALID_MIP_LEVEL :
return (char*)" CL_INVALID_MIP_LEVEL ";
case CL_INVALID_GLOBAL_WORK_SIZE :
return (char*)" CL_INVALID_GLOBAL_WORK_SIZE ";
case CL_INVALID_PROPERTY :
return (char*)" CL_INVALID_PROPERTY ";
default:
return (char*)"UNKNOWN ERROR";
}
}
void error(char *str, cl_int clStatus){
fprintf(stderr,"ERROR: %s\n%s\n",str,err_code(clStatus));
exit(EXIT_FAILURE);
}
int output_device_info(cl_device_id device_id)
{
int err,i; // error code returned from OpenCL calls
cl_device_type device_type; // Parameter defining the type of the compute device
cl_uint comp_units; // the max number of compute units on a device
cl_char vendor_name[1024] = {0}; // string to hold vendor name for compute device
cl_char device_name[1024] = {0}; // string to hold name of compute device
#ifdef VERBOSE
cl_uint max_work_itm_dims;
size_t max_wrkgrp_size;
size_t *max_loc_size;
#endif
err = clGetDeviceInfo(device_id, CL_DEVICE_NAME, sizeof(device_name), &device_name, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to access device name!\n");
return EXIT_FAILURE;
}
printf(" \n Device is %s ",device_name);
err = clGetDeviceInfo(device_id, CL_DEVICE_TYPE, sizeof(device_type), &device_type, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to access device type information!\n");
return EXIT_FAILURE;
}
if(device_type == CL_DEVICE_TYPE_GPU)
printf(" GPU from ");
else if (device_type == CL_DEVICE_TYPE_CPU)
printf("\n CPU from ");
else
printf("\n non CPU or GPU processor from ");
err = clGetDeviceInfo(device_id, CL_DEVICE_VENDOR, sizeof(vendor_name), &vendor_name, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to access device vendor name!\n");
return EXIT_FAILURE;
}
printf(" %s ",vendor_name);
err = clGetDeviceInfo(device_id, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(cl_uint), &comp_units, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to access device number of compute units !\n");
return EXIT_FAILURE;
}
printf(" with a max of %d compute units \n",comp_units);
#ifdef VERBOSE
//
// Optionally print information about work group sizes
//
err = clGetDeviceInfo( device_id, CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS, sizeof(cl_uint),
&max_work_itm_dims, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to get device Info (CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS)\n%s\n", err_code(err));
return EXIT_FAILURE;
}
max_loc_size = (size_t*)malloc(max_work_itm_dims * sizeof(size_t));
if(max_loc_size == NULL){
printf(" malloc failed\n");
return EXIT_FAILURE;
}
err = clGetDeviceInfo( device_id, CL_DEVICE_MAX_WORK_ITEM_SIZES, max_work_itm_dims* sizeof(size_t),
max_loc_size, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to get device Info (CL_DEVICE_MAX_WORK_ITEM_SIZES)\n%s\n",err_code(err));
return EXIT_FAILURE;
}
err = clGetDeviceInfo( device_id, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(size_t),
&max_wrkgrp_size, NULL);
if (err != CL_SUCCESS)
{
printf("Error: Failed to get device Info (CL_DEVICE_MAX_WORK_GROUP_SIZE)\n%s\n",err_code(err));
return EXIT_FAILURE;
}
printf("work group, work item information");
printf("\n max loc dim ");
for(i=0; i< max_work_itm_dims; i++)
printf(" %d ",(int)(*(max_loc_size+i)));
printf("\n");
printf(" Max work group size = %d\n",(int)max_wrkgrp_size);
#endif
return CL_SUCCESS;
}
//------------------------------------------------------------------------------
char * getKernelSource(char *filename)
{
FILE *file = fopen(filename, "r");
if (!file)
{
fprintf(stderr, "Error: Could not open kernel source file\n");
exit(EXIT_FAILURE);
}
fseek(file, 0, SEEK_END);
int len = ftell(file) + 1;
rewind(file);
char *source = (char *)calloc(sizeof(char), len);
if (!source)
{
fprintf(stderr, "Error: Could not allocate memory for source string\n");
exit(EXIT_FAILURE);
}
fread(source, sizeof(char), len, file);
fclose(file);
return source;
}
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
#define MIN(X,Y) ((X) < (Y) ? (X) : (Y))
//#define EGM2008_GM 3986004.418e+08
//#define EGM96_GM 3986004.415e+08
//#define PZ90_2_GM 3986004.418e+08
#define EGM2008_GM 398600.4418
#define EGM96_GM 398600.4418
#define PZ90_2_GM 398600.4418
#define EGM2008_A 6378.1370
#define EGM96_A 6378.1370
#define PZ90_2_A 6378.1360
#define EGM2008_F (1./298.257223563)
#define EGM96_F (1./298.257223563)
#define PZ90_2_F (1./298.25784)
#define TRUE 1
#define FALSE !TRUE
#define MAXHARMONIX 8
enum {
jplMercury,
jplVenus,
jplEarth,
jplMars,
jplJupiter,
jplSaturn,
jplUranus,
jplNeptune,
jplPluto,
jplMoon,
jplSun,
jplSS_Bary,
jplEM_Bary,
jplNutations,
jplLibrations,
jplEmbary
};
enum {
PF_EGM2008,
PF_EGM96,
PF_PZ90_2
};
#define MAIN_PF_MODEL PF_PZ90_2
typedef struct {
int planet;
int count;
double GM; // геоцентрическая гравитационная постоянная, [м^3/с^2]
double A; // экваториальный радиус, м
double GM_div_A; // ;)
double f; // сжатие эллипсоида
double *C;
double *S;
double nC20;
int mem;
char *fname;
}PlanetField;
int initplanetfield(int PF,PlanetField *field);
extern void planetfieldacc(int N,PlanetField *f,double *x,double *grad,double *U);
extern void c_harmo(int N, double *x, double *grad);
extern void c_harmo_(int N, double *x, double *grad);
extern void c_harmo8(double *x, double *grad);
extern void legendre(int N, double sin_fi, double cos_fi, double *P);
extern void legendre8(double sin_fi, double cos_fi, double *P);
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "gravity.h"
#define VECTOR_SIZE 4
#define VECTORS 180
extern double wtime();
double fRand(double fMin, double fMax)
{
double f = (double)rand() / RAND_MAX;
return fMin + f * (fMax - fMin);
}
int main() {
int i,j;
// PlanetField f;
// double u;
double *A = (double*)malloc(sizeof(double)*VECTOR_SIZE*VECTORS);
double *B = (double*)malloc(sizeof(double)*VECTOR_SIZE*VECTORS);
/*
A[0] = 30544.803906240781;
A[1] = 28708.707913342118;
A[2] = 4777.8592942243504;
*
double a[3];
*/
for (i=0; i<VECTORS; i++) {
A[i*VECTOR_SIZE+0] = fRand(-35000,35000);
A[i*VECTOR_SIZE+1] = fRand(-35000,35000);
A[i*VECTOR_SIZE+2] = fRand(-35000,35000);
A[i*VECTOR_SIZE+3] = 0;
}
double rtime = wtime();
/* for (i=0;i<3;i++)
if ( InitPlanetField(i,&f) ) {
printf ("%s %d %e\n",f.fname,f.count,f.nC20);
PlanetFieldAcc(128,&f,v,a,&u);
printf("%e %e %e\n",a[0],a[1],a[2]);
}
*/
for (i=0; i<VECTORS; i++) {
c_harmo(8,&A[i*VECTOR_SIZE],&B[i*VECTOR_SIZE]);
}
rtime = wtime() - rtime;
for(i = 0; i < VECTORS; i++) printf("%+10.3lf %+10.3lf %+10.3lf %+.16e %+.16e %+.16e\n", A[i*VECTOR_SIZE],A[i*VECTOR_SIZE+1],A[i*VECTOR_SIZE+2],B[i*VECTOR_SIZE],B[i*VECTOR_SIZE+1],B[i*VECTOR_SIZE+2]);
printf("%lf\n", rtime);
// }
// double *P = (double *)malloc( sizeof(double[81]) );
/* for ( i=0;i<10000000;i++)
c_harmo_(8,v,a);
printf("%.15e %.15e %.15e\n",a[0],a[1],a[2]);
legendre(8,0.2,sqrt(1-0.2*0.2),P);
for (i=0; i<8; i++) {
for (j=0; j<8; j++)
printf("%+.2e ",P[i*8+j]);
printf("\n");
}
printf("\n");
for ( i=0;i<10000000;i++)
c_harmo8(v,a);
printf("%.15e %.15e %.15e\n",a[0],a[1],a[2]);
legendre8(0.2,sqrt(1-0.2*0.2),P);
for (i=0; i<8; i++) {
for (j=0; j<8; j++)
printf("%+.2e ",P[i*8+j]);
printf("\n");
}
*/
// free (P);
/* printf("\n");
legendre8(0.2,sqrt(1-0.2*0.2),P);
for (i=0; i<8; i++) {
for (j=0; j<8; j++)
printf("%+.5e",P[i*8+j]);
printf("\n");
}
*/
/* v[0]=30554.803906240781;
v[1]=28718.707913342118;
v[2]=4787.8592942243504;
c_harmo(8,v,a);
printf("%e %e %e\n",a[0],a[1],a[2]);
*/
return 0;
}
This diff is collapsed.
module test1
implicit none
contains
subroutine HrmEgm (N,x,f)
implicit none
integer N
double precision x(*),f(*)
interface
subroutine c_harmo (N,x,f) bind(c)
use,intrinsic :: iso_c_binding
integer(c_int),value,intent(in) :: N
real(c_double),intent(in) :: x(*)
real(c_double),intent(out) :: f(*)
end subroutine c_harmo
end interface
call c_harmo(N,x,f)
end subroutine
end module test1
program test
use test1
double precision x(3),f(3)
data x(1) / 30544.803906240781D00 /
data x(2) / 28708.707913342118D00 /
data x(3) / 4777.8592942243504D00 /
call HrmEgm(8,x,f)
write (*,*) f(1),f(2),f(3)
call HrmEgm96(8,x,f)
write (*,*) f(1),f(2),f(3)
end program
\ No newline at end of file
#ifdef _OPENMP
#include <omp.h>
#else
#include <sys/time.h>
#endif
#include <stdlib.h>
double wtime()
{
#ifdef _OPENMP
/* Use omp_get_wtime() if we can */
return omp_get_wtime();
#else
/* Use a generic timer */
static int sec = -1;
struct timeval tv;
gettimeofday(&tv, NULL);
if (sec < 0) sec = tv.tv_sec;
return (tv.tv_sec - sec) + 1.0e-6*tv.tv_usec;
#endif
}
memc_gate
!*.c
!*.h
!*.new
!*.f
!.gitingnore
PLATFORM=$(shell uname -s |awk -F- '{print $$1}')
BIN = ../build/$(PLATFORM)/bin
LIBS = ../build/$(PLATFORM)/lib
BUILD = ../build/$(PLATFORM)/obj/memc_gate
SHARE = ../build/$(PLATFORM)/share
INCLUDE = ../build/$(PLATFORM)/include
MODDIR = ../build/$(PLATFORM)/obj/mod
TARGET=memc_gate
TARGLIB=$(LIBS)/libmemc_gate.a
TARGF=test
F=.f
C=.c
SLV=solvers
OBJ=$(BUILD)/memc_main.o
OBJLIB= $(BUILD)/memc_gate.o
OBJF=$(BUILD)/test.o
GCC=gcc
FC=gfortran
AR=ar
GCCFLAGS += -g -Wall -c -O -I/opt/local/include -I/usr/local/include -I$(INCLUDE)
GCCFLAGS += $(FCEXTRA)
FCFLAGS += -Wall -ffree-line-length-none -c -O -J$(MODDIR) -cpp
FCFLAGS += $(FCEXTRA)
LDFLAGS = -lgfortran -lmemc_gate -lgravity -lstdc++ -lmemcached -lsofa_c -lm -L$(LIBS) -L/usr/local/lib -L/opt/local/lib -L/opt/local/lib/gcc48 -L/usr/local/lib/gcc48
ifeq ($(PLATFORM),FreeBSD)
LDFLAGS += -lexecinfo -rpath /usr/local/lib/gcc48
endif
LDFLAGS += $(FCEXTRA)
RM=rm
all: dirs $(TARGLIB) $(TARGET) $(TARGF)
dirs:
mkdir -p $(BIN) $(LIBS) $(BUILD)
cp memc_gate.h $(INCLUDE)
$(TARGLIB) : $(OBJLIB)
$(AR) rvs $@ $(OBJLIB)
$(TARGET) : $(OBJ) $(OBJLIB)
$(GCC) -O -o $@ $(OBJ) $(LDFLAGS)
$(TARGF) : $(OBJF)
$(FC) -O -o $@ $(OBJF) $(LDFLAGS)
$(BUILD)/%.o : %.c
$(GCC) $(GCCFLAGS) -c $< -o $@
$(BUILD)/%.o : %.f
$(FC) $(FCFLAGS) -c $< -o $@
clean:
rm -rf $(BUILD) $(TARGET) $(TARGLIB) $(TARGF)
#all: memc_gate.c
# gcc -Wall -I/usr/local/include -L/usr/local/lib -L/usr/home/viking/develop/solvers/build/FreeBSD/lib -lmemcached -lsofa_c -lm -lexecinfo -lstdc++ -lc -o memc_gate memc_gate.c
double precision function iau_ee00a (dt1,dt2)
implicit none
double precision dt1,dt2
interface
real(c_double) function memc_iau_ee00a (dt1,dt2) bind(c)
use,intrinsic :: iso_c_binding
real(c_double),value,intent(in) :: dt1
real(c_double),value,intent(in) :: dt2
end function memc_iau_ee00a
end interface
iau_ee00a = memc_iau_ee00a(dt1,dt2)
end function
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <libmemcached/memcached.h>
#include <sofa.h>
#include "memc_gate.h"
#define SZOC sizeof(char)
#define SZOI sizeof(int)
#define SZOD sizeof(double)
#define TRUE 1
#define FALSE !TRUE
// #define STRKEYS // long string keys
// #define DEBUG // as is
memcached_st *memc;
char first_run = TRUE;
void close_memc_connect() {
memcached_free(memc);
}
extern void legendre(int N, double sin_fi, double cos_fi, double *P);
#ifdef DEBUG
void debug_keys(int keylen, char* key) {
#ifndef STRKEYS
int j;
#endif
fprintf(stderr,"KEYLEN:%d\n",(int)keylen);
#ifdef STRKEYS
fprintf(stderr,"KEY:%s\n",key);
#else
fprintf(stderr,"HEXKEY:");
for (j=0;j<keylen;j++)
fprintf(stderr,"%02x ",(unsigned char)key[j]);
fprintf(stderr,"\n");
#endif
}
#endif
void memc_iau_nut00a(double dt1, double dt2, double *dpsi, double *deps ) {
char fn = '1'; // unique function descriptor
#ifdef STRKEYS
size_t keylen = 0;
char key[MEMCACHED_MAX_KEY] = "";
#else
size_t keylen = SZOC+SZOD*2;
char key[SZOC+SZOD*2] = "";
#endif
uint32_t flags = 0;
size_t length = 0;
char *retval = NULL;
memcached_return_t rc;
char val[SZOD*2]; // 16 bytes for values
#ifdef STRKEYS
sprintf(key,"%c%+-.9f%+-.9f",fn,dt1,dt2);
keylen = strlen(key);
#else
memcpy(key , &fn, SZOC);
memcpy(key + SZOC , &dt1, SZOD);
memcpy(key + SZOC + SZOD, &dt2, SZOD);
#endif
#ifdef DEBUG
debug_keys(keylen,key);
#endif
if (first_run) {
const char *config_string= "--SOCKET=\"/tmp/memcached.sock\" --BINARY-PROTOCOL";
// add --NAMESPACE=... later ??? for simple function results isolation
memc = memcached(config_string, strlen(config_string));
first_run = ! first_run;
}
retval = memcached_get (memc, key, keylen, &length, &flags, &rc);
if (rc != MEMCACHED_SUCCESS) {
#ifdef DEBUG
fprintf(stderr,"%d: %s\n", rc, memcached_strerror(memc,rc));
#endif
iauNut00a(dt1,dt2,dpsi,deps);
memcpy(val , dpsi, SZOD);
memcpy(val + SZOD, deps, SZOD);
rc = memcached_set(memc, key, keylen, val, SZOD*2, (time_t)0, (uint32_t)0);
#ifdef DEBUG
if (rc != MEMCACHED_SUCCESS)
fprintf(stderr,"%d: %s\n", rc, memcached_strerror(memc,rc));
#endif
} else {
memcpy(dpsi,retval , SZOD);
memcpy(deps,retval + SZOD, SZOD);
}
free(retval);
}
double memc_iau_ee00a ( double dt1, double dt2 ) {
char fn = '2'; // unique function descriptor
#ifdef STRKEYS
size_t keylen = 0;
char key[MEMCACHED_MAX_KEY] = "";
#else
size_t keylen = SZOC+SZOD*2;
char key[SZOC+SZOD*2] = "";
#endif
uint32_t flags = 0;
size_t length = 0;
char *retval = NULL;
memcached_return_t rc;
double ee00a;
char val[SZOD]; // 8 bytes for values
#ifdef STRKEYS
sprintf(key,"%c%+-.9f%+-.9f",fn,dt1,dt2);
keylen = strlen(key);
#else
memcpy(key , &fn, SZOC);
memcpy(key + SZOC , &dt1, SZOD);
memcpy(key + SZOC + SZOD, &dt2, SZOD);
#endif
#ifdef DEBUG
debug_keys(keylen,key);
#endif
if (first_run) {
const char *config_string= "--SOCKET=\"/tmp/memcached.sock\" --BINARY-PROTOCOL";
memc = memcached(config_string, strlen(config_string));
first_run = ! first_run;
}
retval = memcached_get (memc, key, keylen, &length, &flags, &rc);
if (rc != MEMCACHED_SUCCESS) {
#ifdef DEBUG
fprintf(stderr,"%d: %s\n", rc, memcached_strerror(memc,rc));
#endif
ee00a = iauEe00a(dt1,dt2);
memcpy(val , &ee00a, SZOD);
rc = memcached_set(memc, key, keylen, val, SZOD, (time_t)0, (uint32_t)0);
#ifdef DEBUG
if (rc != MEMCACHED_SUCCESS)
fprintf(stderr,"%d: %s\n", rc, memcached_strerror(memc,rc));
#endif
} else {
memcpy(&ee00a,retval , SZOD);
}
free(retval);
return ee00a;
}
void memc_legendre(int N, double sin_fi, double cos_fi , double *P){
char fn = '3'; // unique function descriptor
#ifdef STRKEYS
size_t keylen = 0;
char key[MEMCACHED_MAX_KEY] = "";
#else
size_t keylen = SZOC+SZOI+SZOD;
char key[SZOC+SZOI+SZOD] = "";
#endif
uint32_t flags = 0;
size_t length = 0;
char *retval = NULL;
memcached_return_t rc;
double r_sin_fi, r_cos_fi;
int n, m, i;
char *val;
i = sizeof(double[N+1][N+1]);
val = (char *)malloc( i );
if ( ! ( ( sin_fi < 0.35 ) && ( sin_fi > -0.35 ) ) ) {
legendre( N, sin_fi, cos_fi, P);
} else {
#define PREC 10000000
r_sin_fi = round( ( sin_fi ) * PREC ) / PREC;
r_cos_fi = round( ( cos_fi ) * PREC ) / PREC;
#ifdef STRKEYS
sprintf(key,"%c%d%+-.8f",fn,N,r_sin_fi);
keylen = strlen(key);
#else
memcpy(key , &fn, SZOC);
memcpy(key + SZOC , &N, SZOI);
memcpy(key + SZOC + SZOI, &r_sin_fi, SZOD);
#endif
#ifdef DEBUG
debug_keys(keylen,key);
#endif
if (first_run) {
const char *config_string= "--SOCKET=\"/tmp/memcached.sock\" --BINARY-PROTOCOL";
memc = memcached(config_string, strlen(config_string));
first_run = ! first_run;
}
retval = memcached_get (memc, key, keylen, &length, &flags, &rc);
#ifdef DEBUG
if (rc != MEMCACHED_SUCCESS)
fprintf(stderr,"%d: %s\n", rc, memcached_strerror(memc,rc));
#endif
if ( rc == 16 ) { // NOT FOUND
legendre( N, r_sin_fi, r_cos_fi, P);
i = 0;
for( m = 0; m <= N; m++ )
for( n = m; n <= N; n++ ) {
memcpy(val+i*SZOD, &P[n*N+m], SZOD);
i++;
}
rc = memcached_set(memc, key, keylen, val, SZOD*i, (time_t)0, (uint32_t)0);
#ifdef DEBUG
if (rc != MEMCACHED_SUCCESS)
fprintf(stderr,"%d: %s\n", rc, memcached_strerror(memc,rc));
#endif
} else {
if ( rc == 0 ) { // FOUND
i=0;
for( m = 0; m <= N; m++ )
for( n = m; n <= N; n++ ) {
memcpy(&P[n*N+m],retval+i*SZOD, SZOD);
i++;
}
} else // NO CONNECTION & OTHER ERRORS
legendre( N, sin_fi, cos_fi, P);
}
free(retval);
}
free(val);
}
void memc_iau_nut00a (double dt1, double dt2, double *dpsi, double *deps);
double memc_iau_ee00a (double dt1, double dt2);
void memc_legendre (int N, double sin_fi, double cos_fi , double *P);
void close_memc_connect ();
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <libmemcached/memcached.h>
#include <sofa.h>
#include "memc_gate.h"
int main( void )
{
int i,n,m;
double dt1 = 2451547.1,dt2 = -1423.1;
double dpsi,deps;
double *P;
int N=8;
P = (double *)malloc( sizeof(double)*(N+1)*(N+1) );
memset( (void *)P, 0, sizeof(double)*(N+1)*(N+1) );
for(i=0;i<1;i++) {
iauNut00a(dt1,dt2,&dpsi,&deps);
fprintf(stderr,"%e\n",dpsi);
fprintf(stderr,"%e\n",deps);
memc_iau_nut00a(dt1,dt2,&dpsi,&deps);
fprintf(stderr,"%e\n",dpsi);
fprintf(stderr,"%e\n",deps);
fprintf(stderr,"%e\n",iauEe00a(dt1,dt2));
fprintf(stderr,"%e\n",memc_iau_ee00a(dt1,dt2));
memc_legendre(N, 0.15,0.98868599666426,P);
for( m = 0; m <= N; m++ ) {
for( n = m; n <= N; n++ ) {
fprintf(stderr,"%.15e ",P[n*N+m]);
}
fprintf(stderr,"\n");
}
}
free(P);
return 0;
}
subroutine iau_nut00a (dt1,dt2,dpsi,deps)
implicit none
double precision dt1,dt2,dpsi,deps
interface
subroutine memc_iau_nut00a (dt1,dt2,dpsi,deps) bind(c)
use,intrinsic :: iso_c_binding
real(c_double),value,intent(in) :: dt1
real(c_double),value,intent(in) :: dt2
real(c_double),intent(out) :: dpsi
real(c_double),intent(out) :: deps
end subroutine memc_iau_nut00a
end interface
call memc_iau_nut00a(dt1,dt2,dpsi,deps)
end subroutine
module test1
implicit none
contains
subroutine iau_nut00a (dt1,dt2,dpsi,deps)
implicit none
double precision dt1,dt2,dpsi,deps
interface
subroutine memc_iau_nut00a (dt1,dt2,dpsi,deps) bind(c)
use,intrinsic :: iso_c_binding
real(c_double),value,intent(in) :: dt1
real(c_double),value,intent(in) :: dt2
real(c_double),intent(out) :: dpsi
real(c_double),intent(out) :: deps
end subroutine memc_iau_nut00a
subroutine close_memc_connect() bind(c)
use, intrinsic :: iso_c_binding
end subroutine close_memc_connect
end interface
call memc_iau_nut00a(dt1,dt2,dpsi,deps)
end subroutine
double precision function iau_ee00a (dt1,dt2)
implicit none
double precision dt1,dt2
interface
real(c_double) function memc_iau_ee00a (dt1,dt2) bind(c)
use,intrinsic :: iso_c_binding
real(c_double),value,intent(in) :: dt1
real(c_double),value,intent(in) :: dt2
end function memc_iau_ee00a
end interface
iau_ee00a = memc_iau_ee00a(dt1,dt2)
end function
end module test1
program test
use test1
double precision adt1,adt2,adpsi,adeps,ee
adt1 = 2451547.1D00
adt2 = -1423.1D00
call iau_nut00a(adt1,adt2,adpsi,adeps)
ee = iau_ee00a (adt1,adt2)
write (*,*) adpsi,adeps,ee
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment