Commit fe7a5980 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC/MR) 15/4/15 update laguerre_numbers.h to set automatically the flaoting...

(JEC/MR) 15/4/15 update laguerre_numbers.h to set automatically the flaoting size; use OpenBLAS in laguerreTransform.cc to compute the Analysis/Synthesis
parent 419f3dc5
......@@ -2,6 +2,12 @@
include Darwin_g++_omp_make.inc
# ===== OpenBLAS ========
BLASDIR = /opt/local/
BLASLIB = $(BLASDIR)/lib
BLASLIBN = -lopenblas
# ===== Libsharp =====
#SHARPDIR = <yourdir>/libsharp-code/auto/
......@@ -64,7 +70,7 @@ CXXHDR = lagsht_exceptions.h \
CPPFLAGS += -I$(SHARPINC)
LDFLAGS += -L$(SHARPLIB) $(SHARPLIBN) -lm
LDFLAGS += -L$(SHARPLIB) $(SHARPLIBN) -L$(BLASLIB) $(BLASLIBN) -lm
......
#ifndef LAGSHTNUMBERS_SEEN
#define LAGSHTNUMBERS_SEEN
typedef signed char int_1;
typedef unsigned char uint_1;
typedef short int_2;
typedef unsigned short uint_2;
typedef int int_4;
typedef unsigned int uint_4;
typedef long long int_8;
typedef unsigned long long uint_8;
typedef float r_4;
typedef double r_8;
typedef long double r_16;
// Template magic to select the proper data types. These templates
// should not be used outside this file.
template <typename T, bool equalSize> struct sizeChooserHelper__
{ typedef void TYPE; };
template <typename T> struct sizeChooserHelper__<T,true>
{ typedef T TYPE; };
template <typename T1, typename T2, typename T3> struct sizeChooserHelper2__
{ typedef T1 TYPE; };
template <typename T2, typename T3> struct sizeChooserHelper2__ <void, T2, T3>
{ typedef T2 TYPE; };
template <typename T3> struct sizeChooserHelper2__ <void, void, T3>
{ typedef T3 TYPE; };
template <> struct sizeChooserHelper2__ <void, void, void>
{ };
template <int sz, typename T1, typename T2=char, typename T3=char>
struct sizeChooser__
{
typedef typename sizeChooserHelper2__
<typename sizeChooserHelper__<T1,sizeof(T1)==sz>::TYPE,
typename sizeChooserHelper__<T2,sizeof(T2)==sz>::TYPE,
typename sizeChooserHelper__<T3,sizeof(T3)==sz>::TYPE >::TYPE TYPE;
};
typedef signed char int_1; // correct by definition
typedef unsigned char uint_1; // correct by definition
typedef sizeChooser__<2, short, int>::TYPE int_2;
typedef sizeChooser__<2, unsigned short, unsigned int>::TYPE uint_2;
typedef sizeChooser__<4, int, long, short>::TYPE int_4;
typedef sizeChooser__<4, unsigned int, unsigned long, unsigned short>::TYPE uint_4;
typedef sizeChooser__<8, long, long long>::TYPE int_8;
typedef sizeChooser__<8, unsigned long, unsigned long long>::TYPE uint_8;
typedef sizeChooser__<4, float, double>::TYPE r_4;
typedef sizeChooser__<8, double, long double>::TYPE r_8;
typedef sizeChooser__<16, double, long double>::TYPE r_16;
#endif //LAGSHTNUMBERS_SEEN
......@@ -2,6 +2,10 @@
#include <iostream>
#include <algorithm> // for copy
#include <iterator> // for ostream_iterator
#idef CBLASS
#include <cblas.h>
//#include <cblas_openblas.h> //Darwin
#endif
using namespace std;
#include "laguerreTransform.h"
......@@ -97,13 +101,16 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
lag.Values(rk,LnAll);
for (int n = 0; n<N_; n++ ){
// LnkMtx[n*N_+k] = LnAll[n]*facts[n];
LnkMtx[n+N_*k] = LnAll[n]*facts[n];
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k];
}
}
// tm.Split("lagTrans Test/End LnkMtx initialization");
#ifdef CBLAS
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fn[0]), 2*stride);
#else
vector<complex<r_8> > vtmp(N_);
for (int l=0; l<stride; l++) {
vtmp.assign(N_,0.);
......@@ -112,13 +119,14 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
r_8 wi = weights_[i];
for (int n=0; n<N_; n++){
// vtmp[n] += fli * wi * LnkMtx[n*N_+i];
vtmp[n] += fli * wi * LnkMtx[n+N_*i];
vtmp[n] += fli * LnkMtx[n+N_*i];
}//loop on k
}//loop on n
for (int n=0; n<N_; n++) {
fn[l+n*stride] += vtmp[n];
}
}//loop on l
#endif
}//MultiAnalysis
......@@ -143,12 +151,8 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
facts[0] = invalphaFact;
for(int n=1; n<N_; n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
// Timer tm("lagTrans synthesis Test");
//TODO : voir de meme si cette initialization sert
LaguerreFuncQuad lag(N_-1);
// vector<r_16> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
for (int k=0; k<N_; k++){
......@@ -159,8 +163,9 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
}
}
// tm.Split("lagTrans Test/End LnkMtx initialization");
#if CBLAS
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fn[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fi[0]), 2*stride);
#else
vector<complex<r_8> > vtmp(N_);
for (int l=0; l<stride; l++) {
vtmp.assign(N_,0.);
......@@ -172,7 +177,7 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
}//loop on n
for (int k=0; k<N_; k++) fi[l+k*stride] += vtmp[k];
}//loop on l
#endif
}//MultiSynthesis
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment