SingAnalyser.hpp

/* -----------------------------------------------------------------------------
 *  Diana process modelling, simulation and analysis software
 *  Copyright (c) 2005, Michael Krasnyk
 *  All rights reserved.

 *  This file is a part of Diana process modelling, simulation and analysis software
 *
 *  Diana 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 (see accompanying file LICENSE)
 * -----------------------------------------------------------------------------
 *  $Id: SingAnalyser.hpp 8160 2009-03-16 18:36:30Z miha $
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef SINGULARITY_ANALYSER__HPP
#define SINGULARITY_ANALYSER__HPP

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <Diana.hpp>
#include <DianaContinuation.hpp>
#include <DianaSparseSolver.hpp>
#include <DianaLSQRSolver.hpp>

#ifdef HAVE_UFSPARSE_CHOLMOD
#include <DianaNormalEquationsSolver.hpp>
#endif

namespace Diana {

  class SingAnalyser : public DianaContinuation {
  public:
    SingAnalyser(Numeric::Solvers::Model::ICapeNumericModel* _pModel);
    virtual ~SingAnalyser();

    Numeric::Solvers::Solver::SolveReturn Solve()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadInvOrder,
             Common::Error::ECapeSolvingError, Common::Error::ECapeOutOfResources);

    ContiReturn Continuate()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeSolvingError);

    Common::Types::CapeString GetVariableName(Common::Types::CapeLong index)
      throw (Common::Error::ECapeUnknown);

  protected:
    int getModelEquations();

    int getModelVariables();

    void packModelBoundaries(Common::Types::CapeArrayDouble& arrMinY,
                             Common::Types::CapeArrayDouble& arrMaxY);

    void packModelVector(Common::Types::CapeArrayDouble& modelVector);

    void unpackModelVector(const Common::Types::CapeArrayDouble& modelVector);

    void packModelResidual(Common::Types::CapeArrayDouble& modelResidual);

    void packModelJacobian(Diana::DianaSparseArray& Jaco);

    void addDerivativeByParameter(Diana::DianaSparseArray& Jaco, int nJacCol, int nParamIndex);

    CapeBoolean calcTestCondition(CapeBoolean bFirstCall=false);

  protected:
    // methods defined in the SingAnalyserEigensolver.cpp
    bool areEigenVectorsConsistent(IDianaLinearAlgebraSolver& dssaA);
    CapeBoolean calcInitialEigenpairs();
    void calcArpackEigenpairs(CapeBoolean bTranspose, IDianaLinearAlgebraSolver& dssA,
                              Common::Types::CapeDouble&      dblEigenvalue,
                              Common::Types::CapeArrayDouble& arrEigenvector);
    void calcLapackEigenpairs(const Diana::DianaSparseArray& dsaA);
  protected:
    // methods defined in the SingAnalyserConditions.cpp
    Common::Types::CapeLong getLambdaParameterIndex();
    void getConditionCheckParameters();
    
    double calcGx();
    void addGxJacobianRow(int nPosRow, Diana::DianaSparseArray& dsa);

    double calcGxx();
    void addGxxJacobianRow(int nPosRow, Diana::DianaSparseArray& dsa);
    void addGxxJacobianRowS(int nPosRow, Diana::DianaSparseArray& dsa);
    void addGxxJacobianRowFD(int nPosRow, Diana::DianaSparseArray& dsa);    

    double calcGp(int nParamIndex);
    double calcGp();
    void addGpJacobianRow(int nPosRow, Diana::DianaSparseArray& dsa);

    double calcGxp();
    void addGxpJacobianRow(int nPosRow, Diana::DianaSparseArray& dsa);
    void calcLm1EPhi(int nParamIndex, Common::Types::CapeArrayDouble& vecLm1EResult);

  protected:
    Common::Types::CapeLong               nConditionCheck;
    Common::Types::CapeLong               nConditionShow;
    Common::Types::CapeLong               nConditionCurrent;
    Common::Types::CapeDouble             dblTestGx;
    Common::Types::CapeDouble             dblTestGxx;
    DianaParameter*                       parTestGp;
    Common::Types::CapeArrayDouble        vecTestGp;
    DianaArray                            arrTestGp;
    Common::Types::CapeDouble             dblTestGxp;

    Common::Types::CapeLong               nSingularityEquations;

    DianaParameter*                       parLambdaParameter;
    DianaParameter*                       parConditionCheckParameters;
    DianaParameter*                       parConditionParameter;
    Common::Types::CapeLong               nLambdaParameter;
    Common::Types::CapeArrayString        arrConditionCheckParametersNames;
    Common::Types::CapeArrayLong          arrConditionCheckParameters;

  private:
    Common::Types::CapeDouble             dblNAN;
    Common::Types::CapeLong               nESOSize;
    
    Common::Types::CapeArrayDouble        arrEigenvector;
    Common::Types::CapeArrayDouble        arrAdjEigenvector;
    Common::Types::CapeDouble             dblBeta;
    Common::Types::CapeDouble             dblGamma;

    DianaSparseSolver                     dssEigJacobian;
    DianaLSQRSolver                       dssLm1E;
    //DianaSparseSolver                     dssLm1E;   // no MNS
    //DianaSparseArray                      dsaLm1E;   // no MNS
  };
};

#endif // SINGULARITY_ANALYSER__HPP