diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..50cd4ae --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +*CMakeFiles +*CMakeCache.txt +.cproject +.project +*.cmake +Makefile +source/main +source/libfepc.so +test/*test +*Testing + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..fba9b7b --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,18 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.6) + +PROJECT(FEPC) + +option(HAS_FFTW3 "Use fftw3" 0) + +if (HAS_FFTW3) + set(CMAKE_C_FLAGS "${CCMAKE_C_FLAGS} -DHAS_FFTW3") +endif(HAS_FFTW3) + +include_directories(include) + +ADD_SUBDIRECTORY(source) + +enable_testing() +ADD_SUBDIRECTORY(test) + + diff --git a/Makefile b/Makefile deleted file mode 100644 index cb8dd88..0000000 --- a/Makefile +++ /dev/null @@ -1,55 +0,0 @@ -CC=gcc -CFLAGS=-DHAS_FFTW3 -I/home/peet/Desktop/KO/include -LFLAGS=-L/home/peet/Desktop/KO/lib -lfftw3 -lm - - - -all: fepcd - - - -fepcd: main.o faltung.o kuerzen.o faltung_hilfe.o koeffizienten.o funktion.o folgen_vektor.o folge.o fft_faltung.o basic.o seconds.o - ${CC} ${CFLAGS} main.o faltung.o kuerzen.o faltung_hilfe.o koeffizienten.o funktion.o folgen_vektor.o folge.o fft_faltung.o basic.o seconds.o ${LFLAGS} -o fepcd - - - -main.o: main.c faltung.h kuerzen.h faltung_hilfe.h koeffizienten.h funktion.h folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} main.c - -faltung.o: faltung.c faltung.h kuerzen.h faltung_hilfe.h koeffizienten.h funktion.h folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} faltung.c - -kuerzen.o: kuerzen.c kuerzen.h faltung_hilfe.h koeffizienten.h funktion.h folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} kuerzen.c - -faltung_hilfe.o: faltung_hilfe.c faltung_hilfe.h koeffizienten.h funktion.h folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} faltung_hilfe.c - -koeffizienten.o: koeffizienten.c koeffizienten.h funktion.h folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} koeffizienten.c - -funktion.o: funktion.c funktion.h folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} funktion.c - -folgen_vektor.o: folgen_vektor.c folgen_vektor.h folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} folgen_vektor.c - -folge.o: folge.c folge.h fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} folge.c - -fft_faltung.o: fft_faltung.c fft_faltung.h basic.h config.h - ${CC} -c ${CFLAGS} fft_faltung.c - -basic.o: basic.c basic.h config.h - ${CC} -c ${CFLAGS} basic.c - -seconds.o: seconds.c seconds.h - ${CC} -c ${CFLAGS} seconds.c - - - -clean: - rm -f *~ *.o fepcd - - - diff --git a/README b/README index cd3feab..165ece4 100644 --- a/README +++ b/README @@ -1,6 +1,6 @@ FEPC library -Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) +Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2010 Stefan Handschuh (handschu@mis.mpg.de) This program is free software: you can redistribute it and/or modify @@ -17,5 +17,9 @@ To read the text of the GNU Lesser General Public License, see -This library provides several classes representing tensors in different formats -and algorithms to deal with them. + +Create the makefiles with "cmake .". If the parameter "HAS_FFTW3" is set, +the fftw3 library will be used (if present). This however forces the +resulting binary to be GPL licensed if published as the fftw3 is released +under the GPL. + diff --git a/basic.h b/include/basic.h similarity index 75% rename from basic.h rename to include/basic.h index a0c14d9..a799a35 100644 --- a/basic.h +++ b/include/basic.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,10 @@ #include "config.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { int dim; /* Dimension des Vektors */ int *array; /* Array mit Werten */ @@ -68,6 +72,10 @@ entry_d2one(vec_p s,vec_p n); vec_p entry_one2d(int pos,vec_p n); +/* Berechnet den zugehoerigen Vektor der Position pos eines multidimensionalen Array, dessen Dimension durch den Vektor n bestimmt ist, ohne dabei streng zu sein */ +vec_p +entry_one2d_sloppy(int pos,vec_p n); + /* Berechnet den Vektor a*s + b*n */ vec_p vec_op(int a, vec_p s, int b, vec_p n); @@ -76,6 +84,14 @@ vec_op(int a, vec_p s, int b, vec_p n); vec_p vec_multi(int a, vec_p n); +/* Multipliziert einen reelwertigen Vektor mit einer Konstante */ +vec_real_p +vec_real_multi(fepc_real_t factor, vec_real_p vector); + +/* Multipliziert einen reelwertigen Vektor mit einer Konstante und speichert das Ergebnis im Eingabevektor */ +void +vec_real_multi2(fepc_real_t factor, vec_real_p vector); + /* Dividiert den Vektor n mit der ganzen Zahl a */ vec_p vec_div(int a, vec_p n); @@ -84,14 +100,25 @@ vec_div(int a, vec_p n); vec_p vec_add(vec_p s, vec_p n); +/* Addiert Vektor s und Vektor n und speichert das Ergebnis in s */ +void +vec_add2(vec_p s, vec_p n); + /* Erzeugt eine inhaltliche Kopie des Vektors n */ vec_p vec_copy(vec_p n); +/* Subtrahiert Vektor s und Vektor n */ +vec_real_p +vec_real_substract(vec_real_p s, vec_real_p n); + /* Berechnet das Skalarprodukt der Vektoren r und s */ int vec_skalar_prod(vec_p r, vec_p s); +fepc_real_t +vec_real_skalar_prod(vec_real_p r, vec_real_p s); + /* Berechnet den groessten Vektor n sodass gilt: n<=r und n<=s */ vec_p vec_min(vec_p r, vec_p s); @@ -121,7 +148,28 @@ besitzt die Groesse 2. Es gilt Array[0] = s , Array[1] = r. */ vec_p* vec_zerlegung(vec_p n); +/** + * Prints out a vector. + */ +void +print_vec(vec_p vector); + +/** + * Prints out a real-valued vector + */ +void +print_vec_real(vec_real_p vector); +/* + * Returns the euklidian norm of the vector. + */ +fepc_real_t +vec_real_norm(vec_real_p vector); +#ifdef __cplusplus +} +#endif #endif + + diff --git a/config.h b/include/config.h similarity index 82% rename from config.h rename to include/config.h index c2eb87a..04a2bb2 100644 --- a/config.h +++ b/include/config.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -34,11 +34,22 @@ #include #include "seconds.h" +#ifdef __cplusplus +extern "C" { +#endif + /* * basic types */ -typedef enum { false, true } bool_t; +/* + * Adds c++ compatibility + */ +#if !defined(__cplusplus) + typedef enum { false, true } bool_t; +#else + typedef bool bool_t; +#endif typedef double fepc_real_t; @@ -69,4 +80,7 @@ typedef double fepc_real_t; #define _INLINE_ static #endif +#ifdef __cplusplus +} +#endif #endif /* __CONFIG_H */ diff --git a/include/discont.h b/include/discont.h new file mode 100644 index 0000000..3bc4578 --- /dev/null +++ b/include/discont.h @@ -0,0 +1,108 @@ +/* + * FEPC + * Copyright (C) 2010, 2011 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#ifndef _DISCONT +#define _DISCONT + +#define ONE_THIRD 0.333333333333333333333333333333333333333333333333333333333333333333333 +#define TWO_THIRD 0.666666666666666666666666666666666666666666666666666666666666666666666 +#define SQRT_12 3.464101615137754587054892683011744733885610507620761256111613958903866 + +#include "fepc_easy.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct { + fepc_real_t y_0; + + fepc_real_t slope; + +} linear_function_t; + +typedef linear_function_t * linear_function_p; + + +typedef struct { + int count; + + linear_function_p * functions; +} linear_function_set_t; + +typedef linear_function_set_t * linear_function_set_p; + +typedef struct { + int stepcount; + + interval_p * intervals; + + linear_function_set_p * function_sets; + +} discont_function_t; + +typedef discont_function_t * discont_function_p; + + +linear_function_p +linear_function_new(fepc_real_t x_0, fepc_real_t slope); + +linear_function_p +linear_function_new_points(fepc_real_t y_0, fepc_real_t y_1, fepc_real_t x_0, fepc_real_t x_1); + +void +discont_function_del(discont_function_p function); + +void +linear_function_set_del(linear_function_set_p function_set); + +linear_function_set_p +linear_function_set_new(int count); + +discont_function_p +discont_function_new(int stepcount); + +func_p +convert_discont_function(discont_function_p function, fepc_real_t stepping); + +discont_function_p +convert_func(func_p function, interval_p * intervals, fepc_real_t stepping); + +void +discont_function_print(discont_function_p function); + +void +linear_function_set_print(linear_function_set_p function_set); + +fepc_real_t +integrate_coeff_discont(discont_function_p function, int v, int position, int p, int step, fepc_real_t stepping); + +void +add_folgenentries_discont(func_p function, discont_function_p discont_function, fepc_real_t stepping); + +void +discont_function_setup(discont_function_p function, int step, fepc_real_t start, fepc_real_t end, linear_function_set_p function_set); + +void +discont_function_setup_points(discont_function_p function, int step, fepc_real_t start, fepc_real_t end, fepc_real_t * y1, fepc_real_t * y2, fepc_real_t stepping); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/faltung.h b/include/faltung.h similarity index 84% rename from faltung.h rename to include/faltung.h index b8ced7e..656cdf7 100644 --- a/faltung.h +++ b/include/faltung.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,10 @@ #include "kuerzen.h" +#ifdef __cplusplus +extern "C" { +#endif + /* Berechnet die Projektion der Faltung der Funktionen f und g auf die Gitterstruktur von w. Von f und g sind die Gitterstrukturen sowie die dazugehoerigen Eintraege gegeben. Von w ist nur die Gitterstruktur gegeben. @@ -37,4 +41,11 @@ faltung_ref(func_p f, func_p g, func_p w, fepc_real_t h); func_p faltung_fepc(func_p f, func_p g, func_p w, fepc_real_t h); +void +faltung_fepc_overwrite(func_t * result, func_p f, func_p g, func_p w, fepc_real_t h); + +#ifdef __cplusplus +} +#endif + #endif diff --git a/faltung_hilfe.h b/include/faltung_hilfe.h similarity index 92% rename from faltung_hilfe.h rename to include/faltung_hilfe.h index dfde42a..7c3bfd1 100644 --- a/faltung_hilfe.h +++ b/include/faltung_hilfe.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,9 @@ #include "koeffizienten.h" +#ifdef __cplusplus +extern "C" { +#endif /* Berechnung der Prolongation siehe Dokumentation */ folgen_vektor_p @@ -57,5 +60,18 @@ faltung_hilfe_Gamma_bauen_1( func_p g, vec_real_p mesh, vec_p maxgrad, matrix3_p func2_p faltung_hilfe_Gamma_bauen_2( func_p g, vec_real_p mesh, vec_p maxgrad, matrix3_p gamma_koef, matrix_p xi_koef ); +#ifdef __cplusplus +} +#endif #endif + + + + + + + + + + diff --git a/include/fepc_cp.h b/include/fepc_cp.h new file mode 100644 index 0000000..1fb191e --- /dev/null +++ b/include/fepc_cp.h @@ -0,0 +1,76 @@ +/* + * FEPC + * Copyright (C) 2010, 2011 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#ifndef __FEPCCP +#define __FEPCCP + + +#include "fepc_easy.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct { + int rank; + int dimension; + func_t * functions; +} func_cp; + +func_cp * +func_cp_new(int rank, int dimension, int * maxlevels); + +void +func_cp_init(func_cp * function, int rank, int dimension, int * maxlevels); + + +/*func_add + * Each summand has a different interval + */ +func_cp * +func_cp_new_blockstructure(int rank, int dimension, Funcimpl * functions, interval_t * intervals, int * maxlevels); + +/* + * Each summand has the same interval + */ +func_cp * +func_cp_new_cp(int rank, int dimension, Funcimpl * functions, interval_t * interval, int * maxlevels); + +func_t * +func_cp_extract(int current_rank, int current_dimension, func_cp * func_cp); + +void +func_cp_del(func_cp * func_cp); + +func_cp * +func_cp_faltung(func_cp * function1, func_cp * function2, func_t * resulting_structure, fepc_real_t h); + +func_cp * +func_cp_multi(func_cp * function1, func_cp * function2, fepc_real_t h); + +void +func_cp_print(func_cp * function); + +int * +int_array_new(int length); + +#ifdef __cplusplus +} +#endif + +#endif // __FEPCCP diff --git a/include/fepc_easy.h b/include/fepc_easy.h new file mode 100644 index 0000000..5580728 --- /dev/null +++ b/include/fepc_easy.h @@ -0,0 +1,168 @@ +/* + * FEPC + * Copyright (C) 2010, 2011 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#ifndef __FEPCEASY +#define __FEPCEASY + +#include +#include +#include +#include "faltung.h" +#include "interval.h" + +#ifdef __cplusplus +extern "C" { +#endif + +#ifndef INT_STEPS +#define INT_STEPS 10 +#endif + +#define SQRT_3 1.7320508075688772935274463415058723669428 +#define DIV_1_6 0.16666666666666666666666666666666666666666666 + +typedef fepc_real_t (*Funcimpl_step) (vec_p, vec_p, int, fepc_real_t); + +typedef vec_real_p (*Funcimpl_vec_step) (vec_p, vec_p, int, fepc_real_t); + +typedef fepc_real_t (*Funcimpl) (vec_real_p); + +typedef struct { + int count; + vec_real_p* elements; + fepc_real_t* element_results; +} vec_real_set_t; + +typedef vec_real_set_t * vec_real_set_p; + +/** + * Converts the given intervals from Omega to A form. + */ +interval_p * convert_interval(int steps, interval_p * intervals, fepc_real_t stepping); + +/** + * Sets the gridstructure in the needed form out of a common list of intervals. + */ +void set_gridstructure(func_p function, interval_p* intervals, fepc_real_t stepping); + +/** + * Creates a new set with allocated space for elements. + */ +vec_real_set_p vec_real_set_new(int count); + +/** + * Prints out all elements of a set. + */ +void print_vec_real_set(vec_real_set_p set); + +/** + * Returns the value of the function at the given point. + */ +vec_real_set_p get_value(func_p function, Funcimpl_vec_step generate_points, fepc_real_t stepping); + +void vec_real_set_del(vec_real_set_p set); + +/** + * Returns the h_l value corrensponding to the step. + */ +fepc_real_t get_h_l(int step, fepc_real_t stepping); + +/** + * Generates the vector of the folge at the given element_position. + */ +vec_p generate_folgenvector(folge_p folge, int element_position); + +/** + * Returns the number of elements of the given folge. + */ +int get_interval_element_count(folge_p folge); + +/** + * Integrates via trapezoidal rule. + */ +fepc_real_t integrate_coeff_st(Funcimpl function_impl, vec_p v, vec_p p, int step, fepc_real_t stepping); + +/** + * Integrates via Gauss-Legendre. + */ +fepc_real_t integrate_coeff(Funcimpl function_impl, vec_p v, vec_p p, int step, fepc_real_t stepping); + +/** + * Generates the x for integration based on the calc_position and the stepcount. + */ +vec_real_p generate_x_for_integration_st(vec_p v, int calc_position, fepc_real_t h, fepc_real_t h_l, int stepcount); + +/** + * Returns true if the vector is in the latter interval of the folge. + */ +bool_t is_in_latter_interval(vec_p v, folge_p folge); + +bool_t is_in_latter_interval_real(vec_real_p v, folge_p folge, int step, fepc_real_t stepping); + +/** + * Returns the value for phi_l (basis function) for the given vector, point x and step. + */ +fepc_real_t phi_l(int step, vec_p v, vec_p p, vec_real_p x, fepc_real_t stepping); + +/** + * Sets the entries of the folge. function_impl xor coeff_function may be NULL. coeff_function is preffered if both are given. + */ +void add_folgenentries(func_p function, Funcimpl function_impl, Funcimpl_step coeff_function, fepc_real_t stepping); + +void +setup_fepc_structure(func_t * func, Funcimpl function, interval_p* intervals, int interval_count, int degree, fepc_real_t stepping); + +func_p +create_fepc_structure(Funcimpl function, interval_p* intervals, int interval_count, int degree, fepc_real_t stepping); + +vec_real_set_p get_value(func_p function, Funcimpl_vec_step generate_points, fepc_real_t stepping); + +fepc_real_t get_value_at_step(func_p function, vec_real_p x, int step, fepc_real_t stepping); + +/** + * Finds the best interval. + */ +fepc_real_t get_value_at(func_p function, vec_real_p x, fepc_real_t stepping); + +/** + * Performs a convolution of the given functions in the given intervals with the given degree basis-functions. + */ +func_p fepc_easy(Funcimpl function1, interval_p* intervals1, int interval_count1, int degree_func1, Funcimpl function2, interval_p* intervals2, int interval_count2, int degree_func2, interval_p* intervals3, int interval_count3, fepc_real_t stepping); + + +fepc_real_t +func_integrate(func_p function, fepc_real_t stepping); + +func_p +func_derive(func_p function, int direction, fepc_real_t stepping); + +func_p +func_laplace(func_p function, fepc_real_t stepping); + +fepc_real_t max(fepc_real_t a, fepc_real_t b); + +int get_degree_count(vec_p degree); + +void set_degree(func_p function, int degree); + +#ifdef __cplusplus +} +#endif + +#endif + diff --git a/include/fepc_easy_helper.h b/include/fepc_easy_helper.h new file mode 100644 index 0000000..95f519a --- /dev/null +++ b/include/fepc_easy_helper.h @@ -0,0 +1,83 @@ +/* + * FEPC + * Copyright (C) 2010, 2011 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#ifndef _FEPC_EASY_HELPER +#define _FEPC_EASY_HELPER + +#include "fepc_easy.h" + +#ifdef __cplusplus +extern "C" { +#endif + +func_p +func_multi(func_p f, func_p g, fepc_real_t stepping); + +func_p +func_div(func_p f, func_p g, fepc_real_t stepping); + +func_p +func_reflect(func_p f); + +func_p +func_reflect_x(func_p f); + +folgen_vektor_p +folgen_vektor_multi(folgen_vektor_p f, folgen_vektor_p g, int step, fepc_real_t stepping); + +folgen_vektor_p +folgen_vektor_div(folgen_vektor_p f, folgen_vektor_p g, int step, fepc_real_t stepping); + +fepc_real_t folge_norm2(folge_p folge, int step, fepc_real_t stepping); + +folge_p +folge_multi(folge_p f, folge_p g, int step, fepc_real_t stepping); + +void +func_multi_overwrite(func_t * result, func_p f, func_p g, fepc_real_t stepping); + +folge_p +folge_div(folge_p f, folge_p g, int step, fepc_real_t stepping); + +vec_real_p +get_mean_points(vec_p v, vec_p grad, int step, fepc_real_t stepping); + +void +func_add2(func_p f, func_p g); + +void +func_add3(func_p f, fepc_real_t factor, func_p g); + + +func_p +func_subtract(func_p f, func_p g); + +void +func_subtract2(func_p f, func_p g); + +/** + * Calculates the square of L2 norm of the given function. + */ +fepc_real_t +func_norm_l2_sqr(func_p function, fepc_real_t stepping); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/include/fepc_easy_sparse.h b/include/fepc_easy_sparse.h new file mode 100644 index 0000000..1c91dd3 --- /dev/null +++ b/include/fepc_easy_sparse.h @@ -0,0 +1,60 @@ +/* + * FEPC + * Copyright (C) 2010, 2011 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#ifndef __FEPCEASYSPARSE +#define __FEPCEASYSPARSE + +#include "fepc_easy.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct { + int summands; + int dimension; + func_p * factors; +} func_sparse_t; + +typedef func_sparse_t * func_sparse_p; + + +func_sparse_p +func_sparse_new(int summands, int dimension); + +void +func_sparse_del(func_sparse_p func_sparse); + +func_p +func_sparse_to_func(func_sparse_p); + +func_sparse_p +func_sparse_multi(func_sparse_p f, func_sparse_p g, fepc_real_t stepping); + +func_sparse_p +func_sparse_reflect(func_sparse_p f); + + +func_sparse_p +func_sparse_convolute(func_sparse_p f, func_sparse_p g, func_p * result_structure, fepc_real_t stepping); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/fft_faltung.h b/include/fft_faltung.h similarity index 84% rename from fft_faltung.h rename to include/fft_faltung.h index 33b116a..6cca06f 100644 --- a/fft_faltung.h +++ b/include/fft_faltung.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,8 +21,16 @@ #include "basic.h" +#ifdef __cplusplus +extern "C" { +#endif + /*Berechnung der Faltung von zwei Arrays ueber die Fouriertransformation*/ fepc_real_t* fft_faltung(fepc_real_t* a, vec_p n_a, fepc_real_t* b, vec_p n_b); +#ifdef __cplusplus +} +#endif + #endif diff --git a/folge.h b/include/folge.h similarity index 88% rename from folge.h rename to include/folge.h index 0b1ca88..e2ce16c 100644 --- a/folge.h +++ b/include/folge.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -19,7 +19,11 @@ #ifndef __FOLGE_H #define __FOLGE_H -#include "fft_faltung.h" +#include "basic.h" + +#ifdef __cplusplus +extern "C" { +#endif typedef struct { vec_p start; /* Startvektor der Folge */ @@ -69,6 +73,10 @@ folge_norm(folge_p f, folge_p g); folge_p folge_add(folge_p f, folge_p g); +/* Subtraktion von Folgen f und g */ +folge_p +folge_subtract(folge_p f, folge_p g); + /* Rueckgabe einer einer Kopie der Folge f */ folge_p folge_copy( folge_p f); @@ -78,7 +86,13 @@ folge_copy( folge_p f); folge_p folge_projekt(folge_p f, folge_p g); +folge_p +folge_multi_factor(folge_p folge, fepc_real_t factor); - +#ifdef __cplusplus +} +#endif #endif + + diff --git a/folgen_vektor.h b/include/folgen_vektor.h similarity index 89% rename from folgen_vektor.h rename to include/folgen_vektor.h index 8606b32..e947b4e 100644 --- a/folgen_vektor.h +++ b/include/folgen_vektor.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,10 @@ #include "folge.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { vec_p grad; /*Polynomgrad des Folgenvektors*/ folge_p *vektor; /*Speicherung von Folgen durch multidimensionales Array dessen Dimension durch Vektor grad repraesentiert wird*/ @@ -96,9 +100,22 @@ folgen_vektor_projekt(folgen_vektor_p f,folgen_vektor_p g); folgen_vektor_p folgen_vektor_add(folgen_vektor_p f, folgen_vektor_p g); +/* Subtraktion von Folgenvektor f und g */ +folgen_vektor_p +folgen_vektor_subtract(folgen_vektor_p f, folgen_vektor_p g); + +/* Multiplikation von Folgenvektor f mit faktor a */ +folgen_vektor_p +folgen_vektor_factor_multi(folgen_vektor_p f, fepc_real_t a); + /* Addition von Folgenmatritzen f und g */ folgen_matrix_p folgen_matrix_add(folgen_matrix_p f, folgen_matrix_p g); +#ifdef __cplusplus +} +#endif #endif + + diff --git a/funktion.h b/include/funktion.h similarity index 88% rename from funktion.h rename to include/funktion.h index 9f0de7c..ed75e09 100644 --- a/funktion.h +++ b/include/funktion.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2010,2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,10 @@ #include "folgen_vektor.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { int dim; /* beschreibt zugrundeliegende Dimension */ int maxlevel; /* hoechste Level in Funktionendarstellung */ @@ -48,9 +52,15 @@ Somit ist die neu initialisierte Funktionen volldefiniert und vollwertig. */ func_p func_new(int maxlevel, int dim); +void +func_init(func_t * function, int maxlevel, int dimension); + void func_del(func_p f); +void +func_clear(func_t * f); + /* Funktion2 wird initialisiert. Beachte: alle Matritzen von Folgen der Funktion werden durch folgen_matrix_new( vec_new(dim), vec_new(dim) ) initialisiert. f = func2_new(maxlevel) bedeutet: f besteht aus maxlevel+1 Matritzen von Folgen f->hierarchie[k] fuer k=0,...,maxlevel. @@ -68,7 +78,8 @@ func2_del(func2_p f); func_p func_projekt(func_p f,func_p g); - +func_p +func_clone(func_p f); /* Gibt die einzelnen Details einer Funktion auf dem Bildschirm wieder. Je groesser der Wert info ist, umso mehr Informationen werden wiedergegeben. 0, 1, 2, 3 sind moegliche Infowerte. 0 wenig Info, 3 viele Infos */ @@ -90,6 +101,10 @@ func_build( int maxlevel, int dim , int grad, int a, int n, int mod, bool_t rand func_p func_add(func_p f, func_p g); +void +func_add_overwrite(func_t * result, func_t * f, func_t * g); + + /* Funktion gibt die Anzahl der Freiheitsgrade der Funktionen f, g und w wieder */ int func_count( func_p f, func_p g, func_p w ); @@ -109,4 +124,17 @@ void func_grid_zero(func_p f); +func_p +func_factor_multi(func_p function, fepc_real_t factor); + +void +funcs_del(func_p * array, int length); + +void +funcs_del_type(func_t * array, int length); + +#ifdef __cplusplus +} +#endif + #endif diff --git a/include/interval.h b/include/interval.h new file mode 100644 index 0000000..607635d --- /dev/null +++ b/include/interval.h @@ -0,0 +1,80 @@ +/* + * FEPC + * Copyright (C) 2010, 2011 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#ifndef __INTERVALS +#define __INTERVALS + +#include "basic.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct { + int dimension; + fepc_real_t* start; + fepc_real_t* end; +} interval_t; + +typedef interval_t * interval_p; + +/** + * Creates a new interval with the given dimension. + */ +interval_p +interval_new(int dimension); + +void +interval_init(interval_p interval, int dimension); + +void +interval_del(interval_p interval); + +void +interval_clear(interval_t * interval); + +void +intervals_del(interval_p* intervals, int interval_count); + +void +intervals_del_type(interval_t * intervals, int interval_count); + +/** + * Prints out the range of the given interval. + */ +void +print_interval(interval_p interval); + +/** + * Prints out the ranges of all given intervals. + */ +void +print_intervals(interval_p * intervals, int count); + +interval_p +interval_clone(interval_p interval); + +interval_p * +intervals_clone(interval_p *interval, int count); + +#ifdef __cplusplus +} +#endif + +#endif + diff --git a/koeffizienten.h b/include/koeffizienten.h similarity index 93% rename from koeffizienten.h rename to include/koeffizienten.h index dbc4eb5..24c77ef 100644 --- a/koeffizienten.h +++ b/include/koeffizienten.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,9 @@ #include "funktion.h" +#ifdef __cplusplus +extern "C" { +#endif typedef struct { int zeilen; /* Anzahl der Zeilen */ @@ -80,5 +83,11 @@ koeffizienten_xi(vec_p alpha, vec_p q, matrix_p xi_koef); fepc_real_t koeffizienten_gamma(int level, vec_p r, vec_real_p mesh, vec_p alpha, vec_p mu, vec_p kappa, matrix3_p gamma_koef); +#ifdef __cplusplus +} +#endif #endif + + + diff --git a/kuerzen.h b/include/kuerzen.h similarity index 95% rename from kuerzen.h rename to include/kuerzen.h index ef516c7..31bf957 100644 --- a/kuerzen.h +++ b/include/kuerzen.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), Stefan Handschuh 2011 (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,6 +21,10 @@ #include "faltung_hilfe.h" +#ifdef __cplusplus +extern "C" { +#endif + /* Diese Struktur dient der Charakterisierung des Traegers mehrerer Folgenvektoren. Die Struktur dient der Umsetzung der Ergebnisse, die in der Dokumentation zum Kapitel "Einkuerzung" gemacht wurden. Ausgehend von einer gegebenen Funktion f des Maxlevel L werden fuer 0<=l<=L mittels der Arrays a und b Traeger fuer die Folgenvektoren f->hierarchie[l] abgespeichert. Das Intervall [ a[l], b[l] ] diff --git a/seconds.h b/include/seconds.h similarity index 82% rename from seconds.h rename to include/seconds.h index 8e1e6af..6047a08 100644 --- a/seconds.h +++ b/include/seconds.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -19,6 +19,10 @@ #ifndef __SECONDS_H #define __SECONDS_H +#ifdef __cplusplus +extern "C" { +#endif + /* * return time in seconds used by program */ @@ -26,4 +30,8 @@ double seconds (); +#ifdef __cplusplus +} +#endif + #endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt new file mode 100644 index 0000000..6ac8583 --- /dev/null +++ b/source/CMakeLists.txt @@ -0,0 +1,30 @@ +set(ALL_SOURCE_FILES + basic.c + faltung_hilfe.c + folge.c + fepc_cp.c + folgen_vektor.c + fepc_easy.c + funktion.c + fepc_easy_helper.c + interval.c + discont.c + fepc_easy_sparse.c + koeffizienten.c + seconds.c + faltung.c + kuerzen.c) + +if (HAS_FFTW3) + set(ALL_SOURCE_FILES fft_faltung.c ${ALL_SOURCE_FILES}) +endif(HAS_FFTW3) + +ADD_LIBRARY(fepc SHARED ${ALL_SOURCE_FILES}) + +add_executable (main main.c) +target_link_libraries(main fepc m) + +if (HAS_FFTW3) + target_link_libraries(main fftw3) +endif(HAS_FFTW3) + diff --git a/basic.c b/source/basic.c similarity index 78% rename from basic.c rename to source/basic.c index 384e516..e9a9586 100644 --- a/basic.c +++ b/source/basic.c @@ -75,8 +75,6 @@ vec_one(int dim) { } - - void vec_del(vec_p s) { free(s->array); @@ -149,6 +147,30 @@ entry_one2d(int pos,vec_p n) { return s; } +vec_p +entry_one2d_sloppy(int pos,vec_p n) { + int k, size, test; + vec_p s; + + size = vec_size( n ); + ASSERT(pos >= 0); + + /*Berechnung*/ + s = vec_new( n->dim ); + + for(k=s->dim-1;k>=1;k--) { + if (n->array[k] == 0) { + s->array[k] = 0; + } else { + s->array[k] = pos % n->array[k]; + pos = (pos - s->array[k]) / n->array[k]; + } + } + s->array[0] = pos; + + return s; +} + vec_p vec_add(vec_p s, vec_p n) { int k; @@ -164,6 +186,33 @@ vec_add(vec_p s, vec_p n) { return back; } +void +vec_add2(vec_p s, vec_p n) { + int k; + + /*Testen auf Konsistenz (siehe Dokumentation)*/ + ASSERT( s->dim == n->dim ); + + for(k=0; kdim;k++) { + s->array[k] += n->array[k]; + } +} + +vec_real_p +vec_real_substract(vec_real_p s, vec_real_p n) { + int k; + vec_real_p back; + + /*Testen auf Konsistenz (siehe Dokumentation)*/ + ASSERT( s->dim == n->dim ); + + back = vec_real_new( s->dim ); + for(k=0; kdim;k++) { + back->array[k] = s->array[k] - n->array[k]; + } + return back; +} + vec_p vec_multi(int a, vec_p n) { @@ -189,6 +238,29 @@ vec_div(int a, vec_p n) { return back; } + +vec_real_p +vec_real_multi(fepc_real_t factor, vec_real_p vector) { + int n; + + vec_real_p back = vec_real_new(vector->dim); + + for (n = 0; n < vector->dim; n++) { + back->array[n] = vector->array[n]*factor; + } + return back; +} + + +void +vec_real_multi2(fepc_real_t factor, vec_real_p vector){ + int n; + + for (n = 0; n < vector->dim; n++) { + vector->array[n] *= factor; + } +} + vec_p vec_copy(vec_p n) { int k; @@ -215,6 +287,21 @@ vec_skalar_prod(vec_p r, vec_p s) { return sum; } +fepc_real_t +vec_real_skalar_prod(vec_real_p r, vec_real_p s) { + int k; + + fepc_real_t sum = 0.; + + /*Testen auf Konsistenz (siehe Dokumentation)*/ + ASSERT( s->dim == r->dim ); + + for(k=0;kdim;k++) { + sum += r->array[k] * s->array[k]; + } + return sum; +} + vec_p vec_min(vec_p r, vec_p s) { int k; @@ -271,6 +358,7 @@ vec_size(vec_p r) { + vec_p vec_r_s_n(vec_p r, vec_p n) { int k, d, dim; @@ -442,3 +530,33 @@ vec_0_s_r(vec_p r, vec_p n) { vec_del( modulo ); return back; } + + + +void +print_vec(vec_p vector) { + int n; + + for (n = 0; n < vector->dim; n++) { + printf("\t%i\n", vector->array[n]); + } +} + + + +void +print_vec_real(vec_real_p vector) { + int n; + + for (n = 0; n < vector->dim; n++) { + printf("\t%f\n", vector->array[n]); + } +} + +fepc_real_t +vec_real_norm(vec_real_p vector) { + return sqrt(vec_real_skalar_prod(vector, vector));; +} + + + diff --git a/source/discont.c b/source/discont.c new file mode 100644 index 0000000..81b7e0c --- /dev/null +++ b/source/discont.c @@ -0,0 +1,242 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "discont.h" + + +fepc_real_t +_round(fepc_real_t value) { + if (value - ((int) value) < 0.5) { + return (fepc_real_t) (int) value; + } else { + return (fepc_real_t) (int) value +1.0; + } +} + +linear_function_p +linear_function_new(fepc_real_t y_0, fepc_real_t slope) { + linear_function_p result = (linear_function_p) malloc(sizeof(linear_function_t)); + + result->y_0 = y_0; + result->slope = slope; + return result; +} + + +linear_function_p +linear_function_new_points(fepc_real_t y_0, fepc_real_t y_1, fepc_real_t x_0, fepc_real_t x_1) { + fepc_real_t slope = (y_1-y_0) / (x_1-x_0); + + return linear_function_new(y_0-slope*x_0, slope); +} + + +linear_function_set_p +linear_function_set_new(int count) { + linear_function_set_p result = (linear_function_set_p) malloc(sizeof(linear_function_set_t)); + + result->count = count; + result->functions = (linear_function_p*) malloc(sizeof(linear_function_p)*count); + return result; +} + + +discont_function_p +discont_function_new(int stepcount) { + discont_function_p result = (discont_function_p) malloc(sizeof(discont_function_t)); + + result->stepcount = stepcount; + result->intervals = (interval_p*) malloc(sizeof(interval_p)*stepcount); + result->function_sets = (linear_function_set_p*) malloc(sizeof(linear_function_set_p)*stepcount); + return result; +} + +void +discont_function_del(discont_function_p function) { + int n; + + for (n = 0; n < function->stepcount; n++) { + interval_del(function->intervals[n]); + linear_function_set_del(function->function_sets[n]); + } + free(function->intervals); + free(function->function_sets); + free(function); +} + +void +linear_function_set_del(linear_function_set_p function_set) { + int n; + + + for (n = 0; n < function_set->count; n++) { + if (function_set->functions[n] != NULL) { + free(function_set->functions[n]); + } + } + free(function_set->functions); + free(function_set); +} + +void +discont_function_setup(discont_function_p function, int step, fepc_real_t start, fepc_real_t end, linear_function_set_p function_set) { + function->intervals[step] = interval_new(1); + function->intervals[step]->start[0] = start; + function->intervals[step]->end[0] = end; + function->function_sets[step] = function_set; +} + +void +discont_function_setup_points(discont_function_p function, int step, fepc_real_t start, fepc_real_t end, fepc_real_t * y1, fepc_real_t * y2, fepc_real_t stepping) { + function->intervals[step] = interval_new(1); + function->intervals[step]->start[0] = start; + function->intervals[step]->end[0] = end; + + fepc_real_t h_l = get_h_l(step, stepping); + + int count, n; + + count = _round((end-start)/h_l); // this has to be equal to the array sizes of y1 and y2 + + function->function_sets[step] = linear_function_set_new(count); + + if (y1 != NULL && y2 != NULL) { + for (n = 0; n < count; n++) { + function->function_sets[step]->functions[n] = linear_function_new_points(y1[n], y2[n], start+n*h_l, start+(n+1)*h_l ); + } + } else { + for (n = 0; n < count; n++) { + function->function_sets[step]->functions[n] = NULL; + } + } +} + + +void +discont_function_print(discont_function_p function) { + int n; + + printf("Discontinuous function with %i levels\n=====================================\n", function->stepcount); + + for (n = 0; n < function->stepcount; n++) { + + printf("Level %i\n", n); + print_interval(function->intervals[n]); + linear_function_set_print(function->function_sets[n]); + } +} + + +void +linear_function_set_print(linear_function_set_p function_set) { + int n; + + printf("Linear function set\n=====\n"); + for (n = 0; n < function_set->count; n++) { + printf(" %f*x + %f\n", function_set->functions[n]->slope, function_set->functions[n]->y_0); + } +} + + +func_p +convert_discont_function(discont_function_p function, fepc_real_t stepping) { + func_p result = func_new(function->stepcount-1, 1); + + set_degree(result, 1); // sets the polynomial degree to 1 + set_gridstructure(result, function->intervals, stepping); // set the correct intervals + add_folgenentries_discont(result, function, stepping); // add folgenentries + return result; +} + + +discont_function_p +convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { + discont_function_p result = discont_function_new(function->maxlevel+1); + + fepc_real_t h_l, temp_x1, temp_x2, temp_y1, temp_y2, slope; + + int n, k, stepcount, start; + + vec_real_p x; + + for (n = 0; n < result->stepcount; n++) { + h_l = get_h_l(n, stepping); + result->intervals[n] = intervals[n]; + stepcount = _round((intervals[n]->end[0] - intervals[n]->start[0])/h_l); // number of interval points + result->function_sets[n] = linear_function_set_new(stepcount); + start = function->hierarchie[n]->vektor[0]->start->array[0]; + for (k = 0; k < stepcount; k++) { + temp_x1 = (ONE_THIRD+start+k)*h_l; + temp_x2 = (TWO_THIRD+start+k)*h_l; + x = vec_real_new(1); + x->array[0] = temp_x1; + temp_y1 = get_value_at_step(function, x, n, stepping); + x->array[0] = temp_x2; + temp_y2 = get_value_at_step(function, x, n, stepping); + vec_real_del(x); + slope = 3*(temp_y2-temp_y1)/h_l; + result->function_sets[n]->functions[k] = linear_function_new(temp_y1 - slope*temp_x1, slope); + } + } + return result; +} + + +fepc_real_t +integrate_coeff_discont(discont_function_p function, int v, int position, int p, int level, fepc_real_t stepping) { + // integrate from v[0]*h_l till (v[0]+1)*h_l + + fepc_real_t h_l, slope, y_0; + h_l = get_h_l(level, stepping); + + slope = function->function_sets[level]->functions[position]->slope; + y_0 = function->function_sets[level]->functions[position]->y_0; + if (p == 0) { // legendre(0) = sqrt(1/h_l) + return sqrt(h_l)*(y_0 + h_l*(slope)*(v+0.5)); + } else { // p == 1 --> legendre(1) = sqrt(12)(x-(v+0.5)*h_l)/(h_l^1.5) + return h_l*sqrt(h_l)*slope/SQRT_12; + } +} + + +void +add_folgenentries_discont(func_p function, discont_function_p discont_function, fepc_real_t stepping) { + int level, n, k, interval_element_count, position, grad_count, pos_2; + + vec_p v, p; + + for (level = 0; level <= function->maxlevel; level++) { + grad_count = get_degree_count(function->hierarchie[level]->grad); + for (k = 0; k < grad_count; k++) { + for (n = 0, interval_element_count = get_interval_element_count(function->hierarchie[level]->vektor[k]); n < interval_element_count; n++) { + v = generate_folgenvector(function->hierarchie[level]->vektor[k], n); + position = entry_d2one(v, function->hierarchie[level]->vektor[k]->lang); + pos_2 = v->array[0]; // position + vec_add2(v, function->hierarchie[level]->vektor[k]->start); // needed bc v starts from 0 + if (level == function->maxlevel || !is_in_latter_interval(v, function->hierarchie[level+1]->vektor[k])) { // only set the vector if it is in the valid interval + p = entry_one2d_sloppy(k, function->hierarchie[level]->grad); // degree - don't care about 0 as entries + function->hierarchie[level]->vektor[k]->glied[position] = integrate_coeff_discont(discont_function, v->array[0], pos_2, p->array[0], level, stepping); + vec_del(p); + } + vec_del(v); + } + } + } +} + + diff --git a/faltung.c b/source/faltung.c similarity index 98% rename from faltung.c rename to source/faltung.c index d0c7d8b..f8df10c 100644 --- a/faltung.c +++ b/source/faltung.c @@ -614,7 +614,16 @@ faltung_ref(func_p f, func_p g, func_p w, fepc_real_t h) { func_p faltung_fepc(func_p f, func_p g, func_p w, fepc_real_t h) { - func_p back, temp1, temp2; + func_p back = func_new(w->maxlevel, w->dim); + + faltung_fepc_overwrite(back, f, g, w, h); + return back; +} + +void +faltung_fepc_overwrite(func_t * result, func_p f, func_p g, func_p w, fepc_real_t h) { + func_p temp1, temp2; + int k, maxgrad_1dim, dim; vec_p maxgrad, p; matrix_p xi_koef; @@ -647,32 +656,30 @@ faltung_fepc(func_p f, func_p g, func_p w, fepc_real_t h) { vec_del( maxgrad ); maxgrad = p; } - maxgrad_1dim = 0; for(k=0;karray[k] ) { maxgrad_1dim = maxgrad->array[k]; } } - /*Berechnung der eindimensionalen Xi-Koeffizienten*/ xi_koef = koeffizienten_xi_1dim( 2 * maxgrad_1dim + 1 ); /*Berechnung der eindimensionalen Gamma-Koeffizienten*/ gamma_koef = koeffizienten_gamma_1dim( 2 * maxgrad_1dim + 1 ); - /*Berechnung der Faltung*/ temp1 = faltung_sum1( f, g, w, mesh, maxgrad, gamma_koef, xi_koef ); temp2 = faltung_sum2( f, g, w, mesh, maxgrad, gamma_koef, xi_koef ); - back = func_add( temp1, temp2 ); - + func_add_overwrite(result, temp1, temp2 ); vec_del( maxgrad ); func_del( temp1 ); func_del( temp2 ); matrix3_del( gamma_koef ); matrix_del( xi_koef ); vec_real_del( mesh ); - /*Folgenwerte, die sich auf Intervallen befinden, welche verfeinert werden, werden gleich Null gesetzt*/ - func_grid_zero( back ); - return back; + func_grid_zero( result ); } + + + + diff --git a/faltung_hilfe.c b/source/faltung_hilfe.c similarity index 99% rename from faltung_hilfe.c rename to source/faltung_hilfe.c index 6059dda..24cf201 100644 --- a/faltung_hilfe.c +++ b/source/faltung_hilfe.c @@ -966,3 +966,4 @@ lambda(folgen_vektor_p g, int level, vec_real_p mesh, matrix3_p gamma_koef, vec_ return back; } */ + diff --git a/source/fepc_cp.c b/source/fepc_cp.c new file mode 100644 index 0000000..3b4d9d9 --- /dev/null +++ b/source/fepc_cp.c @@ -0,0 +1,152 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_cp.h" +#include "fepc_easy_helper.h" + +func_cp * +func_cp_new(int rank, int dimension, int * maxlevels) { + func_cp * result = (func_cp*) malloc(sizeof(func_cp)); + + func_cp_init(result, rank, dimension, maxlevels); + return result; +} + +void +func_cp_init(func_cp * function, int rank, int dimension, int * maxlevels) { + int n, k = rank*dimension; + function->rank = rank; + function->dimension = dimension; + function->functions = (func_t *) malloc(sizeof(func_t)*k); + + for (n = 0; n < k; n++) { + func_init(&(function->functions[n]), maxlevels[n], 1); + } +} + +func_cp * +func_cp_new_blockstructure(int rank, int dimension, Funcimpl * functions, interval_p intervals, int * maxlevels) { + func_cp * result = func_cp_new(rank, dimension, maxlevels); + + // TODO implement + return result; +} + +func_cp * +func_cp_new_cp(int rank, int dimension, Funcimpl * functions, interval_p interval, int * maxlevels) { + func_cp * result = func_cp_new(rank, dimension, maxlevels); + + + return result; +} + +func_t * +func_cp_extract(int current_rank, int current_dimension, func_cp * func_cp) { + ASSERT(current_rank < func_cp->rank); + ASSERT(current_dimension < func_cp->dimension); + return &(func_cp->functions[current_rank*func_cp->dimension+current_dimension]); +} + +func_cp * +func_cp_faltung(func_cp * function1, func_cp * function2, func_t * resulting_structure, fepc_real_t h) { + ASSERT(function1->dimension == function2->dimension); + + int n, levels_count = function1->rank*function2->rank*function1->dimension; + + int * maxlevels = (int*) malloc(sizeof(int)*levels_count); + + for (n = 0; n < levels_count; n++) { + maxlevels[n] = resulting_structure[n % function1->dimension].maxlevel; // since the resulting structure is not BLOCK + } + + func_cp * result = func_cp_new(function1->rank*function2->rank, function1->dimension, maxlevels); + + free(maxlevels); + + int k, d; + + for (n = 0; n < function1->rank; n++) { + for (k = 0; k < function2->rank; k++) { + for (d = 0; d < result->dimension; d++) { + faltung_fepc_overwrite(&(result->functions[(n*function2->rank + k)*result->dimension+d]), func_cp_extract(n, d, function1), func_cp_extract(k, d, function2), &resulting_structure[d], h); + } + } + } + + return result; +} + +void +func_cp_del(func_cp * func_cp) { + funcs_del_type(func_cp->functions, func_cp->dimension*func_cp->rank); + free(func_cp); +} + +func_cp * +func_cp_multi(func_cp * function1, func_cp * function2, fepc_real_t h) { + ASSERT(function1->dimension == function2->dimension); + + int levels_count = function1->rank*function2->rank*function1->dimension; + + int * maxlevels = int_array_new(levels_count); // FIXME: this has to be the maximum over all levels in each entry + + func_cp * result = func_cp_new(function1->rank*function2->rank, function1->dimension, maxlevels); + + int n, k, d; + + for (n = 0; n < function1->rank; n++) { + for (k = 0; k < function2->rank; k++) { + for (d = 0; d < result->dimension; d++) { + func_multi_overwrite(&(result->functions[(n*function2->rank + k)*result->dimension+d]), func_cp_extract(n, d, function1), func_cp_extract(k, d, function2), h); + } + } + } + free(maxlevels); + return result; +} + +void +func_cp_print(func_cp * cp) { + int n, k; + + printf("CP function with rank %i:\n", cp->rank); + for (n = 0; n < cp->rank; n++) { + for (k = 0; k < cp->dimension; k++) { + func_print(func_cp_extract(n, k, cp), 3); + } + if (n +1 < cp->rank) { + printf("\n+"); + } else { + printf("\n\n"); + } + } +} + +int * +int_array_new(int length) { + int n; + + int * result = (int*) malloc(sizeof(int)*length); + + for (n = 0; n < length; n++) { + result[n] = 0; + } + return result; +} + + diff --git a/source/fepc_easy.c b/source/fepc_easy.c new file mode 100644 index 0000000..5722bdd --- /dev/null +++ b/source/fepc_easy.c @@ -0,0 +1,904 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_easy.h" + + +/** + * This method converts standard d-dimensional intervals to the form that is needed by the fepc algorithm. + * Only called for debug. The transformation is also done inside of set_gridstructure. + * + * @param steps defines the number of steps that are represented by the invervals; equals the number of intervals + * @param intervals array of intervals that should be converted + */ +interval_p * convert_interval(int steps, interval_p * intervals, fepc_real_t stepping) { + interval_p* result = (interval_p*) malloc(sizeof(interval_p) * steps); + + fepc_real_t pow, h_l; + + pow = 2.; + + int n, k; + + for (n = 0; n < steps; n++) { + pow *= 0.5; + h_l = pow*stepping; + result[n] = interval_new(intervals[n]->dimension); + for (k = 0; k < result[n]->dimension; k++) { + result[n]->start[k] = intervals[n]->start[k] / h_l; + result[n]->end[k] = intervals[n]->end[k] / h_l - 1; + } + } + + return result; +} + +/** + * Prints out all elements of the set and its results. + * + * @param set set to print out + */ +void print_vec_real_set(vec_real_set_p set) { + int n, k; + printf(" Elements \n"); + for (n = 0; n < set->count; n++) { + for (k = 0; k < set->elements[n]->dim; k++) { + printf("%.5f ", set->elements[n]->array[k]); + } + printf("%.5f\n", set->element_results[n]); + } +} + +/** + * Creates a new set with allocated memory for the entries. + * + * @param count number of entries + * + * @return new set + */ +vec_real_set_p vec_real_set_new(int count) { + vec_real_set_p result = (vec_real_set_p) malloc(sizeof(vec_real_set_t)); + + result->count = count; + result->elements = (vec_real_p*) malloc(sizeof(vec_real_p)*count); + result->element_results = (fepc_real_t*) malloc(sizeof(fepc_real_t)*count); + + int n; + + for (n = 0; n < count; n++) { + result->element_results[n] = 0.0; + } + return result; +} + + +void vec_real_set_del(vec_real_set_p set) { + int n; + + for (n = 0; n < set->count; n++) { + vec_real_del(set->elements[n]); + } + free(set->elements); + free(set->element_results); + free(set); +} + + +fepc_real_t get_h_l(int step, fepc_real_t stepping) { + return pow(.5, step)*stepping; +} + +/** + * Returns the value of the legendre-polynomial with the given degree at the given point. + * + * @param degree degree of the legendre-polynomial + * @param x point where the polynomial should be evaluated + * + * @return P_degree(x) + */ +fepc_real_t legendre(int degree, fepc_real_t x) { + if (degree == 0) { + return 1; + } else if (degree == 1) { + return x; + } else { + return ((2*degree-1)*x*legendre(degree-1, x) - (degree-1)*legendre(degree-2, x)) / degree; + } +} + +/** + * This method calculates the basis function. + * + * @param step current step + * @param v vector that describes the integral-bounds; integral is taken from v[n]*h_l upto (v[n]+1)*h_l + * @param p degree-vector of the basis function + * @param x point at which the basis function is evaluated + * + * @return phi_{v,p}(x) + */ +fepc_real_t phi_l(int step, vec_p v, vec_p p, vec_real_p x, fepc_real_t stepping) { + fepc_real_t h_l, result; + + h_l = get_h_l(step, stepping); + result = 1.; + + int n; + + for (n = 0; n < v->dim; n++) { + if (v->array[n]*h_l > x->array[n] || (v->array[n]+1)*h_l < x->array[n]) { + return 0; + } + result *= sqrt((2*p->array[n] + 1)/h_l)*legendre(p->array[n], 2.0*x->array[n]/h_l - 2.0*v->array[n] - 1.0); + } + return result; +} + +/** + * Retransforms the resulting points from S(M) into R^n. + * This method should be used at the end of the calculation to regain the actual R^n coodinates of the resulting points. + * + * @param function function in S(M); this should be a result of fecp + * @param generate_points a function that generates the points where the R^n values should be taken + * + * @return a set of all valid points + */ +vec_real_set_p get_value(func_p function, Funcimpl_vec_step generate_points, fepc_real_t stepping) { + vec_p r, x; + + int n, k, m, l, intervals_element_count, position, degree_count; + + int* interval_element_counts = (int*) malloc(sizeof(int)*(function->maxlevel+1)); + + // get the total count of interval elements as well as the separate interval lengths + for (n = 0, intervals_element_count = 0; n <= function->maxlevel; n++) { + interval_element_counts[n] = get_interval_element_count(function->hierarchie[n]->vektor[0]); // intervals have all the same length, so we can use any "vektor" + intervals_element_count += interval_element_counts[n]; + } + vec_real_set_p result = vec_real_set_new(intervals_element_count); + for (n = 0, m = 0; n <= function->maxlevel; n++) { + degree_count = get_degree_count(function->hierarchie[n]->grad); + for (k = 0; k < interval_element_counts[n]; k++, m++) { + for (l = 0; l < degree_count; l++) { + r = generate_folgenvector(function->hierarchie[n]->vektor[l], k); + position = entry_d2one(r, function->hierarchie[n]->vektor[l]->lang); + vec_add2(r, function->hierarchie[n]->vektor[l]->start); + // only take the point if it is inside of the correct A-interval; if not in interval -> reduce size of set + if (n == function->maxlevel || !is_in_latter_interval(r, function->hierarchie[n+1]->vektor[l])) { + x = entry_one2d_sloppy(l, function->hierarchie[n]->grad); + result->elements[m] = generate_points(r, x, n, stepping); + result->element_results[m] += function->hierarchie[n]->vektor[l]->glied[position]* phi_l(n, r, x, result->elements[m], stepping); + vec_del(x); + } else { + m--; + result->count--; + // TODO free also memory? + } + vec_del(r); + } + } + } + free(interval_element_counts); + return result; +} + +fepc_real_t get_value_at_step(func_p function, vec_real_p x, int step, fepc_real_t stepping) { + ASSERT(function != NULL && step <= function->maxlevel); + + int n, l, k, degree_count, interval_element_count, position; + + vec_p r, p; + + fepc_real_t result = 0.; + + degree_count = get_degree_count(function->hierarchie[step]->grad); + interval_element_count = get_interval_element_count(function->hierarchie[step]->vektor[0]); // intervals have all the same length, so we can use any "vektor" + for (n = 0; n < interval_element_count; n++) { + for (k = 0; k < degree_count; k++) { + r = generate_folgenvector(function->hierarchie[step]->vektor[k], n); + position = entry_d2one(r, function->hierarchie[step]->vektor[k]->lang); + vec_add2(r, function->hierarchie[step]->vektor[k]->start); + // only take the point if it is inside of the correct A-interval + if (step == function->maxlevel || !is_in_latter_interval(r, function->hierarchie[step+1]->vektor[k])) { + p = entry_one2d_sloppy(k, function->hierarchie[step]->grad); + result += function->hierarchie[step]->vektor[k]->glied[position]* phi_l(step, r, p, x, stepping); + vec_del(p); + } + vec_del(r); + } + } + return result; + +} + +fepc_real_t get_value_at(func_p function, vec_real_p x, fepc_real_t stepping) { + int n; + + for (n = 0; n < function->maxlevel; n++) { + if (!is_in_latter_interval_real(x, function->hierarchie[n+1]->vektor[0], n+1, stepping)) { // the intervals are at the level constant, independent of the degree + return get_value_at_step(function, x, n, stepping);//*(function->maxlevel - n +1); + } + } + return get_value_at_step(function, x, function->maxlevel, stepping); +} + +double round(double x) { + fepc_real_t rest = x - (int) x; + + if (rest < 0.5) { + return (int) x; + } else { + return 1 + (int) x; + } +} + +/** + * Sets the correct grid-structure that is needed by the fepc algorithm out of common human-readable intervals. + * + * @param function function whose grid structure should be set + * @param intervals human-readable intervals; note the the inverval count has to be function->maxlevel+1 + */ +void set_gridstructure(func_p function, interval_p* intervals, fepc_real_t stepping) { + vec_p start, lang; + + int n, dimension, k, steps, grad_count; + + steps = function->maxlevel+1; + dimension = function->dim; + + fepc_real_t pow, h_l; + + pow = 2.; + + folge_p folge; + + for (n = 0; n < steps; n++) { + start = vec_new(dimension); + lang = vec_new(dimension); + + pow *= 0.5; + h_l = pow*stepping; + + if (intervals[n] != NULL) + for (k = 0; k < dimension; k++) { // here the interval is converted to the nessesary form + start->array[k] = intervals[n]->start[k] / h_l; + lang->array[k] = round(intervals[n]->end[k] / h_l) - start->array[k]; // lang = length (integer), see p. 68; round() is needed to correct error + } + + folge = folge_new(start, lang); + folge_del(function->hierarchie[n]->vektor[0]); + function->hierarchie[n]->vektor[0] = folge; + for (k = 1, grad_count = get_degree_count(function->hierarchie[n]->grad); k < grad_count; k++) { + folge_del(function->hierarchie[n]->vektor[k]); + function->hierarchie[n]->vektor[k] = folge_copy(folge); + } + + #ifdef __DEBUG + printf("set grid structure\n start\n"); + print_vec(start); + printf(" lang\n"); + print_vec(lang); + #endif + } +} + +fepc_real_t sqr(fepc_real_t x) { + return x*x; +} + +fepc_real_t integrate_legendre(vec_p v, vec_p degree, int step, fepc_real_t stepping) { + int n; + + fepc_real_t h_l = get_h_l(step, stepping), result = 1.; + + for (n = 0; n < v->dim; n++) { + if (degree->array[n] == 0) { + result *= sqrt(h_l); + } else { + result *= sqrt(h_l/(4*(2*degree->array[n] + 1)))*(legendre(degree->array[n] + 1, -1) - legendre(degree->array[n]-1, -1)); + } + //result *= sqrt((2*degree->array[n] + 1)/h_l)*h_l;//(legendre(degree->array[n], sqr((v->array[n]+1)*h_l)/2.) - legendre(degree->array[n], sqr(v->array[n]*h_l)/2.)); + } +} + + +fepc_real_t func_integrate(func_p function, fepc_real_t stepping) { + int step, n, k, interval_element_count, position, grad_count; + + vec_p v, x; + + fepc_real_t result = 0.; + + for (step = 0; step <= function->maxlevel; step++) { + grad_count = get_degree_count(function->hierarchie[step]->grad); + for (k = 0; k < grad_count; k++) { + for (n = 0, interval_element_count = get_interval_element_count(function->hierarchie[step]->vektor[k]); n < interval_element_count; n++) { + v = generate_folgenvector(function->hierarchie[step]->vektor[k], n); + position = entry_d2one(v, function->hierarchie[step]->vektor[k]->lang); + vec_add2(v, function->hierarchie[step]->vektor[k]->start); // needed bc v starts from 0 + if (step == function->maxlevel || !is_in_latter_interval(v, function->hierarchie[step+1]->vektor[k])) { // only use the vector if it is in the valid interval + x = entry_one2d_sloppy(k, function->hierarchie[step]->grad); // don't care about 0 as entries + result += function->hierarchie[step]->vektor[k]->glied[position]*integrate_legendre(v, x, step, stepping); + vec_del(x); + } + vec_del(v); + } + } + } +} + +/** + * Returns the number of discrete elmenets that fit in the given folge. + * + * @param folge the folge whose elements should be count + * + * @return #folge + */ +int get_interval_element_count(folge_p folge) { + int n, result; + + for (n = 0, result = 1; n < folge->lang->dim; n++) { + result *= folge->lang->array[n]; + } + return result; +} + +/** + * This method generates the vector at the given position inside of the inverval that is defined in the given folge. + * + * @param folge the folge that defines the interval + * @param element_position position of the vector that is requested + * + * @return folge.interval(element_position) + */ +vec_p generate_folgenvector(folge_p folge, int element_position) { + vec_p result = vec_new(folge->lang->dim); + + int n, length; + + for (n = 0; n < result->dim; n++) { + length = folge->lang->array[n]; + result->array[n] = element_position % length; + element_position = (element_position - result->array[n]) / length; + //result->array[n] += folge->start->array[n]; that should be here but in that case the algorithm does not work --> use vec_add2 in add_folgenentries with the start-value + } + return result; +} + +vec_real_p generate_x_for_integration(vec_p v, int calc_position, fepc_real_t h, fepc_real_t h_l, int stepcount) { + vec_real_p result = vec_real_new(v->dim); + + int n, current; + + for (n = 0; n < v->dim; n++) { + current = calc_position % stepcount; + calc_position = (calc_position - current) / stepcount; + result->array[n] = v->array[n]*h_l+current*h; + } + return result; +} + +bool_t has_value_until(int dimension, int calc_position, int stepcount, int start, int end) { + int n, current; + + for (n = 0; n < dimension; n++) { + current = calc_position % stepcount; + if (current >= start && current <= end) { + return true; + } + calc_position = (calc_position - current) / stepcount; + } + return false; +} + +int get_current_changing_dimension(int dimension, int calc_position, int stepcount) { + int n, current; + + for (n = 0; n < dimension; n++) { + current = calc_position % stepcount; + if (current != calc_position -1 % stepcount) { + return n; + } + calc_position = (calc_position - current) / stepcount; + } + return 0; + +} + +fepc_real_t integrate_coeff(Funcimpl function_impl, vec_p v, vec_p p, int step, fepc_real_t stepping) { + int i, n, calc_count, position; + + n = INT_STEPS; + + fepc_real_t h, h_l, phi, f_x, result, k; + + bool_t second_sum = false; + h_l = get_h_l(step, stepping); + + + h = h_l/n; // step size + + result = 0.; + + + // interval length is always h_l! (by design) + + calc_count = pow(n+1, v->dim); + + vec_real_p x; + #ifdef __DEBUG + printf(" calc_stepcount = %i\n", calc_count); + #endif + + for (i = 0; i < calc_count; i++) { + x = generate_x_for_integration(v, i, h, h_l, n+1); + f_x = function_impl(x); + phi = phi_l(step, v, p, x, stepping); + if (has_value_until(v->dim, i, n+1, 0, 0)) { // first component + result += 0.5*f_x*phi; + } + if (has_value_until(v->dim, i, n+1, n, n)) { // last component + result += 0.5*f_x*phi; + } + if (has_value_until(v->dim, i, n+1, 1, n-1)) { // first sum + result += f_x*phi; + second_sum = true; + } + if (second_sum || has_value_until(v->dim, i, n+1, 1, n)) { // second sum + position = i % v->dim; + k = (x->array[position] - v->array[n]*h_l)/h; + x->array[position] = (x->array[position] + v->array[n]*h_l+(k-1)*h)/2.;// only one coordinate has to be changed + result += 2*function_impl(x)*phi_l(step, v, p, x, stepping); + vec_real_del(x); + } + vec_real_del(x); + second_sum = false; + } + #ifdef __DEBUG + printf("sum = %f, result = %f", result, pow(h_l, v->dim)*result / (2.*calc_count)); + #endif + return pow(h/3, v->dim)*result; +} + +/** + * Generates the real-valued vector at the calc_position. + * + * @param v vector that descripes the interval that is currently used + * @param calc_position position in the interval of which the vector should be returned + * @param h step size + * @param h_l h_l + * @param stepcount total number of interval elements + * + * @return vector at calc_position + */ +vec_real_p generate_x_for_integration_st(vec_p v, int calc_position, fepc_real_t h, fepc_real_t h_l, int stepcount) { + vec_real_p result = vec_real_new(v->dim); + + int n, current; + + for (n = 0; n < v->dim; n++) { + current = calc_position % stepcount; + calc_position = (calc_position - current) / stepcount; + result->array[n] = v->array[n]*h_l+(2.0*current +1)*h/2.0; + } + return result; +} + +/** + * Returns the number of start or end values of the given vector in the calculated interval. + * @param v vector that descripes the interval that is currently used + * @param calc_position position in the interval of which the vector should be returned + * @param h step size + * @param h_l h_l + * @param stepcount total number of interval elements + * + * @return startvalues + endvalues + */ +int has_start_or_end_value(vec_p v, int calc_position, fepc_real_t h, fepc_real_t h_l, int stepcount) { + int n, current, result; + + result = 0; + for (n = 0; n < v->dim; n++) { + current = calc_position % stepcount; + + calc_position = (calc_position - current) / stepcount; + if (current == 0 || current == stepcount-1) { + //result = 1 + result++; + } + } + //printf("%i\n", result); + return result; +} + +/** + * Integrates the given function via trapezoidal rule. Uses multiple cores where available! + * + * @param function_impl function to integrate + * @param v vector that defines the integration-bounds (v[n]*h_l till (v[n]+1)*h_l) + * @param p current polynomial degree + * @param step step + * + * @return integration result + */ +fepc_real_t integrate_coeff_st(Funcimpl function_impl, vec_p v, vec_p p, int step, fepc_real_t stepping) { + // uses Sehnen-Trapez + int n, calc_count; + + fepc_real_t h, result, h_l, temp; + h_l = get_h_l(step, stepping); + + h = h_l/(INT_STEPS); // step size + result = 0.0; + // interval length is always h_l! (by design) + + calc_count = pow(INT_STEPS, v->dim); + + vec_real_p x; + + //printf(" calc_stepcount = %i\n", calc_count); + + + #pragma omp parallel for private(temp) reduction(+:result) schedule(static,1) + for (n = 0; n < calc_count; n++) { + x = generate_x_for_integration_st(v, n, h, h_l, INT_STEPS); + //printf("%f\n", phi_l(step, v, p, x, stepping)); + temp = function_impl(x)*phi_l(step, v, p, x, stepping)*pow(.5,has_start_or_end_value(v, n, h, h_l, INT_STEPS)); + result += temp; + vec_real_del(x); + } + + return pow(h_l, v->dim)*result / pow(INT_STEPS-1, v->dim); +} + + +/** + * Returns true if the given vector is in the latter interval that is defined by the given folge. + * + * @param v vector whose position is going to be checked + * @param folge defines the interval + * + * @return v in folge.interval + */ +bool_t is_in_latter_interval(vec_p v, folge_p folge) { + int n; + + for (n = 0; n < v->dim; n++) { + if (v->array[n] < folge->start->array[n] / 2. || v->array[n] >= (folge->start->array[n]+folge->lang->array[n]) /2.) { + return false; + } + } + + return true; +} + +/** + * Returns true if the given vector is in the latter interval that is defined by the given folge. + * + * @param v vector whose position is going to be checked + * @param folge defines the interval + * + * @return v in folge.interval + * + * ACTS DIFFERENTLY TO is_in_latter_interval ! + */ +bool_t is_in_latter_interval_real(vec_real_p v, folge_p folge, int step, fepc_real_t stepping) { + int n; + + fepc_real_t h_l = get_h_l(step, stepping); + for (n = 0; n < v->dim; n++) {//printf("Start = %i Laenge = %i\n",folge->start->array[n], folge->lang->array[n]); + if (v->array[n] < h_l*folge->start->array[n] || v->array[n] >= h_l*(folge->start->array[n]+folge->lang->array[n])) { + return false; + } + } + return true; +} + +/** + * Returns the number of degree-vectors that are possible in [0, degree). + * + * @param degree vector that contains the component-degrees + * + * @return |[0, degree)| + */ +int get_degree_count(vec_p degree) { + int result, n; + + result = 1; + for (n = 0; n < degree->dim; n++) { + result *= degree->array[n]+1; + } + return result; +} + +/** + * Sets all entries of the hp structure of the given function. + * + * @param function the function whose entries are goin to be set + * @param function_impl R^n->R function from the convolution; maybe NULL if coeff_function is specified + * @param coeff_function function that calculate the hp-structure as a closed expression without numerical integration; maybe NULL if function_impl is specified; preferred + */ +void add_folgenentries(func_p function, Funcimpl function_impl, Funcimpl_step coeff_function, fepc_real_t stepping) { + ASSERT(function_impl != NULL || coeff_function != NULL); // we need one of the two functions to generate the folge + + int step, n, k, interval_element_count, position, grad_count; + + vec_p v, p; + + for (step = 0; step <= function->maxlevel; step++) { + grad_count = get_degree_count(function->hierarchie[step]->grad); + for (k = 0; k < grad_count; k++) { + for (n = 0, interval_element_count = get_interval_element_count(function->hierarchie[step]->vektor[k]); n < interval_element_count; n++) { + v = generate_folgenvector(function->hierarchie[step]->vektor[k], n); + position = entry_d2one(v, function->hierarchie[step]->vektor[k]->lang); + vec_add2(v, function->hierarchie[step]->vektor[k]->start); // needed bc v starts from 0 + if (step == function->maxlevel || !is_in_latter_interval(v, function->hierarchie[step+1]->vektor[k])) { // only set the vector if it is in the valid interval + p = entry_one2d_sloppy(k, function->hierarchie[step]->grad); // don't care about 0 as entries + // we prefer the coeff-function due to speed and acuracy: + function->hierarchie[step]->vektor[k]->glied[position] = coeff_function != NULL ? coeff_function(v, p, step, stepping) : integrate_coeff_st(function_impl, v, p, step, stepping); + vec_del(p); + } + vec_del(v); + } + } + } +} + +/** + * Sets the degree to every element in the given function. + * + * @param function the function where to set the degree + * @param degree degree + */ +void set_degree(func_p function, int degree) { + int n, dimension; + + dimension = function->dim; + vec_p p = vec_new(dimension); + + for (n = 0; n < dimension; n++) { + p->array[n] = degree; + } + for (n = 0; n <= function->maxlevel; n++) { + folgen_vektor_del(function->hierarchie[n]); + function->hierarchie[n] = folgen_vektor_new(vec_copy(p)); + } + vec_del(p); +} + + + +/** + * Sets up the needed structure from human readable intervals. + * + * @param func function that should be set up + * @param function R^n -> R function + * @param intervals human readable inverals + * @param interval_count number of intervals + * @param degree degree of the function (0 for piecewise constant) + * @param stepping stepping + * + */ +void +setup_fepc_structure(func_t * func, Funcimpl function, interval_p* intervals, int interval_count, int degree, fepc_real_t stepping) { + ASSERT(func->dim == intervals[interval_count-1]->dimension); + ASSERT(func->maxlevel == interval_count-1); + if (degree > 0) { + set_degree(func, degree); + } + + set_gridstructure(func, intervals, stepping); + if (function != NULL) { + add_folgenentries(func, function, NULL, stepping); + } +} + +/** + * Creates a functions and sets up the correct fepc structure. + * + * @param function R^n -> R function + * @param intervals human readable inverals + * @param interval_count number of intervals + * @param degree degree of the function (0 for piecewise constant) + * @param stepping stepping + * + * @return setted up function + */ + +func_t * +create_fepc_structure(Funcimpl function, interval_p* intervals, int interval_count, int degree, fepc_real_t stepping) { + func_t * result = func_new(interval_count-1, intervals[interval_count-1]->dimension); + + setup_fepc_structure(result, function, intervals, interval_count, degree, stepping); + return result; +} + +fepc_real_t max(fepc_real_t a, fepc_real_t b) { + return a > b ? a : b; +} + +bool_t isSuccessor(vec_p vector1, vec_p vector2, int direction, int delta) { + int n; + + for (n = 0; n < vector1->dim; n++) { + if (n == direction) { + if (vector2->array[n] - vector1->array[n] != delta) { + return false; + } + } else { + if (vector2->array[n] != vector1->array[n]) { + return false; + } + } + } + return true; +} + +vec_p getSuccessor(vec_p v, int position, folge_p folge, int direction, int element_count) { + int n; + + vec_p result; + + for (n = position+1; n < element_count; n++) { + result = generate_folgenvector(folge, n); + if (isSuccessor(v, result, direction, 1)) { + return result; + } else { + vec_del(result); + } + } + return NULL; +} + +vec_p getPredeccessor(vec_p v, int position, folge_p folge, int direction) { + int n; + + vec_p result; + + for (n = position-1; n >-1; n--) { + result = generate_folgenvector(folge, n); + if (isSuccessor(v, result, direction, -1)) { + return result; + } else { + vec_del(result); + } + } + return NULL; +} + +/** + * Derives a function in the given direction. + * + * @param function function to derive + * @param direction direction of the derivative + * @param stepping stepping + * + * @return derived function in the direction + */ +func_p func_derive(func_p function, int direction, fepc_real_t stepping) { + ASSERT(function->dim > direction); + + func_p result = func_new(function->maxlevel, function->dim); + + vec_p v, v2; + + int step, n, n2, i, grad_count, position, position2, interval_element_count, position3; + + for (step = 0; step <= function->maxlevel; step++) { + result->hierarchie[step]->vektor[0] = folge_new(vec_copy(function->hierarchie[step]->vektor[0]->start), vec_copy(function->hierarchie[step]->vektor[0]->lang)); + for (n = 0, interval_element_count = get_interval_element_count(function->hierarchie[step]->vektor[0]); n < interval_element_count; n++) { + v = generate_folgenvector(function->hierarchie[step]->vektor[0], n); + + v2 = getSuccessor(v, n, function->hierarchie[step]->vektor[0], direction, interval_element_count); + position = entry_d2one(v, function->hierarchie[step]->vektor[0]->lang); + if (v2 != NULL) { + position2 = entry_d2one(v2, function->hierarchie[step]->vektor[0]->lang); + result->hierarchie[step]->vektor[0]->glied[position] = (function->hierarchie[step]->vektor[0]->glied[position2] - function->hierarchie[step]->vektor[0]->glied[position])/get_h_l(0, stepping); + vec_del(v2); + vec_del(v); + //printf("%f\n", result->hierarchie[step]->vektor[0]->glied[position]); + } else { + + v2 = getPredeccessor(v, n, function->hierarchie[step]->vektor[0], direction); + vec_del(v); + if (v2 != NULL) { + v = getPredeccessor(v2, n, function->hierarchie[step]->vektor[0], direction); + position2 = entry_d2one(v2, function->hierarchie[step]->vektor[0]->lang); + position3 = entry_d2one(v, function->hierarchie[step]->vektor[0]->lang); + result->hierarchie[step]->vektor[0]->glied[position] = 2*result->hierarchie[step]->vektor[0]->glied[position2] - result->hierarchie[step]->vektor[0]->glied[position3]; + vec_del(v2); + vec_del(v); + } else { // this happens only, if stepping == real_interval_length[direction] + result->hierarchie[step]->vektor[0]->glied[position] = function->hierarchie[step]->vektor[0]->glied[position]; + } + //vec_del(v2); + + } + + + + } + + } + + return result; + +} + +/** + * Computes the laplacian of a given function. + * + * @param function function that will be 'laplaced' + * @param stepping stepping + * + * @return laplacian of the given function + */ +func_p func_laplace(func_p function, fepc_real_t stepping) { + func_p result, temp1, temp2; + + int n; + + for (n = 0; n < function->dim; n++) { + temp1 = func_derive(function, n, stepping); + temp2 = func_derive(temp1, n, stepping); + func_del(temp1); + if (n == 0) { + result = temp2; + } else { + func_add2(result, temp2); + func_del(temp2); + } + } + return result; +} + +/** + * This method simplifies the fecp convolution algorithm in a way that one does not need to take care of the structural expections the fepc algorithm has. + * + * @param function1 first R^n -> R function of the convolution + * @param intervals1 refined grid of the first function where the i-th grid has to be finer than the (i-1)-th grid + * @param interval_count1 number of the intervals + * @param degree_func1 degree of every component variable of the basis function that are associated to the first function + * @param function2 second R^n -> R function of the convolution + * @param intervals2 refined grid of the second function where the i-th grid has to be finer than the (i-1)-th grid + * @param interval_count2 number of the intervals + * @param degree_func2 degree of every component variable of the basis function that are associated to the second function + * @param intervals3 refined grid of the convolution result where the i-th grid has to be finer than the (i-1)-th grid + * @param interval_count3 number of the intervals + * @param stepping stepping of the first level refinement + * + * @return function1 * function2 in S(M) + */ +func_p fepc_easy(Funcimpl function1, interval_p* intervals1, int interval_count1, int degree_func1, Funcimpl function2, interval_p* intervals2, int interval_count2, int degree_func2, interval_p* intervals3, int interval_count3, fepc_real_t stepping) { + // check for coherence of the parameters + ASSERT(interval_count1 > 0 && interval_count2 > 0 && interval_count3 > 0 && degree_func1 > -1 && degree_func2 > -1); + ASSERT(function1 != NULL && intervals1 != NULL && function2 != NULL && intervals2 != NULL && intervals3 != NULL); + + func_p f, g, w, fepc; + + f = create_fepc_structure(function1, intervals1, interval_count1, degree_func1, stepping); + g = create_fepc_structure(function2, intervals2, interval_count2, degree_func2, stepping); + w = create_fepc_structure(NULL, intervals3, interval_count3, 0, stepping); + + fepc = faltung_fepc(f, g, w, stepping); + func_del(f); + func_del(g); + func_del(w); + return fepc; +} + + diff --git a/source/fepc_easy_helper.c b/source/fepc_easy_helper.c new file mode 100644 index 0000000..b7b0817 --- /dev/null +++ b/source/fepc_easy_helper.c @@ -0,0 +1,520 @@ +/* + * FEPC + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_easy_helper.h" +#include "fepc_easy.h" + + +func_p +func_multi(func_p f, func_p g, fepc_real_t stepping) { + func_p result = func_new(f->maxlevel > g->maxlevel ? f->maxlevel : g->maxlevel, f->dim); + + func_multi_overwrite(result, f, g, stepping); + return result; +} + +void +func_multi_overwrite(func_t * result, func_p f, func_p g, fepc_real_t stepping) { + ASSERT(f != NULL && g != NULL); + ASSERT(f->dim == g->dim); + + folgen_vektor_p temp; + + int n, k, l, m, steps; + + steps = (f->maxlevel > g->maxlevel ? f->maxlevel : g->maxlevel)+1; + + ASSERT(result->maxlevel == steps -1); + ASSERT(result->dim == f->dim); + + for (n = 0; n < steps; n++) { + temp = folgen_vektor_multi( f->maxlevel < n ? NULL : f->hierarchie[n], g->maxlevel < n ? NULL : g->hierarchie[n], n, stepping ); + folgen_vektor_del(result->hierarchie[n]); + result->hierarchie[n] = temp; + } +} + +func_p +func_div(func_p f, func_p g, fepc_real_t stepping) { + ASSERT(f != NULL && g != NULL); + ASSERT(f->dim == g->dim); + + folgen_vektor_p temp; + + int n, k, l, m, steps; + + steps = (f->maxlevel > g->maxlevel ? f->maxlevel : g->maxlevel)+1; + + func_p result = func_new(steps-1, f->dim); + + for (n = 0; n < steps; n++) { + temp = folgen_vektor_div( f->maxlevel < n ? NULL : f->hierarchie[n], g->maxlevel < n ? NULL : g->hierarchie[n], n, stepping ); + folgen_vektor_del(result->hierarchie[n]); + result->hierarchie[n] = temp; + } + return result; +} + + +folgen_vektor_p +folgen_vektor_multi(folgen_vektor_p f, folgen_vektor_p g, int step, fepc_real_t stepping) { + folgen_vektor_p back; + int k, d, size, dim, a, b, test_f, test_g; + vec_p max, vec_1, n, n_f, n_g, r; + + ASSERT( f->grad->dim == g->grad->dim ); + dim = f->grad->dim; + max = vec_max( f->grad, g->grad ); + + + vec_1 = vec_one( dim ); + n = vec_add( vec_1, max ); + n_f = vec_add( vec_1, f->grad ); + n_g = vec_add( vec_1, g->grad ); + size = vec_size( n ); + + back = folgen_vektor_new( max ); + for(k=0;karray[d] > f->grad->array[d] ) { + test_f = test_f + 1; + } + if( r->array[d] > g->grad->array[d] ) { + test_g = test_g + 1; + } + } + if( (test_f == 0) && (test_g == 0) ) { + a = entry_d2one( r, n_f ); + b = entry_d2one( r, n_g ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_multi( f->vektor[a], g->vektor[b], step, stepping ); + } + if( (test_f != 0) && (test_g == 0) ) { + b = entry_d2one( r, n_g ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_copy( g->vektor[b] ); + } + if( (test_f == 0) && (test_g != 0) ) { + a = entry_d2one( r, n_f ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_copy( f->vektor[a] ); + } + vec_del( r ); + } + vec_del( n ); + vec_del( n_f ); + vec_del( n_g ); + vec_del( vec_1 ); + + return back; +} + + +folgen_vektor_p +folgen_vektor_div(folgen_vektor_p f, folgen_vektor_p g, int step, fepc_real_t stepping) { + folgen_vektor_p back; + int k, d, size, dim, a, b, test_f, test_g; + vec_p max, vec_1, n, n_f, n_g, r; + + ASSERT( f->grad->dim == g->grad->dim ); + dim = f->grad->dim; + max = vec_max( f->grad, g->grad ); + + + vec_1 = vec_one( dim ); + n = vec_add( vec_1, max ); + n_f = vec_add( vec_1, f->grad ); + n_g = vec_add( vec_1, g->grad ); + size = vec_size( n ); + + back = folgen_vektor_new( max ); + for(k=0;karray[d] > f->grad->array[d] ) { + test_f = test_f + 1; + } + if( r->array[d] > g->grad->array[d] ) { + test_g = test_g + 1; + } + } + if( (test_f == 0) && (test_g == 0) ) { + a = entry_d2one( r, n_f ); + b = entry_d2one( r, n_g ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_div( f->vektor[a], g->vektor[b], step, stepping ); + } + if( (test_f != 0) && (test_g == 0) ) { + b = entry_d2one( r, n_g ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_copy( g->vektor[b] ); + } + if( (test_f == 0) && (test_g != 0) ) { + a = entry_d2one( r, n_f ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_copy( f->vektor[a] ); + } + vec_del( r ); + } + vec_del( n ); + vec_del( n_f ); + vec_del( n_g ); + vec_del( vec_1 ); + + return back; +} + + +fepc_real_t folge_norm2(folge_p folge, int step, fepc_real_t stepping) { + int n, interval_element_count, position, k; + + vec_p v, p; + + vec_real_p x = vec_real_new(folge->start->dim); + + p = vec_new(folge->start->dim); + + fepc_real_t result, h_l; + + result = 0.; + h_l = get_h_l(step, stepping); + + for (n = 0, interval_element_count = get_interval_element_count(folge); n < interval_element_count; n++) { + v = generate_folgenvector(folge, n); + position = entry_d2one(v, folge->lang); + vec_add2(v, folge->start); + for (k = 0; k < v->dim; k++) { + x->array[k] = v->array[k]*h_l; // dummy values because only p = 0 is supported + } + result = phi_l(step, v, p, x, stepping); + //printf("Result %f\n", result); + vec_del(v); + } + vec_del(p); + vec_real_del(x); + return result; +} + + +folge_p +folge_multi(folge_p f, folge_p g, int step, fepc_real_t stepping) { + folge_p back; + int size, k, size_g, size_f; + fepc_real_t x, y; + vec_p temp1, temp2, max, lang, min, r; + + + size_f = vec_size( f->lang ); + size_g = vec_size( g->lang ); + if(size_f == 0) { + back = folge_copy( g ); + return back; + } + + if(size_g == 0) { + back = folge_copy( f ); + return back; + } + + fepc_real_t norm = folge_norm2(g, step, stepping); + if( (size_g!=0) && (size_f!=0) ) { + min = vec_min( f->start, f->start ); + temp1 = vec_add( f->start, f->lang ); + temp2 = vec_add( f->start, f->lang ); + max = vec_max( temp1, temp2 ); + vec_del( temp1 ); + vec_del( temp2 ); + lang = vec_op( 1, max, -1, min); + vec_del( max ); + back = folge_new( min, lang ); + size = vec_size( lang ); + for(k=0;kglied[k] = x*y*norm; + //printf("%i\t%f\t%f\n", step, x, y); + } + return back; + } +} + +folge_p +folge_div(folge_p f, folge_p g, int step, fepc_real_t stepping) { + folge_p back; + int size, k, size_g, size_f; + fepc_real_t x, y; + vec_p temp1, temp2, max, lang, min, r; + + + size_f = vec_size( f->lang ); + size_g = vec_size( g->lang ); + if(size_f == 0) { + back = folge_copy( g ); + return back; + } + + if(size_g == 0) { + back = folge_copy( f ); + return back; + } + + fepc_real_t norm = folge_norm2(g, step, stepping); + + if( (size_g!=0) && (size_f!=0) ) { + min = vec_min( f->start, g->start ); + temp1 = vec_add( f->start, f->lang ); + temp2 = vec_add( f->start, f->lang ); + max = vec_max( temp1, temp2 ); + vec_del( temp1 ); + vec_del( temp2 ); + lang = vec_op( 1, max, -1, min); + vec_del( max ); + back = folge_new( min, lang ); + size = vec_size( lang ); + for(k=0;kglied[k] = x == 0. ? 0. : x/(y*norm); + } + return back; + } +} + + +vec_real_p get_mean_points(vec_p v, vec_p grad, int step, fepc_real_t stepping) { + fepc_real_t h_l = get_h_l(step, stepping); + + vec_real_p result = vec_real_new(v->dim); + + int n; + + for (n = 0; n < v->dim; n++) { + result->array[n] = (v->array[n] + .5)*h_l; + } + return result; +} + +folgen_vektor_p +folgen_vektor_reflect(folgen_vektor_p f, int step) { + folgen_vektor_p result = folgen_vektor_copy_structure(f); + + fepc_real_t temp; + + vec_p r; + + int n, k, degree_count, interval_element_count, position; + + degree_count = get_degree_count(f->grad); + interval_element_count = get_interval_element_count(f->vektor[0]); // since they have all the same element-count + + for (n = 0; n < interval_element_count / 2; n++) { + for (k = 0; k < degree_count; k++) { + r = generate_folgenvector(f->vektor[k], n); + position = entry_d2one(r, f->vektor[k]->lang); + result->vektor[k]->glied[position] = f->vektor[k]->glied[interval_element_count - position - 1]; + result->vektor[k]->glied[interval_element_count - position - 1] = f->vektor[k]->glied[position]; + vec_del(r); + } + } + return result; +} + + +folgen_vektor_p +folgen_vektor_multi_factor(folgen_vektor_p f, fepc_real_t factor, int step) { + folgen_vektor_p result = folgen_vektor_copy_structure(f); + + + vec_p r; + + int n, k, degree_count, interval_element_count, position; + + degree_count = get_degree_count(f->grad); + interval_element_count = get_interval_element_count(f->vektor[0]); // since they have all the same element-count + + for (n = 0; n < interval_element_count; n++) { + for (k = 0; k < degree_count; k++) { + r = generate_folgenvector(f->vektor[k], n); + position = entry_d2one(r, f->vektor[k]->lang); + result->vektor[k]->glied[position] = factor*f->vektor[k]->glied[position]; // TODO maybe this factor needs to be adjusted + vec_del(r); + } + } + return result; +} + +folgen_vektor_p +folgen_vektor_reflect_x(folgen_vektor_p f, int step) { + return folgen_vektor_multi_factor(f, -1., step); +} + + + + +func_p +func_reflect(func_p f) { + func_p result = func_new(f->maxlevel, f->dim); + + int n; + + folgen_vektor_p temp; + + for (n = 0; n <= result->maxlevel; n++) { + temp = folgen_vektor_reflect(f->hierarchie[n], n); + folgen_vektor_del(result->hierarchie[n]); + result->hierarchie[n] = temp; + } + return result; +} + +func_p +func_reflect_x(func_p f) { + func_p result = func_new(f->maxlevel, f->dim); + + int n; + + folgen_vektor_p temp; + + for (n = 0; n <= result->maxlevel; n++) { + temp = folgen_vektor_reflect_x(f->hierarchie[n], n); + folgen_vektor_del(result->hierarchie[n]); + result->hierarchie[n] = temp; + } + return result; +} + +void +func_add2(func_p f, func_p g) { + int k, maxlevel, dim; + folgen_vektor_p temp; + + ASSERT( f->dim == g->dim ); + + if (f->maxlevel > g->maxlevel) { + maxlevel = g->maxlevel; + } + else { + maxlevel = f->maxlevel; + } + + for (k=0;k<=maxlevel;k++) { + temp = folgen_vektor_add( f->hierarchie[k], g->hierarchie[k] ); + folgen_vektor_del( f->hierarchie[k] ); + f->hierarchie[k] = temp; + } +} + +void +func_add3(func_p f, fepc_real_t factor, func_p g) { + int k, maxlevel, dim; + folgen_vektor_p temp, temp2; + + ASSERT( f->dim == g->dim ); + + if (f->maxlevel > g->maxlevel) { + maxlevel = g->maxlevel; + } + else { + maxlevel = f->maxlevel; + } + + for (k=0;k<=maxlevel;k++) { + temp = folgen_vektor_multi_factor(g->hierarchie[k], factor, k); + temp2 = folgen_vektor_add( f->hierarchie[k], temp); + folgen_vektor_del(temp); + folgen_vektor_del( f->hierarchie[k] ); + f->hierarchie[k] = temp2; + } +} + +func_p +func_subtract(func_p f, func_p g) { + func_p back; + int k, maxlevel, dim; + folgen_vektor_p temp, temp2; + + ASSERT( f->dim == g->dim ); + dim = f->dim; + if (f->maxlevel > g->maxlevel) { + maxlevel = g->maxlevel; + } + else { + maxlevel = f->maxlevel; + } + + + back = func_new( maxlevel, dim ); + for (k=0;k<=maxlevel;k++) { + temp = folgen_vektor_multi_factor(g->hierarchie[k], -1., k); + temp2 = folgen_vektor_add( f->hierarchie[k], temp); + folgen_vektor_del(temp); + folgen_vektor_del( f->hierarchie[k] ); + back->hierarchie[k] = temp2; + } + + return back; +} + + +void +func_subtract2(func_p f, func_p g) { + int k, maxlevel, dim; + folgen_vektor_p temp; + + ASSERT( f->dim == g->dim ); + + if (f->maxlevel > g->maxlevel) { + maxlevel = g->maxlevel; + } + else { + maxlevel = f->maxlevel; + } + + for (k=0;k<=maxlevel;k++) { + temp = folgen_vektor_subtract( f->hierarchie[k], g->hierarchie[k] ); + folgen_vektor_del( f->hierarchie[k] ); + f->hierarchie[k] = temp; + } +} + +fepc_real_t func_norm_l2_sqr(func_p function, fepc_real_t stepping) { + func_p temp = func_multi(function, function, stepping); + + fepc_real_t result = func_integrate(temp, stepping); + func_del(temp); + + return result; +} + + + diff --git a/source/fepc_easy_sparse.c b/source/fepc_easy_sparse.c new file mode 100644 index 0000000..5bb43ba --- /dev/null +++ b/source/fepc_easy_sparse.c @@ -0,0 +1,97 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_easy_sparse.h" +#include "fepc_easy_helper.h" + +func_sparse_p +func_sparse_new(int summands, int dimension) { + func_sparse_p result = (func_sparse_p) malloc(sizeof(func_sparse_t)); + + result->summands = summands; + result->dimension = dimension; + result->factors = (func_p*) malloc(sizeof(func_p)*summands*dimension); + return result; +} + + +void +func_sparse_del(func_sparse_p func_sparse) { + int n , i = func_sparse->summands*func_sparse->dimension; + + for (n = 0; n < i; n++) { + func_del(func_sparse->factors[n]); + } + free(func_sparse->factors); + free(func_sparse); +} + + +func_p +func_sparse_to_func(func_sparse_p func_sparse) { + func_p result = func_new(1, func_sparse->dimension); + + return result; +} + + +func_sparse_p +func_sparse_multi(func_sparse_p f, func_sparse_p g, fepc_real_t stepping) { + int n, i, k, dimension = f->dimension; + + func_sparse_p result = func_sparse_new(f->summands * g->summands, dimension); + + for (n = 0; n < f->summands; n++) { + for (i = 0; i < g->summands; i++) { + for (k = 0; k < dimension; k++) { + result->factors[(n*g->summands+i)*dimension+k] = func_multi(f->factors[n*dimension+k], g->factors[i*dimension+k], stepping); + } + } + } + return result; +} + +func_sparse_p +func_sparse_reflect(func_sparse_p f) { + int n, i; + + func_sparse_p result = func_sparse_new(f->summands, f->dimension); + + for (n = 0, i = f->summands*f->dimension; n < i; n++) { + result->factors[n] = func_reflect(f->factors[n]); + } + return result; +} + +func_sparse_p +func_sparse_convolute(func_sparse_p f, func_sparse_p g, func_p * result_structure, fepc_real_t stepping) { + int n, i, k, dimension = f->dimension; + + func_sparse_p result = func_sparse_new(f->summands * g->summands, dimension); + + for (n = 0; n < f->summands; n++) { + for (i = 0; i < g->summands; i++) { + for (k = 0; k < dimension; k++) { + result->factors[(n*g->summands+i)*dimension+k] = faltung_fepc(f->factors[n*dimension+k], g->factors[i*dimension+k], result_structure[k], stepping); + } + } + } + return result; +} + + diff --git a/fft_faltung.c b/source/fft_faltung.c similarity index 77% rename from fft_faltung.c rename to source/fft_faltung.c index 5d48651..46e2cc1 100644 --- a/fft_faltung.c +++ b/source/fft_faltung.c @@ -1,6 +1,7 @@ /* * FEPC * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * 2011 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -16,11 +17,7 @@ * along with this program. If not, see . */ -#if defined(HAS_FFTW3) #include - -#endif - #include "fft_faltung.h" @@ -41,14 +38,14 @@ *******************************************************/ fepc_real_t* fft_faltung(fepc_real_t* a, vec_p n_a, fepc_real_t* b, vec_p n_b) { -#if defined(HAS_FFTW3) - int size_a, size_b, size_c, dim; int k, i, wert, test; int *n; vec_p temp, n_c; fepc_real_t *c; - fftw_complex *in, *out,*in_a, *out_a, *in_b, *out_b; + fftw_complex *in, *out_a, *out_b; + double *out, *in_a, *in_b; + fftw_plan p; @@ -72,7 +69,7 @@ fft_faltung(fepc_real_t* a, vec_p n_a, fepc_real_t* b, vec_p n_b) { /*Berechnen der Fouriertrafo von in_a*/ - in_a = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size_c); + in_a = (double*) fftw_malloc(sizeof(double) * size_c); out_a = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size_c); for (k=0;kglied = w_glied; return back; +#else + return folge_slow_faltung(f, g); +#endif } @@ -243,7 +251,7 @@ folge_print(folge_p f, int info) { if(info == 1) { printf("\n\t#:glied"); for(k=0;kglied[k]); + printf("\t%.5lf",f->glied[k]); } } printf("\n------------------------------------------------------------\n"); @@ -408,6 +416,51 @@ folge_add(folge_p f, folge_p g) { } } +folge_p +folge_subtract(folge_p f, folge_p g) { + folge_p back; + int size, k, size_g, size_f; + fepc_real_t x, y; + vec_p temp1, temp2, max, lang, min, r; + + + size_f = vec_size( f->lang ); + size_g = vec_size( g->lang ); + if(size_f == 0) { + back = folge_copy( g ); + folge_multi_factor(back, -1); + return back; + } + + if(size_g == 0) { + back = folge_copy( f ); + return back; + } + + if( (size_g!=0) && (size_f!=0) ) { + min = vec_min( f->start, g->start ); + temp1 = vec_add( f->start, f->lang ); + temp2 = vec_add( g->start, g->lang ); + max = vec_max( temp1, temp2 ); + vec_del( temp1 ); + vec_del( temp2 ); + lang = vec_op( 1, max, -1, min); + vec_del( max ); + back = folge_new( min, lang ); + size = vec_size( lang ); + for(k=0;kglied[k] = x - y; + } + return back; + } +} + folge_p folge_copy( folge_p f) { @@ -450,3 +503,27 @@ folge_projekt(folge_p f, folge_p g) { return back; } + +folge_p +folge_multi_factor(folge_p folge, fepc_real_t factor) { + folge_p back; + vec_p lang, start; + int k, size; + + start = vec_copy( folge->start ); + lang = vec_copy( folge->lang ); + + back = folge_new( start, lang ); + size = vec_size( lang ); + + for(k=0;kglied[k] = folge->glied[k]*factor; + } + + return back; +} + + + + diff --git a/folgen_vektor.c b/source/folgen_vektor.c similarity index 87% rename from folgen_vektor.c rename to source/folgen_vektor.c index 6fb0a35..8ef6660 100644 --- a/folgen_vektor.c +++ b/source/folgen_vektor.c @@ -95,8 +95,6 @@ folgen_vektor_del(folgen_vektor_p f) { - - folgen_matrix_p folgen_matrix_new(vec_p grad1, vec_p grad2) { folgen_matrix_p back; @@ -514,6 +512,87 @@ folgen_vektor_add(folgen_vektor_p f, folgen_vektor_p g) { return back; } +folgen_vektor_p +folgen_vektor_subtract(folgen_vektor_p f, folgen_vektor_p g) { + folgen_vektor_p back; + int k, d, size, dim, a, b, test_f, test_g; + vec_p max, vec_1, n, n_f, n_g, r; + + ASSERT( f->grad->dim == g->grad->dim ); + dim = f->grad->dim; + max = vec_max( f->grad, g->grad ); + + + vec_1 = vec_one( dim ); + n = vec_add( vec_1, max ); + n_f = vec_add( vec_1, f->grad ); + n_g = vec_add( vec_1, g->grad ); + size = vec_size( n ); + + back = folgen_vektor_new( max ); + + folge_p temp; + + for(k=0;karray[d] > f->grad->array[d] ) { + test_f = test_f + 1; + } + if( r->array[d] > g->grad->array[d] ) { + test_g = test_g + 1; + } + } + if( (test_f == 0) && (test_g == 0) ) { + a = entry_d2one( r, n_f ); + b = entry_d2one( r, n_g ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_subtract( f->vektor[a], g->vektor[b] ); + } + if( (test_f != 0) && (test_g == 0) ) { + b = entry_d2one( r, n_g ); + folge_del( back->vektor[k] ); + temp = folge_copy( g->vektor[b] ); + back->vektor[k] = folge_multi_factor(temp, -1.); + folge_del(temp); + } + if( (test_f == 0) && (test_g != 0) ) { + a = entry_d2one( r, n_f ); + folge_del( back->vektor[k] ); + back->vektor[k] = folge_copy( f->vektor[a] ); + } + vec_del( r ); + } + vec_del( n ); + vec_del( n_f ); + vec_del( n_g ); + vec_del( vec_1 ); + + return back; +} + +folgen_vektor_p +folgen_vektor_factor_multi(folgen_vektor_p f, fepc_real_t factor) { + folgen_vektor_p back; + int k, size; + vec_p vec_1, n; + + vec_1 = vec_one( f->grad->dim ); + n = vec_add( vec_1, f->grad ); + vec_del(vec_1); + size = vec_size( n ); + vec_del(n); + back = folgen_vektor_new( f->grad ); + + + for(k=0;kvektor[k] = folge_multi_factor(f->vektor[k], factor); + } + return back; + +} @@ -607,3 +686,26 @@ folgen_matrix_add(folgen_matrix_p f, folgen_matrix_p g) { return back; } + + + + + + + + + + + + + + + + + + + + + + + diff --git a/funktion.c b/source/funktion.c similarity index 79% rename from funktion.c rename to source/funktion.c index 713a5d8..f29e3d4 100644 --- a/funktion.c +++ b/source/funktion.c @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de) + * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2010 Stefan Handschuh (handschu@mis.mpg.de) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,15 +38,20 @@ func_p func_new(int maxlevel, int dim) { - func_p back; + func_p back = (func_p) malloc(sizeof(func_t)); + func_init(back, maxlevel, dim); + return back; +} + +void +func_init(func_t * function, int maxlevel, int dim) { folgen_vektor_p *hierarchie; int k; vec_p grad; ASSERT(maxlevel >= 0); ASSERT(dim > 0); - back = (func_p) malloc(sizeof(func_t)); - ASSERT(back != NULL); + ASSERT(function != NULL); hierarchie = (folgen_vektor_p*) malloc(sizeof(folgen_vektor_p)*(maxlevel+1)); ASSERT(hierarchie != NULL); @@ -56,23 +61,27 @@ func_new(int maxlevel, int dim) { grad = vec_new(dim); hierarchie[k] = folgen_vektor_new(grad); } - back->hierarchie = hierarchie; - back->maxlevel = maxlevel; - back->dim = dim; - return back; + function->hierarchie = hierarchie; + function->maxlevel = maxlevel; + function->dim = dim; } - - void func_del(func_p f) { + if (f && f != NULL) { + func_clear(f); + free(f); + } +} + +void +func_clear(func_t * f) { int k; for(k=0;k<=f->maxlevel;k++) { folgen_vektor_del(f->hierarchie[k]); } free(f->hierarchie); - free(f); } @@ -140,6 +149,10 @@ func_projekt(func_p f,func_p g) { return back; } +func_p +func_clone(func_p f) { + return func_projekt(f, f); +} @@ -189,41 +202,40 @@ func_build( int maxlevel, int dim , int grad, int a, int n, int mod, bool_t rand } - - func_p func_add(func_p f, func_p g) { - func_p back; - int k, maxlevel, dim; + ASSERT(f->dim == g->dim); + + func_p back = func_new( f->maxlevel < g->maxlevel ? g->maxlevel : f->maxlevel, f->dim); + func_add_overwrite(back, f, g); + return back; +} + +void +func_add_overwrite(func_t * result, func_t * f, func_t * g) { + int k, maxlevel; + folgen_vektor_p temp; ASSERT( f->dim == g->dim ); - dim = f->dim; - if (f->maxlevel < g->maxlevel) { - maxlevel = g->maxlevel; - } - else { - maxlevel = f->maxlevel; - } - + maxlevel = f->maxlevel < g->maxlevel ? g->maxlevel : f->maxlevel; + ASSERT(maxlevel == result->maxlevel); - back = func_new( maxlevel, dim ); for (k=0;k<=maxlevel;k++) { if (k <= f->maxlevel) { - temp = folgen_vektor_add( back->hierarchie[k], f->hierarchie[k] ); - folgen_vektor_del( back->hierarchie[k] ); - back->hierarchie[k] = temp; + temp = folgen_vektor_add( result->hierarchie[k], f->hierarchie[k] ); + folgen_vektor_del( result->hierarchie[k] ); + result->hierarchie[k] = temp; } if (k <= g->maxlevel) { - temp = folgen_vektor_add( back->hierarchie[k], g->hierarchie[k] ); - folgen_vektor_del( back->hierarchie[k] ); - back->hierarchie[k] = temp; + temp = folgen_vektor_add( result->hierarchie[k], g->hierarchie[k] ); + folgen_vektor_del( result->hierarchie[k] ); + result->hierarchie[k] = temp; } } - - return back; } + int func_count( func_p f, func_p g, func_p w ) { int back; @@ -383,3 +395,40 @@ func_grid_zero(func_p f) { } vec_del(vec_1); } + +func_p +func_factor_multi(func_p function, fepc_real_t factor) { + func_p back = func_new(function->maxlevel, function->dim); + + folgen_vektor_p temp; + + int k; + + for (k=0;k<=function->maxlevel;k++) { + temp = folgen_vektor_factor_multi( function->hierarchie[k], factor ); + folgen_vektor_del( back->hierarchie[k] ); + back->hierarchie[k] = temp; + } + + return back; +} + +void +funcs_del(func_p * array, int length) { + int n; + + for (n = 0; n < length; n++) { + func_del(array[n]); + } + free(array); +} + +void +funcs_del_type(func_t * array, int length) { + int n; + + for (n = 0; n < length; n++) { + func_clear(&array[n]); + } + free(array); +} diff --git a/source/interval.c b/source/interval.c new file mode 100644 index 0000000..9951130 --- /dev/null +++ b/source/interval.c @@ -0,0 +1,112 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "interval.h" + + +interval_p +interval_new(int dimension) { + interval_p result = (interval_p) malloc(sizeof(interval_t)); + + interval_init(result, dimension); + return result; +} + +void +interval_init(interval_p interval, int dimension) { + interval->dimension = dimension; + interval->start = (fepc_real_t*) malloc(sizeof(fepc_real_t)*dimension); + interval->end = (fepc_real_t*) malloc(sizeof(fepc_real_t)*dimension); + + int n; + + for (n = 0; n < dimension; n++) { + interval->start[n] = 0.; + } +} + +void interval_del(interval_p interval) { + interval_clear(interval); + free(interval); +} + +void +interval_clear(interval_t * interval) { + free(interval->start); + free(interval->end); +} + +void intervals_del(interval_p* intervals, int interval_count) { + int n; + + for (n = 0; n < interval_count; n++) { + interval_del(intervals[n]); + } + free(intervals); +} + +void intervals_del_type(interval_t * intervals, int interval_count) { + int n; + + for (n = 0; n < interval_count; n++) { + interval_clear(&intervals[n]); + } + free(intervals); +} + +void print_interval(interval_p interval) { + int n; + + for (n = 0; n < interval->dimension; n++) { + printf("%f\t<\t%f\n", interval->start[n], interval->end[n]); + } +} + +void print_intervals(interval_p * intervals, int count) { + int n; + + for (n = 0; n < count; n++) { + printf("\nInterval %i\n", n+1); + print_interval(intervals[n]); + } +} + +interval_p interval_clone(interval_p interval) { + int n, dimension; + + dimension = interval->dimension; + interval_p result = interval_new(dimension); + for (n = 0; n < dimension; n++) { + result->start[n] = interval->start[n]; + result->end[n] = interval->end[n]; + } + return result; +} + +interval_p * intervals_clone(interval_p* intervals, int count) { + interval_p* result = (interval_p*) malloc(sizeof(interval_p)*count); + + int n; + + for (n = 0; n < count; n++) { + result[n] = interval_clone(intervals[n]); + } + return result; +} + + diff --git a/koeffizienten.c b/source/koeffizienten.c similarity index 99% rename from koeffizienten.c rename to source/koeffizienten.c index 72e4327..622d922 100644 --- a/koeffizienten.c +++ b/source/koeffizienten.c @@ -392,3 +392,8 @@ matrix3_del(matrix3_p matrix) { free(matrix->a); free(matrix); } + + + + + diff --git a/kuerzen.c b/source/kuerzen.c similarity index 99% rename from kuerzen.c rename to source/kuerzen.c index e74b839..4b3832d 100644 --- a/kuerzen.c +++ b/source/kuerzen.c @@ -605,3 +605,6 @@ kuerzen_F_bauen_c2( func_p f, func_p g, func_p w, func2_p Gamma, matrix_p xi_koe return back; } + + + diff --git a/main.c b/source/main.c similarity index 100% rename from main.c rename to source/main.c diff --git a/seconds.c b/source/seconds.c similarity index 99% rename from seconds.c rename to source/seconds.c index 75fef72..3d4bbd7 100644 --- a/seconds.c +++ b/source/seconds.c @@ -36,3 +36,4 @@ seconds () return sec + (usec * 1e-6); } + diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt new file mode 100644 index 0000000..b935582 --- /dev/null +++ b/test/CMakeLists.txt @@ -0,0 +1,29 @@ +add_executable(fepc_easy_test fepc_easy_test.c) +target_link_libraries(fepc_easy_test fepc m) + +add_executable(stest stest.c) +target_link_libraries(stest fepc m) + +add_executable(discont_test discont_test.c) +target_link_libraries(discont_test fepc m) + +add_executable(sparsetest sparsetest.c) +target_link_libraries(sparsetest fepc m) + +add_executable(fepc_cp_test fepc_cp_test.c) +target_link_libraries(fepc_cp_test fepc m) + +if (HAS_FFTW3) + target_link_libraries(fepc_easy_test fftw3) + target_link_libraries(discont_test fftw3) + target_link_libraries(sparsetest fftw3) + target_link_libraries(fepc_cp_test fftw3) + target_link_libraries(stest fftw3) +endif(HAS_FFTW3) + +add_test(fepc_easy_test fepc_easy_test) +add_test(stest stest) +add_test(discont_test discont_test) +add_test(sparsetest sparsetest) +add_test(fepc_cp_test fepc_cp_test) + diff --git a/test/discont_test.c b/test/discont_test.c new file mode 100644 index 0000000..9b39088 --- /dev/null +++ b/test/discont_test.c @@ -0,0 +1,124 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "discont.h" +#include "fepc_easy_helper.h" + + +fepc_real_t x(vec_real_p vec) { + return vec->array[0]+1.0; +} + +int +main(void) { + /* Example */ + + fepc_real_t stepping = 0.1; + + int steps, n; + + func_p f1, f2, w, convolution_result; + + steps = 1; + + discont_function_p function, result; + + function = discont_function_new(steps); + + fepc_real_t y11[] = {1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9}; // the stepcount is 10 in the first level because h_0 is 0.1 + fepc_real_t y21[] = {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0}; + + + // first step + discont_function_setup_points(function, 0, 0.0, 1.0, y11, y21, stepping); // interval [0, 1] at level 0 + + fepc_real_t y12[] = {0.1, 2.6, 0.3, 4.7, 0.0, 7.5}; // the stepcount is 6 in the second level because h_1 is 0.05 + fepc_real_t y22[] = {9.4, 2.5, 3.0, 1.9, 4.2, 0.5}; + + // second step + //discont_function_setup_points(function, 1, 0.2, 0.5, y12, y22, stepping); // interval [0.2, 0.5] at level 1 + + discont_function_print(function); + + f1 = convert_discont_function(function, stepping); + + + func_print(f1, 3); + func_print(create_fepc_structure(x, function->intervals, 1, 1, stepping), 3); + + vec_real_set_p points = get_value(create_fepc_structure(x, function->intervals, 1, 1, stepping), get_mean_points, stepping); + print_vec_real_set(points); + + discont_function_del(function); + + // Create f2 in the same way! + + function = discont_function_new(steps); + + fepc_real_t y13[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + fepc_real_t y23[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + + // first step + discont_function_setup_points(function, 0, 0.0, 1.0, y13, y23, stepping); // interval [0, 1] at level 0 + + fepc_real_t y14[] = {0.1, 2.6, 0.3, 4.7, 0.0, 7.5}; // the stepcount is 6 in the second level because h_1 is 0.05 + fepc_real_t y24[] = {9.4, 2.5, 3.0, 1.9, 4.2, 0.5}; + + // second step + //discont_function_setup_points(function, 1, 0.2, 0.5, y14, y24, stepping); // interval [0.2, 0.5] at level 1 + + discont_function_print(function); + + f2 = convert_discont_function(function, stepping); + discont_function_del(function); + + func_print(f2, 3); + + // set up the structure w of the result + + function = discont_function_new(steps); + + discont_function_setup_points(function, 0, 0.0, 1.0, NULL, NULL, stepping); // interval [0, 1] at level 0 + + //discont_function_setup_points(function, 1, 0.3, 0.6, NULL, NULL, stepping); // interval [0.3, 0.6] at level 1 + + w = func_new(steps-1, 1); + + set_degree(w, 1); + set_gridstructure(w, function->intervals, stepping); + + convolution_result = faltung_fepc(f1, f2, w, stepping); + //func_del(f1); + //func_del(f2); + //func_del(w); + + func_print(convolution_result, 3); + + result = convert_func(convolution_result, intervals_clone(function->intervals, function->stepcount), stepping); + + //discont_function_del(function); + // points = get_value(convolution_result, get_mean_points, stepping); + // print_vec_real_set(points); + + discont_function_print(convert_func(f1, intervals_clone(function->intervals, function->stepcount), stepping)); + func_del(convolution_result); + + // the result is stored in "result" as a discont-function + discont_function_print(result); + return 0; +} diff --git a/test/fepc_cp_test.c b/test/fepc_cp_test.c new file mode 100644 index 0000000..543e3cf --- /dev/null +++ b/test/fepc_cp_test.c @@ -0,0 +1,84 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_cp.h" + +fepc_real_t +id(vec_real_p vector) { + return vector->array[0]; +} + +int +main() { + fepc_real_t h = 0.2; + + interval_t * intervals = (interval_t*) malloc(sizeof(interval_t)*1); + + interval_init(&intervals[0], 1); + intervals[0].start[0] = -1.; + intervals[0].end[0] = 1.; + + int rank1 = 3, rank2 = 4, dimension = 2; + + int n, k, interval_count = 1; + + int * maxlevels_cp1 = int_array_new(rank1*dimension); + + func_cp * cp1 = func_cp_new(rank1, dimension, maxlevels_cp1); + + free(maxlevels_cp1); + + for (n = 0; n < rank1; n++) { + for (k = 0; k < dimension; k++) { + setup_fepc_structure(&(cp1->functions[n*dimension+k]), id, &intervals, 1, 0, h); + } + } + + int * maxlevels_cp2 = int_array_new(rank2*dimension); + + func_cp * cp2 = func_cp_new(rank2, dimension, maxlevels_cp2); + free(maxlevels_cp2); + for (n = 0; n < rank2; n++) { + for (k = 0; k < dimension; k++) { + setup_fepc_structure(&(cp2->functions[n*dimension+k]), id, &intervals, 1, 0, h); + } + } + + func_t * result_intervals = (func_t*) malloc(sizeof(func_t)*dimension); + + for (k = 0; k < dimension; k++) { + func_init(&result_intervals[k], 0, 1); + set_gridstructure(&result_intervals[k], &intervals, h); + } + + func_cp * result = func_cp_faltung(cp1, cp2, result_intervals, h); + + + func_cp_del(cp1); + func_cp_del(cp2); + + func_cp_print(result); + + funcs_del_type(result_intervals, dimension); + intervals_del_type(intervals, 1); + func_cp_del(result); +#ifdef HAS_FFTW3 + fftw_cleanup(); +#endif + return 0; +} diff --git a/test/fepc_easy_test.c b/test/fepc_easy_test.c new file mode 100644 index 0000000..9985e8b --- /dev/null +++ b/test/fepc_easy_test.c @@ -0,0 +1,533 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + + + + +#define __DEBUG + +#define A 0.15 +#define B 4. +#define C -3.5 + +#include "fepc_easy.h" +#include "fepc_easy_helper.h" +#include "seconds.h" +#include + +#define STEPPING 0.5 + + +fepc_real_t const_1(vec_real_p x) { + return 1.; +} + +fepc_real_t f_coeff_peter(vec_p v, vec_p grad, int step) { + fepc_real_t h_l = get_h_l(step, STEPPING); + + return log((v->array[0] +1)*h_l + A) - log(v->array[0]*h_l + A); +} + +fepc_real_t g_coeff_peter(vec_p v, vec_p grad, int step) { + fepc_real_t h_l = get_h_l(step, STEPPING); + + return B*h_l*h_l*(v->array[0] + .5) + C *h_l; +} + +fepc_real_t f_peter(vec_real_p vector) { + return 1./(vector->array[0]+A); +} + +fepc_real_t g_peter(vec_real_p vector) { + return B*vector->array[0] + C; +} + +void set_result_peter_real(vec_real_set_p set) { + int n; + + for (n = 0; n < set->count; n++) { + set->element_results[n] -= set->elements[n]->array[1]*((B*(set->elements[n]->array[0]+A)+C)*(log(set->elements[n]->array[0]+A)-log(A))-B*set->elements[n]->array[0]); + } +} + +void test_example_peter() { + func_p fepc; + + int dimension, steps_f, steps_g, steps_w, grad_f, grad_g; + + dimension = 2; + + steps_f = 4; + steps_g = 1; + steps_w = 3; + + grad_f = 0; + grad_g = 0; + + interval_p * intervals_f = (interval_p*) malloc(sizeof(interval_p)*steps_f); + + + intervals_f[0] = interval_new(dimension); + + intervals_f[0]->start[0] = 0; + intervals_f[0]->start[1] = 0; + intervals_f[0]->end[0] = 1; + intervals_f[0]->end[1] = 1; + + + intervals_f[1] = interval_new(dimension); + + intervals_f[1]->start[0] = 0.; + intervals_f[1]->start[1] = 0.; + intervals_f[1]->end[0] = 0.2; + intervals_f[1]->end[1] = 1.; + + intervals_f[2] = interval_new(dimension); + + intervals_f[2]->start[0] = 0.; + intervals_f[2]->start[1] = 0.; + intervals_f[2]->end[0] = 0.2; + intervals_f[2]->end[1] = 1.0; + + intervals_f[3] = interval_new(dimension); + + intervals_f[3]->start[0] = 0.; + intervals_f[3]->start[1] = 0.; + intervals_f[3]->end[0] = 0.2; + intervals_f[3]->end[1] = 1.0; + + + interval_p * intervals_g = (interval_p*) malloc(sizeof(interval_p)*steps_g); + + intervals_g[0] = interval_new(dimension); + + intervals_g[0]->start[0] = 0; + intervals_g[0]->start[1] = 0; + intervals_g[0]->end[0] = 1; + intervals_g[0]->end[1] = 1; + + interval_p * intervals_w = (interval_p*) malloc(sizeof(interval_p)*steps_w); + + intervals_w[0] = interval_new(dimension); + + intervals_w[0]->start[0] = 0; + intervals_w[0]->start[1] = 0; + intervals_w[0]->end[0] = 1; + intervals_w[0]->end[1] = 1; + + intervals_w[1] = interval_new(dimension); + + intervals_w[1]->start[0] = 0.2; + intervals_w[1]->start[1] = 0.4; + intervals_w[1]->end[0] = 0.8; + intervals_w[1]->end[1] = 1; + + intervals_w[2] = interval_new(dimension); + + intervals_w[2]->start[0] = 0.4; + intervals_w[2]->start[1] = 0.8; + intervals_w[2]->end[0] = 0.6; + intervals_w[2]->end[1] = 1; + + fepc = fepc_easy(f_peter, intervals_f, steps_f, grad_f, g_peter, intervals_g, steps_g, grad_g, intervals_w, steps_w, STEPPING); // this will be our result + + #ifdef __DEBUG + + printf("==========================================\n"); + func_print(fepc, 3); + #endif + vec_real_set_p points = get_value(fepc, get_mean_points, STEPPING); + print_vec_real_set(points); + set_result_peter_real(points); + print_vec_real_set(points); +} + +void test_example_peter_orig() { + func_p f, g, w, fepc, ref; + fepc_real_t norm; + fepc_real_t a,b,c, mesh, temp1, temp2; + int dim, i,j , l, pos; + vec_p grad, start, lang, r, s; + vec_real_p h; + char dateiname1[20]; + + + + dim = 2; + a = 0.15; + b = 4.; + c = -3.5; + mesh = 0.2; + + + h = vec_real_new(5); + for(l=0;l<5;l++) { + h->array[l] = mesh*pow(2,-l); // 5-stufiges h-array anlegen + } + + + + + +/*Die Funktion f erstellen*/ + f = func_new(2,dim); + + // hp-Funktion erstellen + for(l=0;l<=f->maxlevel;l++) { + grad = vec_new(dim); + grad->array[0] = 0; + grad->array[1] = 0; + f->hierarchie[l] = folgen_vektor_new( grad ); + } + + // Gitterstruktur Stufe 0 + start = vec_new(dim); + lang = vec_new(dim); + start->array[0] = 0; + start->array[1] = 0; + lang->array[0] = 5; + lang->array[1] = 5; + f->hierarchie[0]->vektor[0] = folge_new(start,lang); + + // Gitterstruktur Stufe 1 und 2 + for(l=1;l<=f->maxlevel;l++) { + start = vec_new(dim); + lang = vec_new(dim); + start->array[0] = 0; + start->array[1] = 0; + lang->array[0] = 4; + lang->array[1] = pow(2,l)*5; + f->hierarchie[l]->vektor[0] = folge_new(start,lang); + } + + + l=0; + for(i=2;i<5;i++) { + for(j=0;j<5;j++) { + r = vec_new(dim); + r->array[0] = i; + r->array[1] = j; + pos = entry_d2one(r,f->hierarchie[l]->vektor[0]->lang); + vec_del(r); + f->hierarchie[l]->vektor[0]->glied[pos] = log((i+1)*h->array[l]+a) - log(i*h->array[l]+a); + } + } + + + for(l=1;lmaxlevel;l++) { + for(i=2;i<4;i++) { + for(j=0;j<(5*pow(2,l));j++) { + r = vec_new(dim); + r->array[0] = i; + r->array[1] = j; + pos = entry_d2one(r,f->hierarchie[l]->vektor[0]->lang); + vec_del(r); + f->hierarchie[l]->vektor[0]->glied[pos] = log((i+1)*h->array[l]+a) - log(i*h->array[l]+a); + } + } + } + + + l=f->maxlevel; + for(i=0;i<4;i++) { + for(j=0;j<(5*pow(2,l));j++) { + r = vec_new(dim); + r->array[0] = i; + r->array[1] = j; + pos = entry_d2one(r,f->hierarchie[l]->vektor[0]->lang); + vec_del(r); + f->hierarchie[l]->vektor[0]->glied[pos] = log((i+1)*h->array[l]+a) - log(i*h->array[l]+a); + } + } + + + + + /*Die Funktion g erstellen*/ + + g = func_new(0,dim); + + for(l=0;l<=g->maxlevel;l++) { + grad = vec_new(dim); + grad->array[0] = 0; + grad->array[1] = 0; + g->hierarchie[l] = folgen_vektor_new( grad ); + } + + for(l=0;l<=g->maxlevel;l++) { + start = vec_new(dim); + lang = vec_new(dim); + start->array[0] = 0; + lang->array[0] = 5*pow(2,l); + start->array[1] = 0; + lang->array[1] = 5*pow(2,l); + g->hierarchie[l]->vektor[0] = folge_new(start,lang); + } + + l=g->maxlevel; + for(i=0;i<(5*pow(2,l));i++) { + for(j=0;j<(5*pow(2,l));j++) { + r = vec_new(dim); + r->array[0] = i; + r->array[1] = j; + pos = entry_d2one(r,g->hierarchie[l]->vektor[0]->lang); + vec_del(r); + temp1 = b/2.*pow(h->array[l],2); + temp2 = 2*i + 1; + g->hierarchie[l]->vektor[0]->glied[pos] = temp1*temp2 + c*h->array[l]; + } + } + + + /*Die Funktion w erstellen*/ + w = func_new(2,dim); + + for(l=0;l<=w->maxlevel;l++) { + grad = vec_new(dim); + grad->array[0] = 0; + grad->array[1] = 0; + w->hierarchie[l] = folgen_vektor_new( grad ); + } + + start = vec_new(dim); + lang = vec_new(dim); + start->array[0] = 0; + lang->array[0] = 5; + start->array[1] = 0; + lang->array[1] = 5; + w->hierarchie[0]->vektor[0] = folge_new(start,lang); + + start = vec_new(dim); + lang = vec_new(dim); + start->array[0] = 2; + lang->array[0] = 6; + start->array[1] = 4; + lang->array[1] = 6; + w->hierarchie[1]->vektor[0] = folge_new(start,lang); + + start = vec_new(dim); + lang = vec_new(dim); + start->array[0] = 8; + lang->array[0] = 4; + start->array[1] = 16; + lang->array[1] = 4; + w->hierarchie[2]->vektor[0] = folge_new(start,lang); +printf("==========================================\n");printf("==========================================\n"); + #ifdef __DEBUG + //print_func(f); + //print_func(g); + //print_func(w); + #endif + + + /*Durchführen der Berechnungen*/ + fepc = faltung_fepc(f,g,w,mesh); + ref = faltung_ref(f,g,w,mesh); + #ifdef __DEBUG + printf("==========================================\n"); + //print_func(fepc); + #endif + vec_real_set_p points = get_value(fepc, get_mean_points, STEPPING); + norm = faltung_hilfe_norm(fepc,ref); + printf("\n Norm %lf \n",norm); + print_vec_real_set(points); + +} + + +void some_other_test() { + func_p f, g, w; + + int dimension = 2; + + int steps_f = 3; + + interval_p * intervals = (interval_p*) malloc(sizeof(interval_p)*steps_f); + + intervals[0] = interval_new(dimension); + + intervals[0]->start[0] = 0; + intervals[0]->start[1] = 0; + intervals[0]->end[0] = 2; + intervals[0]->end[1] = 2; + + intervals[1] = interval_new(dimension); + + intervals[1]->start[0] = 0.8; + intervals[1]->start[1] = 0.8; + intervals[1]->end[0] = 2; + intervals[1]->end[1] = 2; + + intervals[2] = interval_new(dimension); + + intervals[2]->start[0] = 1.2; + intervals[2]->start[1] = 1.2; + intervals[2]->end[0] = 2; + intervals[2]->end[1] = 2; + + f = func_new(steps_f-1, dimension); + + #ifdef __DEBUG + print_intervals(intervals, steps_f); + printf("\nwill be converted to\n"); + print_intervals(convert_interval(steps_f, intervals, STEPPING), steps_f); + #endif + //print_intervals(intervals, 3); + set_gridstructure(f, intervals, STEPPING); + + int n; + printf("Hier2\n"); + for (n = 0; n < steps_f; n++) {// add the hp-structure of f + //add_folgenentries(f, n); + } + + + + + + + + g = func_new(0, dimension); + + + w = func_new(2, dimension); + + // generate the grid structure for w + + //func_p fepc = faltung_fepc(f, g, w, STEPPING); // this will be our result + + +} + +void test_integration() { + int step = 2; + + + + vec_p v = vec_new(2); + v->array[0] = 0.; + v->array[1] = 0.; + + vec_p p = vec_new(2); + + fepc_real_t result, h_l, check_result, result_2; + + h_l = get_h_l(step, STEPPING); + + + check_result = g_coeff_peter(v, p, step); + result = integrate_coeff(g_peter, v, p, step, STEPPING); + result_2 = integrate_coeff_st(g_peter, v, p, step, STEPPING); + double seconds_used = seconds(); + printf("\nIntegration %.10f\n", integrate_coeff(f_peter, v, p, step, STEPPING)); + printf("seconds: %f\n", seconds()-seconds_used); + seconds_used = seconds(); + printf("\nIntegration_st %.10lf\n", integrate_coeff_st(f_peter, v, p, step, STEPPING)); + printf("seconds: %f\n", seconds()-seconds_used); + printf("real result \n======\n %.10lf\n", f_coeff_peter(v, p, step)); + + printf("\nIntegrate %.10f\nExact %.10lf\nIntegrate_st %.10lf\n", result, check_result, result_2); + vec_del(v); + vec_del(p); +} + +void test_multi() { + func_p f, g, multi; + + int dimension, steps_f, steps_g, grad_f, grad_g; + + dimension = 2; + steps_f = 2; + steps_g = 2; + + grad_f = 1; + grad_g = 0; + + interval_p * intervals = (interval_p*) malloc(sizeof(interval_p)*steps_f); + + intervals[0] = interval_new(dimension); + + intervals[0]->start[0] = 0; + intervals[0]->start[1] = 0; + intervals[0]->end[0] = 1; + intervals[0]->end[1] = 1; + + intervals[0] = interval_new(dimension); + + intervals[0]->start[0] = 0; + intervals[0]->start[1] = 0; + intervals[0]->end[0] = 1; + intervals[0]->end[1] = 1; + + intervals[1] = interval_new(dimension); + + intervals[1]->start[0] = 0.; + intervals[1]->start[1] = 0.; + intervals[1]->end[0] = 0.4; + intervals[1]->end[1] = 1.; + + f = create_fepc_structure(f_peter, intervals, steps_f, grad_f, STEPPING); + g = create_fepc_structure(g_peter, intervals, steps_g, grad_g, STEPPING); + multi = func_multi(f, g, STEPPING); + + vec_real_set_p points = get_value(multi, get_mean_points, STEPPING); + print_vec_real_set(points); + points = get_value(f, get_mean_points, STEPPING); + print_vec_real_set(points); + points = get_value(g, get_mean_points, STEPPING); + print_vec_real_set(points); +} + +void test_norm() { + int dimension = 2; + + interval_p * intervals = (interval_p*) malloc(sizeof(interval_p)*1); + + + intervals[0] = interval_new(dimension); + + intervals[0]->start[0] = 0.; + intervals[0]->start[1] = 0.; + intervals[0]->end[0] = 1.; + intervals[0]->end[1] = 1.; + + func_p function = create_fepc_structure(const_1, intervals, 1, 0, STEPPING); + + //func_print(function, 3); + + fepc_real_t integral = func_integrate(function, STEPPING); + + //vec_real_set_p set = get_value(function, get_mean_points, STEPPING); + + //print_vec_real_set(set); + printf("Integral = %.10lf\n", integral); + +} + +int main(void) { + test_norm(); + test_integration(); + //test_example_peter_orig(); + //test_example_peter(); + //test_multi(); + return 0; +} + + + + diff --git a/test/sparsetest.c b/test/sparsetest.c new file mode 100644 index 0000000..96e378c --- /dev/null +++ b/test/sparsetest.c @@ -0,0 +1,28 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_easy_sparse.h" + + +int main(void) { + func_sparse_p func_sparse = func_sparse_new(2,3); + + func_sparse_del(func_sparse); + + return 0; +} \ No newline at end of file diff --git a/test/stest.c b/test/stest.c new file mode 100644 index 0000000..8884bf2 --- /dev/null +++ b/test/stest.c @@ -0,0 +1,111 @@ +/* + * FEPC + * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +#include "fepc_easy.h" +#include "fepc_easy_helper.h" + +fepc_real_t norm_x_sqr(vec_real_p x) { + int n; + + fepc_real_t result = 0.; + + for (n = 0; n < x->dim; n++) { + result += x->array[n]*x->array[n]; + } + return result; +} + +fepc_real_t const_2dim(vec_real_p x) { + return x->dim*2.; +} + +fepc_real_t two_xone(vec_real_p x) { + return x->array[1]*2; +} + +fepc_real_t nine_xone_sqr(vec_real_p x) { + return 9*x->array[0]*x->array[0]; +} + +fepc_real_t function(vec_real_p x) { + return 3*x->array[0]*x->array[0]*x->array[0] + x->array[1]*x->array[1]*x->array[1]*x->array[1]; + +} + +fepc_real_t function_laplace(vec_real_p x) { + return 18*x->array[0]+12*x->array[1]*x->array[1]; +} + + + +int main() { + fepc_real_t stepping = 0.001; + + interval_p * intervals = (interval_p*) malloc(sizeof(interval_p)*1); + + int n, dimension = 2; + + intervals[0] = interval_new(dimension); + + for (n = 0; n < dimension; n++) { + intervals[0]->start[n] = -.01; + intervals[0]->end[n] = .01; + } + + //func_p func_const_2dim = setup_fepc_structure(const_2dim, intervals, 1, 0, stepping); + + + //func_p func_norm_x_sqr = setup_fepc_structure(norm_x_sqr, intervals, 1, 0, stepping); + + //func_p func_2x1 = setup_fepc_structure(two_xone, intervals, 1, 0, stepping); + + //func_p func_9x1_sqr = setup_fepc_structure(nine_xone_sqr, intervals, 1, 0, stepping); + + func_p func_function = create_fepc_structure(function, intervals, 1, 0, stepping); + + printf("Funktion 1 fertig\n"); + func_p func_function_laplace = create_fepc_structure(function_laplace, intervals, 1, 0, stepping); + printf("Funktion 2 fertig\n"); + func_p func_function_laplace_calc = func_laplace(func_function, stepping); + printf("Laplace fertig\n"); + + //func_p func_laplace_norm_x_sqr = func_laplace(func_norm_x_sqr, stepping); + + //func_print(func_function, 3); + //func_print(func_derive(func_function, 0, stepping), 3); + //func_print(func_9x1_sqr, 3); + + //func_print(func_function_laplace_calc, 3); + //func_print(func_function_laplace, 3); + + func_p temp = func_subtract(func_function_laplace_calc, func_function_laplace); + + + fepc_real_t norm = func_norm_l2_sqr(temp, stepping); + + + + /*func_del(func_const_6); + func_del(func_norm_x_sqr); + func_del(func_laplace_norm_x_sqr); + func_del(temp); + */ + printf("Norm^2 = %1.15f\n", norm); + + return 0; +}