SU2Wrappers.h

Go to the documentation of this file.

#ifndef SU2WRAPPERS_H_
#define SU2WRAPPERS_H_
 
// As the default libraries for the 3nj-symbols, we use the GSL (Gnu scientific library):
#if !defined USE_WIG_SU2_COEFFS && !defined USE_WIG_SU2_COEFFS && !defined USE_FAST_WIG_SU2_COEFFS
#define USE_GSL_SU2_COEFFS 1
#endif
 
#ifdef USE_GSL_SU2_COEFFS
#include <gsl/gsl_sf_coupling.h>
#pragma message("Using GSL library for 3nj-symbols.")
#endif
 
#ifdef USE_WIG_SU2_COEFFS
#include "wigxjpf.h"
#pragma message("Using WIGXJPF library for 3nj-symbols.")
#endif
 
#ifdef USE_FAST_WIG_SU2_COEFFS
#include "fastwigxj.h"
#include "wigxjpf.h"
#pragma message("Using FASTWIGXJ library for 3nj-symbols.")
#endif
 
#include "DmrgExternal.h"
 
#ifdef USE_GSL_SU2_COEFFS
inline double coupl_9j_base(const int q1, const int q2, const int q3, 
                            const int q4, const int q5, const int q6, 
                            const int q7, const int q8, const int q9)
{
    return gsl_sf_coupling_9j(q1-1, q2-1, q3-1,
                              q4-1, q5-1, q6-1,
                              q7-1, q8-1, q9-1);
}
 
inline double coupl_6j_base (const int q1, const int q2, const int q3, 
                             const int q4, const int q5, const int q6)
{
    return gsl_sf_coupling_6j(q1-1, q2-1, q3-1,
                              q4-1, q5-1, q6-1);
}
 
inline double coupl_3j_base (const int q1  , const int q2  , const int q3,
                             const int q1_z, const int q2_z, const int q3_z)
{
    return gsl_sf_coupling_3j(q1-1, q2-1, q3-1,
                              q1_z, q2_z, q3_z);
}
#endif //USE_GSL_SU2_COEFFS
 
#ifdef USE_WIG_SU2_COEFFS
inline double coupl_9j_base(const int q1, const int q2, const int q3, 
                            const int q4, const int q5, const int q6, 
                            const int q7, const int q8, const int q9)
{
    return wig9jj(q1-1, q2-1, q3-1,
                  q4-1, q5-1, q6-1,
                  q7-1, q8-1, q9-1);
}
 
inline double coupl_6j_base (const int q1, const int q2, const int q3, 
                             const int q4, const int q5, const int q6)
{
    return wig6jj(q1-1, q2-1, q3-1,
                  q4-1, q5-1, q6-1);
}
 
inline double coupl_3j_base (const int q1  , const int q2  , const int q3,
                             const int q1_z, const int q2_z, const int q3_z)
{
    return wig3jj(q1-1, q2-1, q3-1,
                  q1_z, q2_z, q3_z);
}
#endif //USE_WIG_SU2_COEFFS
 
#ifdef USE_FAST_WIG_SU2_COEFFS
inline double coupl_9j_base(const int q1, const int q2, const int q3, 
                            const int q4, const int q5, const int q6, 
                            const int q7, const int q8, const int q9)
{
    return fw9jja(q1-1, q2-1, q3-1,
                  q4-1, q5-1, q6-1,
                  q7-1, q8-1, q9-1);
}
 
inline double coupl_6j_base (const int q1, const int q2, const int q3, 
                             const int q4, const int q5, const int q6)
{
    return fw6jja(q1-1, q2-1, q3-1,
                  q4-1, q5-1, q6-1);
}
 
inline double coupl_3j_base (const int q1  , const int q2  , const int q3,
                             const int q1_z, const int q2_z, const int q3_z)
{
    return fw3jja(q1-1, q2-1, q3-1,
                  q1_z, q2_z);
}
#endif //USE_FAST_WIG_SU2_COEFFS
 
#ifdef OWN_HASH_CGC
#ifdef _OPENMP
assert(omp_get_max_threads() == 1 and "Hashing the cgcs is not treadsafe!");
#endif
 
std::unordered_map<std::array<int,9>,double > Table9j;
std::unordered_map<std::array<int,6>,double > Table6j;
 
double coupling_9j (const int q1, const int q2, const int q3, 
                    const int q4, const int q5, const int q6, 
                    const int q7, const int q8, const int q9)
{
    auto it = Table9j.find(std::array<int,9>{q1,q2,q3,q4,q5,q6,q7,q8,q9});
    
    if (it != Table9j.end())
    {
        return Table9j[std::array<int,9>{q1,q2,q3,q4,q5,q6,q7,q8,q9}];
    }
    else
    {
        double out = coupl_9j_base(q1,q2,q3,
                                   q4,q5,q6,
                                   q7,q8,q9);
        Table9j[std::array<int,9>{q1,q2,q3,q4,q5,q6,q7,q8,q9}] = out;
        return out;
    }
}
 
double coupling_6j (const int q1, const int q2, const int q3, 
                    const int q4, const int q5, const int q6)
{
    auto it = Table6j.find(std::array<int,6>{q1,q2,q3,q4,q5,q6});   
    if (it != Table6j.end())
    {
        return Table6j[std::array<int,6>{q1,q2,q3,q4,q5,q6}];
    }
    else
    {
        double out = coupl_6j_base(q1,q2,q3,
                                   q4,q5,q6);
        Table6j[std::array<int,6>{q1,q2,q3,q4,q5,q6}] = out;
        return out;
    }
}
 
inline double coupling_3j (const int q1  , const int q2  , const int q3,
                           const int q1_z, const int q2_z, const int q3_z)
{
    return coupl_3j_base(q1  ,q2  ,q3,
                         q1_z,q2_z,q3_z);
}
 
#else //no OWN_HASH_CGC
inline double coupling_9j (const int q1, const int q2, const int q3, 
                           const int q4, const int q5, const int q6, 
                           const int q7, const int q8, const int q9)
{
    return coupl_9j_base(q1,q2,q3,
                         q4,q5,q6,
                         q7,q8,q9);
}
 
inline double coupling_6j (const int q1, const int q2, const int q3, 
                           const int q4, const int q5, const int q6)
{
    return coupl_6j_base(q1,q2,q3,
                         q4,q5,q6);
}
 
inline double coupling_3j (const int q1  , const int q2  , const int q3,
                           const int q1_z, const int q2_z, const int q3_z)
{
    return coupl_3j_base(q1  ,q2  ,q3,
                         q1_z,q2_z,q3_z);
}
 
#endif //OWN_HASH_CGC
 
#endif