Operators.h

Go to the documentation of this file.

/*-------------------------------------------------------------------
Copyright 2011 Ravishankar Sundararaman

This file is part of JDFTx.

JDFTx 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 3 of the License, or
(at your option) any later version.

JDFTx 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 JDFTx.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------*/

#ifndef JDFTX_CORE_OPERATORS_H
#define JDFTX_CORE_OPERATORS_H


#include <core/ScalarField.h>
#include <core/BlasExtra.h>
#include <core/GridInfo.h>
#include <core/LoopMacros.h>

//----------------- Real / complex conversion ------------------
ScalarField Real(const complexScalarField&); 
ScalarFieldTilde Real(const complexScalarFieldTilde&); 
ScalarField Imag(const complexScalarField&); 
ScalarFieldTilde Imag(const complexScalarFieldTilde&); 
complexScalarField Complex(const ScalarField&); 
complexScalarField Complex(const ScalarField& re, const ScalarField& im); 
complexScalarFieldTilde Complex(const ScalarFieldTilde&); 

//------------------------------ Linear Unary operators ------------------------------

ScalarFieldTilde O(const ScalarFieldTilde&); 
ScalarFieldTilde O(ScalarFieldTilde&&); 
complexScalarFieldTilde O(const complexScalarFieldTilde&); 
complexScalarFieldTilde O(complexScalarFieldTilde&&); 

//Transform operators:
//  compat: force GPU transform output to be FFTW compatible (affects only how redundant Nyquist frequency components in the C->R transforms are handled)
//  nThreads: maximum number of threads to use; use to prevent thread nesting when performing transforms of vector fields etc.

ScalarField I(const ScalarFieldTilde&, bool compat=false, int nThreads=0); 
ScalarField I(ScalarFieldTilde&&, bool compat=false, int nThreads=0); 
complexScalarField I(const complexScalarFieldTilde&, int nThreads=0); 
complexScalarField I(complexScalarFieldTilde&&, int nThreads=0); 

ScalarFieldTilde J(const ScalarField&, int nThreads=0); 
complexScalarFieldTilde J(const complexScalarField&, int nThreads=0); 
complexScalarFieldTilde J(complexScalarField&&, int nThreads=0); 

ScalarFieldTilde Idag(const ScalarField&, int nThreads=0); 
complexScalarFieldTilde Idag(const complexScalarField&, int nThreads=0); 
complexScalarFieldTilde Idag(complexScalarField&&, int nThreads=0); 

ScalarField Jdag(const ScalarFieldTilde&, bool compat=false, int nThreads=0); 
ScalarField Jdag(ScalarFieldTilde&&, bool compat=false, int nThreads=0); 
complexScalarField Jdag(const complexScalarFieldTilde&, int nThreads=0); 
complexScalarField Jdag(complexScalarFieldTilde&&, int nThreads=0); 

ScalarField JdagOJ(const ScalarField&); 
ScalarField JdagOJ(ScalarField&&); 
complexScalarField JdagOJ(const complexScalarField&); 
complexScalarField JdagOJ(complexScalarField&&); 

ScalarFieldTilde L(const ScalarFieldTilde&); 
ScalarFieldTilde L(ScalarFieldTilde&&); 
complexScalarFieldTilde L(const complexScalarFieldTilde&); 
complexScalarFieldTilde L(complexScalarFieldTilde&&); 

ScalarFieldTilde Linv(const ScalarFieldTilde&); 
ScalarFieldTilde Linv(ScalarFieldTilde&&); 
complexScalarFieldTilde Linv(const complexScalarFieldTilde&); 
complexScalarFieldTilde Linv(complexScalarFieldTilde&&); 

void zeroNyquist(RealKernel& Gdata); 
void zeroNyquist(ScalarFieldTilde& Gptr); 
void zeroNyquist(ScalarField& Rptr); 

//------------------------------ Nonlinear Unary operators ------------------------------

ScalarField exp(const ScalarField&); 
ScalarField exp(ScalarField&&); 

ScalarField log(const ScalarField&); 
ScalarField log(ScalarField&&); 

ScalarField sqrt(const ScalarField&); 
ScalarField sqrt(ScalarField&&); 

ScalarField inv(const ScalarField&); 
ScalarField inv(ScalarField&&); 

ScalarField pow(const ScalarField&, double alpha); 
ScalarField pow(ScalarField&&, double alpha); 

#define Tptr std::shared_ptr<T> 

template<class T> Tptr clone(const Tptr& X) { if(X) return X->clone(); else return 0; } 

//------------------------------ Multiplication operators ------------------------------

template<class T> Tptr& operator*=(Tptr& in, double scaleFac) { if(in) in->scale *= scaleFac; return in; } 
template<class T> Tptr operator*(const Tptr& in, double scaleFac) { Tptr out(in->clone()); return out *= scaleFac; } 
template<class T> Tptr operator*(double scaleFac, const Tptr& in) { Tptr out(in->clone()); return out *= scaleFac; } 
template<class T> Tptr operator*(Tptr&& in, double scaleFac) { return in *= scaleFac; } 
template<class T> Tptr operator*(double scaleFac, Tptr&& in) { return in *= scaleFac; } 

template<class T> Tptr conj(Tptr&& in)
{       callPref(eblas_dscal)(in->nElem, -1., ((double*)in->dataPref(false))+1, 2); //negate the imaginary parts
        return in;
}
template<class T> Tptr conj(const Tptr& in) { return conj(clone(in)); }

template<class T> Tptr& operator*=(Tptr& in, const Tptr& other)
{       in->scale *= other->scale;
        callPref(eblas_zmul)(in->nElem, other->dataPref(false), 1, in->dataPref(false), 1);
        return in;
}
ScalarField& operator*=(ScalarField& in, const ScalarField& other); 
template<class T> Tptr operator*(const Tptr& in1, const Tptr& in2) { Tptr out(in1->clone()); return out *= in2; } 
template<class T> Tptr operator*(const Tptr& in1, Tptr&& in2) { return in2 *= in1; } 
template<class T> Tptr operator*(Tptr&& in1, const Tptr& in2) { return in1 *= in2; } 
template<class T> Tptr operator*(Tptr&& in1, Tptr&& in2) { return in1 *= in2; } 

//Extra operators in R-space alone for mixed complex-real elementwise multiplications:
complexScalarField& operator*=(complexScalarField&, const ScalarField&); 
complexScalarField operator*(const complexScalarField&, const ScalarField&);
complexScalarField operator*(const ScalarField&, const complexScalarField&);
complexScalarField operator*(complexScalarField&&, const ScalarField&);
complexScalarField operator*(const ScalarField&, complexScalarField&&);

//Extra operators in G-space alone for real kernel multiplication:
ScalarFieldTilde& operator*=(ScalarFieldTilde&, const RealKernel&); 
ScalarFieldTilde operator*(const RealKernel&, const ScalarFieldTilde&); 
ScalarFieldTilde operator*(const ScalarFieldTilde&, const RealKernel&); 
ScalarFieldTilde operator*(const RealKernel&, ScalarFieldTilde&&); 
ScalarFieldTilde operator*(ScalarFieldTilde&&, const RealKernel&); 


//------------------------------ Linear combine operators ------------------------------

template<typename T> void axpy(double alpha, const Tptr& X, Tptr& Y)
{       if(X)
        {       if(Y)
                {       if(Y->scale == 0.0) { Y = X * alpha; }
                        else callPref(eblas_zaxpy)(X->nElem, alpha*X->scale/Y->scale, X->dataPref(false), 1, Y->dataPref(false), 1);
                }
                else Y = X * alpha;
        }
        //if X is null, nothing needs to be done, Y remains unchanged
}
void axpy(double alpha, const ScalarField& X, ScalarField& Y); 
template<class T> Tptr& operator+=(Tptr& in, const Tptr& other) { axpy(+1.0, other, in); return in; } 
template<class T> Tptr& operator-=(Tptr& in, const Tptr& other) { axpy(-1.0, other, in); return in; } 
template<class T> Tptr operator+(const Tptr& in1, const Tptr& in2) { Tptr out(in1->clone()); return out += in2; } 
template<class T> Tptr operator+(const Tptr& in1, Tptr&& in2) { return in2 += in1; } 
template<class T> Tptr operator+(Tptr&& in1, const Tptr& in2) { return in1 += in2; } 
template<class T> Tptr operator+(Tptr&& in1, Tptr&& in2) { return in1 += in2; } 
template<class T> Tptr operator-(const Tptr& in1, const Tptr& in2) { Tptr out(in1->clone()); return out -= in2; } 
template<class T> Tptr operator-(const Tptr& in1, Tptr&& in2) { return (in2 -= in1) *= -1.0; } 
template<class T> Tptr operator-(Tptr&& in1, const Tptr& in2) { return in1 -= in2; } 
template<class T> Tptr operator-(Tptr&& in1, Tptr&& in2) { return in1 -= in2; } 
template<class T> Tptr operator-(const Tptr& in) { return (-1.0)*in; } 
template<class T> Tptr operator-(Tptr&& in) { return in*=(-1.0); } 
//Extra operators in R-space alone for scalar additions:
ScalarField& operator+=(ScalarField&, double); 
ScalarField operator+(double, const ScalarField&); 
ScalarField operator+(const ScalarField&, double); 
ScalarField operator+(double, ScalarField&&); 
ScalarField operator+(ScalarField&&, double); 
ScalarField& operator-=(ScalarField&, double); 
ScalarField operator-(double, const ScalarField&); 
ScalarField operator-(const ScalarField&, double); 
ScalarField operator-(double, ScalarField&&); 
ScalarField operator-(ScalarField&&, double); 


//------------------------------ Norms and dot products ------------------------------

template<typename T> complex dot(const Tptr& X, const Tptr& Y) 
{       return callPref(eblas_zdotc)(X->nElem, X->dataPref(), 1, Y->dataPref(), 1); 
}
template<typename T> double nrm2(const Tptr& X) 
{       return callPref(eblas_dznrm2)(X->nElem, X->dataPref(), 1); 
}
template<typename T> complex sum(const Tptr& X) 
{       FieldData dataOne(X->gInfo, "complexScalarField", 1, 2, false); *((complex*)dataOne.data()) = 1.0;
        return callPref(eblas_zdotc)(X->nElem, (complex*)dataOne.dataPref(), 0, X->dataPref(), 1); 
}
//Special handling for real scalar fields:
double dot(const ScalarField&, const ScalarField&); 
double dot(const ScalarFieldTilde&, const ScalarFieldTilde&); 
double nrm2(const ScalarField&); 
double nrm2(const ScalarFieldTilde&); 
double sum(const ScalarField&); 
double sum(const ScalarFieldTilde&); 

double integral(const ScalarField&); 
double integral(const ScalarFieldTilde&); 
complex integral(const complexScalarField&); 
complex integral(const complexScalarFieldTilde&); 

//------------------------------ Grid conversion utilities ------------------------------
ScalarFieldTilde changeGrid(const ScalarFieldTilde&, const GridInfo& gInfoNew); //Fourier up/down-sample to get to the new grid
ScalarField changeGrid(const ScalarField&, const GridInfo& gInfoNew); //Fourier up/down-sample to get to the new grid (real-space version)
complexScalarFieldTilde changeGrid(const complexScalarFieldTilde&, const GridInfo& gInfoNew); //Fourier up/down-sample to get to the new grid
complexScalarField changeGrid(const complexScalarField&, const GridInfo& gInfoNew); //Fourier up/down-sample to get to the new grid (real-space version)

//------------------------------ Initialization utilities ------------------------------

#include <string.h>

template<typename T> void initZero(Tptr& X) { X->zero(); }
template<typename T> void initZero(Tptr& X, const GridInfo& gInfo) { if(X) X->zero(); else nullToZero(X, gInfo); }
template<typename T> void nullToZero(Tptr& X, const GridInfo& gInfo) { if(!X) { X=T::alloc(gInfo,isGpuEnabled()); initZero(X); } } 
void initRandom(ScalarField&, double cap=0.0); 
void initRandomFlat(ScalarField&); 
void initGaussianKernel(RealKernel&, double x0); 
void initTranslation(ScalarFieldTilde&, const vector3<>& r);  
ScalarFieldTilde gaussConvolve(const ScalarFieldTilde&, double sigma); 
ScalarFieldTilde gaussConvolve(ScalarFieldTilde&&, double sigma); 

template<typename Func, typename... Args> void applyFuncGsq(const GridInfo& gInfo, const Func& f, Args... args);

template<typename Func, typename... Args> void applyFunc_r(const GridInfo& gInfo, const Func& f, Args... args);

//------------------------------ Debug utilities ------------------------------

void printStats(const ScalarField& X, const char* name, FILE* fp=stdout); 


template<typename Callable, typename Vec> void checkSymmetry(Callable* func, const Vec& v1, const Vec& v2, const char* funcName)
{       double dot1 = dot(v1, (*func)(v2));
        double dot2 = dot(v2, (*func)(v1));
        double dotDiff = fabs(dot1-dot2);
        logPrintf("Relative error in symmetry of %s: %le\n", funcName, dotDiff/sqrt(fabs(dot1)*fabs(dot2)));
}


#undef Tptr



template<typename Func, typename... Args>
void applyFuncGsq_sub(size_t iStart, size_t iStop, const vector3<int> S, const matrix3<> GGT, const Func* f, Args... args)
{       THREAD_halfGspaceLoop( (*f)(i, GGT.metric_length_squared(iG), args...); )
}
template<typename Func, typename... Args> void applyFuncGsq(const GridInfo& gInfo, const Func& f, Args... args)
{       threadLaunch(applyFuncGsq_sub<Func,Args...>, gInfo.nG, gInfo.S, gInfo.GGT, &f, args...);
}

template<typename Func, typename... Args>
void applyFunc_r_sub(size_t iStart, size_t iStop, const vector3<int> S, const vector3<> h[3], const Func* f, Args... args)
{       THREAD_rLoop
        (       vector3<> ri = iv[0]*h[0] + iv[1]*h[1] + iv[2]*h[2];
                (*f)(i, ri, args...);
        )
}
template<typename Func, typename... Args> void applyFunc_r(const GridInfo& gInfo, const Func& f, Args... args)
{       threadLaunch(applyFunc_r_sub<Func,Args...>, gInfo.nr, gInfo.S, gInfo.h, f, args...);
}
#endif //JDFTX_CORE_OPERATORS_H