VM2D::GmresSolver Class Reference

#include <Gmres2D.h>

Collaboration diagram for VM2D::GmresSolver:

Public Member Functions
	GmresSolver (const World2D &W_)
	~GmresSolver ()
bool	IterRot (const double nrmRhs, double &gs, int m, bool residualShow)
	Контроль невязки после выполнения очередной итерации GMRES.
void	PreCalculateCoef (int aflIndex)
void	SolCircleRundirect (const std::vector< double > &A, const std::vector< double > &rhs, size_t startRow, size_t startRowReg, size_t np, std::vector< double > &res)
void	SolMdirect (const std::vector< double > &A, const std::vector< double > &rhs, size_t startRow, size_t startRowReg, size_t np, std::vector< double > &res, bool lin)
void	GMRES_Direct (int nAllVars, int nafl, const std::vector< double > &mtrDir, const std::vector< double > &rhsDir, const std::vector< int > &pos, const std::vector< int > &vsize, std::vector< std::vector< double > > &gam, std::vector< double > &R)
void	SolM (std::vector< double > &AX, const std::vector< double > &rhs, int p)
void	SolCircleRun (std::vector< double > &AX, const std::vector< double > &rhs, int p)

Private Attributes
const World2D &	W
std::vector< double >	w
std::vector< double >	c
std::vector< double >	s
std::vector< double >	Vflat
std::vector< std::vector< double > >	H
std::vector< std::vector< double > >	diag
std::vector< std::vector< Point2D > >	prea
std::vector< std::vector< Point2D > >	prec
std::vector< std::vector< Point2D > >	prea1
std::vector< std::vector< Point2D > >	prec1
std::vector< sweepVectors >	allSW
std::vector< std::vector< double > >	allBuf2
std::vector< size_t >	np
std::vector< double >	bufnewSol
std::vector< double >	bufcurrentSol
std::vector< double >	g
std::vector< double >	Y
size_t	totalVsize
size_t	nTotPan
const size_t	iterSize = 100

Detailed Description

Definition at line 145 of file Gmres2D.h.

Constructor & Destructor Documentation

GmresSolver()

GmresSolver::GmresSolver

(

const World2D &

W_

)

Definition at line 602 of file Gmres2D.cpp.

    : W(W_)
{
#ifdef USE_CUDA
    cublasCreate(&cublas_handle);
#endif
 
    bool linScheme = (W.getPassport().numericalSchemes.boundaryCondition.second == 2);
    const size_t nAfl = W.getNumberOfAirfoil();
 
    nTotPan = 0;
    for (size_t s = 0; s < nAfl; ++s)
        nTotPan += W.getAirfoil(s).getNumberOfPanels();
 
    totalVsize = linScheme ? nTotPan * 2 : nTotPan;
 
#ifdef USE_CUDA
    cudaMalloc((void**)&devViBund, (totalVsize + nAfl) * sizeof(double) * iterSize);
    cudaMalloc((void**)&devw, (totalVsize + nAfl) * sizeof(double));
 
    //cudaMalloc(&mulptr, iterSize * sizeof(double));
    //std::vector<double> mul(iterSize);
#endif
 
    w.resize(totalVsize + nAfl);
    c.reserve(iterSize);
    s.reserve(iterSize);
 
    Vflat.reserve((totalVsize + nAfl) * iterSize);
 
    diag.resize(nAfl);
    H.resize(iterSize + 1);
    for (int i = 0; i < iterSize + 1; ++i)
        H[i].resize(iterSize);
 
    for (int p = 0; p < (int)nAfl; ++p)
    {
        size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
        if (!linScheme)
            diag[p].resize(nPanelsP);
        else
            diag[p].resize(2 * nPanelsP);       
    }
 
    // PRECONDITIONER
    prea.resize(nAfl);
    prec.resize(nAfl);
 
    if (linScheme)
    {
        prea1.resize(nAfl); 
        prec1.resize(nAfl);
    }
 
#pragma omp parallel for
    for (int p = 0; p < nAfl; ++p)
    {
        const auto& aflP = W.getAirfoil(p);
        size_t nPanelsP = aflP.getNumberOfPanels();
 
        prea[p].resize(nPanelsP);
        prec[p].resize(nPanelsP);
        if (linScheme)
        {
            prea1[p].resize(nPanelsP);
            prec1[p].resize(nPanelsP);
        }       
    }
 
    allSW.resize(nAfl);
    allBuf2.resize(nAfl);
    np.resize(nAfl, 0);
 
    for (size_t p = 0; p < nAfl; ++p)
    {
        size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
 
        allSW[p].resize(nPanelsP);
 
        if (!linScheme)
            allBuf2[p].resize(nPanelsP + 1);
        else
            allBuf2[p].resize(2 * nPanelsP + 1);
 
        if (p == 0)
            continue;
 
        np[p] = np[p - 1];
        if (!linScheme)
            np[p] += (/*(p == 0) ? 0 :*/ W.getAirfoil(p - 1).getNumberOfPanels());
        else
            np[p] += (/*(p == 0) ? 0 :*/ 2 * W.getAirfoil(p - 1).getNumberOfPanels());
    }
 
    bufnewSol.resize(totalVsize);
    bufcurrentSol.resize(totalVsize, 0.0);
 
    g.reserve(iterSize);
    Y.reserve(iterSize);
 
#ifdef USE_CUDA 
    W.getNonConstCuda().AllocateSolution(W.getNonConstCuda().dev_sol, nTotPan);
    if (linScheme)
        W.getNonConstCuda().AllocateSolution(W.getNonConstCuda().dev_solLin, nTotPan);
    else
        W.getNonConstCuda().dev_solLin = nullptr;
#endif
};

Here is the call graph for this function:

~GmresSolver()

GmresSolver::~GmresSolver

(

)

Definition at line 711 of file Gmres2D.cpp.

{
#ifdef USE_CUDA
    cublasDestroy(cublas_handle);
 
    bool linScheme = (W.getPassport().numericalSchemes.boundaryCondition.second == 2);
    W.getNonConstCuda().ReleaseSolution(W.getNonConstCuda().dev_sol);   
    if (linScheme)
        W.getNonConstCuda().ReleaseSolution(W.getNonConstCuda().dev_solLin);
 
    //W.getNonConstCuda().inflTreePnlVortex->MemoryFreeForGMRES();
    cudaFree(devViBund);
    cudaFree(devw);
    //cudaFree(mulptr);
 
#endif
}

Here is the call graph for this function:

Member Function Documentation

GMRES_Direct()

void GmresSolver::GMRES_Direct	(	int	nAllVars,
		int	nafl,
		const std::vector< double > &	mtrDir,
		const std::vector< double > &	rhsDir,
		const std::vector< int > &	pos,
		const std::vector< int > &	vsize,
		std::vector< std::vector< double > > &	gam,
		std::vector< double > &	R
	)

Definition at line 204 of file Gmres2D.cpp.

{
    const size_t maxIter = nAllVars + 1;
    std::vector<double> x(nAllVars, 0.0);
    std::vector<double> y(nAllVars, 0.0);
    std::vector<double> residual(nAllVars, 0.0);
    std::vector<double> r(nAllVars, 0.0);
    std::vector<double> w(nAllVars, 0.0);
    std::vector<double> wOld(nAllVars, 0.0);
    std::vector<double> g(nAllVars, 0.0);
    std::vector<double> c(nAllVars, 0.0);
    std::vector<double> s(nAllVars, 0.0);
    std::vector<double> hisGs(nAllVars + 1, 0.0);
    std::vector<double> v(nAllVars * (maxIter + 1), 0.0);
 
    std::vector<double> vecFromV(nAllVars, 0.0);
 
    std::vector<double> h((nAllVars + 1) * maxIter, 0.0);
    double beta = 0.0, gs, normb = 0.0;
    size_t m = 0;
 
    for (size_t i = 0; i < nAllVars; ++i)
        normb += rhsDir[i] * rhsDir[i];
    normb = sqrt(normb);
 
    for (size_t i = 0; i < nAllVars; ++i)   
        residual[i] = rhsDir[i];
 
#pragma omp parallel for 
    for (int aflI = 0; aflI < nafl; ++aflI)
    {
        // SolMdirect(mtrDir, residual, pos[aflI], nAllVars - nafl + aflI, W.getAirfoil(aflI).getNumberOfPanels(), r, (vsize[aflI] == 2*n[aflI]));
        r = residual;
    }
 
    for (size_t i = 0; i < nAllVars; ++i)
        beta += r[i] * r[i];
    beta = sqrt(beta);
    gs = beta;
    g[0] = beta;
 
    //��������� �����, ���� ������ ����� �������
    if (beta == 0)
    {
        for (size_t aflI = 0; aflI < nafl; ++aflI)
        {
            for (size_t i = 0; i < vsize[aflI]; ++i)
                gam[aflI][i] = 0.0;
            R[aflI] = 0.0;
        }   
        return;
    }
 
    for (size_t aflI = 0; aflI < nafl; ++aflI)
        R[aflI] = x[x.size() - nafl + aflI];
 
    for (size_t i = 0; i < nAllVars; ++i)
        v[i * (maxIter + 1) + 0] = r[i] / beta;
 
    for (size_t j = 0; j < nAllVars; ++j)
    {
        double timer1 = omp_get_wtime();
 
        for (size_t p = 0; p < nAllVars; ++p)
            vecFromV[p] = v[p * (maxIter + 1) + j];
 
        //cblas_dgemv(CblasRowMajor, CblasNoTrans, nAllVars, nAllVars, 1.0, mtrDir.data(), nAllVars, &v[0 * (maxIter + 1) + j], (nAllVars + 1), 0.0, wOld.data(), 1);
        //cblas_dgemv(CblasRowMajor, CblasNoTrans, nAllVars, nAllVars, 1.0, mtrDir.data(), nAllVars, vecFromV.data(), 1, 0.0, wOld.data(), 1);
        for (int i = 0; i < nAllVars; ++i)
        {
            wOld[i] = 0.0;
            for (int jj = 0; jj < nAllVars; ++jj)
            {
                wOld[i] += mtrDir[i * nAllVars + jj] * vecFromV[jj];
            }
        }
 
        double timer2 = omp_get_wtime();
 
#pragma omp parallel for
        for (int aflI = 0; aflI < nafl; ++aflI)
        {
            // SolMdirect(mtrDir, wOld, pos[aflI], nAllVars - nafl + aflI, W.getAirfoil(aflI).getNumberOfPanels(), w, (vsize[aflI] == 2 * n[aflI]));
            w = wOld;
        }
 
        for (size_t i = 0; i <= j; ++i)
        {
            h[i * maxIter + j] = 0.0;
            for (size_t q = 0; q < nAllVars; ++q)
                h[i * maxIter + j] += w[q] * v[q * (maxIter + 1) + i];
 
            for (size_t q = 0; q < nAllVars; ++q)
                w[q] -= h[i * maxIter + j] * v[q * (maxIter + 1) + i];
        }//for i
 
        h[(j + 1) * maxIter + j] = 0.0;
        for (size_t p = 0; p < nAllVars; ++p)
            h[(j + 1) * maxIter + j] += w[p] * w[p];
        h[(j + 1) * maxIter + j] = sqrt(h[(j + 1) * maxIter + j]);
 
        for (size_t p = 0; p < nAllVars; ++p)
            v[p * (maxIter + 1) + (j + 1)] = w[p] / h[(j + 1) * maxIter + j];
 
        for (size_t i = 0; i < j; ++i)
        {
            double buf = h[i * maxIter + j];
            h[i * maxIter + j] = c[i] * buf + s[i] * h[(i + 1) * maxIter + j];
            h[(i + 1) * maxIter + j] = -s[i] * buf + c[i] * h[(i + 1) * maxIter + j];
        }
 
        double zn = sqrt(h[j * maxIter + j] * h[j * maxIter + j] + h[(j + 1) * maxIter + j] * h[(j + 1) * maxIter + j]);
        c[j] = fabs(h[j * maxIter + j] / zn);
        s[j] = c[j] * h[(j + 1) * maxIter + j] / h[j * maxIter + j];
 
        gs *= -s[j];
 
        h[j * maxIter + j] = c[j] * h[j * maxIter + j] + s[j] * h[(j + 1) * maxIter + j];
        h[(j + 1) * maxIter + j] = 0.0;
        hisGs[j] = fabs(gs) / normb;
 
        if ((fabs(gs) / normb) < W.getPassport().numericalSchemes.gmresEps)
        {
            m = j + 1;
            break;
        }
    }//for j
 
    //��������������� �������
    for (size_t i = 0; i < m; ++i)
    {
        double buf = g[i];
        g[i] *= c[i];
        if (i < (m - 1))
            g[i + 1] = -s[i] * buf;
 
    }
    y[m - 1] = g[m - 1] / h[(m - 1) * maxIter + (m - 1)];
 
    for (size_t i = m - 2; i + 1 > 0; --i)
    {
        double sum = 0;
        for (size_t j = i + 1; j < m; ++j)
            sum += h[i * maxIter + j] * y[j];
 
        y[i] = (g[i] - sum) / h[i * maxIter + i];
    }
 
    for (size_t p = 0; p < nAllVars; ++p)
    {
        for (size_t q = 0; q < m; ++q)
            x[p] += v[p * (maxIter + 1) + q] * y[q];
    }
 
    std::cout << "#iterations:     " << m << std::endl;
 
    for (size_t aflI = 0; aflI < nafl; ++aflI)
    {
        for (size_t i = 0; i < vsize[aflI]; ++i)
            gam[aflI][i] = x[pos[aflI] + i];
    }
 
    for (size_t aflI = 0; aflI < nafl; ++aflI)
        R[aflI] = x[x.size() - nafl + aflI];
}

Here is the call graph for this function:

IterRot()

bool GmresSolver::IterRot	(	const double	nrmRhs,
		double &	gs,
		int	m,
		bool	residualShow
	)

Контроль невязки после выполнения очередной итерации GMRES.

Parameters

[in,out]	H	ссылка на матрицу вращений Гивенса
[in]	rhs	константная ссылка на вектор правой части решаемой СЛАУ
[in,out]	gs	ссылка на правую часть решаемой СЛАУ

Definition at line 59 of file Gmres2D.cpp.

{
    bool fl;
    double buf;
 
    double epsGMRES = W.getPassport().numericalSchemes.gmresEps;
 
    for (int i = 0; i < m - 1; ++i)
    {
        buf = H[i][m - 1];
 
        H[i][m - 1] = c[i] * buf + s[i] * H[i + 1][m - 1];
        H[i + 1][m - 1] = -s[i] * buf + c[i] * H[i + 1][m - 1];
    }
 
    double izn = 1.0 / sqrt(H[m - 1][m - 1] * H[m - 1][m - 1] + H[m][m - 1] * H[m][m - 1]);
 
    c.push_back(H[m - 1][m - 1] * izn);
    s.push_back(H[m][m - 1] * izn);
    gs *= -s[m - 1];
 
    H[m - 1][m - 1] = c[m - 1] * H[m - 1][m - 1] + s[m - 1] * H[m][m - 1];
    H[m][m - 1] = 0.0;
 
    if (residualShow)
        std::cout << "Iteration: " << m << ", residual = " << (fabs(gs) / nrmRhs) << std::endl;
 
    fl = ((fabs(gs) / nrmRhs) < epsGMRES);
    return fl;
}

Here is the call graph for this function:

PreCalculateCoef()

void GmresSolver::PreCalculateCoef

(

int

aflIndex

)

Definition at line 532 of file Gmres2D.cpp.

{
    Point2D p1, s1, p2, s2, di, dj;
    numvector<double, 3> alpha, lambda;
    numvector<Point2D, 3> v00, v11;
    Point2D i00, i11;
 
    const Airfoil& aflP = W.getAirfoil(aflIndex);
 
    size_t nPanelsP = aflP.getNumberOfPanels();
    bool linScheme = W.getPassport().numericalSchemes.boundaryCondition.second == 2;
 
    for (size_t i = 0; i < nPanelsP; ++i)
    {
        std::array<int, 2> vecV = { ((int)i - 1 >= 0) ? int(i - 1) : int(nPanelsP - 1), (i + 1 < nPanelsP) ? int(i + 1) : 0 };
        for (int j : vecV)
        {
            p1 = aflP.getR(i + 1) - aflP.getR(j + 1);
            s1 = aflP.getR(i + 1) - aflP.getR(j);
            p2 = aflP.getR(i) - aflP.getR(j + 1);
            s2 = aflP.getR(i) - aflP.getR(j);
            di = aflP.getR(i + 1) - aflP.getR(i);
            dj = aflP.getR(j + 1) - aflP.getR(j);
 
            alpha = { \
                (aflP.isAfter(j, i)) ? 0.0 : VMlib::Alpha(s2, s1), \
                VMlib::Alpha(s2, p1), \
                (aflP.isAfter(i, j)) ? 0.0 : VMlib::Alpha(p1, p2) \
            };
 
            lambda = { \
                (aflP.isAfter(j, i)) ? 0.0 : VMlib::Lambda(s2, s1), \
                VMlib::Lambda(s2, p1), \
                (aflP.isAfter(i, j)) ? 0.0 : VMlib::Lambda(p1, p2) \
            };
 
            v00 = {
                VMlib::Omega(s1, di.unit(), dj.unit()),
                -VMlib::Omega(di, di.unit(), dj.unit()),
                VMlib::Omega(p2, di.unit(), dj.unit())
            };
 
            i00 = IDPI / di.length() * (-(alpha[0] * v00[0] + alpha[1] * v00[1] + alpha[2] * v00[2]).kcross() \
                + (lambda[0] * v00[0] + lambda[1] * v00[1] + lambda[2] * v00[2]));
 
            if (aflP.isAfter(j, i))
                prec[aflIndex][i] = (i00.kcross()) * (1.0 / dj.length());
            if (aflP.isAfter(i, j))
                prea[aflIndex][i] = (i00.kcross()) * (1.0 / dj.length());
 
            if (linScheme) {
 
                v11 = { 1.0 / (12.0 * di.length() * dj.length()) * \
                (2.0 * (s1 & Omega(s1 - 3.0 * p2, di.unit(), dj.unit())) * VMlib::Omega(s1, di.unit(), dj.unit()) - \
                    s1.length2() * (s1 - 3.0 * p2)) - 0.25 * VMlib::Omega(s1, di.unit(), dj.unit()),
                -di.length() / (12.0 * dj.length()) * Omega(di, dj.unit(), dj.unit()),
                -dj.length() / (12.0 * di.length()) * Omega(dj, di.unit(), di.unit()) };
 
                i11 = (IDPI / di.length()) * ((alpha[0] + alpha[2]) * v11[0] + (alpha[1] + alpha[2]) * v11[1] + alpha[2] * v11[2]\
                    + ((lambda[0] + lambda[2]) * v11[0] + (lambda[1] + lambda[2]) * v11[1] + lambda[2] * v11[2] \
                        + 1.0 / 12.0 * (dj.length() * di.unit() + di.length() * dj.unit() - 2.0 * VMlib::Omega(s1, di.unit(), dj.unit()))).kcross());
                if (aflP.isAfter(j, i))
                    prec1[aflIndex][i] = (i11) * (1.0 / dj.length());
                if (aflP.isAfter(i, j))
                    prea1[aflIndex][i] = (i11) * (1.0 / dj.length());
            }
        }
    }
}

Here is the call graph for this function:

SolCircleRun()

void GmresSolver::SolCircleRun	(	std::vector< double > &	AX,
		const std::vector< double > &	rhs,
		int	p
	)

Definition at line 390 of file Gmres2D.cpp.

{
    const Airfoil& afl = W.getAirfoil(p);
    int n = (int)afl.getNumberOfPanels();
    sweepVectors& sw = allSW[p];
    bool linScheme = (W.getPassport().numericalSchemes.boundaryCondition.second == 2);
    
    std::vector<double>& alpha = sw.alpha;
    std::vector<double>& beta = sw.beta;    
    std::vector<double>& delta = sw.delta;
    std::vector<double>& phi = sw.phi;
    std::vector<double>& xi = sw.xi;
    
    double yn, a;
 
    double len24 = afl.len[0] * 24.0;
 
    alpha[1] = (afl.tau[0] & prec1[p][0]) * len24;
    beta[1] = (afl.tau[0] & prea1[p][0]) * len24;
    delta[1] = -len24 * rhs[n + 0];
    double zn;
 
    for (int i = 1; i < n - 1; ++i)
    {
        a = (afl.tau[i] & prea1[p][i]);
        zn = a * alpha[i] - 1.0 / (24.0 * afl.len[i]);
        alpha[i + 1] = -(afl.tau[i] & prec1[p][i]) / zn;
        beta[i + 1] = -a * beta[i] / zn;
        delta[i + 1] = (rhs[n + i] - a * delta[i]) / zn;
    }
 
    a = afl.tau[n - 2] & prea1[p][n - 2];
    zn = alpha[n - 2] * a - 1.0 / (24.0 * afl.len[n - 2]);
    phi[n - 1] = -((afl.tau[n - 2] & prec1[p][n - 2]) + beta[n - 2] * a) / zn;
    xi[n - 1] = (rhs[n + n - 2] - delta[n - 2] * a) / zn;
 
    for (int i = n - 2; i > 0; --i)
    {
        phi[i] = alpha[i] * phi[i + 1] + beta[i];
        xi[i] = alpha[i] * xi[i + 1] + delta[i];
    }
 
 
    double e = (afl.tau[n - 1] & prec1[p][n - 1]);
    a = (afl.tau[n - 1] & prea1[p][n - 1]);
 
    yn = (rhs[n + n - 1] - e * xi[1] - a * xi[n - 1]) / (e * phi[1] + a * phi[n - 1] - 1.0 / (24.0 * afl.len[n - 1]));
 
    AX[n + n - 1] = yn;
    for (int i = n - 2; i >= 0; --i)
        AX[n + i] = phi[i + 1] * yn + xi[i + 1];
}

Here is the call graph for this function:

SolCircleRundirect()

void GmresSolver::SolCircleRundirect	(	const std::vector< double > &	A,
		const std::vector< double > &	rhs,
		size_t	startRow,
		size_t	startRowReg,
		size_t	np,
		std::vector< double > &	res
	)

Definition at line 91 of file Gmres2D.cpp.

{
    size_t nAllVars = rhs.size();
    std::vector<double> alpha(np + 2), beta(np + 2), gamma(np + 2);
    std::vector<double> u(np), v(np);
 
    double zn = A[(startRow + np + 1) * nAllVars + (startRow + np + 1)];
    alpha[2] = -A[(startRow + np + 1) * nAllVars + (startRow + np + 2)] / zn;
    beta[2] = rhs[startRow + np + 1] / zn;
    gamma[2] = -A[(startRow + np + 1) * nAllVars + (startRow + np + 0)] / zn;
 
    for (size_t i = 2; i <= np; ++i)
    {
        zn = A[(startRow + np + (i % np)) * nAllVars + (startRow + np + (i % np))] + alpha[i] * A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i - 1) % np))];
        alpha[i + 1] = -A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i + 1) % np))] / zn;
        beta[i + 1] = (rhs[startRow + np + (i % np)] - beta[i] * A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i - 1) % np))]) / zn;
        gamma[i + 1] = -(gamma[i] * A[(startRow + np + (i % np)) * nAllVars + (startRow + np + ((i - 1) % np))]) / zn;
    }
 
    u[np - 1] = beta[np];
    v[np - 1] = alpha[np] + gamma[np];
    for (size_t i = np - 2; i >= 1; --i)
    {
        u[i] = alpha[i + 1] * u[i + 1] + beta[i + 1];
        v[i] = alpha[i + 1] * v[i + 1] + gamma[i + 1];
    }
 
    res[startRow + np] = (beta[np + 1] + alpha[np + 1] * u[1]) / (1.0 - gamma[np + 1] - alpha[np + 1] * v[1]);
    for (size_t i = 1; i < np; ++i)
        res[startRow + np + i] = u[i] + res[startRow + np] * v[i];
 
    return;
}

Here is the caller graph for this function:

SolM()

void GmresSolver::SolM	(	std::vector< double > &	AX,
		const std::vector< double > &	rhs,
		int	p
	)

Definition at line 445 of file Gmres2D.cpp.

{   
    const Airfoil& afl = W.getAirfoil(p);
    int n = (int)afl.getNumberOfPanels();
    sweepVectors& sw = allSW[p];
    bool linScheme = (W.getPassport().numericalSchemes.boundaryCondition.second == 2);  
 
    std::vector<double>& alpha = sw.alpha;
    std::vector<double>& beta = sw.beta;
    std::vector<double>& gamma = sw.gamma;
    std::vector<double>& delta = sw.delta;
    std::vector<double>& phi = sw.phi;
    std::vector<double>& xi = sw.xi;
    std::vector<double>& psi = sw.psi;
 
    double yn, ynn;
    double lam1, lam2, mu1, mu2, xi1, xi2;
    double a;
 
    double len2 = afl.len[0] * 2.0;
 
    alpha[1] = (afl.tau[0] & prec[p][0]) * len2;
    beta[1] = (afl.tau[0] & prea[p][0]) * len2;
    gamma[1] = len2;
    delta[1] = -len2 * rhs[0];
 
    double zn;
 
    for (int i = 1; i < n - 1; ++i)
    {
        a = (afl.tau[i] & prea[p][i]);
        zn = a * alpha[i] - 0.5 / afl.len[i];
        alpha[i + 1] = -(afl.tau[i] & prec[p][i]) / zn;
        beta[i + 1] = -a * beta[i] / zn;
        gamma[i + 1] = -(1.0 + a * gamma[i]) / zn;
        delta[i + 1] = (rhs[i] - a * delta[i]) / zn;
    }
 
    a = afl.tau[n - 2] & prea[p][n - 2];
    zn = alpha[n - 2] * a - 0.5 / afl.len[n - 2];
    phi[n - 1] = -((afl.tau[n - 2] & prec[p][n - 2]) + beta[n - 2] * a) / zn;
    psi[n - 1] = -(1.0 + gamma[n - 2] * a) / zn;
    xi[n - 1] = (rhs[n - 2] - delta[n - 2] * a) / zn;
 
    for (int i = n - 2; i > 0; --i)
    {
        phi[i] = alpha[i] * phi[i + 1] + beta[i];
        psi[i] = alpha[i] * psi[i + 1] + gamma[i];
        xi[i] = alpha[i] * xi[i + 1] + delta[i];
    }
 
    double e = (afl.tau[n - 1] & prec[p][n - 1]);
    a = (afl.tau[n - 1] & prea[p][n - 1]);
    lam1 = e * phi[1] + a * phi[n - 1] - 0.5 / afl.len[n - 1];
    mu1 = e * psi[1] + a * psi[n - 1] + 1.0;
    xi1 = rhs[n - 1] - e * xi[1] - a * xi[n - 1];
    lam2 = mu2 = xi2 = 0.0;
    for (int j = 0; j < n - 1; ++j)
    {
        lam2 += phi[j + 1];
        mu2 += psi[j + 1];
        xi2 -= xi[j + 1];
    }
    lam2 += 1.0;
 
    if (!linScheme)
        xi2 += rhs[n];
    else if (linScheme)
        xi2 += rhs[2 * n];
 
    zn = lam1 * mu2 - lam2 * mu1;
    yn = (xi1 * mu2 - xi2 * mu1) / zn;
    ynn = -(xi1 * lam2 - xi2 * lam1) / zn;
 
    if (!linScheme)
        AX[n] = ynn;
    else if (linScheme)
        AX[2 * n] = ynn;
 
    AX[n - 1] = yn;
    for (int i = n - 2; i >= 0; --i)
        AX[i] = phi[i + 1] * yn + psi[i + 1] * ynn + xi[i + 1];
}

Here is the call graph for this function:

SolMdirect()

void GmresSolver::SolMdirect	(	const std::vector< double > &	A,
		const std::vector< double > &	rhs,
		size_t	startRow,
		size_t	startRowReg,
		size_t	np,
		std::vector< double > &	res,
		bool	lin
	)

Definition at line 128 of file Gmres2D.cpp.

{
    std::vector<double> alpha(np), beta(np), gamma(np), delta(np), phi(np), xi(np), psi(np);
    double yn, ynn;
    double lam1, lam2, mu1, mu2, xi1, xi2;
    double a;
 
    size_t nAllVars = rhs.size();
 
    double izn = 1.0 / A[(startRow + 0) * nAllVars + (startRow + 0)];
    alpha[1] = -A[(startRow + 0) * nAllVars + (startRow + 1)] * izn;
    beta[1] = -A[(startRow + 0) * nAllVars + (startRow + np - 1)] * izn;
    gamma[1] = -A[(startRow + 0) * nAllVars + (startRowReg)] * izn;
    delta[1] = rhs[startRow] * izn;
    double zn;
 
 
    for (size_t i = 1; i < np - 1; ++i)
    {
        a = A[(startRow + i) * nAllVars + (startRow + i - 1)];
        zn = a * alpha[i] + A[(startRow + i) * nAllVars + (startRow + i)];
        alpha[i + 1] = -A[(startRow + i) * nAllVars + (startRow + i + 1)] / zn;
        beta[i + 1] = -a * beta[i] / zn;
        gamma[i + 1] = -(A[(startRow + i) * nAllVars + (startRowReg)] + a * gamma[i]) / zn;
        delta[i + 1] = (rhs[startRow + i] - a * delta[i]) / zn;
    }
 
    a = A[(startRow + (np - 2)) * nAllVars + (startRow + (np - 3))];
    zn = alpha[np - 2] * a + A[(startRow + (np - 2)) * nAllVars + (startRow + (np - 2))];
 
    phi[np - 1] = -(A[(startRow + np - 2) * nAllVars + (startRow + np - 1)] + a * beta[np - 2]) / zn;
    psi[np - 1] = -(A[(startRow + np - 2) * nAllVars + (startRowReg)] + gamma[np - 2] * a) / zn;
    xi[np - 1] = (rhs[startRow + np - 2] - delta[np - 2] * a) / zn;
 
    for (size_t i = np - 2; i > 0; --i)
    {
        phi[i] = alpha[i] * phi[i + 1] + beta[i];
        psi[i] = alpha[i] * psi[i + 1] + gamma[i];
        xi[i] = alpha[i] * xi[i + 1] + delta[i];
    }
 
 
    double e = A[(startRow + np - 1) * nAllVars + (startRow + 0)];
    a = A[(startRow + np - 1) * nAllVars + (startRow + np - 2)];
    lam1 = e * phi[1] + a * phi[np - 1] + A[(startRow + np - 1) * nAllVars + (startRow + np - 1)];
    mu1 = e * psi[1] + a * psi[np - 1] + A[(startRow + np - 1) * nAllVars + (startRowReg)];
    xi1 = rhs[startRow + np - 1] - e * xi[1] - a * xi[np - 1];
 
    lam2 = mu2 = xi2 = 0.0;
    for (size_t j = 0; j < np - 1; ++j)
    {
        lam2 += phi[j + 1];
        mu2 += psi[j + 1];
        xi2 -= xi[j + 1];
    }
    lam2 += 1.0;
 
    xi2 += rhs[startRowReg];
 
    zn = lam1 * mu2 - lam2 * mu1;
    yn = (xi1 * mu2 - xi2 * mu1) / zn;
    ynn = -(xi1 * lam2 - xi2 * lam1) / zn;
 
    res[startRowReg] = ynn;
 
    res[startRow + np - 1] = yn;
    for (size_t i = np - 2; i + 1 > 0; --i)
        res[startRow + i] = phi[i + 1] * yn + psi[i + 1] * ynn + xi[i + 1];
 
//����������� �������� ��� �������� �����
    if (lin)
        SolCircleRundirect(A, rhs, startRow, startRowReg, np, res);
}

Here is the call graph for this function:

Member Data Documentation

allBuf2

std::vector<std::vector<double> > VM2D::GmresSolver::allBuf2

private

Definition at line 163 of file Gmres2D.h.

allSW

std::vector<sweepVectors> VM2D::GmresSolver::allSW

private

Definition at line 162 of file Gmres2D.h.

bufcurrentSol

std::vector<double> VM2D::GmresSolver::bufcurrentSol

private

Definition at line 166 of file Gmres2D.h.

bufnewSol

std::vector<double> VM2D::GmresSolver::bufnewSol

private

Definition at line 166 of file Gmres2D.h.

c

std::vector<double> VM2D::GmresSolver::c

private

Definition at line 156 of file Gmres2D.h.

diag

std::vector<std::vector<double> > VM2D::GmresSolver::diag

private

Definition at line 158 of file Gmres2D.h.

g

std::vector<double> VM2D::GmresSolver::g

private

Definition at line 167 of file Gmres2D.h.

H

std::vector<std::vector<double> > VM2D::GmresSolver::H

private

Definition at line 158 of file Gmres2D.h.

iterSize

const size_t VM2D::GmresSolver::iterSize = 100

private

Definition at line 172 of file Gmres2D.h.

np

std::vector<size_t> VM2D::GmresSolver::np

private

Definition at line 164 of file Gmres2D.h.

nTotPan

size_t VM2D::GmresSolver::nTotPan

private

Definition at line 170 of file Gmres2D.h.

prea

std::vector<std::vector<Point2D> > VM2D::GmresSolver::prea

private

Definition at line 160 of file Gmres2D.h.

prea1

std::vector<std::vector<Point2D> > VM2D::GmresSolver::prea1

private

Definition at line 160 of file Gmres2D.h.

prec

std::vector<std::vector<Point2D> > VM2D::GmresSolver::prec

private

Definition at line 160 of file Gmres2D.h.

prec1

std::vector<std::vector<Point2D> > VM2D::GmresSolver::prec1

private

Definition at line 160 of file Gmres2D.h.

s

std::vector<double> VM2D::GmresSolver::s

private

Definition at line 156 of file Gmres2D.h.

totalVsize

size_t VM2D::GmresSolver::totalVsize

private

Definition at line 169 of file Gmres2D.h.

Vflat

std::vector<double> VM2D::GmresSolver::Vflat

private

Definition at line 157 of file Gmres2D.h.

W

const World2D& VM2D::GmresSolver::W

private

Definition at line 148 of file Gmres2D.h.

w

std::vector<double> VM2D::GmresSolver::w

private

Definition at line 156 of file Gmres2D.h.

Y

std::vector<double> VM2D::GmresSolver::Y

private

Definition at line 167 of file Gmres2D.h.

---

The documentation for this class was generated from the following files:

• VM2D/Gmres2D/Gmres2D.h
• VM2D/Gmres2D/Gmres2D.cpp