libftsh.h

Go to the documentation of this file.

/*
 * Copyright (c) 2000, 2012 Martin J. Mohlenkamp
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 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 General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


/* make sure we haven't already included this library **/
#ifndef LIBFTSH_PRESENT
#define LIBFTSH_PRESENT


#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include <time.h>


#ifndef PI
#define PI 3.14159265358979323846264338328
#endif
 
#ifndef USE_double
#define USE_double 1
#endif

#if USE_double
#define REAL double
#else
#define REAL float
#endif

/* Deal with possible ways fftw was compiled */
#if 0

#if USE_double
  /* double precision header */
#include "dfftw.h"
#include "drfftw.h"
#else
  /* single precision header */
#include "sfftw.h"
#include "srfftw.h"
#endif
#else /*fftw in only one precision. pray it is the right precision */
/*  Make sure it is compiled with fftw_real == REAL */
#include "fftw.h"
#include "rfftw.h"
#endif

/*************************************************************************/
/**** data structures ****************************************************/
/*************************************************************************/


typedef struct{
  int numpts;
  fftw_complex *complex_work;
  fftw_plan w_plan;
  rfftw_plan rw_plan_forward;
  rfftw_plan rw_plan_reverse;
  REAL *trig_save; 
}Ftrigtw_Plan; 

typedef struct 
{
  int max_size;
  int n_row;
  int n_col;
  REAL *matrix;
} Dense_Matrix;

typedef struct 
{
  int max_size;
  int n_row;
  REAL *vector;
} Dense_Vector;

typedef struct 
{
  int max_size;
  int max_rows;
  int max_bands;
  int n_row;
  int n_col;
  int non_zero;
  int *band_p;
  int *band_list;
  REAL *matrix;
} Banded_Matrix;

typedef struct{
  unsigned int max_size; 
  int *int_work;
  REAL *real_work;
} Bdmat_Op_Workspace;

typedef struct{
  int best_level; 
  int num_levels;
  unsigned short int **branching;
  unsigned int **cost_full;
  unsigned int **cost_asym;
  unsigned int *order_level;
  unsigned int *order_within;
}Dyadic_Gsearch_Save;

typedef struct{
  int bell_length;
  int *bell_offsets;
  REAL *bell_handle;
}Fast1_Save;

typedef struct{
  int numpts;
  int left_extra;
  int right_extra;
  int fcn_offset;
  int coef_offset;
  REAL *bell;
  REAL *dualbell;
  int which_transform;
  Ftrigtw_Plan transform_plan;
}Local_Trig_Interval;

typedef struct{
  int max_intervals;
  int n_interval[2];
  Banded_Matrix matrix[2];
  Local_Trig_Interval *ilist[2];
}Pm_1D;

typedef struct{
  Banded_Matrix matrix[2];
}Pm_Direct_bd;

typedef struct{
  Dense_Matrix matrix[2];
}Pm_Direct_d;

typedef struct{
  int numpts;
  int band_limit;
  Dense_Vector dvect_coef[2];
  Dense_Vector dvect_fcn[2];
  REAL *eval_fcn[2];
  REAL *over_rootsin;
  REAL *mod_weight;
}Pm_1D_Workspace;


typedef struct{
  int numpts;
  REAL **sofar;
  int *size;
  REAL *ex;
  REAL *why;
}Pmncos_Recurrence_Save;

/*************************************************************************/
/******* END Data Structures *********************************************/
/*************************************************************************/


/*************************************************************************/
/******* BEGIN function prototypes ***************************************/
/*************************************************************************/

/* dvect.c */
Dense_Vector init_dvect_memory(unsigned  int n
                               );
void give_dvect_memory(Dense_Vector *current,
                       Dense_Vector *previous
                       );
/* dmat.c */
Dense_Matrix init_dmat_memory(unsigned  int n
                               );
void give_dmat_memory(Dense_Matrix *current,
                       Dense_Matrix *previous
                       );
void free_dmat_memory(Dense_Matrix *in
                      );
void dump_bin_dmat(Dense_Matrix *x,
                   FILE *out
                   );
void load_bin_dmat(Dense_Matrix *x,
                   FILE *in
                   );
/* dmat_ops.c */
void dmat_dvect_multiply(Dense_Vector *out,
                         Dense_Matrix *inmat,
                         Dense_Vector *invect
                         );
void dmat_t_dvect_multiply(Dense_Vector *out,
                           Dense_Matrix *inmat,
                           Dense_Vector *invect
                           );
/* bdmat.c */
Banded_Matrix init_bdmat_memory(int max_size,
                                int max_rows,
                                int max_bands
                                );
void give_bdmat_memory(Banded_Matrix *current,
                       Banded_Matrix *previous
                       );
void free_bdmat_memory(Banded_Matrix *in
                       );
void bdmat_copy(Banded_Matrix *out,
                Banded_Matrix *in
                );
void dump_bin_bdmat(Banded_Matrix *x,
                    FILE *out
                    );
void load_bin_bdmat(Banded_Matrix *x,
                    FILE *in
                    );
/* bdmat_ops.c */
void bdmat_dvect_multiply(Dense_Vector *out,
                          Banded_Matrix *inmat,
                          Dense_Vector *invect
                          );
void bdmat_t_dvect_multiply(Dense_Vector *out,
                            Banded_Matrix *inmat,
                            Dense_Vector *invect
                            );
void dmat_to_bdmat(Banded_Matrix *out,
                   Dense_Matrix *in,
                   REAL cutoff,
                   int min_gap
                   );
/* dyadic_gsearch.c */
int dyadic_gsearch(Dyadic_Gsearch_Save *dy_gs_save
                   );
void dyadic_gsearch_init(Dyadic_Gsearch_Save *dy_gs_save,
                         int num_levels
                         );
void dyadic_gsearch_free(Dyadic_Gsearch_Save *dy_gs_save
                         );
/* fast1_fromfile.c*/
void fast1_fromfile(int *numpts_ptr,
                    int *band_limit_ptr,
                    REAL **node_ptr,
                    REAL **weight_ptr,
                    Pm_1D **pmn_matrix_ptr,
                    Pm_1D_Workspace *pmnwork,
                    FILE *fp
                    );
void load_fast1(int *numpts_ptr,
                int *band_limit_ptr,
                int *num_levels_ptr,
                REAL **node_ptr,
                REAL **weight_ptr,
                Fast1_Save *saves_fordump,
                Pm_1D **pmn_matrix_ptr,
                FILE *fp
                );
/* direct1_fromfile.c*/
void direct1_fromfile(int *numpts_ptr,
                      int *band_limit_ptr,
                      REAL **node_ptr,
                      REAL **weight_ptr,
                      Pm_Direct_d **pmn_matrix_ptr,
                      Pm_1D_Workspace *pmnwork,
                      FILE *fp
                      );
void load_direct1(int *numpts_ptr,
                  int *band_limit_ptr,
                  REAL **node_ptr,
                  REAL **weight_ptr,
                  Pm_Direct_d **pmn_matrix_ptr,
                  FILE *fp
                  );
/* dctw.c */
void dctw(REAL *fcn,
          Ftrigtw_Plan *main_plan,
          int dir
          );
void dctw_init(Ftrigtw_Plan *main_plan,
               int flags,
               int numpts
               );
void rfftw_to_NR_format(REAL *data_nr,
                        REAL *data_rfftw,
                        int numpts,
                        int dir
                        );
/* lctw.c */
void lctw(REAL *fcn,
          Ftrigtw_Plan *main_plan
          );
void lctw_init(Ftrigtw_Plan *main_plan,
               int flags,
               int numpts
               );
/* bell_utils.c */
REAL dualabell(REAL (*origbell) (int,int),
                 int laplength,
                 int t
                 );
void bell_1size_init(REAL *bellout,
                     REAL (*whichbell) (int,int),
                     int lapleft,
                     int center,
                     int lapright,
                     int dual_flag
                     );
void fold_e_e_1s(REAL *coef,
                 REAL *fcn,
                 int numpts,
                 REAL *left,
                 int n_left,
                 REAL *right,
                 int n_right,
                 REAL *bell,
                 int dir
                 );
void fold_e_o_1s(REAL *coef,
                 REAL *fcn,
                 int numpts,
                 REAL *left,
                 int n_left,
                 REAL *right,
                 int n_right,
                 REAL *bell,
                 int dir
                 );
/* Pmn_construct.c */
REAL Pmncos(int order,
            int index,
            REAL phi
            );
void Pmncos_resave_init(Pmncos_Recurrence_Save *pmn_save,
                        REAL *node,
                        int numpts
                        );
void next_rtsinPmncos(REAL *outfcn,
                      Pmncos_Recurrence_Save *pmn_save,
                      int halfpts,
                      int order,
                      int index
                      );
void Pmncos_resave_destroy(Pmncos_Recurrence_Save *pmn_save
                           );
/* direct1_tofile.c */
void direct1_tofile(int  numpts,
                    int band_limit,
                    REAL *node,
                    REAL *weight,
                    FILE *fp
                    );
void init_direct1(Pm_Direct_d **out,
                  int halfpts,
                  int band_limit,
                  REAL *node
                  );
void dump_direct1(int numpts,
                  int band_limit,
                  REAL *node,
                  REAL *weight,
                  Pm_Direct_d *pmn_matrix,
                  FILE *fp
                  );
/* fast1_tofile.c */
REAL fast1_tofile(int  numpts,
                  int band_limit,
                  REAL *node,
                  REAL *weight,
                  REAL epsilon,
                  FILE *fp
                  );
void init_fast1(Pm_1D **out,
                Fast1_Save *saves_fordump,
                int halfpts,
                int band_limit,
                REAL *node,
                int *numlevels,
                REAL epsilon,
                REAL (*bell_choice) (int,int)
                );
void dump_fast1(int numpts,
                int band_limit,
                int num_levels,
                REAL *node,
                REAL *weight,
                Fast1_Save *saves_fordump,
                Pm_1D *pmn_matrix,
                FILE *fp
                );
/* direct1_transform.c */
void Pm_direct1_analyse(REAL *outcoef,
                        Pm_Direct_d *inmat,
                        REAL *infcn,
                        Pm_1D_Workspace *pmnwork
                        );
void Pm_direct1_synthesize(REAL *outfcn,
                           Pm_Direct_d *inmat,
                           REAL *incoef,
                           Pm_1D_Workspace *pmnwork
                           );
/* matbells.c  */
REAL bellmat1(int laplength,
                int t
                );
REAL bellmat2(int laplength,
                int t
                );
REAL bellmat3(int laplength,
                int t
                );
REAL bellmat4(int laplength,
                int t
                );
REAL bellmat5(int laplength,
                int t
                );
REAL bellmat6(int laplength,
                int t
                );
REAL bellmat7(int laplength,
                int t
                );
REAL bellmat8(int laplength,
                int t
                );
REAL bellmat9(int laplength,
                int t
                );
REAL bellmat10(int laplength,
                int t
                );
REAL bellmat11(int laplength,
                int t
                );
REAL bellmat12(int laplength,
                int t
                );
REAL bellmat13(int laplength,
                int t
                );
REAL bellmat14(int laplength,
                int t
                );
REAL bellmat15(int laplength,
                int t
                );
REAL bellmat16(int laplength,
                int t
                );
REAL bellmat17(int laplength,
                int t
                );
REAL bellmat18(int laplength,
                int t
                );
REAL bellmat19(int laplength,
                int t
                );
REAL bellmat20(int laplength,
                int t
                );
/* fast1_transform.c */
void Pm_1d_analyse(REAL *outcoef,
                   Pm_1D *inmat,
                   REAL *infcn,
                   Pm_1D_Workspace *pmnwork
                   );
void Pm_1d_synthesize(REAL *outfcn,
                      Pm_1D *inmat,
                      REAL *incoef,
                      Pm_1D_Workspace *pmnwork
                      );
/* mk_gnode_gweight.c */
void mk_gnode_gweight(REAL *gnode,
                      REAL *gweight,
                      REAL epsilon,
                      int numpts
                      );
void mk_cosi_gnode(REAL *gnode,
                   REAL epsilon,
                   int numpts
                   );
void cosi_gnode_to_gweight(REAL *gweight,
                           REAL *gnode,
                           int numpts
                           );
/* second.c */
double second();
/* getoption.c */
char* getoption(const char* flag,
                int argc,
                char **argv
                );
/* pointwise_ops.c */
void pointwise_multiply(REAL *out,
                        REAL *in_1,
                        REAL *in_2,
                        int numpts
                        );
void reflect_evenodd(REAL *even,
                     REAL *odd,
                     REAL *whole,
                     int numpts,
                     int dir
                     );
void remove_rootsin(REAL *out,
                    REAL *in,
                    REAL *node,
                    int numpts
                    );
/* dyadic_trig_gsearch.c */
void dyadic_g_ilist_init(Local_Trig_Interval *ilist[2],
                         int *ints_used,
                         int *coefs_used,
                         REAL **bell_handle,
                         int *bell_length,
                         int **bell_offsets,
                         int halfpts,
                         int num_levels,
                         REAL (*bell_choice) (int,int),
                         int flags
                         );
void ltrigint_init(Local_Trig_Interval *out,
                   int numpts,
                   int left_extra,
                   int right_extra,
                   int fcn_offset,
                   int coef_offset,
                   REAL (*bell_choice)(int,int),
                   int which_transform,
                   int fftw_flags
                   );
void gsearch_to_ilist(Local_Trig_Interval *ilist_out,
                      int *intervals_kept,
                      int max_intervals,
                      Local_Trig_Interval *ilist_in,
                      Dyadic_Gsearch_Save *dy_gs_save,
                      int halfpts
                      );
/* trig_int_ops.c */
void evaluate_intlist(REAL *outfcn,
                      REAL *incoef,
                      Local_Trig_Interval *interval_list,
                      int num_intervals
                      );
void expand_d_intlist(REAL *outcoef,
                      REAL *infcn,
                      Local_Trig_Interval *interval_list,
                      int num_intervals
                      );
void expand_intlist(REAL *outcoef,
                   REAL *infcn,
                   Local_Trig_Interval *interval_list,
                   int num_intervals
                   );
/* count_thresh.c */
int count_thresh(REAL *lambda,
                 int howmany,
                 REAL thresh
                 );
/* init_Pm_1d_workspace */
void init_Pm_1d_workspace(Pm_1D_Workspace *out,
                          int numpts,
                          int band_limit,
                          REAL *node,
                          REAL *weight
                          );

/* */
void ftgaqd(int *nlat,
            REAL *theta,
            REAL *wts,
            REAL *work,
            int *lwork,
            int *ierror
            );
void gnode_to_Pmninit(REAL ***sofar,
                      REAL **ex,
                      REAL **why,
                      int **size,
                      REAL *gnode,
                      int hemipts
                      );

/*************************************************************************/
/******* END function prototypes *****************************************/
/*************************************************************************/

#endif /* LIBFTSH_PRESENT */