| 1 | #ifndef mymath_h
|
|---|
| 2 | #define mymath_h
|
|---|
| 3 | /** \file mymath.h
|
|---|
| 4 | * Header file for \ref mymath.c
|
|---|
| 5 | *
|
|---|
| 6 | * Contains declarations of the functions implemented in \ref mymath.c, shorter
|
|---|
| 7 | * definitions of mathematical constants such as PI, square root of 2 SQRT2
|
|---|
| 8 | * and hard-coded constructs for calculating maximum MAX(), row majors CalcRowMajor2D(),
|
|---|
| 9 | * CalcRowMajor3D(), determinants RDET2(), RDET3(), scalar products real RSP3()
|
|---|
| 10 | * and complex CSP3re(), CSP3im(), euclidian norm real RNORMSQ3() and complex CNORMSQ3(),
|
|---|
| 11 | * complex multiplication CCMULTre(), CCMULTim(), complex multiplication with scalar
|
|---|
| 12 | * RCMULTre(), RCMULTim().
|
|---|
| 13 | *
|
|---|
| 14 | Project: ParallelCarParrinello
|
|---|
| 15 | Jan Hamaekers
|
|---|
| 16 | 2000
|
|---|
| 17 |
|
|---|
| 18 | File: mymath.h
|
|---|
| 19 | $Id: mymath.h,v 1.15 2007-03-29 13:35:51 foo Exp $
|
|---|
| 20 | */
|
|---|
| 21 |
|
|---|
| 22 | // use double precision fft when we have it
|
|---|
| 23 | #ifdef HAVE_CONFIG_H
|
|---|
| 24 | #include <config.h>
|
|---|
| 25 | #endif
|
|---|
| 26 |
|
|---|
| 27 | #if defined _BSD_SOURCE || defined _XOPEN_SOURCE
|
|---|
| 28 | //! short form for pi from math.h
|
|---|
| 29 | # define PI M_PI
|
|---|
| 30 | //! short form for square root of 2 from math.h
|
|---|
| 31 | # define SQRT2 M_SQRT2
|
|---|
| 32 | #else /* generische Form */
|
|---|
| 33 | //! short form for pi
|
|---|
| 34 | # define PI (acos(-1.0))
|
|---|
| 35 | //! short form for square root of 2
|
|---|
| 36 | # define SQRT2 (sqrt(2.))
|
|---|
| 37 | #endif
|
|---|
| 38 |
|
|---|
| 39 | //! Bohr radius in Angstrï¿œm
|
|---|
| 40 | //#define BOHRRADIUS 0.5291772108
|
|---|
| 41 |
|
|---|
| 42 | #include "defs.h"
|
|---|
| 43 |
|
|---|
| 44 | #ifdef HAVE_DFFTW_H
|
|---|
| 45 | #include "dfftw.h"
|
|---|
| 46 | #else
|
|---|
| 47 | #include "fftw.h"
|
|---|
| 48 | #endif
|
|---|
| 49 |
|
|---|
| 50 | #define MAX(a,b) ((a) > (b) ? (a) : (b)) //!< returns maximum of a or b
|
|---|
| 51 | #define CalcRowMajor3D(R0,R1,R2,N0,N1,N2) ((R2)+(N2)*((R1)+(N1)*(R0)))//!< calculates row major of 3x3 matrix
|
|---|
| 52 | #define CalcRowMajor2D(R0,R1,N0,N1) ((R1)+(N1)*(R0)) //!< calculates row major of 2x2 matrix
|
|---|
| 53 | #define RSP3(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2]) //!< scalar product of two 3-dim vectors
|
|---|
| 54 | #define RNORMSQ3(a) ((a)[0]*(a)[0] + (a)[1]*(a)[1] + (a)[2]*(a)[2]) //!< squared euclidian norm
|
|---|
| 55 | #define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3]) //!< hard-coded determinant of a 3x3 matrix
|
|---|
| 56 | #define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2)) //!< hard-coded determinant of a 2x2 matrix
|
|---|
| 57 | #define CCMULTre(a,b) ((a).re*(b).re - (a).im*(b).im) //!< real part of a complex multiplication
|
|---|
| 58 | #define CCMULTim(a,b) ((a).re*(b).im + (a).im*(b).re) //!< imaginary part of a complex multiplication
|
|---|
| 59 | #define RCMULTre(a,b) ((a).re*(b)) //!< real part of a complex number scaled by a real number
|
|---|
| 60 | #define RCMULTim(a,b) ((a).im*(b)) //!< imaginary part of a complex number scaled by a real number
|
|---|
| 61 | #define CSP3re(a,b) (CCMULTre((a)[0],(b)[0]) + CCMULTre((a)[1],(b)[1]) + CCMULTre((a)[2],(b)[2])) //!< real part of a scalar product of two 3x3 complex vectors
|
|---|
| 62 | #define CSP3im(a,b) (CCMULTim((a)[0],(b)[0]) + CCMULTim((a)[1],(b)[1]) + CCMULTim((a)[2],(b)[2])) //!< imaginary part of a scalar product of two 3x3 complex vectors
|
|---|
| 63 | #define CNORMSQ3(a) ((a).re[0]*(a).re[0] + (a).re[1]*(a).re[1] + (a).re[2]*(a).re[2] + (a).im[0]*(a).im[0] + (a).im[1]*(a).im[1] + (a).im[2]*(a).im[2]) //!< square of complex euclidian norm
|
|---|
| 64 |
|
|---|
| 65 | inline double tpow(double, int);
|
|---|
| 66 | inline int Rest(int n, int m);
|
|---|
| 67 | inline void RTranspose3(double *A);
|
|---|
| 68 | inline void RMatMat33(double C[NDIM*NDIM], const double A[NDIM*NDIM], const double B[NDIM*NDIM]);
|
|---|
| 69 | inline void RMat33Vec3(double C[NDIM], const double M[NDIM*NDIM], const double V[NDIM]);
|
|---|
| 70 | inline int RMatReci3(double B[NDIM_NDIM], const double A[NDIM_NDIM]);
|
|---|
| 71 | inline void RVec3Mat33(double C[NDIM], const double V[NDIM], const double M[NDIM*NDIM]);
|
|---|
| 72 | inline void VP3(double V[NDIM], double A[NDIM], double B[NDIM]);
|
|---|
| 73 | /* Skalarprodukt */
|
|---|
| 74 | inline double SP(const double *a, const double *b, const int n);
|
|---|
| 75 | /* Multiplikation mit Skalar */
|
|---|
| 76 | inline void SM(double *a, const double c, const int n);
|
|---|
| 77 | /* Nullvektor erzeugen */
|
|---|
| 78 | inline void NV(double *a, int n);
|
|---|
| 79 | inline double dSum(int n, double *dx, int incx);
|
|---|
| 80 | inline double Simps(int n, double *f, double h);
|
|---|
| 81 | inline double derf(double x);
|
|---|
| 82 | /* Initialisiere a array[3] mit b - c Orte mit periodisch */
|
|---|
| 83 | double Dist(const double *a, const double *b, const int n);
|
|---|
| 84 | inline void SetArrayToDouble0(double *a, int n);
|
|---|
| 85 | void PrintCMat330(fftw_complex M[NDIM_NDIM]);
|
|---|
| 86 | void PrintRMat330(fftw_real M[NDIM_NDIM]);
|
|---|
| 87 | void PrintCVec30(fftw_complex M[NDIM]);
|
|---|
| 88 | void PrintRVec30(fftw_real M[NDIM]);
|
|---|
| 89 | void RotateToAlign(fftw_real Q[NDIM_NDIM], fftw_real matrix[NDIM_NDIM], fftw_real vector[NDIM]);
|
|---|
| 90 | #endif
|
|---|
