From b652c096dcc0690f8fdcd5d6500c391ece9e845a Mon Sep 17 00:00:00 2001 From: Stefan Date: Tue, 11 May 2010 13:01:30 +0200 Subject: [PATCH 01/24] Adds CMake as a build tool.\nAdds easy invokation. Adds discontinous function support --- Makefile | 55 -- basic.h => include/basic.h | 40 + config.h => include/config.h | 9 +- include/discont.h | 101 +++ faltung.h => include/faltung.h | 0 faltung_hilfe.h => include/faltung_hilfe.h | 10 + include/fepc_easy.h | 156 ++++ include/fepc_easy_helper.h | 73 ++ include/fepc_easy_sparse.h | 65 ++ fft_faltung.h => include/fft_faltung.h | 0 folge.h => include/folge.h | 8 + folgen_vektor.h => include/folgen_vektor.h | 11 +- funktion.h => include/funktion.h | 8 +- include/interval.h | 45 ++ koeffizienten.h => include/koeffizienten.h | 3 + kuerzen.h => include/kuerzen.h | 0 seconds.h => include/seconds.h | 0 source/CMakeLists.txt | 6 + basic.c => source/basic.c | 122 ++- source/discont.c | 243 ++++++ faltung.c => source/faltung.c | 6 + faltung_hilfe.c => source/faltung_hilfe.c | 1 + source/fepc_easy.c | 861 +++++++++++++++++++++ source/fepc_easy_helper.c | 508 ++++++++++++ source/fepc_easy_sparse.c | 97 +++ fft_faltung.c => source/fft_faltung.c | 1 + folge.c => source/folge.c | 68 ++ folgen_vektor.c => source/folgen_vektor.c | 106 ++- funktion.c => source/funktion.c | 41 +- source/interval.c | 93 +++ koeffizienten.c => source/koeffizienten.c | 5 + kuerzen.c => source/kuerzen.c | 3 + main.c => source/main.c | 0 seconds.c => source/seconds.c | 1 + test/CMakeLists.txt | 23 + test/discont_test.c | 86 ++ test/fepc_easy_test.c | 533 +++++++++++++ test/sparsetest.c | 28 + test/stest.c | 111 +++ 39 files changed, 3456 insertions(+), 71 deletions(-) delete mode 100644 Makefile rename basic.h => include/basic.h (78%) mode change 100644 => 100755 rename config.h => include/config.h (91%) mode change 100644 => 100755 create mode 100755 include/discont.h rename faltung.h => include/faltung.h (100%) mode change 100644 => 100755 rename faltung_hilfe.h => include/faltung_hilfe.h (99%) mode change 100644 => 100755 create mode 100755 include/fepc_easy.h create mode 100755 include/fepc_easy_helper.h create mode 100644 include/fepc_easy_sparse.h rename fft_faltung.h => include/fft_faltung.h (100%) mode change 100644 => 100755 rename folge.h => include/folge.h (94%) mode change 100644 => 100755 rename folgen_vektor.h => include/folgen_vektor.h (93%) mode change 100644 => 100755 rename funktion.h => include/funktion.h (97%) mode change 100644 => 100755 create mode 100755 include/interval.h rename koeffizienten.h => include/koeffizienten.h (99%) mode change 100644 => 100755 rename kuerzen.h => include/kuerzen.h (100%) mode change 100644 => 100755 rename seconds.h => include/seconds.h (100%) mode change 100644 => 100755 create mode 100755 source/CMakeLists.txt rename basic.c => source/basic.c (78%) mode change 100644 => 100755 create mode 100755 source/discont.c rename faltung.c => source/faltung.c (99%) mode change 100644 => 100755 rename faltung_hilfe.c => source/faltung_hilfe.c (99%) mode change 100644 => 100755 create mode 100755 source/fepc_easy.c create mode 100755 source/fepc_easy_helper.c create mode 100644 source/fepc_easy_sparse.c rename fft_faltung.c => source/fft_faltung.c (99%) mode change 100644 => 100755 rename folge.c => source/folge.c (87%) mode change 100644 => 100755 rename folgen_vektor.c => source/folgen_vektor.c (87%) mode change 100644 => 100755 rename funktion.c => source/funktion.c (92%) mode change 100644 => 100755 create mode 100755 source/interval.c rename koeffizienten.c => source/koeffizienten.c (99%) mode change 100644 => 100755 rename kuerzen.c => source/kuerzen.c (99%) mode change 100644 => 100755 rename main.c => source/main.c (100%) rename seconds.c => source/seconds.c (99%) mode change 100644 => 100755 create mode 100755 test/CMakeLists.txt create mode 100755 test/discont_test.c create mode 100755 test/fepc_easy_test.c create mode 100644 test/sparsetest.c create mode 100755 test/stest.c 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/basic.h b/include/basic.h old mode 100644 new mode 100755 similarity index 78% rename from basic.h rename to include/basic.h index a0c14d9..00c8ca8 --- a/basic.h +++ b/include/basic.h @@ -68,6 +68,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 +80,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 +96,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 +144,24 @@ 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); #endif + + diff --git a/config.h b/include/config.h old mode 100644 new mode 100755 similarity index 91% rename from config.h rename to include/config.h index c2eb87a..a3408cd --- a/config.h +++ b/include/config.h @@ -38,7 +38,14 @@ * 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; diff --git a/include/discont.h b/include/discont.h new file mode 100755 index 0000000..3cb38ea --- /dev/null +++ b/include/discont.h @@ -0,0 +1,101 @@ +/* + * 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 . + */ + +#ifndef _DISCONT +#define _DISCONT + +#define ONE_THIRD 0.333333333333333333333333333333333333333333333333333333333333333333333 +#define TWO_THIRD 0.666666666666666666666666666666666666666666666666666666666666666666666 +#define SQRT_12 3.464101615137754587054892683011744733885610507620761256111613958903866 + +#include "fepc_easy.h" + +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 x_0, fepc_real_t x_1, fepc_real_t dx); + +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 position, int v, 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); + + +#endif diff --git a/faltung.h b/include/faltung.h old mode 100644 new mode 100755 similarity index 100% rename from faltung.h rename to include/faltung.h diff --git a/faltung_hilfe.h b/include/faltung_hilfe.h old mode 100644 new mode 100755 similarity index 99% rename from faltung_hilfe.h rename to include/faltung_hilfe.h index dfde42a..24663c5 --- a/faltung_hilfe.h +++ b/include/faltung_hilfe.h @@ -59,3 +59,13 @@ faltung_hilfe_Gamma_bauen_2( func_p g, vec_real_p mesh, vec_p maxgrad, matrix3_p #endif + + + + + + + + + + diff --git a/include/fepc_easy.h b/include/fepc_easy.h new file mode 100755 index 0000000..d3af7ac --- /dev/null +++ b/include/fepc_easy.h @@ -0,0 +1,156 @@ +/* + * 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 . + */ + +#ifndef __FEPCEASY +#define __FEPCEASY + +#include +#include +#include +#include "faltung.h" +#include "interval.h" + +#ifndef INT_STEPS +#define INT_STEPS 20 +#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); + +func_p setup_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); + +#endif + diff --git a/include/fepc_easy_helper.h b/include/fepc_easy_helper.h new file mode 100755 index 0000000..a27cca1 --- /dev/null +++ b/include/fepc_easy_helper.h @@ -0,0 +1,73 @@ +/* + * 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 . + */ + +#ifndef _FEPC_EASY_HELPER +#define _FEPC_EASY_HELPER + +#include "fepc_easy.h" + +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); + +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); + + +#endif diff --git a/include/fepc_easy_sparse.h b/include/fepc_easy_sparse.h new file mode 100644 index 0000000..9a65c12 --- /dev/null +++ b/include/fepc_easy_sparse.h @@ -0,0 +1,65 @@ +/* + * 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 . + */ + +#ifndef __FEPCEASYSPARSE +#define __FEPCEASYSPARSE + +#include "fepc_easy.h" + + + +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); + + + + + + + + + + + + +#endif \ No newline at end of file diff --git a/fft_faltung.h b/include/fft_faltung.h old mode 100644 new mode 100755 similarity index 100% rename from fft_faltung.h rename to include/fft_faltung.h diff --git a/folge.h b/include/folge.h old mode 100644 new mode 100755 similarity index 94% rename from folge.h rename to include/folge.h index 0b1ca88..b67f355 --- a/folge.h +++ b/include/folge.h @@ -69,6 +69,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 +82,11 @@ 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); #endif + + diff --git a/folgen_vektor.h b/include/folgen_vektor.h old mode 100644 new mode 100755 similarity index 93% rename from folgen_vektor.h rename to include/folgen_vektor.h index 8606b32..5f0e70c --- a/folgen_vektor.h +++ b/include/folgen_vektor.h @@ -96,9 +96,18 @@ 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); - #endif + + diff --git a/funktion.h b/include/funktion.h old mode 100644 new mode 100755 similarity index 97% rename from funktion.h rename to include/funktion.h index 9f0de7c..669174d --- a/funktion.h +++ b/include/funktion.h @@ -68,7 +68,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 +91,7 @@ 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); + /* 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 +111,8 @@ void func_grid_zero(func_p f); +func_p +func_factor_multi(func_p function, fepc_real_t factor); + + #endif diff --git a/include/interval.h b/include/interval.h new file mode 100755 index 0000000..6b166ab --- /dev/null +++ b/include/interval.h @@ -0,0 +1,45 @@ +#ifndef __INTERVALS +#define __INTERVALS + +#include "basic.h" + +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_del(interval_p interval); + +void +intervals_del(interval_p* 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); + +#endif + diff --git a/koeffizienten.h b/include/koeffizienten.h old mode 100644 new mode 100755 similarity index 99% rename from koeffizienten.h rename to include/koeffizienten.h index dbc4eb5..70428f6 --- a/koeffizienten.h +++ b/include/koeffizienten.h @@ -82,3 +82,6 @@ koeffizienten_gamma(int level, vec_p r, vec_real_p mesh, vec_p alpha, vec_p mu, #endif + + + diff --git a/kuerzen.h b/include/kuerzen.h old mode 100644 new mode 100755 similarity index 100% rename from kuerzen.h rename to include/kuerzen.h diff --git a/seconds.h b/include/seconds.h old mode 100644 new mode 100755 similarity index 100% rename from seconds.h rename to include/seconds.h diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt new file mode 100755 index 0000000..90fed2e --- /dev/null +++ b/source/CMakeLists.txt @@ -0,0 +1,6 @@ +file(GLOB ALL_SOURCE_FILES *.c) + +# Find all header files in the include directory +file(GLOB ALL_HEADER_FILES ${CMAKE_SOURCE_DIR}/include/*.h) + +ADD_LIBRARY(fepc SHARED ${ALL_SOURCE_FILES} ${ALL_HEADER_FILES}) diff --git a/basic.c b/source/basic.c old mode 100644 new mode 100755 similarity index 78% rename from basic.c rename to source/basic.c index 384e516..e9a9586 --- 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 100755 index 0000000..2a09017 --- /dev/null +++ b/source/discont.c @@ -0,0 +1,243 @@ +/* + * 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_p)); + + 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 dx) { + linear_function_p result = (linear_function_p) malloc(sizeof(linear_function_p)); + + result->y_0 = y_0; + result->slope = (y_1 - y_0)/dx; + return result; +} + + +linear_function_set_p +linear_function_set_new(int count) { + linear_function_set_p result = (linear_function_set_p) malloc(sizeof(linear_function_set_p)); + + result->count = count; + result->functions = (linear_function_p*) malloc(sizeof(linear_function_t)*count); + return result; +} + + +discont_function_p +discont_function_new(int stepcount) { + discont_function_p result = (discont_function_p) malloc(sizeof(discont_function_p)); + + 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++) { + 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], h_l); + } + } +} + + +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; + + 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); + result->function_sets[n] = linear_function_set_new(stepcount); + for (k = 0; k < stepcount; k++) { + /* + to make sure we are are using the correct linear function, we use the points 1/3 and 2/3 instead of 0 and 1 to + calculate the linear function + */ + temp_x1 = (k + ONE_THIRD)*h_l; + temp_x2 = (k + TWO_THIRD)*h_l; + x = vec_real_new(1); + x->array[0] = temp_x2; + temp_y1 = get_value_at_step(function, x, n, stepping); + vec_real_del(x); + x = vec_real_new(1); + 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); + result->function_sets[n]->functions[k] = linear_function_new(slope*k*h_l+temp_y1, slope); + } + } + return result; +} + + +fepc_real_t +integrate_coeff_discont(discont_function_p function, int position, int v, int p, int step, 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(step, stepping); + + slope = function->function_sets[step]->functions[position]->slope; + y_0 = function->function_sets[step]->functions[position]->y_0; + + if (p == 0) { // legendre(0) = sqrt(1/h_l) + return sqrt(h_l)*(y_0 + h_l*(slope/2.0)*pow(v+1, 2)- pow(v, 2)); + } else { // p == 1 --> legendre(1) = sqrt(12)(x-(v+0.5)*h_l)/(h_l^1.5) + return SQRT_12*h_l*((slope/3.0)*(pow(v+1, 3)- pow(v, 3))*sqrt(h_l) + ((y_0 / h_l -1*slope*(v+0.5))/2.0)*(pow(v+1, 2)- pow(v, 2)) - y_0*(v+0.5)*h_l); + } +} + + +void +add_folgenentries_discont(func_p function, discont_function_p discont_function, fepc_real_t stepping) { + + int step, n, k, interval_element_count, position, grad_count, pos_2; + + vec_p v, x; + + 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); + pos_2 = v->array[0]; + 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 + x = entry_one2d_sloppy(k, function->hierarchie[step]->grad); // don't care about 0 as entries + function->hierarchie[step]->vektor[k]->glied[position] = integrate_coeff_discont(discont_function, pos_2 ,v->array[0], x->array[0], step, stepping); + vec_del(x); + } + vec_del(v); + } + } + } +} + + diff --git a/faltung.c b/source/faltung.c old mode 100644 new mode 100755 similarity index 99% rename from faltung.c rename to source/faltung.c index d0c7d8b..0744024 --- a/faltung.c +++ b/source/faltung.c @@ -676,3 +676,9 @@ faltung_fepc(func_p f, func_p g, func_p w, fepc_real_t h) { func_grid_zero( back ); return back; } + + + + + + diff --git a/faltung_hilfe.c b/source/faltung_hilfe.c old mode 100644 new mode 100755 similarity index 99% rename from faltung_hilfe.c rename to source/faltung_hilfe.c index 6059dda..24cf201 --- 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_easy.c b/source/fepc_easy.c new file mode 100755 index 0000000..17b662e --- /dev/null +++ b/source/fepc_easy.c @@ -0,0 +1,861 @@ +/* + * 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*x->array[n]/h_l - 2*v->array[n] - 1); + } + 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); + for (k = 0, 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 + folge_del(folge); + } +} + +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+current*h; + } + 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; + + + // 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); + } + + //printf("sum = %f, result = %f ", result, pow(h_l, v->dim)*result/calc_count); + + 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, x; + + 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 + x = 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, x, step, stepping) : integrate_coeff_st(function_impl, v, x, step, stepping); + vec_del(x); + } + 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); +} + + + +func_p setup_fepc_structure(Funcimpl function, interval_p* intervals, int interval_count, int degree, fepc_real_t stepping) { + func_p result = func_new(interval_count-1, intervals[interval_count-1]->dimension); + + if (degree > 0) { + set_degree(result, degree); + } + + set_gridstructure(result, intervals, stepping); + if (function != NULL) { + add_folgenentries(result, function, NULL, 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; +} + +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; + +} + +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 = setup_fepc_structure(function1, intervals1, interval_count1, degree_func1, stepping); + g = setup_fepc_structure(function2, intervals2, interval_count2, degree_func2, stepping); + w = setup_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 100755 index 0000000..5cd4fdd --- /dev/null +++ b/source/fepc_easy_helper.c @@ -0,0 +1,508 @@ +/* + * 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_helper.h" +#include "fepc_easy.h" + + +func_p +func_multi(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_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; + } + return result; +} + +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; + } + + 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*folge_norm2(g, step, stepping); + //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; + } + + 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*folge_norm2(g, step, stepping)); + } + 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 old mode 100644 new mode 100755 similarity index 99% rename from fft_faltung.c rename to source/fft_faltung.c index 5d48651..baaad2a --- a/fft_faltung.c +++ b/source/fft_faltung.c @@ -160,3 +160,4 @@ fft_faltung(fepc_real_t* a, vec_p n_a, fepc_real_t* b, vec_p n_b) { exit( 1 ); #endif } + diff --git a/folge.c b/source/folge.c old mode 100644 new mode 100755 similarity index 87% rename from folge.c rename to source/folge.c index 47c4d9d..e81bd95 --- a/folge.c +++ b/source/folge.c @@ -408,6 +408,50 @@ 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 ); // TODO multiply with -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 +494,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 old mode 100644 new mode 100755 similarity index 87% rename from folgen_vektor.c rename to source/folgen_vektor.c index 6fb0a35..8ef6660 --- 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 old mode 100644 new mode 100755 similarity index 92% rename from funktion.c rename to source/funktion.c index 713a5d8..a889822 --- a/funktion.c +++ b/source/funktion.c @@ -66,13 +66,15 @@ func_new(int maxlevel, int dim) { void func_del(func_p f) { - int k; - - for(k=0;k<=f->maxlevel;k++) { - folgen_vektor_del(f->hierarchie[k]); - } - free(f->hierarchie); - free(f); + if (f && f != NULL) { + int k; + + for(k=0;k<=f->maxlevel;k++) { + folgen_vektor_del(f->hierarchie[k]); + } + free(f->hierarchie); + free(f); + } } @@ -140,6 +142,10 @@ func_projekt(func_p f,func_p g) { return back; } +func_p +func_clone(func_p f) { + return func_projekt(f, f); +} @@ -189,8 +195,6 @@ 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; @@ -224,6 +228,7 @@ func_add(func_p f, func_p g) { return back; } + int func_count( func_p f, func_p g, func_p w ) { int back; @@ -383,3 +388,21 @@ 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; +} + diff --git a/source/interval.c b/source/interval.c new file mode 100755 index 0000000..91111ae --- /dev/null +++ b/source/interval.c @@ -0,0 +1,93 @@ +/* + * 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)); + + result->dimension = dimension; + result->start = (fepc_real_t*) malloc(sizeof(fepc_real_t)*dimension); + result->end = (fepc_real_t*) malloc(sizeof(fepc_real_t)*dimension); + + int n; + + for (n = 0; n < dimension; n++) { + result->start[n] = 0.; + } + return result; +} + +void interval_del(interval_p interval) { + free(interval->start); + free(interval->end); + free(interval); +} + +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 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 old mode 100644 new mode 100755 similarity index 99% rename from koeffizienten.c rename to source/koeffizienten.c index 72e4327..622d922 --- 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 old mode 100644 new mode 100755 similarity index 99% rename from kuerzen.c rename to source/kuerzen.c index e74b839..4b3832d --- 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 old mode 100644 new mode 100755 similarity index 99% rename from seconds.c rename to source/seconds.c index 75fef72..3d4bbd7 --- 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 100755 index 0000000..f9fb84b --- /dev/null +++ b/test/CMakeLists.txt @@ -0,0 +1,23 @@ +include_directories (${${PROJECT_NAME}_SOURCE_DIR}/src) + +link_directories (${${PROJECT_NAME}_BINARY_DIR}/src) + + +add_executable (fepc_easy_test fepc_easy_test.c) +target_link_libraries (fepc_easy_test fepc fftw3 m) + +add_executable (stest stest.c) +target_link_libraries (stest fepc fftw3 m) + +add_executable (discont_test discont_test.c) +target_link_libraries (discont_test fepc fftw3 m) + + +add_executable (sparsetest sparsetest.c) +target_link_libraries (sparsetest fepc fftw3 m) + + +add_test(fepc_easy_test fepc_easy_test) +add_test(stest stest) +add_test(discont_test discont_test) +add_test(sparsetest sparsetest) \ No newline at end of file diff --git a/test/discont_test.c b/test/discont_test.c new file mode 100755 index 0000000..e6de36e --- /dev/null +++ b/test/discont_test.c @@ -0,0 +1,86 @@ +/* + * 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" + + +int +main(void) { + /* Example */ + + fepc_real_t stepping = 0.00001; + + int steps, n; + + func_p f1, f2, w, convolution_result; + + steps = 2; + + discont_function_p function, result; + + function = discont_function_new(steps); + + fepc_real_t y11[] = {0.1, 0.0, 0.3, 3.0, 4.2, 0.3, 5.8, 7.3, 1.0, 2.4}; // the stepcount is 10 in the first level because h_0 is 0.1 + fepc_real_t y21[] = {0.2, 1.3, 0.1, 6.5, 4.2, 3.3, 4.1, 0.3, 9.3, 2.2}; + + + // 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[] = {3.4, 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); + discont_function_del(function); + + func_print(f1, 3); + + + // Create f2 in the same way! + + + // 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_gridstructure(w, function->intervals, stepping); + + convolution_result = faltung_fepc(f1, f2, w, stepping); + + func_del(f1); + func_del(f2); + func_del(w); + + result = convert_func(convolution_result, function->intervals, stepping); + discont_function_del(function); + + func_del(convolution_result); + + + // the result is stored in "result" as a discont-function + + return 0; +} diff --git a/test/fepc_easy_test.c b/test/fepc_easy_test.c new file mode 100755 index 0000000..bc99334 --- /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 = setup_fepc_structure(f_peter, intervals, steps_f, grad_f, STEPPING); + g = setup_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 = setup_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 100755 index 0000000..99d56b0 --- /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 = setup_fepc_structure(function, intervals, 1, 0, stepping); + + printf("Funktion 1 fertig\n"); + func_p func_function_laplace = setup_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; +} From f92129a85de910649c380d768e7cb0b8c3e52cf3 Mon Sep 17 00:00:00 2001 From: Stefan Date: Tue, 11 May 2010 13:02:55 +0200 Subject: [PATCH 02/24] Changes copyright holders. --- README | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README b/README index cd3feab..5dab4c0 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 From 7b96447b65352292b1950edc10ca89f7cd605c74 Mon Sep 17 00:00:00 2001 From: Stefan Date: Tue, 11 May 2010 13:03:42 +0200 Subject: [PATCH 03/24] Adds forgotten cmake instructions. --- CMakeLists.txt | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100755 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100755 index 0000000..1766197 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,12 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.6) + +PROJECT(FEPC) + +include_directories(include) + +ADD_SUBDIRECTORY(source) + +enable_testing() +ADD_SUBDIRECTORY(test) + + From 6866ceb4cdb22233bcb5b01172e0352e3565ce67 Mon Sep 17 00:00:00 2001 From: Stefan Date: Tue, 11 May 2010 14:59:53 +0200 Subject: [PATCH 04/24] Memory allocation fixes. --- source/discont.c | 15 ++++++++------- test/discont_test.c | 7 ++++--- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/source/discont.c b/source/discont.c index 2a09017..b304fb7 100755 --- a/source/discont.c +++ b/source/discont.c @@ -21,7 +21,7 @@ fepc_real_t _round(fepc_real_t value) { - if (value - ((int) value) >= 0.5) { + if (value - ((int) value) < 0.5) { return (fepc_real_t) (int) value; } else { return (fepc_real_t) (int) value +1.0; @@ -30,7 +30,7 @@ _round(fepc_real_t value) { 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_p)); + linear_function_p result = (linear_function_p) malloc(sizeof(linear_function_t)); result->y_0 = y_0; result->slope = slope; @@ -40,7 +40,7 @@ linear_function_new(fepc_real_t y_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 dx) { - linear_function_p result = (linear_function_p) malloc(sizeof(linear_function_p)); + linear_function_p result = (linear_function_p) malloc(sizeof(linear_function_t)); result->y_0 = y_0; result->slope = (y_1 - y_0)/dx; @@ -50,17 +50,17 @@ linear_function_new_points(fepc_real_t y_0, fepc_real_t y_1, fepc_real_t dx) { linear_function_set_p linear_function_set_new(int count) { - linear_function_set_p result = (linear_function_set_p) malloc(sizeof(linear_function_set_p)); + 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_t)*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_p)); + discont_function_p result = (discont_function_p) malloc(sizeof(discont_function_t)); result->stepcount = stepcount; result->intervals = (interval_p*) malloc(sizeof(interval_p)*stepcount); @@ -111,8 +111,9 @@ discont_function_setup_points(discont_function_p function, int step, fepc_real_t 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], h_l); diff --git a/test/discont_test.c b/test/discont_test.c index e6de36e..4018134 100755 --- a/test/discont_test.c +++ b/test/discont_test.c @@ -22,8 +22,8 @@ int main(void) { /* Example */ - - fepc_real_t stepping = 0.00001; + + fepc_real_t stepping = 0.1; int steps, n; @@ -55,12 +55,13 @@ main(void) { discont_function_del(function); func_print(f1, 3); - + // Create f2 in the same way! // 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 From 57b08374413928a7a9c2e7d420a95780b0a8475f Mon Sep 17 00:00:00 2001 From: Stefan Date: Tue, 11 May 2010 16:12:26 +0200 Subject: [PATCH 05/24] Fixes integration bug at level 0. --- source/discont.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/discont.c b/source/discont.c index b304fb7..24829e8 100755 --- a/source/discont.c +++ b/source/discont.c @@ -116,7 +116,7 @@ discont_function_setup_points(discont_function_p function, int step, fepc_real_t 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], h_l); + function->function_sets[step]->functions[n] = linear_function_new_points(y1[0], y1[0], h_l); // eigentlich y2[n] als 2ter parameter } } } @@ -200,15 +200,15 @@ convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { fepc_real_t integrate_coeff_discont(discont_function_p function, int position, int v, int p, int step, 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(step, stepping); - + slope = function->function_sets[step]->functions[position]->slope; y_0 = function->function_sets[step]->functions[position]->y_0; if (p == 0) { // legendre(0) = sqrt(1/h_l) - return sqrt(h_l)*(y_0 + h_l*(slope/2.0)*pow(v+1, 2)- pow(v, 2)); + return sqrt(h_l)*(y_0 + h_l*(slope/2.0)*(pow(v+1, 2)- pow(v, 2))); } else { // p == 1 --> legendre(1) = sqrt(12)(x-(v+0.5)*h_l)/(h_l^1.5) return SQRT_12*h_l*((slope/3.0)*(pow(v+1, 3)- pow(v, 3))*sqrt(h_l) + ((y_0 / h_l -1*slope*(v+0.5))/2.0)*(pow(v+1, 2)- pow(v, 2)) - y_0*(v+0.5)*h_l); } From 73fda48318408eff579aaa1b511a57c4d17647b9 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Wed, 12 May 2010 01:55:29 +0200 Subject: [PATCH 06/24] Fixes calculation of basis coefficients. --- source/discont.c | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/source/discont.c b/source/discont.c index 24829e8..14c2fc9 100755 --- a/source/discont.c +++ b/source/discont.c @@ -116,7 +116,7 @@ discont_function_setup_points(discont_function_p function, int step, fepc_real_t if (y1 != NULL && y2 != NULL) { for (n = 0; n < count; n++) { - function->function_sets[step]->functions[n] = linear_function_new_points(y1[0], y1[0], h_l); // eigentlich y2[n] als 2ter parameter + function->function_sets[step]->functions[n] = linear_function_new_points(y1[n], y2[n], h_l); } } } @@ -201,16 +201,17 @@ fepc_real_t integrate_coeff_discont(discont_function_p function, int position, int v, int p, int step, fepc_real_t stepping) { // integrate from v[0]*h_l till (v[0]+1)*h_l - fepc_real_t h_l, slope, y_0; + fepc_real_t h_l, slope, y_0, sqrt_h_l; h_l = get_h_l(step, stepping); slope = function->function_sets[step]->functions[position]->slope; y_0 = function->function_sets[step]->functions[position]->y_0; + sqrt_h_l = sqrt(h_l); if (p == 0) { // legendre(0) = sqrt(1/h_l) - return sqrt(h_l)*(y_0 + h_l*(slope/2.0)*(pow(v+1, 2)- pow(v, 2))); + return sqrt_h_l*(y_0 + h_l*(slope/2.0)*(pow(v+1, 2)- pow(v, 2))); } else { // p == 1 --> legendre(1) = sqrt(12)(x-(v+0.5)*h_l)/(h_l^1.5) - return SQRT_12*h_l*((slope/3.0)*(pow(v+1, 3)- pow(v, 3))*sqrt(h_l) + ((y_0 / h_l -1*slope*(v+0.5))/2.0)*(pow(v+1, 2)- pow(v, 2)) - y_0*(v+0.5)*h_l); + return SQRT_12*h_l*((slope/3.0)*(pow(v+1, 3)- pow(v, 3))*sqrt_h_l + sqrt_h_l*((y_0 / h_l -1*slope*(v+0.5))/2.0)*(pow(v+1, 2)- pow(v, 2)) - y_0*(v+0.5)/sqrt_h_l); } } From e32a2e1750020c3105f005028ad13a702b0d1828 Mon Sep 17 00:00:00 2001 From: Stefan Date: Wed, 12 May 2010 15:13:00 +0200 Subject: [PATCH 07/24] Simplification of coefficient integration. --- source/discont.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/discont.c b/source/discont.c index 14c2fc9..bdf5df3 100755 --- a/source/discont.c +++ b/source/discont.c @@ -209,9 +209,9 @@ integrate_coeff_discont(discont_function_p function, int position, int v, int p, sqrt_h_l = sqrt(h_l); if (p == 0) { // legendre(0) = sqrt(1/h_l) - return sqrt_h_l*(y_0 + h_l*(slope/2.0)*(pow(v+1, 2)- pow(v, 2))); + return sqrt_h_l*(y_0 + h_l*(slope/2.0)); } else { // p == 1 --> legendre(1) = sqrt(12)(x-(v+0.5)*h_l)/(h_l^1.5) - return SQRT_12*h_l*((slope/3.0)*(pow(v+1, 3)- pow(v, 3))*sqrt_h_l + sqrt_h_l*((y_0 / h_l -1*slope*(v+0.5))/2.0)*(pow(v+1, 2)- pow(v, 2)) - y_0*(v+0.5)/sqrt_h_l); + return h_l*sqrt_h_l*slope/SQRT_12; } } From ae1da7677777c3ba400c87e298d555273c016882 Mon Sep 17 00:00:00 2001 From: Stefan Date: Wed, 12 May 2010 15:48:24 +0200 Subject: [PATCH 08/24] Fixes deallocation and completes example. --- source/discont.c | 8 +++++++- source/fft_faltung.c | 14 -------------- test/discont_test.c | 35 +++++++++++++++++++++++++++++------ 3 files changed, 36 insertions(+), 21 deletions(-) diff --git a/source/discont.c b/source/discont.c index bdf5df3..7ad0742 100755 --- a/source/discont.c +++ b/source/discont.c @@ -86,7 +86,9 @@ linear_function_set_del(linear_function_set_p function_set) { int n; for (n = 0; n < function_set->count; n++) { - free(function_set->functions[n]); + if (function_set->functions[n] != NULL) { + free(function_set->functions[n]); + } } free(function_set->functions); free(function_set); @@ -118,6 +120,10 @@ discont_function_setup_points(discont_function_p function, int step, fepc_real_t for (n = 0; n < count; n++) { function->function_sets[step]->functions[n] = linear_function_new_points(y1[n], y2[n], h_l); } + } else { + for (n = 0; n < count; n++) { + function->function_sets[step]->functions[n] = NULL; + } } } diff --git a/source/fft_faltung.c b/source/fft_faltung.c index baaad2a..c481075 100755 --- a/source/fft_faltung.c +++ b/source/fft_faltung.c @@ -16,11 +16,7 @@ * along with this program. If not, see . */ -#if defined(HAS_FFTW3) #include - -#endif - #include "fft_faltung.h" @@ -41,8 +37,6 @@ *******************************************************/ 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; @@ -151,13 +145,5 @@ fft_faltung(fepc_real_t* a, vec_p n_a, fepc_real_t* b, vec_p n_b) { fftw_free(out_a); fftw_free(out_b); return c; - - - -#else - printf( "\n (fft_faltung) FEHLER : keine FFT Bibliothek verfuegbar\n" ); - - exit( 1 ); -#endif } diff --git a/test/discont_test.c b/test/discont_test.c index 4018134..48f5655 100755 --- a/test/discont_test.c +++ b/test/discont_test.c @@ -41,14 +41,13 @@ main(void) { // 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[] = {3.4, 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 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); @@ -58,12 +57,34 @@ main(void) { // Create f2 in the same way! + + function = discont_function_new(steps); + + fepc_real_t y13[] = {0.1, 0.0, 0.3, 3.0, 4.2, 0.3, 5.8, 7.3, 1.0, 2.4}; + fepc_real_t y23[] = {0.2, 1.3, 0.1, 6.5, 4.2, 3.3, 4.1, 0.3, 9.3, 2.2}; + + // 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); @@ -75,13 +96,15 @@ main(void) { func_del(f2); func_del(w); + func_print(convolution_result, 3); + result = convert_func(convolution_result, function->intervals, stepping); - discont_function_del(function); + discont_function_del(function); func_del(convolution_result); // the result is stored in "result" as a discont-function - + //discont_function_print(result); return 0; } From 20c688e1b6e89cdfcbcd82795d15795d7f71d526 Mon Sep 17 00:00:00 2001 From: Stefan Date: Wed, 12 May 2010 16:54:59 +0200 Subject: [PATCH 09/24] Fixes example to not throw a segmentation fault. --- test/discont_test.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/discont_test.c b/test/discont_test.c index 48f5655..89c8660 100755 --- a/test/discont_test.c +++ b/test/discont_test.c @@ -86,6 +86,7 @@ main(void) { 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_gridstructure(w, function->intervals, stepping); @@ -98,13 +99,12 @@ main(void) { func_print(convolution_result, 3); - result = convert_func(convolution_result, function->intervals, stepping); + result = convert_func(convolution_result, intervals_clone(function->intervals, function->stepcount), stepping); discont_function_del(function); func_del(convolution_result); - // the result is stored in "result" as a discont-function - //discont_function_print(result); + discont_function_print(result); return 0; } From f4302c506d483f5b0d626d5a5621239cba105dce Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Thu, 13 May 2010 18:18:35 +0200 Subject: [PATCH 10/24] Fixes selection of x-values. --- source/discont.c | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/source/discont.c b/source/discont.c index 7ad0742..f20facc 100755 --- a/source/discont.c +++ b/source/discont.c @@ -171,7 +171,7 @@ convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { fepc_real_t h_l, temp_x1, temp_x2, temp_y1, temp_y2, slope; - int n, k, stepcount; + int n, k, stepcount, start; vec_real_p x; @@ -180,15 +180,12 @@ convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { result->intervals[n] = intervals[n]; stepcount = _round((intervals[n]->end[0] - intervals[n]->start[0])/h_l); result->function_sets[n] = linear_function_set_new(stepcount); + start = function->hierarchie[n]->vektor[0]->start->array[0]; for (k = 0; k < stepcount; k++) { - /* - to make sure we are are using the correct linear function, we use the points 1/3 and 2/3 instead of 0 and 1 to - calculate the linear function - */ - temp_x1 = (k + ONE_THIRD)*h_l; - temp_x2 = (k + TWO_THIRD)*h_l; + temp_x1 = (start+k)*h_l; + temp_x2 = (1+start+k)*h_l; x = vec_real_new(1); - x->array[0] = temp_x2; + x->array[0] = temp_x1; temp_y1 = get_value_at_step(function, x, n, stepping); vec_real_del(x); x = vec_real_new(1); From c6b0e70d864c0cfd82eb191571c1c862d13e03d6 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Fri, 14 May 2010 02:18:06 +0200 Subject: [PATCH 11/24] Changes function handling. --- include/discont.h | 2 +- source/discont.c | 23 ++++++++++------------- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/include/discont.h b/include/discont.h index 3cb38ea..5a03ab3 100755 --- a/include/discont.h +++ b/include/discont.h @@ -59,7 +59,7 @@ 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 x_0, fepc_real_t x_1, fepc_real_t dx); +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); diff --git a/source/discont.c b/source/discont.c index f20facc..7624f8f 100755 --- a/source/discont.c +++ b/source/discont.c @@ -39,12 +39,10 @@ linear_function_new(fepc_real_t y_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 dx) { - linear_function_p result = (linear_function_p) malloc(sizeof(linear_function_t)); - - result->y_0 = y_0; - result->slope = (y_1 - y_0)/dx; - return result; +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_1-slope*x_0, slope); } @@ -85,6 +83,7 @@ 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]); @@ -118,7 +117,7 @@ discont_function_setup_points(discont_function_p function, int step, fepc_real_t 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], h_l); + 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++) { @@ -182,18 +181,16 @@ convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { 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 = (start+k)*h_l; - temp_x2 = (1+start+k)*h_l; + 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); - vec_real_del(x); - x = vec_real_new(1); 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); - result->function_sets[n]->functions[k] = linear_function_new(slope*k*h_l+temp_y1, slope); + 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; From 61cb35c77fcc1ec2947f29f5807e1c9b101a787f Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Fri, 14 May 2010 02:29:41 +0200 Subject: [PATCH 12/24] Fixes typo. --- source/discont.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/discont.c b/source/discont.c index 7624f8f..212e8c3 100755 --- a/source/discont.c +++ b/source/discont.c @@ -42,7 +42,7 @@ 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_1-slope*x_0, slope); + return linear_function_new(y_0-slope*x_0, slope); } From ab8c059315da66d12a24002020e7def6e3588b42 Mon Sep 17 00:00:00 2001 From: Stefan Date: Fri, 14 May 2010 19:41:27 +0200 Subject: [PATCH 13/24] Fixes coeff calculation (again). Fixes numerical integration. --- include/discont.h | 2 +- include/fepc_easy.h | 2 +- source/discont.c | 45 +++++++++++++++++++++------------------------ source/fepc_easy.c | 23 ++++++++--------------- test/discont_test.c | 41 ++++++++++++++++++++++++++++------------- 5 files changed, 59 insertions(+), 54 deletions(-) diff --git a/include/discont.h b/include/discont.h index 5a03ab3..374b56c 100755 --- a/include/discont.h +++ b/include/discont.h @@ -86,7 +86,7 @@ void linear_function_set_print(linear_function_set_p function_set); fepc_real_t -integrate_coeff_discont(discont_function_p function, int position, int v, int p, int step, fepc_real_t stepping); +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); diff --git a/include/fepc_easy.h b/include/fepc_easy.h index d3af7ac..edf4a01 100755 --- a/include/fepc_easy.h +++ b/include/fepc_easy.h @@ -26,7 +26,7 @@ #include "interval.h" #ifndef INT_STEPS -#define INT_STEPS 20 +#define INT_STEPS 10 #endif #define SQRT_3 1.7320508075688772935274463415058723669428 diff --git a/source/discont.c b/source/discont.c index 212e8c3..81b7e0c 100755 --- a/source/discont.c +++ b/source/discont.c @@ -177,7 +177,7 @@ convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { 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); + 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++) { @@ -198,43 +198,40 @@ convert_func(func_p function, interval_p * intervals, fepc_real_t stepping) { fepc_real_t -integrate_coeff_discont(discont_function_p function, int position, int v, int p, int step, fepc_real_t stepping) { +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, sqrt_h_l; - h_l = get_h_l(step, stepping); + fepc_real_t h_l, slope, y_0; + h_l = get_h_l(level, stepping); - slope = function->function_sets[step]->functions[position]->slope; - y_0 = function->function_sets[step]->functions[position]->y_0; - - sqrt_h_l = sqrt(h_l); + 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/2.0)); + 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; + 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; - int step, n, k, interval_element_count, position, grad_count, pos_2; - - vec_p v, x; + vec_p v, p; - for (step = 0; step <= function->maxlevel; step++) { - grad_count = get_degree_count(function->hierarchie[step]->grad); + 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[step]->vektor[k]); n < interval_element_count; n++) { - v = generate_folgenvector(function->hierarchie[step]->vektor[k], n); - pos_2 = v->array[0]; - 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 - x = entry_one2d_sloppy(k, function->hierarchie[step]->grad); // don't care about 0 as entries - function->hierarchie[step]->vektor[k]->glied[position] = integrate_coeff_discont(discont_function, pos_2 ,v->array[0], x->array[0], step, stepping); - vec_del(x); + 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/source/fepc_easy.c b/source/fepc_easy.c index 17b662e..bcf9eb3 100755 --- a/source/fepc_easy.c +++ b/source/fepc_easy.c @@ -143,7 +143,7 @@ fepc_real_t phi_l(int step, vec_p v, vec_p p, vec_real_p x, fepc_real_t stepping 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*x->array[n]/h_l - 2*v->array[n] - 1); + 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; } @@ -177,7 +177,6 @@ vec_real_set_p get_value(func_p function, Funcimpl_vec_step generate_points, fep 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); @@ -213,13 +212,12 @@ fepc_real_t get_value_at_step(func_p function, vec_real_p x, int step, fepc_real 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); } } @@ -493,7 +491,7 @@ vec_real_p generate_x_for_integration_st(vec_p v, int calc_position, fepc_real_t 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; + result->array[n] = v->array[n]*h_l+(2.0*current +1)*h/2.0; } return result; } @@ -543,9 +541,7 @@ fepc_real_t integrate_coeff_st(Funcimpl function_impl, vec_p v, vec_p p, int ste h_l = get_h_l(step, stepping); h = h_l/(INT_STEPS); // step size - result = 0; - - + result = 0.0; // interval length is always h_l! (by design) calc_count = pow(INT_STEPS, v->dim); @@ -564,8 +560,6 @@ fepc_real_t integrate_coeff_st(Funcimpl function_impl, vec_p v, vec_p p, int ste vec_real_del(x); } - //printf("sum = %f, result = %f ", result, pow(h_l, v->dim)*result/calc_count); - return pow(h_l, v->dim)*result / pow(INT_STEPS-1, v->dim); } @@ -641,7 +635,7 @@ void add_folgenentries(func_p function, Funcimpl function_impl, Funcimpl_step co int step, n, k, interval_element_count, position, grad_count; - vec_p v, x; + vec_p v, p; for (step = 0; step <= function->maxlevel; step++) { grad_count = get_degree_count(function->hierarchie[step]->grad); @@ -651,11 +645,10 @@ void add_folgenentries(func_p function, Funcimpl function_impl, Funcimpl_step co 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 - x = entry_one2d_sloppy(k, function->hierarchie[step]->grad); // don't care about 0 as entries - + 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, x, step, stepping) : integrate_coeff_st(function_impl, v, x, step, stepping); - vec_del(x); + 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); } diff --git a/test/discont_test.c b/test/discont_test.c index 89c8660..e36020b 100755 --- a/test/discont_test.c +++ b/test/discont_test.c @@ -17,8 +17,13 @@ */ #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 */ @@ -29,14 +34,14 @@ main(void) { func_p f1, f2, w, convolution_result; - steps = 2; + steps = 1; discont_function_p function, result; function = discont_function_new(steps); - fepc_real_t y11[] = {0.1, 0.0, 0.3, 3.0, 4.2, 0.3, 5.8, 7.3, 1.0, 2.4}; // the stepcount is 10 in the first level because h_0 is 0.1 - fepc_real_t y21[] = {0.2, 1.3, 0.1, 6.5, 4.2, 3.3, 4.1, 0.3, 9.3, 2.2}; + 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 @@ -46,22 +51,27 @@ main(void) { 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_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); - discont_function_del(function); + func_print(f1, 3); + func_print(setup_fepc_structure(x, function->intervals, 1, 1, stepping), 3); + + vec_real_set_p points = get_value(setup_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[] = {0.1, 0.0, 0.3, 3.0, 4.2, 0.3, 5.8, 7.3, 1.0, 2.4}; - fepc_real_t y23[] = {0.2, 1.3, 0.1, 6.5, 4.2, 3.3, 4.1, 0.3, 9.3, 2.2}; + 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 @@ -70,7 +80,7 @@ main(void) { 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_setup_points(function, 1, 0.2, 0.5, y14, y24, stepping); // interval [0.2, 0.5] at level 1 discont_function_print(function); @@ -85,23 +95,28 @@ main(void) { 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 + //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_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); + //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 From 787c6578714aea2b08c248e8d0e77dd97f808a6a Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Wed, 15 Dec 2010 11:21:07 +0100 Subject: [PATCH 14/24] Adds cp structure and its convolution. --- include/fepc_cp.h | 53 ++++++++++++++++++++++++ include/funktion.h | 5 ++- source/CMakeLists.txt | 8 ++++ source/fepc_cp.c | 96 +++++++++++++++++++++++++++++++++++++++++++ source/folge.c | 2 +- source/funktion.c | 9 ++++ test/CMakeLists.txt | 16 +++++++- test/discont_test.c | 1 - test/fepc_cp_test.c | 77 ++++++++++++++++++++++++++++++++++ 9 files changed, 263 insertions(+), 4 deletions(-) create mode 100644 include/fepc_cp.h create mode 100644 source/fepc_cp.c create mode 100644 test/fepc_cp_test.c diff --git a/include/fepc_cp.h b/include/fepc_cp.h new file mode 100644 index 0000000..2cdef3b --- /dev/null +++ b/include/fepc_cp.h @@ -0,0 +1,53 @@ +/* + * 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 . + */ + +#ifndef __FEPCCP +#define __FEPCCP + + +#include "fepc_easy.h" + + +typedef struct { + int rank; + int dimension; + func_p * functions; +} func_cp_t; + +typedef func_cp_t * func_cp_p; + +func_cp_p +func_cp_new(int rank, int dimension); + +func_p +func_cp_extract(int current_rank, int current_dimension, func_cp_p func_cp); + +void +func_cp_del(func_cp_p func_cp); + +func_cp_p +func_cp_faltung(func_cp_p function1, func_cp_p function2, func_p * resulting_structure, fepc_real_t h); + +func_cp_p +func_cp_multi(func_cp_p function1, func_cp_p function2, fepc_real_t h); + +void +func_cp_print(func_cp_p function); + + +#endif // __FEPCCP diff --git a/include/funktion.h b/include/funktion.h index 669174d..7dae755 100755 --- a/include/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 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 @@ -114,5 +114,8 @@ 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); + #endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 90fed2e..2583465 100755 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -3,4 +3,12 @@ file(GLOB ALL_SOURCE_FILES *.c) # Find all header files in the include directory file(GLOB ALL_HEADER_FILES ${CMAKE_SOURCE_DIR}/include/*.h) +LIST(REMOVE_ITEM ALL_SOURCE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/gauss.c) +LIST(REMOVE_ITEM ALL_HEADER_FILES ${CMAKE_SOURCE_DIR}/include/gauss.h) + ADD_LIBRARY(fepc SHARED ${ALL_SOURCE_FILES} ${ALL_HEADER_FILES}) + +if (EXISTS ${CMAKE_SOURCE_DIR}/include/gauss.h) + ADD_LIBRARY(gauss SHARED gauss.c ${CMAKE_SOURCE_DIR}/include/gauss.h) +endif() + diff --git a/source/fepc_cp.c b/source/fepc_cp.c new file mode 100644 index 0000000..1d94ba9 --- /dev/null +++ b/source/fepc_cp.c @@ -0,0 +1,96 @@ +/* + * 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_p +func_cp_new(int rank, int dimension) { + func_cp_p result = (func_cp_p) malloc(sizeof(func_cp_t)); + + result->rank = rank; + result->dimension = dimension; + result->functions = (func_p*) malloc(sizeof(func_p)*rank*dimension); + return result; +} + +func_p +func_cp_extract(int current_rank, int current_dimension, func_cp_p 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_p +func_cp_faltung(func_cp_p function1, func_cp_p function2, func_p * resulting_structure, fepc_real_t h) { + ASSERT(function1->dimension == function2->dimension); + + func_cp_p result = func_cp_new(function1->rank*function2->rank, function1->dimension); + + 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++) { + result->functions[(n*function2->rank + k)*result->dimension+d] = faltung_fepc(func_cp_extract(n, d, function1), func_cp_extract(k, d, function2), resulting_structure[d], h); + } + } + } + return result; +} + +void +func_cp_del(func_cp_p func_cp) { + funcs_del(func_cp->functions, func_cp->dimension*func_cp->rank); + free(func_cp); +} + +func_cp_p +func_cp_multi(func_cp_p function1, func_cp_p function2, fepc_real_t h) { + ASSERT(function1->dimension == function2->dimension); + + func_cp_p result = func_cp_new(function1->rank*function2->rank, function1->dimension); + + 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++) { + result->functions[(n*function2->rank + k)*result->dimension+d] = func_multi(func_cp_extract(n, d, function1), func_cp_extract(k, d, function2), h); + } + } + } + return result; +} + +void +func_cp_print(func_cp_p 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"); + } + } +} diff --git a/source/folge.c b/source/folge.c index e81bd95..f22265a 100755 --- a/source/folge.c +++ b/source/folge.c @@ -243,7 +243,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"); diff --git a/source/funktion.c b/source/funktion.c index a889822..b5aafcd 100755 --- a/source/funktion.c +++ b/source/funktion.c @@ -406,3 +406,12 @@ func_factor_multi(func_p function, fepc_real_t factor) { return back; } +void +funcs_del(func_p * array, int length) { + int n; + + for (n = 0; n < length; n++) { + func_del(array[n]); + } + free(array); +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f9fb84b..07df6ad 100755 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,7 +17,21 @@ add_executable (sparsetest sparsetest.c) target_link_libraries (sparsetest fepc fftw3 m) +if (EXISTS ${CMAKE_SOURCE_DIR}/include/gauss.h) + add_executable (gauss_test gauss_test.c) + target_link_libraries (gauss_test gauss blas lapack gsl m) + add_test(gauss_test gauss_test) +endif() + + +add_executable(fepc_cp_test fepc_cp_test.c) +target_link_libraries(fepc_cp_test fepc fftw3 m) + + add_test(fepc_easy_test fepc_easy_test) add_test(stest stest) add_test(discont_test discont_test) -add_test(sparsetest sparsetest) \ No newline at end of file +add_test(sparsetest sparsetest) + +add_test(fepc_cp_test fepc_cp_test) + diff --git a/test/discont_test.c b/test/discont_test.c index e36020b..2c43042 100755 --- a/test/discont_test.c +++ b/test/discont_test.c @@ -103,7 +103,6 @@ main(void) { set_gridstructure(w, function->intervals, stepping); convolution_result = faltung_fepc(f1, f2, w, stepping); - //func_del(f1); //func_del(f2); //func_del(w); diff --git a/test/fepc_cp_test.c b/test/fepc_cp_test.c new file mode 100644 index 0000000..e632ab9 --- /dev/null +++ b/test/fepc_cp_test.c @@ -0,0 +1,77 @@ +/* + * 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_p * intervals = (interval_p*) malloc(sizeof(interval_p)*1); + + intervals[0] = interval_new(1); + intervals[0]->start[0] = -1.; + intervals[0]->end[0] = 1.; + + int rank1 = 3, rank2 = 4, dimension = 2; + + int n, k; + + func_cp_p cp1 = func_cp_new(rank1, dimension); + + for (n = 0; n < rank1; n++) { + for (k = 0; k < dimension; k++) { + cp1->functions[n*dimension+k] = setup_fepc_structure(id, intervals, 1, 0, h); + } + } + + func_cp_p cp2 = func_cp_new(rank2, dimension); + + for (n = 0; n < rank2; n++) { + for (k = 0; k < dimension; k++) { + cp2->functions[n*dimension+k] = setup_fepc_structure(id, intervals, 1, 0, h); + } + } + + func_p * result_intervals = (func_p*) malloc(sizeof(func_p)*dimension); + + for (k = 0; k < dimension; k++) { + result_intervals[k] = func_new(0, 1); + set_gridstructure(result_intervals[k], intervals, h); + } + + func_cp_p result = func_cp_faltung(cp1, cp2, result_intervals, h); + + + func_cp_del(cp1); + func_cp_del(cp2); + + func_cp_print(result); + + funcs_del(result_intervals, dimension); + intervals_del(intervals, 1); + func_cp_del(result); + + + return 0; +} From 20b89c7f57965e6c54935aa7fe2b0572252d9095 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Thu, 16 Dec 2010 12:32:36 +0100 Subject: [PATCH 15/24] Increases memory efficiency. --- include/faltung.h | 3 ++ include/fepc_cp.h | 43 +++++++++++++------ include/fepc_easy.h | 6 ++- include/fepc_easy_helper.h | 3 ++ include/funktion.h | 12 ++++++ include/interval.h | 9 ++++ source/faltung.c | 23 +++++----- source/fepc_cp.c | 86 ++++++++++++++++++++++++++++--------- source/fepc_easy.c | 38 +++++++++++------ source/fepc_easy_helper.c | 40 +++++++++++------- source/funktion.c | 87 +++++++++++++++++++++++--------------- source/interval.c | 41 +++++++++++++----- test/discont_test.c | 4 +- test/fepc_cp_test.c | 37 ++++++++-------- test/fepc_easy_test.c | 6 +-- test/stest.c | 4 +- 16 files changed, 296 insertions(+), 146 deletions(-) diff --git a/include/faltung.h b/include/faltung.h index b8ced7e..852c262 100755 --- a/include/faltung.h +++ b/include/faltung.h @@ -37,4 +37,7 @@ 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); + #endif diff --git a/include/fepc_cp.h b/include/fepc_cp.h index 2cdef3b..1aca062 100644 --- a/include/fepc_cp.h +++ b/include/fepc_cp.h @@ -26,28 +26,45 @@ typedef struct { int rank; int dimension; - func_p * functions; -} func_cp_t; + func_t * functions; +} func_cp; -typedef func_cp_t * func_cp_p; +func_cp * +func_cp_new(int rank, int dimension, int * maxlevels); -func_cp_p -func_cp_new(int rank, int dimension); +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); -func_p -func_cp_extract(int current_rank, int current_dimension, func_cp_p func_cp); +/* + * 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_p func_cp); +func_cp_del(func_cp * func_cp); -func_cp_p -func_cp_faltung(func_cp_p function1, func_cp_p function2, func_p * resulting_structure, fepc_real_t h); +func_cp * +func_cp_faltung(func_cp * function1, func_cp * function2, func_t * resulting_structure, fepc_real_t h); -func_cp_p -func_cp_multi(func_cp_p function1, func_cp_p function2, 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_p function); +func_cp_print(func_cp * function); + +int * +int_array_new(int length); #endif // __FEPCCP diff --git a/include/fepc_easy.h b/include/fepc_easy.h index edf4a01..ee1c113 100755 --- a/include/fepc_easy.h +++ b/include/fepc_easy.h @@ -120,7 +120,11 @@ fepc_real_t phi_l(int step, vec_p v, vec_p p, vec_real_p x, fepc_real_t stepping */ void add_folgenentries(func_p function, Funcimpl function_impl, Funcimpl_step coeff_function, fepc_real_t stepping); -func_p setup_fepc_structure(Funcimpl function, interval_p* intervals, int interval_count, int degree, 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); diff --git a/include/fepc_easy_helper.h b/include/fepc_easy_helper.h index a27cca1..37c24d5 100755 --- a/include/fepc_easy_helper.h +++ b/include/fepc_easy_helper.h @@ -44,6 +44,9 @@ 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); diff --git a/include/funktion.h b/include/funktion.h index 7dae755..ff3d745 100755 --- a/include/funktion.h +++ b/include/funktion.h @@ -48,9 +48,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. @@ -91,6 +97,9 @@ 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 @@ -117,5 +126,8 @@ 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); + #endif diff --git a/include/interval.h b/include/interval.h index 6b166ab..e47a3ff 100755 --- a/include/interval.h +++ b/include/interval.h @@ -17,12 +17,21 @@ typedef interval_t * interval_p; 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. */ diff --git a/source/faltung.c b/source/faltung.c index 0744024..f8df10c 100755 --- a/source/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,38 +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/source/fepc_cp.c b/source/fepc_cp.c index 1d94ba9..acfa639 100644 --- a/source/fepc_cp.c +++ b/source/fepc_cp.c @@ -19,67 +19,99 @@ #include "fepc_cp.h" #include "fepc_easy_helper.h" -func_cp_p -func_cp_new(int rank, int dimension) { - func_cp_p result = (func_cp_p) malloc(sizeof(func_cp_t)); +func_cp * +func_cp_new(int rank, int dimension, int * maxlevels) { + func_cp * result = (func_cp*) malloc(sizeof(func_cp)); - result->rank = rank; - result->dimension = dimension; - result->functions = (func_p*) malloc(sizeof(func_p)*rank*dimension); + func_cp_init(result, rank, dimension, maxlevels); return result; } -func_p -func_cp_extract(int current_rank, int current_dimension, func_cp_p func_cp) { +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); + + + 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]; + return &(func_cp->functions[current_rank*func_cp->dimension+current_dimension]); } -func_cp_p -func_cp_faltung(func_cp_p function1, func_cp_p function2, func_p * resulting_structure, fepc_real_t h) { +func_cp * +func_cp_faltung(func_cp * function1, func_cp * function2, func_t * resulting_structure, fepc_real_t h) { ASSERT(function1->dimension == function2->dimension); - func_cp_p result = func_cp_new(function1->rank*function2->rank, function1->dimension); + int * maxlevels = int_array_new(100); // FIXME + + 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++) { - result->functions[(n*function2->rank + k)*result->dimension+d] = faltung_fepc(func_cp_extract(n, d, function1), func_cp_extract(k, d, function2), resulting_structure[d], h); + 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); } } } + free(maxlevels); return result; } void -func_cp_del(func_cp_p func_cp) { - funcs_del(func_cp->functions, func_cp->dimension*func_cp->rank); +func_cp_del(func_cp * func_cp) { + funcs_del_type(func_cp->functions, func_cp->dimension*func_cp->rank); free(func_cp); } -func_cp_p -func_cp_multi(func_cp_p function1, func_cp_p function2, fepc_real_t h) { +func_cp * +func_cp_multi(func_cp * function1, func_cp * function2, fepc_real_t h) { ASSERT(function1->dimension == function2->dimension); - func_cp_p result = func_cp_new(function1->rank*function2->rank, function1->dimension); + int * maxlevels = int_array_new(100); // FIXME + + 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++) { - result->functions[(n*function2->rank + k)*result->dimension+d] = func_multi(func_cp_extract(n, d, function1), func_cp_extract(k, d, function2), h); + 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_p cp) { +func_cp_print(func_cp * cp) { int n, k; printf("CP function with rank %i:\n", cp->rank); @@ -94,3 +126,17 @@ func_cp_print(func_cp_p cp) { } } } + +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 index bcf9eb3..a4f4586 100755 --- a/source/fepc_easy.c +++ b/source/fepc_easy.c @@ -280,17 +280,19 @@ void set_gridstructure(func_p function, interval_p* intervals, fepc_real_t stepp } folge = folge_new(start, lang); - for (k = 0, grad_count = get_degree_count(function->hierarchie[n]->grad); k < grad_count; k++) { + 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 - folge_del(folge); } } @@ -680,18 +682,26 @@ void set_degree(func_p function, int degree) { -func_p setup_fepc_structure(Funcimpl function, interval_p* intervals, int interval_count, int degree, fepc_real_t stepping) { - func_p result = func_new(interval_count-1, intervals[interval_count-1]->dimension); - - if (degree > 0) { - set_degree(result, degree); +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(result, intervals, stepping); + set_gridstructure(func, intervals, stepping); if (function != NULL) { - add_folgenentries(result, function, NULL, stepping); + add_folgenentries(func, function, NULL, stepping); } - return result; +} + +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) { @@ -840,9 +850,9 @@ func_p fepc_easy(Funcimpl function1, interval_p* intervals1, int interval_count1 func_p f, g, w, fepc; - f = setup_fepc_structure(function1, intervals1, interval_count1, degree_func1, stepping); - g = setup_fepc_structure(function2, intervals2, interval_count2, degree_func2, stepping); - w = setup_fepc_structure(NULL, intervals3, interval_count3, 0, stepping); + 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); diff --git a/source/fepc_easy_helper.c b/source/fepc_easy_helper.c index 5cd4fdd..95421a5 100755 --- a/source/fepc_easy_helper.c +++ b/source/fepc_easy_helper.c @@ -22,25 +22,33 @@ func_p func_multi(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_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 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); diff --git a/source/funktion.c b/source/funktion.c index b5aafcd..f29e3d4 100755 --- a/source/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,27 +61,29 @@ 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) { - int k; - - for(k=0;k<=f->maxlevel;k++) { - folgen_vektor_del(f->hierarchie[k]); - } - free(f->hierarchie); + 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); +} + @@ -197,35 +204,35 @@ 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; } @@ -415,3 +422,13 @@ funcs_del(func_p * array, int length) { } 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 index 91111ae..9951130 100755 --- a/source/interval.c +++ b/source/interval.c @@ -23,24 +23,34 @@ interval_p interval_new(int dimension) { interval_p result = (interval_p) malloc(sizeof(interval_t)); - result->dimension = dimension; - result->start = (fepc_real_t*) malloc(sizeof(fepc_real_t)*dimension); - result->end = (fepc_real_t*) malloc(sizeof(fepc_real_t)*dimension); - - int n; - - for (n = 0; n < dimension; n++) { - result->start[n] = 0.; - } + 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) { - free(interval->start); - free(interval->end); + 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; @@ -50,6 +60,15 @@ void intervals_del(interval_p* intervals, int interval_count) { 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; diff --git a/test/discont_test.c b/test/discont_test.c index 2c43042..9b39088 100755 --- a/test/discont_test.c +++ b/test/discont_test.c @@ -59,9 +59,9 @@ main(void) { func_print(f1, 3); - func_print(setup_fepc_structure(x, function->intervals, 1, 1, stepping), 3); + func_print(create_fepc_structure(x, function->intervals, 1, 1, stepping), 3); - vec_real_set_p points = get_value(setup_fepc_structure(x, function->intervals, 1, 1, stepping), get_mean_points, stepping); + 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); diff --git a/test/fepc_cp_test.c b/test/fepc_cp_test.c index e632ab9..076ed60 100644 --- a/test/fepc_cp_test.c +++ b/test/fepc_cp_test.c @@ -27,40 +27,42 @@ int main() { fepc_real_t h = 0.2; - interval_p * intervals = (interval_p*) malloc(sizeof(interval_p)*1); + interval_t * intervals = (interval_t*) malloc(sizeof(interval_t)*1); - intervals[0] = interval_new(1); - intervals[0]->start[0] = -1.; - intervals[0]->end[0] = 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; + int n, k, interval_count = 1; - func_cp_p cp1 = func_cp_new(rank1, dimension); + int * maxlevels = int_array_new(10); + + func_cp * cp1 = func_cp_new(rank1, dimension, maxlevels); for (n = 0; n < rank1; n++) { for (k = 0; k < dimension; k++) { - cp1->functions[n*dimension+k] = setup_fepc_structure(id, intervals, 1, 0, h); + setup_fepc_structure(&(cp1->functions[n*dimension+k]), id, &intervals, 1, 0, h); } } - func_cp_p cp2 = func_cp_new(rank2, dimension); - + func_cp * cp2 = func_cp_new(rank2, dimension, maxlevels); + free(maxlevels); for (n = 0; n < rank2; n++) { for (k = 0; k < dimension; k++) { - cp2->functions[n*dimension+k] = setup_fepc_structure(id, intervals, 1, 0, h); + setup_fepc_structure(&(cp2->functions[n*dimension+k]), id, &intervals, 1, 0, h); } } - func_p * result_intervals = (func_p*) malloc(sizeof(func_p)*dimension); + func_t * result_intervals = (func_t*) malloc(sizeof(func_t)*dimension); for (k = 0; k < dimension; k++) { - result_intervals[k] = func_new(0, 1); - set_gridstructure(result_intervals[k], intervals, h); + func_init(&result_intervals[k], 0, 1); + set_gridstructure(&result_intervals[k], &intervals, h); } - func_cp_p result = func_cp_faltung(cp1, cp2, result_intervals, h); + func_cp * result = func_cp_faltung(cp1, cp2, result_intervals, h); func_cp_del(cp1); @@ -68,10 +70,9 @@ main() { func_cp_print(result); - funcs_del(result_intervals, dimension); - intervals_del(intervals, 1); + funcs_del_type(result_intervals, dimension); + intervals_del_type(intervals, 1); func_cp_del(result); - - + fftw_cleanup(); return 0; } diff --git a/test/fepc_easy_test.c b/test/fepc_easy_test.c index bc99334..9985e8b 100755 --- a/test/fepc_easy_test.c +++ b/test/fepc_easy_test.c @@ -481,8 +481,8 @@ void test_multi() { intervals[1]->end[0] = 0.4; intervals[1]->end[1] = 1.; - f = setup_fepc_structure(f_peter, intervals, steps_f, grad_f, STEPPING); - g = setup_fepc_structure(g_peter, intervals, steps_g, grad_g, STEPPING); + 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); @@ -506,7 +506,7 @@ void test_norm() { intervals[0]->end[0] = 1.; intervals[0]->end[1] = 1.; - func_p function = setup_fepc_structure(const_1, intervals, 1, 0, STEPPING); + func_p function = create_fepc_structure(const_1, intervals, 1, 0, STEPPING); //func_print(function, 3); diff --git a/test/stest.c b/test/stest.c index 99d56b0..8884bf2 100755 --- a/test/stest.c +++ b/test/stest.c @@ -76,10 +76,10 @@ int main() { //func_p func_9x1_sqr = setup_fepc_structure(nine_xone_sqr, intervals, 1, 0, stepping); - func_p func_function = setup_fepc_structure(function, 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 = setup_fepc_structure(function_laplace, intervals, 1, 0, stepping); + 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"); From 54def9de1e2c3488e16826035ad5767d34e288ac Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Thu, 16 Dec 2010 14:06:20 +0100 Subject: [PATCH 16/24] Sets correct array size --- source/fepc_cp.c | 20 +++++++++++++++----- test/fepc_cp_test.c | 12 ++++++++---- 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/source/fepc_cp.c b/source/fepc_cp.c index acfa639..3b4d9d9 100644 --- a/source/fepc_cp.c +++ b/source/fepc_cp.c @@ -43,7 +43,7 @@ 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; } @@ -66,11 +66,19 @@ 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 * maxlevels = int_array_new(100); // FIXME + 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); - int n, k, d; + free(maxlevels); + + int k, d; for (n = 0; n < function1->rank; n++) { for (k = 0; k < function2->rank; k++) { @@ -79,7 +87,7 @@ func_cp_faltung(func_cp * function1, func_cp * function2, func_t * resulting_str } } } - free(maxlevels); + return result; } @@ -93,7 +101,9 @@ func_cp * func_cp_multi(func_cp * function1, func_cp * function2, fepc_real_t h) { ASSERT(function1->dimension == function2->dimension); - int * maxlevels = int_array_new(100); // FIXME + 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); diff --git a/test/fepc_cp_test.c b/test/fepc_cp_test.c index 076ed60..c16dbe9 100644 --- a/test/fepc_cp_test.c +++ b/test/fepc_cp_test.c @@ -37,9 +37,11 @@ main() { int n, k, interval_count = 1; - int * maxlevels = int_array_new(10); + int * maxlevels_cp1 = int_array_new(rank1*dimension); - func_cp * cp1 = func_cp_new(rank1, dimension, maxlevels); + 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++) { @@ -47,8 +49,10 @@ main() { } } - func_cp * cp2 = func_cp_new(rank2, dimension, maxlevels); - free(maxlevels); + 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); From eda1f8ce9765b14a6781ba7c410b8862c2cbe2cc Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Mon, 21 Feb 2011 00:31:11 +0100 Subject: [PATCH 17/24] Adds c++ compatibility to headers. --- include/basic.h | 10 +++++++++- include/config.h | 9 ++++++++- include/discont.h | 9 ++++++++- include/faltung.h | 10 +++++++++- include/faltung_hilfe.h | 8 +++++++- include/fepc_cp.h | 8 +++++++- include/fepc_easy.h | 10 +++++++++- include/fepc_easy_helper.h | 9 ++++++++- include/fepc_easy_sparse.h | 21 ++++++++----------- include/fft_faltung.h | 10 +++++++++- include/folge.h | 10 ++++++++-- include/folgen_vektor.h | 10 +++++++++- include/funktion.h | 9 ++++++++- include/gauss.h | 41 ++++++++++++++++++++++++++++++++++++++ include/interval.h | 26 ++++++++++++++++++++++++ include/koeffizienten.h | 8 +++++++- include/kuerzen.h | 6 +++++- include/seconds.h | 10 +++++++++- 18 files changed, 195 insertions(+), 29 deletions(-) create mode 100644 include/gauss.h diff --git a/include/basic.h b/include/basic.h index 00c8ca8..a799a35 100755 --- a/include/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 */ @@ -162,6 +166,10 @@ print_vec_real(vec_real_p vector); fepc_real_t vec_real_norm(vec_real_p vector); +#ifdef __cplusplus +} +#endif + #endif diff --git a/include/config.h b/include/config.h index a3408cd..04a2bb2 100755 --- a/include/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,6 +34,10 @@ #include #include "seconds.h" +#ifdef __cplusplus +extern "C" { +#endif + /* * basic types */ @@ -76,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 index 374b56c..3bc4578 100755 --- a/include/discont.h +++ b/include/discont.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * 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 @@ -25,6 +25,10 @@ #include "fepc_easy.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { fepc_real_t y_0; @@ -97,5 +101,8 @@ discont_function_setup(discont_function_p function, int step, fepc_real_t start, 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/include/faltung.h b/include/faltung.h index 852c262..656cdf7 100755 --- a/include/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. @@ -40,4 +44,8 @@ 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/include/faltung_hilfe.h b/include/faltung_hilfe.h index 24663c5..7c3bfd1 100755 --- a/include/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,6 +60,9 @@ 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 index 1aca062..1fb191e 100644 --- a/include/fepc_cp.h +++ b/include/fepc_cp.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * 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 @@ -22,6 +22,9 @@ #include "fepc_easy.h" +#ifdef __cplusplus +extern "C" { +#endif typedef struct { int rank; @@ -66,5 +69,8 @@ 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 index ee1c113..5580728 100755 --- a/include/fepc_easy.h +++ b/include/fepc_easy.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * 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 @@ -25,6 +25,10 @@ #include "faltung.h" #include "interval.h" +#ifdef __cplusplus +extern "C" { +#endif + #ifndef INT_STEPS #define INT_STEPS 10 #endif @@ -156,5 +160,9 @@ 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 index 37c24d5..95f519a 100755 --- a/include/fepc_easy_helper.h +++ b/include/fepc_easy_helper.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * 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 @@ -21,6 +21,10 @@ #include "fepc_easy.h" +#ifdef __cplusplus +extern "C" { +#endif + func_p func_multi(func_p f, func_p g, fepc_real_t stepping); @@ -72,5 +76,8 @@ func_subtract2(func_p f, func_p g); 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 index 9a65c12..1c91dd3 100644 --- a/include/fepc_easy_sparse.h +++ b/include/fepc_easy_sparse.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2010 Stefan Handschuh (handschu@mis.mpg.de) + * 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 @@ -21,7 +21,9 @@ #include "fepc_easy.h" - +#ifdef __cplusplus +extern "C" { +#endif typedef struct { int summands; @@ -51,15 +53,8 @@ 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 \ No newline at end of file +#endif diff --git a/include/fft_faltung.h b/include/fft_faltung.h index 33b116a..6cca06f 100755 --- a/include/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/include/folge.h b/include/folge.h index b67f355..f1e62f2 100755 --- a/include/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 @@ -21,6 +21,10 @@ #include "fft_faltung.h" +#ifdef __cplusplus +extern "C" { +#endif + typedef struct { vec_p start; /* Startvektor der Folge */ vec_p lang; /* Laenge der Folge */ @@ -85,7 +89,9 @@ 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/include/folgen_vektor.h b/include/folgen_vektor.h index 5f0e70c..e947b4e 100755 --- a/include/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*/ @@ -108,6 +112,10 @@ folgen_vektor_factor_multi(folgen_vektor_p f, fepc_real_t a); folgen_matrix_p folgen_matrix_add(folgen_matrix_p f, folgen_matrix_p g); +#ifdef __cplusplus +} +#endif + #endif diff --git a/include/funktion.h b/include/funktion.h index ff3d745..ed75e09 100755 --- a/include/funktion.h +++ b/include/funktion.h @@ -1,6 +1,6 @@ /* * FEPC - * Copyright (C) 2009 Peter Gerds (gerds@mis.mpg.de), 2010 Stefan Handschuh (handschu@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 */ @@ -129,5 +133,8 @@ funcs_del(func_p * array, int length); void funcs_del_type(func_t * array, int length); +#ifdef __cplusplus +} +#endif #endif diff --git a/include/gauss.h b/include/gauss.h new file mode 100644 index 0000000..7aac736 --- /dev/null +++ b/include/gauss.h @@ -0,0 +1,41 @@ +/* + * Copyright (C) 2010 Philipp Waehnert (waehnert@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 + * 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 __gauss_h +#define __gauss_h + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct node_ node; +struct node_ { + double x; + double w; +}; + +node * create_gauss_nodes(int n); +int node_compare(const void* a, const void* b); +void sort_gauss_nodes(node* nodes, unsigned int n); +void shift_gauss_nodes(node* nodes, unsigned int n, double a, double b); +void free_gauss_nodes(node* nodes); + +#ifdef __cplusplus +} +#endif + +#endif // __gauss_h diff --git a/include/interval.h b/include/interval.h index e47a3ff..607635d 100755 --- a/include/interval.h +++ b/include/interval.h @@ -1,8 +1,30 @@ +/* + * 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; @@ -50,5 +72,9 @@ interval_clone(interval_p interval); interval_p * intervals_clone(interval_p *interval, int count); +#ifdef __cplusplus +} +#endif + #endif diff --git a/include/koeffizienten.h b/include/koeffizienten.h index 70428f6..24c77ef 100755 --- a/include/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,6 +83,9 @@ 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/include/kuerzen.h b/include/kuerzen.h index ef516c7..31bf957 100755 --- a/include/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/include/seconds.h b/include/seconds.h index 8e1e6af..6047a08 100755 --- a/include/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 From 5ebfbdd75c9802ebd8d74176f380d66cc8cf3086 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Mon, 21 Feb 2011 00:33:55 +0100 Subject: [PATCH 18/24] Removes invalid header file. --- include/gauss.h | 41 ----------------------------------------- 1 file changed, 41 deletions(-) delete mode 100644 include/gauss.h diff --git a/include/gauss.h b/include/gauss.h deleted file mode 100644 index 7aac736..0000000 --- a/include/gauss.h +++ /dev/null @@ -1,41 +0,0 @@ -/* - * Copyright (C) 2010 Philipp Waehnert (waehnert@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 - * 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 __gauss_h -#define __gauss_h - -#ifdef __cplusplus -extern "C" { -#endif - -typedef struct node_ node; -struct node_ { - double x; - double w; -}; - -node * create_gauss_nodes(int n); -int node_compare(const void* a, const void* b); -void sort_gauss_nodes(node* nodes, unsigned int n); -void shift_gauss_nodes(node* nodes, unsigned int n, double a, double b); -void free_gauss_nodes(node* nodes); - -#ifdef __cplusplus -} -#endif - -#endif // __gauss_h From e0c51b1ac70e1dd51d80b264a0ef2d70dfb2673f Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Fri, 4 Mar 2011 09:53:02 +0100 Subject: [PATCH 19/24] Corrects readme. --- README | 3 --- 1 file changed, 3 deletions(-) diff --git a/README b/README index 5dab4c0..ced6f12 100644 --- a/README +++ b/README @@ -16,6 +16,3 @@ GNU Lesser General Public License for more details. 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. From b27a7d76cb632c51ad433c367617ef613662b778 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Mon, 2 May 2011 10:47:34 +0200 Subject: [PATCH 20/24] Increases speed of folge multiplication and division. --- source/fepc_easy.c | 40 +++++++++++++++++++++++++++++++++++++++ source/fepc_easy_helper.c | 7 +++++-- 2 files changed, 45 insertions(+), 2 deletions(-) diff --git a/source/fepc_easy.c b/source/fepc_easy.c index a4f4586..5722bdd 100755 --- a/source/fepc_easy.c +++ b/source/fepc_easy.c @@ -682,6 +682,17 @@ void set_degree(func_p function, int degree) { +/** + * 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); @@ -696,6 +707,18 @@ setup_fepc_structure(func_t * func, Funcimpl function, interval_p* intervals, in } } +/** + * 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); @@ -757,6 +780,15 @@ vec_p getPredeccessor(vec_p v, int position, folge_p folge, int direction) { 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); @@ -807,6 +839,14 @@ func_p func_derive(func_p function, int direction, fepc_real_t stepping) { } +/** + * 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; diff --git a/source/fepc_easy_helper.c b/source/fepc_easy_helper.c index 95421a5..92cdafc 100755 --- a/source/fepc_easy_helper.c +++ b/source/fepc_easy_helper.c @@ -236,6 +236,7 @@ folge_multi(folge_p f, folge_p g, int step, fepc_real_t stepping) { 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 ); @@ -255,7 +256,7 @@ folge_multi(folge_p f, folge_p g, int step, fepc_real_t stepping) { y = folge_glied( r, g ); //print_vec(r); vec_del( r ); - back->glied[k] = x*y*folge_norm2(g, step, stepping); + back->glied[k] = x*y*norm; //printf("%i\t%f\t%f\n", step, x, y); } return back; @@ -282,6 +283,8 @@ folge_div(folge_p f, folge_p g, int step, fepc_real_t stepping) { 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 ); @@ -300,7 +303,7 @@ folge_div(folge_p f, folge_p g, int step, fepc_real_t stepping) { x = folge_glied( r, f ); y = folge_glied( r, g ); vec_del( r ); - back->glied[k] = x == 0. ? 0. : x/(y*folge_norm2(g, step, stepping)); + back->glied[k] = x == 0. ? 0. : x/(y*norm); } return back; } From 5152b58f760144a19ce6fa2cdf87a6b9d4281b64 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Thu, 5 May 2011 11:29:35 +0200 Subject: [PATCH 21/24] Implements subtraction if f is 0. --- source/folge.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/folge.c b/source/folge.c index f22265a..67543ca 100755 --- a/source/folge.c +++ b/source/folge.c @@ -419,7 +419,8 @@ folge_subtract(folge_p f, folge_p g) { size_f = vec_size( f->lang ); size_g = vec_size( g->lang ); if(size_f == 0) { - back = folge_copy( g ); // TODO multiply with -1 + back = folge_copy( g ); + folge_multi_factor(back, -1); return back; } From a142a8522595d926e86a075a09a4c1bbe6036c20 Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Wed, 22 Jun 2011 18:56:13 +0200 Subject: [PATCH 22/24] Improves fft speed by using r2c and c2r codes. --- source/CMakeLists.txt | 4 ++++ source/fft_faltung.c | 37 ++++++++++++++++++------------------- 2 files changed, 22 insertions(+), 19 deletions(-) diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 2583465..c8a3151 100755 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -8,6 +8,10 @@ LIST(REMOVE_ITEM ALL_HEADER_FILES ${CMAKE_SOURCE_DIR}/include/gauss.h) ADD_LIBRARY(fepc SHARED ${ALL_SOURCE_FILES} ${ALL_HEADER_FILES}) +add_executable (main main) +target_link_libraries (main fepc fftw3 m) + + if (EXISTS ${CMAKE_SOURCE_DIR}/include/gauss.h) ADD_LIBRARY(gauss SHARED gauss.c ${CMAKE_SOURCE_DIR}/include/gauss.h) endif() diff --git a/source/fft_faltung.c b/source/fft_faltung.c index c481075..46e2cc1 100755 --- a/source/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 @@ -42,7 +43,9 @@ fft_faltung(fepc_real_t* a, vec_p n_a, fepc_real_t* b, vec_p n_b) { 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; @@ -66,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;k Date: Tue, 5 Jun 2012 17:56:11 +0200 Subject: [PATCH 23/24] Corrects file modes. Adds .gitignore. Adds cmake option to use fftw3 w/ hint to resulting GPL-binary. --- .gitignore | 11 +++++++++++ CMakeLists.txt | 6 ++++++ README | 7 +++++++ include/basic.h | 0 include/config.h | 0 include/discont.h | 0 include/faltung.h | 0 include/faltung_hilfe.h | 0 include/fepc_easy.h | 0 include/fepc_easy_helper.h | 0 include/fft_faltung.h | 0 include/folge.h | 2 +- include/folgen_vektor.h | 0 include/funktion.h | 0 include/interval.h | 0 include/koeffizienten.h | 0 include/kuerzen.h | 0 include/seconds.h | 0 source/CMakeLists.txt | 38 +++++++++++++++++++++++------------- source/basic.c | 0 source/discont.c | 0 source/faltung.c | 0 source/faltung_hilfe.c | 0 source/fepc_easy.c | 0 source/fepc_easy_helper.c | 0 source/fft_faltung.c | 0 source/folge.c | 8 ++++++++ source/folgen_vektor.c | 0 source/funktion.c | 0 source/interval.c | 0 source/koeffizienten.c | 0 source/kuerzen.c | 0 source/seconds.c | 0 test/CMakeLists.txt | 40 +++++++++++++++----------------------- test/discont_test.c | 0 test/fepc_cp_test.c | 2 ++ test/fepc_easy_test.c | 0 test/stest.c | 0 38 files changed, 76 insertions(+), 38 deletions(-) create mode 100644 .gitignore mode change 100755 => 100644 CMakeLists.txt mode change 100755 => 100644 include/basic.h mode change 100755 => 100644 include/config.h mode change 100755 => 100644 include/discont.h mode change 100755 => 100644 include/faltung.h mode change 100755 => 100644 include/faltung_hilfe.h mode change 100755 => 100644 include/fepc_easy.h mode change 100755 => 100644 include/fepc_easy_helper.h mode change 100755 => 100644 include/fft_faltung.h mode change 100755 => 100644 include/folge.h mode change 100755 => 100644 include/folgen_vektor.h mode change 100755 => 100644 include/funktion.h mode change 100755 => 100644 include/interval.h mode change 100755 => 100644 include/koeffizienten.h mode change 100755 => 100644 include/kuerzen.h mode change 100755 => 100644 include/seconds.h mode change 100755 => 100644 source/CMakeLists.txt mode change 100755 => 100644 source/basic.c mode change 100755 => 100644 source/discont.c mode change 100755 => 100644 source/faltung.c mode change 100755 => 100644 source/faltung_hilfe.c mode change 100755 => 100644 source/fepc_easy.c mode change 100755 => 100644 source/fepc_easy_helper.c mode change 100755 => 100644 source/fft_faltung.c mode change 100755 => 100644 source/folge.c mode change 100755 => 100644 source/folgen_vektor.c mode change 100755 => 100644 source/funktion.c mode change 100755 => 100644 source/interval.c mode change 100755 => 100644 source/koeffizienten.c mode change 100755 => 100644 source/kuerzen.c mode change 100755 => 100644 source/seconds.c mode change 100755 => 100644 test/CMakeLists.txt mode change 100755 => 100644 test/discont_test.c mode change 100755 => 100644 test/fepc_easy_test.c mode change 100755 => 100644 test/stest.c 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 old mode 100755 new mode 100644 index 1766197..fba9b7b --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,12 @@ 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) diff --git a/README b/README index ced6f12..165ece4 100644 --- a/README +++ b/README @@ -16,3 +16,10 @@ GNU Lesser General Public License for more details. To read the text of the GNU Lesser General Public License, see + + +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/include/basic.h b/include/basic.h old mode 100755 new mode 100644 diff --git a/include/config.h b/include/config.h old mode 100755 new mode 100644 diff --git a/include/discont.h b/include/discont.h old mode 100755 new mode 100644 diff --git a/include/faltung.h b/include/faltung.h old mode 100755 new mode 100644 diff --git a/include/faltung_hilfe.h b/include/faltung_hilfe.h old mode 100755 new mode 100644 diff --git a/include/fepc_easy.h b/include/fepc_easy.h old mode 100755 new mode 100644 diff --git a/include/fepc_easy_helper.h b/include/fepc_easy_helper.h old mode 100755 new mode 100644 diff --git a/include/fft_faltung.h b/include/fft_faltung.h old mode 100755 new mode 100644 diff --git a/include/folge.h b/include/folge.h old mode 100755 new mode 100644 index f1e62f2..e2ce16c --- a/include/folge.h +++ b/include/folge.h @@ -19,7 +19,7 @@ #ifndef __FOLGE_H #define __FOLGE_H -#include "fft_faltung.h" +#include "basic.h" #ifdef __cplusplus extern "C" { diff --git a/include/folgen_vektor.h b/include/folgen_vektor.h old mode 100755 new mode 100644 diff --git a/include/funktion.h b/include/funktion.h old mode 100755 new mode 100644 diff --git a/include/interval.h b/include/interval.h old mode 100755 new mode 100644 diff --git a/include/koeffizienten.h b/include/koeffizienten.h old mode 100755 new mode 100644 diff --git a/include/kuerzen.h b/include/kuerzen.h old mode 100755 new mode 100644 diff --git a/include/seconds.h b/include/seconds.h old mode 100755 new mode 100644 diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt old mode 100755 new mode 100644 index c8a3151..6ac8583 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -1,18 +1,30 @@ -file(GLOB ALL_SOURCE_FILES *.c) +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) -# Find all header files in the include directory -file(GLOB ALL_HEADER_FILES ${CMAKE_SOURCE_DIR}/include/*.h) +if (HAS_FFTW3) + set(ALL_SOURCE_FILES fft_faltung.c ${ALL_SOURCE_FILES}) +endif(HAS_FFTW3) -LIST(REMOVE_ITEM ALL_SOURCE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/gauss.c) -LIST(REMOVE_ITEM ALL_HEADER_FILES ${CMAKE_SOURCE_DIR}/include/gauss.h) - -ADD_LIBRARY(fepc SHARED ${ALL_SOURCE_FILES} ${ALL_HEADER_FILES}) +ADD_LIBRARY(fepc SHARED ${ALL_SOURCE_FILES}) -add_executable (main main) -target_link_libraries (main fepc fftw3 m) +add_executable (main main.c) +target_link_libraries(main fepc m) - -if (EXISTS ${CMAKE_SOURCE_DIR}/include/gauss.h) - ADD_LIBRARY(gauss SHARED gauss.c ${CMAKE_SOURCE_DIR}/include/gauss.h) -endif() +if (HAS_FFTW3) + target_link_libraries(main fftw3) +endif(HAS_FFTW3) diff --git a/source/basic.c b/source/basic.c old mode 100755 new mode 100644 diff --git a/source/discont.c b/source/discont.c old mode 100755 new mode 100644 diff --git a/source/faltung.c b/source/faltung.c old mode 100755 new mode 100644 diff --git a/source/faltung_hilfe.c b/source/faltung_hilfe.c old mode 100755 new mode 100644 diff --git a/source/fepc_easy.c b/source/fepc_easy.c old mode 100755 new mode 100644 diff --git a/source/fepc_easy_helper.c b/source/fepc_easy_helper.c old mode 100755 new mode 100644 diff --git a/source/fft_faltung.c b/source/fft_faltung.c old mode 100755 new mode 100644 diff --git a/source/folge.c b/source/folge.c old mode 100755 new mode 100644 index 67543ca..d3d8fb7 --- a/source/folge.c +++ b/source/folge.c @@ -18,6 +18,10 @@ #include "folge.h" +#ifdef HAS_FFTW3 +#include "fft_faltung.h" +#endif + /******************************************************* * @@ -160,6 +164,7 @@ folge_slow_faltung(folge_p f,folge_p g) { folge_p folge_faltung(folge_p f,folge_p g) { +#ifdef HAS_FFTW3 folge_p back; int k; int size_f, size_g, dim; @@ -219,6 +224,9 @@ folge_faltung(folge_p f,folge_p g) { back->glied = w_glied; return back; +#else + return folge_slow_faltung(f, g); +#endif } diff --git a/source/folgen_vektor.c b/source/folgen_vektor.c old mode 100755 new mode 100644 diff --git a/source/funktion.c b/source/funktion.c old mode 100755 new mode 100644 diff --git a/source/interval.c b/source/interval.c old mode 100755 new mode 100644 diff --git a/source/koeffizienten.c b/source/koeffizienten.c old mode 100755 new mode 100644 diff --git a/source/kuerzen.c b/source/kuerzen.c old mode 100755 new mode 100644 diff --git a/source/seconds.c b/source/seconds.c old mode 100755 new mode 100644 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt old mode 100755 new mode 100644 index 07df6ad..b935582 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,37 +1,29 @@ -include_directories (${${PROJECT_NAME}_SOURCE_DIR}/src) +add_executable(fepc_easy_test fepc_easy_test.c) +target_link_libraries(fepc_easy_test fepc m) -link_directories (${${PROJECT_NAME}_BINARY_DIR}/src) +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 (fepc_easy_test fepc_easy_test.c) -target_link_libraries (fepc_easy_test fepc fftw3 m) - -add_executable (stest stest.c) -target_link_libraries (stest fepc fftw3 m) - -add_executable (discont_test discont_test.c) -target_link_libraries (discont_test fepc fftw3 m) - - -add_executable (sparsetest sparsetest.c) -target_link_libraries (sparsetest fepc fftw3 m) - - -if (EXISTS ${CMAKE_SOURCE_DIR}/include/gauss.h) - add_executable (gauss_test gauss_test.c) - target_link_libraries (gauss_test gauss blas lapack gsl m) - add_test(gauss_test gauss_test) -endif() - +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 fftw3 m) +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 old mode 100755 new mode 100644 diff --git a/test/fepc_cp_test.c b/test/fepc_cp_test.c index c16dbe9..543e3cf 100644 --- a/test/fepc_cp_test.c +++ b/test/fepc_cp_test.c @@ -77,6 +77,8 @@ main() { 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 old mode 100755 new mode 100644 diff --git a/test/stest.c b/test/stest.c old mode 100755 new mode 100644 From 1aa9531b93c53073aff1f68b4190401a2be084ff Mon Sep 17 00:00:00 2001 From: Stefan Handschuh Date: Thu, 28 Jun 2012 13:50:21 +0200 Subject: [PATCH 24/24] Adds missing copyright. --- source/fepc_easy_helper.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/fepc_easy_helper.c b/source/fepc_easy_helper.c index 92cdafc..b7b0817 100644 --- a/source/fepc_easy_helper.c +++ b/source/fepc_easy_helper.c @@ -1,6 +1,7 @@ /* * FEPC - * Copyright (C) 2010 Stefan Handschuh (handschu@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