Commit 7e9b887a authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

CorrectionVectorAction: supports complex Correction Vector

parent 35fea50f
......@@ -90,6 +90,7 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
#include "NoPthreadsNg.h"
#include "TridiagRixsStatic.h"
#include "KrylovHelper.h"
#include "CorrectionVectorAction.h"
namespace Dmrg {
......@@ -146,84 +147,19 @@ public:
public:
class Action {
typedef CorrectionVectorAction<ComplexOrRealType, TargetParamsType> ActionType;
public:
typedef typename ThisType::ModelType::SolverParamsType SolverParamsType;
typedef typename ThisType::MatrixComplexOrRealType MatrixComplexOrRealType;
typedef typename ThisType::VectorWithOffsetType VectorWithOffsetType;
typedef typename PsimagLite::Vector<RealType>::Type VectorRealType;
enum ActionEnum {ACTION_IMAG, ACTION_REAL};
Action(const TargetParamsType& tstStruct,
RealType E0,
const VectorRealType& eigs)
: tstStruct_(tstStruct),E0_(E0),eigs_(eigs)
{}
RealType operator()(SizeType k) const
{
if (tstStruct_.omega().first == PsimagLite::FREQ_REAL)
return actionWhenReal(k);
return actionWhenMatsubara(k);
}
void setReal() const
{
action_ = ACTION_REAL;
}
void setImag() const
{
action_ = ACTION_IMAG;
}
private:
RealType actionWhenReal(SizeType k) const
{
RealType sign = (tstStruct_.type() == 0) ? -1.0 : 1.0;
RealType part1 = (eigs_[k] - E0_)*sign + tstStruct_.omega().second;
RealType denom = part1*part1 + tstStruct_.eta()*tstStruct_.eta();
return (action_ == ACTION_IMAG) ? tstStruct_.eta()/denom :
-part1/denom;
}
RealType actionWhenMatsubara(SizeType k) const
{
RealType sign = (tstStruct_.type() == 0) ? -1.0 : 1.0;
RealType wn = tstStruct_.omega().second;
RealType part1 = (eigs_[k] - E0_)*sign;
RealType denom = part1*part1 + wn*wn;
return (action_ == ACTION_IMAG) ? wn/denom : -part1 / denom;
}
const TargetParamsType& tstStruct_;
RealType E0_;
const VectorRealType& eigs_;
mutable ActionEnum action_;
};
public:
typedef Action ActionType;
CalcR(const TargetParamsType& tstStruct,
RealType E0,
const VectorRealType& eigs)
CalcR(const TargetParamsType& tstStruct, RealType E0, const VectorRealType& eigs)
: action_(tstStruct,E0,eigs)
{}
const Action& imag() const
const ActionType& imag() const
{
action_.setImag();
return action_;
}
const Action& real() const
const ActionType& real() const
{
action_.setReal();
return action_;
......@@ -231,10 +167,17 @@ public:
private:
Action action_;
ActionType action_;
};
struct TypeWrapper {
typedef typename ThisType::MatrixComplexOrRealType MatrixComplexOrRealType;
typedef typename ThisType::VectorWithOffsetType VectorWithOffsetType;
typedef typename ModelType::SolverParamsType SolverParamsType;
};
typedef KrylovHelper<typename CalcR::ActionType> KrylovHelperType;
typedef KrylovHelper<typename CalcR::ActionType, TypeWrapper> KrylovHelperType;
CorrectionVectorSkeleton(InputValidatorType& ioIn,
const TargetParamsType& tstStruct,
......@@ -363,6 +306,8 @@ private:
xi.resize(n);
psimag::BLAS::GEMV('N',n,n2,zone,&(V(0,0)),n,&(tmp[0]),1,zzero,&(xi[0]),1);
if (CalcR::ActionType::isValueComplex()) return;
krylovHelper_.calcR(r, what.real(), T, V, phi, n2, i0);
psimag::BLAS::GEMV('N',n2,n2,zone,&(T(0,0)),n2,&(r[0]),1,zzero,&(tmp[0]),1);
......
......@@ -5,15 +5,15 @@
namespace Dmrg {
template<typename ActionType>
template<typename ActionType, typename TypeWrapperType>
class KrylovHelper {
public:
typedef typename ActionType::MatrixComplexOrRealType MatrixComplexOrRealType;
typedef typename ActionType::VectorWithOffsetType VectorWithOffsetType;
typedef typename TypeWrapperType::MatrixComplexOrRealType MatrixComplexOrRealType;
typedef typename TypeWrapperType::VectorWithOffsetType VectorWithOffsetType;
typedef typename ActionType::VectorRealType VectorRealType;
typedef typename ActionType::SolverParamsType SolverParamsType;
typedef typename TypeWrapperType::SolverParamsType SolverParamsType;
typedef typename VectorWithOffsetType::VectorType VectorType;
typedef typename VectorWithOffsetType::value_type ComplexOrRealType;
......
......@@ -133,13 +133,17 @@ public:
struct Action {
typedef typename ThisType::VectorRealType VectorRealType;
};
struct TypeWrapper {
typedef typename ThisType::MatrixComplexOrRealType MatrixComplexOrRealType;
typedef typename ThisType::VectorWithOffsetType VectorWithOffsetType;
typedef typename ThisType::VectorRealType VectorRealType;
typedef typename ModelType::SolverParamsType SolverParamsType;
};
typedef KrylovHelper<Action> KrylovHelperType;
typedef KrylovHelper<Action, TypeWrapper> KrylovHelperType;
TimeVectorsKrylov(const SizeType& currentTimeStep,
const TargetParamsType& tstStruct,
......
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