Gmres2D.cpp

Go to the documentation of this file.

/*--------------------------------*- VM2D -*-----------------*---------------*\
| ##  ## ##   ##  ####  #####   |                            | Version 1.14   |
| ##  ## ### ### ##  ## ##  ##  |  VM2D: Vortex Method       | 2026/03/06     |
| ##  ## ## # ##    ##  ##  ##  |  for 2D Flow Simulation    *----------------*
|  ####  ##   ##   ##   ##  ##  |  Open Source Code                           |
|   ##   ##   ## ###### #####   |  https://www.github.com/vortexmethods/VM2D  |
|                                                                             |
| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova   |
*-----------------------------------------------------------------------------*
| File name: Gmres2D.cpp                                                      |
| Info: Source code of VM2D                                                   |
|                                                                             |
| This file is part of VM2D.                                                  |
| VM2D 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.                                         |
|                                                                             |
| VM2D 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 VM2D.  If not, see <http://www.gnu.org/licenses/>.               |
\*---------------------------------------------------------------------------*/
 
 
#include "Gmres2D.h"
 
#include "Airfoil2D.h"
#include "Boundary2D.h"
#include "MeasureVP2D.h"
#include "Mechanics2D.h"
#include "Preprocessor.h"
#include "StreamParser.h"
#include "Velocity2D.h"
#include "Wake2D.h"
#include "World2D.h"
 
 
using namespace VM2D;
 
typedef Eigen::Map<const Eigen::VectorXd> MapVec;
 
bool GmresSolver::IterRot(const double nrmRhs, double& gs, int m, bool residualShow)
{
    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;
}
 
 
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)
{
    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;
}
 
 
 
 
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)
{
    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);
}
 
 
 
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)
{
    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];
}
 
 
 
 
 
 
 
 
 
 
 
 
void GmresSolver::SolCircleRun(std::vector<double>& AX, const std::vector<double>& rhs, int p)
{
    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];
}
 
 
 
void GmresSolver::SolM(std::vector<double>& AX, const std::vector<double>& rhs, int p)
{   
    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];
}
 
 
 
// ���������� ������������� ��� ���������� ������������������� �� ������ ��������
void GmresSolver::PreCalculateCoef(int aflIndex) 
{
    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());
            }
        }
    }
}
 
GmresSolver::GmresSolver(const World2D& W_)
    : 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
};
 
GmresSolver::~GmresSolver()
{
#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
}
 
 
#ifdef USE_CUDA
void GmresSolver::GMRES(
    std::vector<std::vector<double>>& X,
    std::vector<double>& R,
    const std::vector<std::vector<double>>& rhs,
    const std::vector<double>& rhsReg,
    int& niter
)
{
    VMlib::vmTimer tAll("All"), tPre("Pre"), tIter("Iter"), tPost("Post"), tIter1("Iter1"), tIterOther("IterOther"), 
        tCG("tCG"), tGC("tGC"), tWrapper("tWrapper"), tSplit("tSplit"), tPrecond("tPrecond"), tRot("tRot"), 
        tRotA("tRotA"), tRotB("tRotB"), tRotC("tRotC"), tRotD("tRotD");
    
    tAll.start();
    tPre.start();
 
    size_t m;   
    double beta;
    const size_t nAfl = W.getNumberOfAirfoil();
    bool linScheme = (W.getPassport().numericalSchemes.boundaryCondition.second == 2);
    int order = W.getPassport().numericalSchemes.gmresMultipoleOrder;
    double theta = W.getPassport().numericalSchemes.gmresTheta;
    double itheta2 = 1.0 / (theta * theta);
    auto& treePnlInfl = *W.getNonConstCuda().inflTreePnlVortex;
 
    if (W.isAnyMovableOrDeformable() || W.getCurrentStep() == 0)
        treePnlInfl.MemoryAllocateForGMRES((W.getCurrentStep() == 0), itheta2);
 
    if (W.getCurrentStep() == 0 || W.isAnyMovableOrDeformable())
#pragma omp parallel for
    for (int p = 0; p < (int)nAfl; ++p)
    {
        size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
        for (int i = 0; i < nPanelsP; ++i)
        {
            double lenI = W.getAirfoil(p).len[i];
            diag[p][i] = 0.5 / lenI;
 
            if (linScheme)
                diag[p][i + nPanelsP] = (1.0 / 24.0) / lenI;
        }
    }
 
    double nrmRhs = 0.0;
#pragma omp parallel for reduction(+:nrmRhs)
    for (int p = 0; p < (int)nAfl; ++p)
        nrmRhs += norm2(rhs[p]);
    nrmRhs = sqrt(nrmRhs);
 
    c.resize(0);
    s.resize(0);
    Vflat.resize(0);
    for (size_t i = 0; i < nAfl; ++i)
        Vflat.insert(Vflat.end(), rhs[i].begin(), rhs[i].end());
 
    for (size_t i = 0; i < nAfl; ++i)
        Vflat.push_back(rhsReg[i]);
 
    
    if (W.getCurrentStep() == 0 || W.isAnyMovableOrDeformable())
    {
#pragma omp parallel for
        for (int p = 0; p < nAfl; ++p)
            PreCalculateCoef(p);
    }
    
 
#pragma omp parallel for
    for (int p = 0; p < (int)nAfl; ++p)
    {
        std::vector<double>& vbuf1 = allBuf2[p];
        size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
 
        for (size_t i = 0; i < nPanelsP; ++i)
        {
            vbuf1[i] = Vflat[np[p] + i];
            if (linScheme)
                vbuf1[i + nPanelsP] = Vflat[np[p] + nPanelsP + i];
        }
 
        if (!linScheme)
            vbuf1[nPanelsP] = Vflat[(totalVsize + nAfl) - (nAfl - p)];
        else
            vbuf1[2 * nPanelsP] = Vflat[(totalVsize + nAfl) - (nAfl - p)];
 
        
        //�������� 
        SolM(vbuf1, vbuf1, p);
        if (linScheme)
            SolCircleRun(vbuf1, vbuf1, p);
 
        for (size_t j = np[p]; j < np[p] + nPanelsP; ++j)
        {
            Vflat[j] = vbuf1[j - np[p]];
            if (linScheme)
                Vflat[j + nPanelsP] = vbuf1[j - np[p] + nPanelsP];
        }
        Vflat[(totalVsize + nAfl) - (nAfl - p)] = vbuf1[vbuf1.size() - 1];
    }
 
    beta = 0;
#pragma omp parallel for reduction(+:beta)
    for (int t = 0; t < (totalVsize + nAfl); ++t)
        beta += Vflat[t] * Vflat[t];
 
    beta = sqrt(beta);  
    if (beta > 0)
#pragma omp parallel for
        for (int t = 0; t < (totalVsize + nAfl); ++t)
            Vflat[t] /= beta;
    else
    {
        for (auto& xp : X)
            for (auto& x : xp)
                x = 0.0;    
        return;
    }
 
    double gs = beta;
    
    double* dev_ptr_rhs = W.getAirfoil(0).devRhsPtr;
    double* dev_ptr_rhsLin = (linScheme ? W.getAirfoil(0).devRhsLinPtr : nullptr);
 
 
 
    const double zero = 0.0, unit = 1.0;
    
    tPre.stop();
    tIter.start();
 
 
 
    //Iterations
    for (int j = 0; j < totalVsize - 1; ++j) //+ n.size() 
    {
        if (j == 0)
            tIter1.start();
        else if (j == 1)
            tIterOther.start();
 
        tCG.start();
 
        size_t npred = 0;
        if (!linScheme)
            for (size_t i = 0; i < W.getNumberOfAirfoil(); ++i)
            {
#pragma omp parallel for
                for (int p = 0; p < (int)W.getAirfoil(i).getNumberOfPanels(); ++p)
                    bufcurrentSol[npred + p] = Vflat[(totalVsize + nAfl)*(j) + npred + p] / W.getAirfoil(i).len[p];
                
                npred += W.getAirfoil(i).getNumberOfPanels();
            }
        else
        {
            for (size_t i = 0; i < W.getNumberOfAirfoil(); ++i)
            {
#pragma omp parallel for
                for (int p = 0; p < (int)W.getAirfoil(i).getNumberOfPanels(); ++p)
                {
                    bufcurrentSol[npred + p] = Vflat[(totalVsize + nAfl) * j + (2 * npred + p)] / W.getAirfoil(i).len[p];
                    bufcurrentSol[nTotPan + npred + p] = Vflat[(totalVsize + nAfl) * j + (2 * npred + W.getAirfoil(i).getNumberOfPanels() + p)] / W.getAirfoil(i).len[p];
                }
                npred += W.getAirfoil(i).getNumberOfPanels();
            }
            W.getNonConstCuda().SetSolution(bufcurrentSol.data() + nTotPan, W.getNonConstCuda().dev_solLin, nTotPan);
        }
        W.getNonConstCuda().SetSolution(bufcurrentSol.data(), W.getNonConstCuda().dev_sol, nTotPan);
        
        tCG.stop();
 
        if (j>0)
            tWrapper.start();
        
        treePnlInfl.UpdatePanelFreeVortexIntensity(W.getNonConstCuda().dev_sol, W.getNonConstCuda().dev_solLin);
        treePnlInfl.UpwardTraversal(order);
        treePnlInfl.DownwardTraversalGMRES(dev_ptr_rhs, dev_ptr_rhsLin, itheta2, order, j, (W.isAnyMovableOrDeformable() || W.getCurrentStep()==0));
        
        if (j>0)
            tWrapper.stop();
 
        tGC.start();
 
        if (!linScheme)
            W.getCuda().CopyMemFromDev<double, 1>(nTotPan, dev_ptr_rhs, w.data(), 22);
 
        if (linScheme)
        {
            W.getCuda().CopyMemFromDev<double, 1>(nTotPan, dev_ptr_rhs, bufnewSol.data(), 22);
            W.getCuda().CopyMemFromDev<double, 1>(nTotPan, dev_ptr_rhsLin, bufnewSol.data() + nTotPan, 22);
        }
 
        tGC.stop();
 
        tSplit.start();
 
        npred = 0;
        if (linScheme)
            for (size_t i = 0; i < W.getNumberOfAirfoil(); ++i)
            {
#pragma omp parallel for
                for (int j = 0; j < (int)W.getAirfoil(i).getNumberOfPanels(); ++j)
                {
                    w[2 * npred + j] = bufnewSol[npred + j];
                    w[2 * npred + W.getAirfoil(i).getNumberOfPanels() + j] = bufnewSol[nTotPan + npred + j];
                }
                npred += W.getAirfoil(i).getNumberOfPanels();
            }
 
        size_t cntr = 0;
        for (size_t pi = 0; pi < nAfl; ++pi)
        {
            size_t nPanelsPi = W.getAirfoil(pi).getNumberOfPanels();
#pragma omp parallel for
            for (int i = 0; i < (int)nPanelsPi; ++i)
            {
                w[cntr + i] += Vflat[(totalVsize + nAfl) * j + (totalVsize + pi)] - Vflat[(totalVsize + nAfl) * j + (cntr + i)] * diag[pi][i];
                if (linScheme)
                    w[cntr + i + nPanelsPi] += -Vflat[(totalVsize + nAfl) * j + (cntr + i + nPanelsPi)] * diag[pi][i + nPanelsPi];
            }
 
            cntr += nPanelsPi;
            if (linScheme)
                cntr += nPanelsPi;
        }
 
        for (size_t i = 0; i < nAfl; ++i)
            w[totalVsize + i] = 0.0;
 
        cntr = 0;
        for (size_t i = 0; i < nAfl; ++i)
        {
            size_t nPanelsI = W.getAirfoil(i).getNumberOfPanels();
            for (size_t k = 0; k < nPanelsI; ++k)
                w[totalVsize + i] += Vflat[(totalVsize + nAfl) * j + (cntr + k)];
            cntr += nPanelsI;
            if (linScheme)
                cntr += nPanelsI;
        }
        tSplit.stop();
 
        tPrecond.start();
 
        // PRECONDITIONER   
         
#pragma omp parallel for
        for (int p = 0; p < (int)nAfl; ++p)
        {
            std::vector<double>& vbuf2 = allBuf2[p];
 
            size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
            
            for (size_t i = 0; i < nPanelsP; ++i)
            {
                vbuf2[i] = w[np[p] + i];
                if (linScheme)
                    vbuf2[i + nPanelsP] = w[np[p] + nPanelsP + i];
            }
 
            if (!linScheme)
                vbuf2[nPanelsP] = w[w.size() - (nAfl - p)];
            else
                vbuf2[2 * nPanelsP] = w[w.size() - (nAfl - p)];
                        
 
            //��������
            SolM(vbuf2, vbuf2, p);
            if (linScheme)
                SolCircleRun(vbuf2, vbuf2, p);          
 
            for (size_t j = np[p]; j < np[p] + nPanelsP; ++j)
            {
                w[j] = vbuf2[j - np[p]];
                if (linScheme)
                    w[j + nPanelsP] = vbuf2[j - np[p] + nPanelsP];
            }
 
            w[w.size() - (nAfl - p)] = vbuf2[vbuf2.size() - 1];
        }
        tPrecond.stop();
 
 
 
        tRot.start();
 
        tRotA.start();
 
        if (j == H.size() - 1)
        {
            std::cout << "Reallocation for GMRES" << std::endl;
            H.resize(j * 2 + 2);
            for (int i = 0; i < j * 2 + 2; ++i)
                H[i].resize(j * 2 + 1);
        
            Vflat.reserve((totalVsize + nAfl) * j * 2);
            c.reserve(j * 2);
            s.reserve(j * 2);
 
            //mul.reserve(j * 2);
            //cudaFree(mulptr);
            //cudaMalloc((void**)&mulptr, j * 2 * sizeof(double));
 
            double* devViBund_new;
            cudaMalloc((void**)&devViBund_new, (totalVsize + nAfl) * j * 2 * sizeof(double));
            cudaMemcpy(devViBund_new, devViBund, (totalVsize + nAfl) * j * sizeof(double), cudaMemcpyDeviceToDevice);
            cudaFree(devViBund);
            devViBund = devViBund_new;
        }
 
        tRotA.stop();       
        tRotB.start();
                
        //blas
        /*
        cudaMemcpy(devViBund + w.size() * j, Vflat.data() + w.size() * j, (1) * w.size() * sizeof(double), cudaMemcpyHostToDevice); //was       
        cudaMemcpy(devw, w.data(), w.size() * sizeof(double), cudaMemcpyHostToDevice); //was
        cublasDgemv(cublas_handle, CUBLAS_OP_T, (int)w.size(), (int)Vflat.size() / (int)w.size(), &unit, devViBund, (int)w.size(), devw, 1, &zero, mulptr, 1);
        cudaMemcpy(mul.data(), mulptr, (j+1) * sizeof(double), cudaMemcpyDeviceToHost);
 
        for (int i = 0; i <= j; ++i)
        {   
            double scal = 0.0;  
            scal = mul[i];
            H[i][j] = scal; 
            scal *= -1;         
            
            cublasDaxpy(cublas_handle, \
            (int)w.size(), \
            &scal, \
            devViBund + w.size() * i, 1, \
            devw, 1);
        }
        cudaMemcpy(w.data(), devw, w.size() * sizeof(double), cudaMemcpyDeviceToHost);
        //*/
 
        //half-blas
        /*
        cudaMemcpy(devViBund + w.size() * j, Vflat.data() + w.size() * j, (1) * w.size() * sizeof(double), cudaMemcpyHostToDevice); //was       
        cudaMemcpy(devw, w.data(), w.size() * sizeof(double), cudaMemcpyHostToDevice); //was
         
        cublasDgemv(cublas_handle, CUBLAS_OP_T, (int)w.size(), (int)Vflat.size() / (int)w.size(), &unit, devViBund, (int)w.size(), devw, 1, &zero, mulptr, 1);
        cudaMemcpy(mul.data(), mulptr, (j+1) * sizeof(double), cudaMemcpyDeviceToHost);
 
        for (int i = 0; i <= j; ++i)
        {
            double scal = 0.0;
            scal = mul[i];
            H[i][j] = scal;
            scal *= -1;
 
#pragma omp for
            for (int q = 0; q < (int)w.size(); ++q)
                w[q] += scal * Vflat[(totalVsize + nAfl) * i + q];//V[i][q];
 
            //cublasDaxpy(cublas_handle, \
            (int)w.size(), \
            &scal, \
            devViBund + w.size() * i, 1, \
            devw, 1);
        }
        //cudaMemcpy(w.data(), devw, w.size() * sizeof(double), cudaMemcpyDeviceToHost);
        //*/
 
 
        //non-blas
        //*
        for (int i = 0; i <= j; ++i)
        {
            double scal = 0.0;
#pragma omp parallel
            {
#pragma omp for reduction(+:scal)
                for (int q = 0; q < (int)w.size(); ++q)
                    scal += w[q] * Vflat[(totalVsize + nAfl) * i + q];//V[i][q];
 
                H[i][j] = scal;
#pragma omp for
                for (int q = 0; q < (int)w.size(); ++q)
                    w[q] -= scal * Vflat[(totalVsize + nAfl) * i + q];//V[i][q];
            }           
        }//*/
 
 
        tRotB.stop();
 
        tRotC.start();
    
        m = j + 1;
        
        double nrmw = 0.0;
        
        //cublasDdot(cublas_handle, \
        //  w.size(), \
        //  devw, 1,  /* ��������� �� w, ��� 1 */ \
        //  devw, 1,  /* ��������� �� w, ��� 1 */ \
        //  &nrmw);
        //nrmw = sqrt(nrmw);        
        //H[m][j] = nrmw;
 
        nrmw = norm(w);
        H[m][j] = nrmw;
 
        Vflat.insert(Vflat.end(), w.begin(), w.end());
        for (int t = 0; t < (totalVsize + nAfl); ++t)
            Vflat[Vflat.size() - t - 1] /= nrmw;    
 
        tRotC.stop();
 
        tRotD.start();
 
        if (IterRot(nrmRhs, gs, (int)m, false))
        {
            //std::cout << "iterations in GMRES = " << j + 1 << std::endl;
            if (j != 0) 
                tIterOther.stop();
            tRotD.stop();
            tRot.stop();
 
            break;
        }
 
        tRotD.stop();
        tRot.stop();
 
        if (j == 0)
            tIter1.stop();
    } // end of iterations
    
    tIter.stop();
    
    tPost.start();
 
 
    niter = (int)m;
 
    g.resize(m + 1);
    g[0] = beta;
 
    //GivensRotations
    double oldValue;
    for (int i = 0; i < m; i++)
    {
        oldValue = g[i];
        g[i] = c[i] * oldValue;
        g[i + 1] = -s[i] * oldValue;
    }
    //end of GivensRotations
 
    // Solve HY=g
    Y.resize(m);
    Y[m - 1] = g[m - 1] / H[m - 1][m - 1];
 
    double sum;
    for (int k = (int)m - 2; k >= 0; --k)
    {
        sum = 0.0;
        for (int s = k + 1; s < m; ++s)
            sum += H[k][s] * Y[s];
        Y[k] = (g[k] - sum) / H[k][k];
    }
    // end of Solve HY=g
 
    size_t cntr = 0;
    for (size_t p = 0; p < nAfl; p++) {
 
        size_t nPanelsP = W.getAirfoil(p).getNumberOfPanels();
 
        if (!linScheme)
            for (size_t i = 0; i < nPanelsP; i++)
            {
                sum = 0.0;
                for (size_t j = 0; j < m; j++)
                    sum += Vflat[(totalVsize + nAfl) * j + (i + cntr)] * Y[j];
                X[p][i] += sum;
            }
        else {
            for (size_t i = 0; i < 2 * nPanelsP; i++)
            {
                sum = 0.0;
                for (size_t j = 0; j < m; j++)
                    sum += Vflat[(totalVsize + nAfl) * j + (i + cntr)] * Y[j];
                X[p][i] += sum;
            }
        }
        cntr += nPanelsP;
 
        if (linScheme)
            cntr += nPanelsP;
 
    }
 
    sum = 0.0;
    for (size_t p = 0; p < nAfl; p++)
    {
        for (size_t j = 0; j < m; j++)
            sum += Vflat[(totalVsize + nAfl) * j + (totalVsize + p)] * Y[j];
        R[p] += sum;
    }
    
    tPost.stop();
 
    tAll.stop();
 
    /*
    std::cout << "tAll   = " << tAll.duration() << std::endl;
    std::cout << "tPre   = " << tPre.duration() << std::endl;
    std::cout << "tIter  = " << tIter.duration() << std::endl;
    std::cout << "tIter1 = " << tIter1.duration() << std::endl;
    if(m > 1)
        std::cout << "tIterN = " << tIterOther.duration() / (m-1.0) << std::endl;
    else 
        std::cout << "tIterN = " << 0.0 << std::endl;
    std::cout << "tPost  = " << tPost.duration() << std::endl;
    
    std::cout << "tCG      = " << tCG.duration() / (double)(m) << std::endl;
    std::cout << "tGC      = " << tGC.duration() / (double)(m) << std::endl;
    std::cout << "tWrap(>0)= " << tWrapper.duration() / (double)(m - 1) << std::endl;
    std::cout << "tSplit   = " << tSplit.duration() / (double)(m) << std::endl;
    std::cout << "tPrecond = " << tPrecond.duration() / (double)(m) << std::endl;
    std::cout << "tRot     = " << (tRotA.duration() + tRotB.duration() + tRotC.duration() +tRotD.duration()) / (double)(m) << std::endl;
    //*/
 
    //std::cout << "   tRotA = " << tRotA.duration() / (double)(m) << std::endl;
    //std::cout << "   tRotB = " << tRotB.duration() / (double)(m) << std::endl;
    //std::cout << "   tRotC = " << tRotC.duration() / (double)(m) << std::endl;
    //std::cout << "   tRotD = " << tRotD.duration() / (double)(m) << std::endl;
}
 
#endif