src/eigsolv/JDBSYM.cpp

Go to the documentation of this file.

//**************************************************************************
//
//                                 NOTICE
//
// This software is a result of the research described in the report
//
// " A comparison of algorithms for modal analysis in the absence
//   of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,
//  Sandia National Laboratories, Technical report SAND2003-1028J.
//
// It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
// framework ( http://software.sandia.gov/trilinos/ ).
//
// The distribution of this software follows also the rules defined in Trilinos.
// This notice shall be marked on any reproduction of this software, in whole or
// in part.
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This program 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.
//
// Code Authors: U. Hetmaniuk (ulhetma@sandia.gov), R. Lehoucq (rblehou@sandia.gov)
//
//**************************************************************************
/**************************************************************************
 *                                                                         *
 *               Swiss Federal Institute of Technology (ETH)               *
 *                       CH-8092 Zuerich, Switzerland                      *
 *                                                                         *
 *                       (C) 1999 All Rights Reserved                      *
 *                                                                         *
 *                                NOTICE                                   *
 *                                                                         *
 *  Permission to use, copy, modify, and distribute this software and      *
 *  its documentation for any purpose and without fee is hereby granted    *
 *  provided that the above copyright notice appear in all copies and      *
 *  that both the copyright notice and this permission notice appear in    *
 *  supporting documentation.                                              *
 *                                                                         *
 *  Neither the Swiss Federal Institute of Technology nor the author make  *
 *  any representations about the suitability of this software for any     *
 *  purpose.  This software is provided ``as is'' without express or       *
 *  implied warranty.                                                      *
 *                                                                         *
 ***************************************************************************/

#include <iomanip>
#include "rlog/rlog.h"
#include "Teuchos_ParameterList.hpp"
#include "JDBSYM.h"

JDBSYM::JDBSYM(const Epetra_Comm &_Comm, const Epetra_Operator *KK,
               const Epetra_Operator *MM, const Epetra_Operator *PP,
               const Epetra_Operator *YHCP,
               Teuchos::ParameterList& params) :
    myVerify(_Comm),
    callBLAS(),
    callLAPACK(),
    ortho(),
    modalTool(_Comm),
    mySort(),
    MyComm(_Comm),
    K(KK),
    M(MM),
    Prec(PP),
    MyWatch(_Comm),
    YHC(YHCP),
    linSolver(0),
    projK(0),
    projP(0),
    historyCount(0),
    resHistory(0),
    maxSpaceSize(0),
    sumSpaceSize(0),
    spaceSizeHistory(0),
    maxIterPCG(0),
    sumIterPCG(0),
    iterPCGHistory(0),
    memRequested(0.0),
    highMem(0.0),
    numCorrectionSolve(0),
    numOrtho(0),
    outerIter(0),
    precOp(0),
    residual(0),
    stifOp(0),
    timeBuildG(0.0),
    timeBuildH(0.0),
    timeCorrectionRHS(0.0),
    timeCorrectionSolve(0.0),
    timeLocalProj(0.0),
    timeLocalSolve(0.0),
    timeLocalUpdate(0.0),
    timeNorm(0.0),
    timeOrtho(0.0),
    timeOuterLoop(0.0),
    timePostProce(0.0),
    timePrecOp(0.0),
    timeResidual(0.0),
    timeRestart(0.0),
    timeStifOp(0.0),
    params_(params)
{}

JDBSYM::~JDBSYM() {
    delete[] resHistory;
    resHistory = 0;
    delete[] spaceSizeHistory;
    spaceSizeHistory = 0;
    delete[] iterPCGHistory;
    iterPCGHistory = 0;
    delete linSolver;
    linSolver = 0;
    delete projK;
    projK = 0;
    delete projP;
    projP = 0;
}


/****************************************************************************
 *                                                                          *
 * Main eigensolver routine                                                 *
 *                                                                          *
 ****************************************************************************/


int JDBSYM::jdbsym(double tau, double jdtol,
                   int kmax, int jmax, int jmin, int blockSize, int blkwise,
                   int V0dim, double *V0,
                   opType _optype, double eps_tr, double toldecay,
                   int strategy, int max_it, int max_inner_it,
                   int clvl,
                   int &k_conv, double *Q, double *lambda, int &it) {

    //***************************************************************************
    //                                                                          *
    // Local variables                                                          *
    //                                                                          *
    //***************************************************************************/

    highMem = (highMem > currentSize()) ? highMem : currentSize();

    double* normWeight = 0;         // disable M-inverse norm (for now)

    // allocatables:
    // initialize with 0, so we can free even unallocated ptrs */
    double *V = 0, *Vtmp = 0, *U = 0, *Qb = 0, *Y = 0;
    double *s = 0, *Res = 0, *resnrm = 0;
    double *matK = 0, *H = 0, *Hlu = 0, *G = 0, *Glu = 0;
    double *temp1 = 0, *temp3 = 0;

    int *Hpiv = 0, *Gpiv = 0, *idx1 = 0, *idx2 = 0;
    int *convind = 0, *keepind = 0, *solvestep = 0;
    int *actcorrits = 0;

    // non-allocated ptrs */
    double *matdummy;

    // scalar vars */
    double alpha, it_tol;

    int k, j, actblockSize, found, conv, keep, Gisused;
    int Misempty, Pisempty, act, cnt, idummy, info;

    // Define Epetra_Map and the global and local sizes */
    Epetra_Map MyMap = K->OperatorDomainMap();

    int nGlobal = MyMap.NumGlobalElements();
    int nLocal = MyMap.NumMyElements();

    int myPid = MyComm.MyPID();

    //***************************************************************************
    //                                                                          *
    // Execution starts here...                                                 *
    //                                                                          *
    //***************************************************************************

    // Set "derived" parameters
    Misempty = (M == 0);
    Pisempty = (Prec == 0);

    // print info header
    if (clvl > 0) {
        ostringstream buf;
        buf << "JDBSYM parameters:" << endl;
        buf << "\tSolving  ";
        if (Misempty)
            buf << "K*x = lambda*x  ";
        else
            buf << "K*x = lambda*M*x  ";
        if (Pisempty)
            buf << "without ";
        else
            buf << "with ";
        buf << "preconditioning" << endl;
        buf << "\tN=";
        buf.width(10);
        buf << nGlobal << "  ITMAX=";
        buf.width(4);
        buf << max_it << endl;
        buf << "\tKMAX=";
        buf.width(3);
        buf << kmax << "  JMIN=";
        buf.width(3);
        buf << jmin << "  JMAX=";
        buf.width(3);
        buf << jmax << "  V0DIM=";
        buf.width(3);
        buf << V0dim << endl;
        buf << "\tBLKSIZE=";
        buf.width(5);
        buf << blockSize << "  BLKWISE=";
        if (blkwise)
            buf << "TRUE\n";
        else
            buf << "FALSE\n";
        buf << "\tTAU=";
        buf.precision(4);
        buf.setf(ios::scientific, ios::floatfield);
        buf << tau << "  JDTOL= " << jdtol << " STRATEGY= " << strategy;
        buf << "  LINSOLVER=    QMRS  OPTYPE=";
        switch (_optype) {
        case OP_UNSYM1:
            buf << "UNSYM1" << endl;
            break;
        case OP_UNSYM2:
            buf << "UNSYM2" << endl;
            break;
        case OP_SYM:
            buf << "SYM" << endl;
            break;
        }
        buf << "\tLINITMAX=";
        buf.width(5);
        buf << max_inner_it << "  EPS_TR=" << eps_tr << "  TOLDECAY=" << toldecay;
        rInfo(buf.str().c_str());
    }

    // Validate input parameters
    assert(0 < jdtol);
    assert(0 < kmax && kmax <= nGlobal);
    assert(0 < jmax && jmax <= nGlobal);
    assert(0 < jmin && jmin < jmax);
    assert(0 <= max_it);
    assert(blockSize <= jmin);
    assert(blockSize <= jmax - jmin);
    assert(0 < blockSize && blockSize <= kmax);
    assert(blkwise == 0 || blkwise == 1);
    assert(0 <= V0dim && V0dim < jmax);
    assert(_optype == OP_UNSYM1 || _optype == OP_UNSYM2 || _optype == OP_SYM);
    assert(0 <= max_inner_it);
    assert(0.0 <= eps_tr);
    assert(1.0 < toldecay);

    // Allocating memory for matrices & vectors
    V = new (nothrow) double[nLocal*jmax];
    U = new (nothrow) double[jmax*jmax];
    s = new (nothrow) double[jmax];
    Res = new (nothrow) double[nLocal*blockSize];
    resnrm = new (nothrow) double[blockSize];
    matK = new (nothrow) double[jmax*jmax];
    Vtmp = new (nothrow) double[jmax*jmax];
    assert(V && U && s && Res && resnrm && matK && Vtmp);
    memRequested += sizeof(double)*(nLocal*jmax + 3*jmax*jmax + jmax + nLocal*blockSize +
                                    blockSize)/(1024.0*1024.0);

    // Define Epetra_MultiVector objects
    Epetra_MultiVector *VV = new Epetra_MultiVector(View, MyMap, V, nLocal, jmax);
    Epetra_MultiVector *QQ = new Epetra_MultiVector(View, MyMap, Q, nLocal, kmax);
    Epetra_MultiVector *RR = new Epetra_MultiVector(View, MyMap, Res, nLocal, blockSize);

    // Allocate matrices H, Y and G only if necessary
    Epetra_MultiVector *YY = 0;
    if (!Pisempty) {
        H = new (nothrow) double[kmax*kmax];
        Hlu = new (nothrow) double[kmax*kmax];
        Hpiv = new (nothrow) int[kmax];
        Y = new (nothrow) double[nLocal*kmax];
        assert(H && Hlu && Hpiv && Y);
        memRequested += sizeof(double)*(nLocal*kmax + 2*kmax*kmax)/(1024.0*1024.0);
        memRequested += sizeof(int)*(kmax)/(1024.0*1024.0);
        YY = new Epetra_MultiVector(View, MyMap, Y, nLocal, kmax);
    }

    Gisused = (_optype == OP_UNSYM2 || _optype == OP_SYM) && (!Misempty);
    if (Gisused) {
        G   = new (nothrow) double[kmax*kmax];
        Glu = new (nothrow) double[kmax*kmax];
        Gpiv = new (nothrow) int[kmax];
        assert(G && Glu && Gpiv);
        memRequested += sizeof(double)*(2*kmax*kmax)/(1024.0*1024.0);
        memRequested += sizeof(int)*(kmax)/(1024.0*1024.0);
    }

    // Analogous for Qb
    Epetra_MultiVector *QbQb = 0;
    if (!Misempty) {
        Qb = new (nothrow) double[nLocal*kmax];
        assert(Qb != 0);
        memRequested += sizeof(double)*(nLocal*kmax)/(1024.0*1024.0);
        QbQb = new Epetra_MultiVector(View, MyMap, Qb, nLocal, kmax);
    }

    idx1 = new (nothrow) int[jmax];
    idx2 = new (nothrow) int[jmax];
    assert(idx1 && idx2);
    memRequested += sizeof(int)*(2*jmax)/(1024.0*1024.0);

    // Idices for (non-)converged approximations
    convind = new (nothrow) int[blockSize];
    keepind = new (nothrow) int[blockSize];
    solvestep = new (nothrow) int[blockSize];
    actcorrits = new (nothrow) int[blockSize];
    assert(convind && keepind && solvestep && actcorrits);
    memRequested += sizeof(int)*(4*blockSize)/(1024.0*1024.0);

    // temp1 is used for the orthogonalization, the linear solves
    int ltemp1 = 8*nLocal;
    temp1 = new (nothrow) double[ltemp1];

    // temp3 is used for computing inner products
    int ltemp3 = (jmax > kmax) ? jmax : kmax;
    temp3 = new (nothrow) double[ltemp3];

    assert(temp1 && temp3);
    memRequested += sizeof(double)*(nLocal+jmax)/(1024.0*1024.0);

    // Allocating work-space for projection operations
    Epetra_MultiVector *copyQ = 0;
    Epetra_MultiVector *copyQb = 0;
    Epetra_MultiVector *copyY = 0;

    if (projK)
        delete projK;
    projK = new Projectors(MyComm, K, M, Prec, _optype, 0, kmax, Hpiv, Gpiv, 0.0, Hlu, Glu,
                           copyQ, copyQb, copyY, temp1, temp1 + nLocal);

    if (projP)
        delete projP;
    projP = new Projectors(MyComm, K, M, Prec, _optype, 0, kmax, Hpiv, Gpiv, 0.0, Hlu, Glu,
                           copyQ, copyQb, copyY, temp1, temp1 + nLocal);

    // Define the linear solver
    if (linSolver)
        delete linSolver;
    linSolver = new LinSolvers(projK, projP, temp1 + 2*nLocal);

    // Get the weight for approximating the M-inverse norm
    Epetra_Vector *vectWeight = 0;
    if (normWeight) {
        vectWeight = new Epetra_Vector(View, MyMap, normWeight);
    }

    //*************************************************************************
    //                                                                        *
    // Generate initial search subspace V. Vectors are taken from V0 and if   *
    // necessary randomly generated.                                          *
    //                                                                        *
    //*************************************************************************

    highMem = (highMem > currentSize()) ? highMem : currentSize();
    timeOuterLoop -= MyWatch.WallTime();

    // copy V0 to V
    memcpy(V, V0, nLocal*V0dim*sizeof(double));
    j = V0dim;

    // if V0dim < blockSize: generate additional random vectors //
    if (V0dim < blockSize) {
        Epetra_MultiVector Vinit(View, *VV, V0dim, blockSize - V0dim);
        Vinit.Random();
        j = blockSize;
    }

    // Applying projector to vectors of V:
    // v = (I - Y * H^-1 * C^T) * v
    if (YHC) {
        Epetra_MultiVector Vinit(View, *VV, 0, blockSize);
        YHC->Apply(Vinit, Vinit);
    }

    // (M-)orthogonalize columns of V //
    timeOrtho -= MyWatch.WallTime();
    if (Misempty) {
        for (cnt = 0; cnt < j; cnt ++) {
            Epetra_Vector Vcnt(View, *VV, cnt);
            if (cnt > 0) {
                Epetra_MultiVector Vprevious(View, *VV, 0, cnt);
                ortho.ModifiedGS(&Vcnt, &Vprevious);
            }
            Vcnt.Norm2(&alpha);
            Vcnt.Scale(1.0/alpha);
        }
    } // if (Misempty)
    else {
        for (cnt = 0; cnt < j; cnt ++) {
            Epetra_Vector Vcnt(View, *VV, cnt);
            if (cnt == 0) {
                ortho.IteratedClassicalGS(&Vcnt, &alpha, 0, temp1, M);
            } else {
                Epetra_MultiVector Vprevious(View, *VV, 0, cnt);
                ortho.IteratedClassicalGS(&Vcnt, &alpha, &Vprevious, temp1, M);
            }
            alpha = 1.0/alpha;
            assert(alpha > 0.0);
            Vcnt.Scale(alpha);
        }
    } // if (Misempty)
    timeOrtho += MyWatch.WallTime();
    numOrtho += j;

    // Generate interaction matrix matK = V'*K*V.
    // Note: Only the upper triangle is computed
    //       Use Res as a workspace
    for (cnt = 0; cnt < j; cnt += blockSize) {
        int numVec = (cnt + blockSize < j) ? blockSize : j - cnt;
        Epetra_MultiVector Vcnt(View, *VV, cnt, numVec);
        timeStifOp -= MyWatch.WallTime();
        K->Apply(Vcnt, *RR);
        timeStifOp += MyWatch.WallTime();
        stifOp += numVec;
        timeLocalProj -= MyWatch.WallTime();
        callBLAS.GEMM('T', 'N', j, numVec, nLocal, 1.0, V, nLocal, Res, nLocal, 0.0, Vtmp, j);
        int iC;
        for (iC = 0; iC < numVec; ++iC) {
            MyComm.SumAll(Vtmp + iC*j, matK + (cnt+iC)*jmax, j);
        }
        timeLocalProj += MyWatch.WallTime();
    } // for (cnt = 0; cnt < j; cnt += blockSize)

    // Other initializations
    k = 0; it = 0; actblockSize = blockSize;
    for(act = 0; act < blockSize; act ++)
        solvestep[act] = 1;

    //***************************************************************************
    //                                                                          *
    // Main JD-iteration loop                                                   *
    //                                                                          *
    //***************************************************************************/

    while (it < max_it) {

        highMem = (highMem > currentSize()) ? highMem : currentSize();

        //***************************************************************************
        //                                                                          *
        // Solving the projected eigenproblem                                       *
        //                                                                          *
        // matK*u = V'*K*V*u = s*u                                                  *
        // matK is symmetric, only the upper triangle is stored                     *
        //                                                                          *
        //***************************************************************************/

        timeLocalSolve -= MyWatch.WallTime();
        info = modalTool.directSolver(j, matK, jmax, 0, 0, j, U, jmax, s, 0, 10);
        timeLocalSolve += MyWatch.WallTime();

        if (info != 0)
            rError("Error solving the projected eigenproblem. dsyev: info = %d", info);
        rAssert(info == 0);

        // sort eigenpairs, such that |S(i)-tau| <= |S(i+1)-tau| for i=1..j-1
        sorteig(j, s, U, jmax, tau, temp3, idx1, idx2, strategy);

        //***************************************************************************
        //                                                                          *
        // Convergence/Restart Check                                                *
        //                                                                          *
        // In case of convergence, strip off a whole block or just the converged    *
        // ones and put 'em into Q.  Update the matrices Q, V, U, s and if          *
        // necessary Qb, Y, H and G.                                                *
        //                                                                          *
        // In case of a restart update the V, U and M matrices and recompute the    *
        // Eigenvectors                                                             *
        //                                                                          *
        //***************************************************************************

        found = 1;
        while(found) {

            // conv/keep = Number of converged/non-converged Approximations
            conv = 0; keep = 0;

            timeLocalUpdate -= MyWatch.WallTime();
            callBLAS.GEMM('N', 'N', nLocal, actblockSize, j, 1.0, V, nLocal, U, jmax,
                          0.0, Q+k*nLocal, nLocal);
            timeLocalUpdate += MyWatch.WallTime();

            // Compute the residual and update the matrix Qb if necessary.
            timeResidual -= MyWatch.WallTime();
            if (Misempty==1) {
                // M is Identity
                Epetra_MultiVector QVec(View, *QQ, k, actblockSize);
                K->Apply(QVec, *RR);
                for (act = 0; act < actblockSize; ++act)
                    callBLAS.AXPY(nLocal, -s[act], QVec.Values() + act*nLocal, Res + act*nLocal);
            } else {
                // M is NOT Identity
                Epetra_MultiVector QVec(View, *QQ, k, actblockSize);
                Epetra_MultiVector QbVec(View, *QbQb, k, actblockSize);
                M->Apply(QVec, QbVec);
                K->Apply(QVec, *RR);
                for (act = 0; act < actblockSize; ++act)
                    callBLAS.AXPY(nLocal, -s[act], QbVec.Values() + act*nLocal, Res + act*nLocal);
            }
            timeResidual += MyWatch.WallTime();
            residual += actblockSize;

            // Update G, if necessary
            timeBuildG -= MyWatch.WallTime();
            if (Gisused) {
                // Compute G(:,k+act) = Q' * q
                // Use Glu as workspace
                callBLAS.GEMM('T', 'N', k+actblockSize, actblockSize, nLocal, 1.0, Q, nLocal,
                              Q+k*nLocal, nLocal, 0.0, Glu, k+actblockSize);
                for (act = 0; act < actblockSize; ++act) {
                    int idummy = k + act;
                    MyComm.SumAll(Glu + act*(k+actblockSize), G + idummy*kmax, idummy+1);
                    int iC;
                    for (iC = 0; iC < idummy; ++iC)
                        G[k + act + iC*kmax] = G[iC + idummy*kmax];
                }
            }
            timeBuildG += MyWatch.WallTime();

            // Finally update matrices H, Y  if necessary */
            if (!Pisempty) {
                if (!Misempty) {
                    // If B exists, then also Qb does
                    matdummy = Qb;
                } else {
                    // Without B, no Qb
                    matdummy = Q;
                }
                // Calculation of Y depends on optype (see header of Projectors.c)
                timePrecOp -= MyWatch.WallTime();
                if (_optype == OP_UNSYM1) {
                    Epetra_MultiVector QVec(View, *QQ, k, actblockSize);
                    Epetra_MultiVector YVec(View, *QQ, k, actblockSize);
                    Prec->ApplyInverse(QVec, YVec);
                } else {
                    Epetra_MultiVector DVec(View, MyMap, matdummy + k*nLocal, nLocal, actblockSize);
                    Epetra_MultiVector YVec(View, MyMap, Y + k*nLocal, nLocal, actblockSize);
                    Prec->ApplyInverse(DVec, YVec);
                }
                timePrecOp += MyWatch.WallTime();
                precOp += actblockSize;

                timeBuildH -= MyWatch.WallTime();
                // update H, starting with the column ... */
                // Note: Use Hlu as workspace             */
                callBLAS.GEMM('T', 'N', k+actblockSize, actblockSize, nLocal, 1.0, matdummy, nLocal,
                              Y+k*nLocal, nLocal, 0.0, Hlu, k+actblockSize);
                for (act = 0; act < actblockSize; ++act) {
                    int idummy = k + act;
                    MyComm.SumAll(Hlu + act*(k+actblockSize), H + idummy*kmax, idummy+1);
                }
                if (kmax > 1) {
                    // ... and then the row
                    // Note: Use Hlu as workspace
                    callBLAS.GEMM('T', 'N', k+actblockSize, actblockSize, nLocal, 1.0, Y, nLocal,
                                  matdummy+k*nLocal, nLocal, 0.0, Hlu, k+actblockSize);
                    for (act = 0; act < actblockSize; ++act) {
                        int idummy = k + act;
                        MyComm.SumAll(Hlu + act*(k+actblockSize), temp3, idummy);
                        int iC;
                        for (iC = 0; iC < idummy; ++iC)
                            H[k + act + iC*kmax] = temp3[iC];
                    }
                }
                timeBuildH += MyWatch.WallTime();
            } // if (!Pisempty)

            // Compute norm of the residual and update arrays convind/keepind*/
            timeNorm -= MyWatch.WallTime();
            if (vectWeight)
                RR->NormWeighted(*vectWeight, resnrm);
            else
                RR->Norm2(resnrm);
            for (act = 0; act < actblockSize; ++act) {
                //resnrm[act] = resnrm[act]/s[act];
                if (resnrm[act] < jdtol) { convind[conv] = act; conv = conv + 1; } else { keepind[keep] = act; keep = keep + 1; }
            }
            timeNorm += MyWatch.WallTime();

            // Perform some (expensive runtime validation checks) */
            if (clvl > 2) {
                Epetra_MultiVector copyQ(View, *QQ, 0, k + actblockSize);
                Epetra_MultiVector copyQb(View, *QbQb, 0, k + actblockSize);
                Epetra_MultiVector copyY(View, *YY, 0, k + actblockSize);
                Epetra_MultiVector copyV(View, *VV, 0, j);
                validate(&copyQ, &copyQb, &copyY, &copyV, H, kmax, _optype, temp1);
            }

            // Check whether the blkwise-mode is chosen and ALL the
            // approximations converged, or whether the strip-off mode is
            // active and SOME of the approximations converged

            found = ((blkwise==1 && conv==actblockSize) || (blkwise==0 && conv!=0))
                && (j > actblockSize || k == kmax - actblockSize);

            //**************************************************************************
            //*                                                                        *
            //* Convergence Case                                                       *
            //*                                                                        *
            //* In case of convergence, strip off a whole block or just the converged  *
            //* ones and put 'em into Q.  Update the matrices Q, V, U, s and if        *
            //* necessary Qb, Y, H and G.                                              *
            //*                                                                        *
            //**************************************************************************

            if (found) {

                // Store Eigenvalues */
                for(act = 0; act < conv; act++)
                    lambda[k+act] = s[convind[act]];

                // Re-use non approximated Ritz-Values */
                for(act = 0; act < keep; act++)
                    s[act] = s[keepind[act]];

                // Shift the others in the right position */
                idummy = j - actblockSize;
                for(act = 0; act < idummy; act ++)
                    s[act+keep] = s[act+actblockSize];

                // Update V. Re-use the V-Vectors not looked at yet. */
                timeLocalUpdate -= MyWatch.WallTime();
                for (act = 0; act < nLocal; act += jmax) {
                    cnt = (act + jmax < nLocal) ? jmax : nLocal - act;
                    int iC;
                    for (iC = 0; iC < j; ++iC)
                        memcpy(Vtmp + iC*cnt, V + act + iC*nLocal, cnt*sizeof(double));
                    callBLAS.GEMM('N', 'N', cnt, idummy, j, 1.0, Vtmp, cnt, U+actblockSize*jmax, jmax,
                                  0.0, V + act + keep*nLocal, nLocal);
                }
                timeLocalUpdate += MyWatch.WallTime();

                // Insert the not converged approximations as first columns in V */
                for(act = 0; act < keep; act++)
                    memcpy(V + act*nLocal, Q + (k+keepind[act])*nLocal, nLocal*sizeof(double));

                // Store Eigenvectors */
                for(act = 0; act < conv; act++)
                    if (act != convind[act])
                        memcpy(Q + (k+act)*nLocal, Q + (k+convind[act])*nLocal, nLocal*sizeof(double));
                // Update Qb if necessary */
                if (Misempty==0) {
                    for(act = 0; act < conv; act++)
                        if (act != convind[act])
                            memcpy(Qb + (k+act)*nLocal, Qb + (k+convind[act])*nLocal, nLocal*sizeof(double));
                }

                // Update H and Y if necessary */
                if (Pisempty==0) {
                    for (act = 0; act < conv; act++)
                        if (act != convind[act])
                            memcpy(Y + (k+act)*nLocal, Y + (k+convind[act])*nLocal, nLocal*sizeof(double));

                    for(act=0; act < conv; act++) {
                        if (act != convind[act]) {
                            idummy = k + act;
                            // Copy column ... */
                            memcpy(H + idummy*kmax, H + (k+convind[act])*kmax, idummy*sizeof(double));
                            // ... diagonal element ... */
                            H[idummy*(kmax+1)] = H[(k+convind[act])*(kmax+1)];
                            // ... and row */
                            int iD;
                            for (iD = 0; iD < idummy; ++iD)
                                H[idummy + iD*kmax] = H[k+convind[act] + iD*kmax];
                        }
                    }
                }

                // Update SearchSpaceSize j */
                j = j - conv;

                // Let matK become a diagonal matrix with the Ritz values as entries ... */
                for (act = 0; act < j; act++) {
                    memset(matK + act*jmax, 0, act*sizeof(double));
                    matK[act*jmax + act] = s[act];
                }

                // ... and U the Identity(jnew,jnew) */
                memset(U, 0, j*jmax*sizeof(double));
                for (act = 0; act < j; act++)
                    U[act*jmax + act] = 1.0;

                // Avoid computation of zero eigenvalues:
                //
                //   If STRATEGY == 1: set tau to the largest of the now
                //   converged eigenvalues.
                //
                //   Warning: This may not work well if BLKSIZE > 1.
                if (strategy == 1)
                    for(act = 0; act < conv; act ++)
                        if (lambda[k+act] > tau)
                            tau = lambda[k+act];

                // Update Converged-Eigenpair-counter and Pro_k */
                k = k + conv;

                // Update the new blocksize */
                actblockSize = (blockSize > kmax-k) ? kmax-k : blockSize;

                // Exit this loop when kmax eigenpairs have converged */
                if (k == kmax)
                    break;

                // Counter for the linear-solver-accuracy */
                for(act = 0; act < keep; act++)
                    solvestep[act] = solvestep[keepind[act]];

                for(act = keep; act < blockSize; act++)
                    solvestep[act] = 1;

            } // if(found)

            //*************************************************************************
            //                                                                        *
            // Restart                                                                *
            //                                                                        *
            // The Eigenvector-Aproximations corresponding to the first jmin          *
            // Petrov-Vectors are kept.                                               *
            //                                                                        *
            //*************************************************************************

            timeRestart -= MyWatch.WallTime();
            if (j+actblockSize > jmax) {

                idummy = j; j = jmin;

                for (act = 0; act < nLocal; act += jmax) {
                    cnt = (act + jmax < nLocal) ? jmax : nLocal - act;
                    int iC;
                    for (iC = 0; iC < idummy; ++iC)
                        memcpy(Vtmp + iC*cnt, V + act + iC*nLocal, cnt*sizeof(double));
                    callBLAS.GEMM('N', 'N', cnt, j, idummy, 1.0, Vtmp, cnt, U, jmax, 0.0, V + act, nLocal);
                }

                for (act = 0; act < j; act++) {
                    memset(matK + act*jmax, 0, act*sizeof(double));
                    matK[act*jmax + act] = s[act];
                    memset(U + act*jmax, 0, j*sizeof(double));
                    U[act*jmax + act] = 1.0;
                }

            }
            timeRestart += MyWatch.WallTime();

        } // while(found)

        // Exit main iteration loop when kmax eigenpairs have converged */
        if (k == kmax)
            break;

        //***************************************************************************
        //                                                                          *
        // Solving the correction equations                                         *
        //                                                                          *
        // Depending on the input-arguments we choose an appropriate lin.solver.    *
        //                                                                          *
        //***************************************************************************

        // calculate Hlu (LU-factorization), if necessary
        if (!Pisempty) {
            idummy = k + actblockSize;
            int iD;
            for (iD = 0; iD < idummy; ++iD) {
                memcpy(Hlu + iD*kmax, H + iD*kmax, idummy*sizeof(double));
            }
            callLAPACK.GETRF(idummy, idummy, Hlu, kmax, Hpiv, &info);
            if (info != 0)
                rError("Factorization of H failed: info=%d", info);
            rAssert(info == 0);
        }

        // Factorize G if necessary
        if (Gisused) {
            idummy = k + actblockSize;
            int iD;
            for (iD = 0; iD < idummy; ++iD) {
                memcpy(Glu + iD*kmax, G + iD*kmax, idummy*sizeof(double));
            }
            callLAPACK.GETRF(idummy, idummy, Glu, kmax, Gpiv, &info);
            if (info != 0)
                rError("Factorization of G failed: info=%d", info);
            rAssert(info == 0);
        }

        // Solve actblockSize times the correction equation ...
        for (act = 0; act < actblockSize; act ++) {

            // Setting start-value for vector v as zeros(n,1).
            // Guarantees (M-)orthogonality
            Epetra_Vector vVec(View, *VV, j);
            vVec.PutScalar(0.0);

            // Adaptive accuracy and shift for the lin.solver. In case the
            // residual is big, we don't need a too precise solution for the
            // correction equation, since even in exact arithmetic the
            // solution wouldn't be too usefull for the Eigenproblem.
            if (resnrm[act] < eps_tr) {
                projK->resetPro_theta(s[act]);
                projP->resetPro_theta(s[act]);
            } else {
                projK->resetPro_theta(tau);
                projP->resetPro_theta(tau);
            }

            it_tol = pow(toldecay, (double)(-solvestep[act]));
            solvestep[act] = solvestep[act] + 1;

            // Tells the Projectors-Lib what size Q, Qb, Y etc. has */
            int prok = k + actblockSize;
            projK->resetPro_k(prok);
            projP->resetPro_k(prok);

            if (copyQ)
                delete copyQ;
            copyQ = new Epetra_MultiVector(View, *QQ, 0, prok);
            projK->resetPro_Q(copyQ);
            projP->resetPro_Q(copyQ);

            if (M) {
                if (copyQb)
                    delete copyQb;
                copyQb = new Epetra_MultiVector(View, *QbQb, 0, prok);
                projK->resetPro_Qb(copyQb);
                projP->resetPro_Qb(copyQb);
            }

            if (Prec) {
                if (copyY)
                    delete copyY;
                copyY = new Epetra_MultiVector(View, *YY, 0, prok);
                projK->resetPro_Y(copyY);
                projP->resetPro_Y(copyY);
            }

            // Apply preconditioner to the right hand side of the correction
            // equation and project if necessary
            timeCorrectionRHS -= MyWatch.WallTime();
            Epetra_Vector rVec(View, MyMap, Res + act*nLocal);
            projK->MakeRHS(rVec);
            timeCorrectionRHS += MyWatch.WallTime();

            // Solve the correction equation
            timeCorrectionSolve -= MyWatch.WallTime();
            info = linSolver->solve(rVec, vVec, it_tol, max_inner_it);
            timeCorrectionSolve += MyWatch.WallTime();
            numCorrectionSolve += 1;

            // Actualizing profiling data
            if (info >= 0) {
                sumIterPCG += info;
                maxIterPCG = (info > maxIterPCG) ? info : maxIterPCG;
                actcorrits[act] = info;
            }
            if (info == -1) {
                sumIterPCG += max_inner_it;
                maxIterPCG = max_inner_it;
                actcorrits[act] = -1;
                info = 0;
            }
            if (info < 0)
                rError("Error detected in correction equation solver (info=%d).", info);
            rAssert(info >= 0);
            info = 0;

            // Applying projector :
            // v = (I - Y * H^-1 * C^T) * v

            if (YHC) {
                YHC->Apply(vVec, vVec);
            }

            // (M-)orthonormalize v to Q, cause the implicit
            // orthogonalization in the solvers may be too inaccurate. Then
            // apply "IteratedCGS" to prevent numerical breakdown
            timeOrtho -= MyWatch.WallTime();
            if (!Misempty) {
                ortho.ModifiedGSB(&vVec, copyQ, copyQb);
                Epetra_MultiVector copyV(View, *VV, 0, j);
                ortho.IteratedClassicalGS(&vVec, &alpha, &copyV, temp1, M);
            } else {
                ortho.ModifiedGS(&vVec, copyQ);
                Epetra_MultiVector copyV(View, *VV, 0, j);
                ortho.IteratedClassicalGS(&vVec, &alpha, &copyV, temp1);
            }
            alpha = 1.0 / alpha;
            vVec.Scale(alpha);
            timeOrtho += MyWatch.WallTime();
            numOrtho += 1;

            // update interaction matrix matK
            Epetra_Vector t1(View, MyMap, temp1);
            timeStifOp -= MyWatch.WallTime();
            K->Apply(vVec, t1);
            timeStifOp += MyWatch.WallTime();
            stifOp += 1;
            idummy = j+1;
            timeLocalProj -= MyWatch.WallTime();
            callBLAS.GEMV('T', nLocal, idummy, 1.0, V, nLocal, temp1, 0.0, temp3);
            MyComm.SumAll(temp3, matK + j*jmax, idummy);
            timeLocalProj += MyWatch.WallTime();

            // Increasing SearchSpaceSize j
            j ++;

        }   // for (act = 0; act < actblockSize; act ++)

        // Increase iteration-counter for outer loop
        it += 1;

        // Print information line
        print_status(clvl*(myPid == 0), it, k, j - blockSize, kmax, blockSize, actblockSize,
                     s, resnrm, actcorrits);

    } // while (it < maxIterEigenSolve)
    timeOuterLoop += MyWatch.WallTime();

    highMem = (highMem > currentSize()) ? highMem : currentSize();

    k_conv = k;

    // Free the memory
    if (copyQ)
        delete copyQ;
    if (copyQb)
        delete copyQb;
    if (copyY)
        delete copyY;
    if (VV)
        delete VV;
    if (QQ)
        delete QQ;
    if (RR)
        delete RR;
    if (YY)
        delete YY;
    if (QbQb)
        delete QbQb;
    if (V)
        delete[] V;
    if (Vtmp)
        delete[] Vtmp;
    if (U)
        delete[] U;
    if (Qb)
        delete[] Qb;
    if (Y)
        delete[] Y;
    if (s)
        delete[] s;
    if (Res)
        delete[] Res;
    if (resnrm)
        delete[] resnrm;
    if (matK)
        delete[] matK;
    if (H)
        delete[] H;
    if (Hlu)
        delete[] Hlu;
    if (Hpiv)
        delete[] Hpiv;
    if (G)
        delete[] G;
    if (Glu)
        delete[] Glu;
    if (Gpiv)
        delete[] Gpiv;
    if (temp1)
        delete[] temp1;
    if (temp3)
        delete[] temp3;
    if (idx1)
        delete[] idx1;
    if (idx2)
        delete[] idx2;
    if (convind)
        delete[] convind;
    if (keepind)
        delete[] keepind;
    if (solvestep)
        delete[] solvestep;
    if (actcorrits)
        delete[] actcorrits;
    if (vectWeight)
        delete vectWeight;

    return info;

} /* JDBSYM::jdbsym(.....) */


/****************************************************************************
 *                                                                          *
 * Supporting functions                                                     *
 *                                                                          *
 ****************************************************************************/

// VALIDATE - perform runtime check to validate JD-routine //
void JDBSYM::validate(Epetra_MultiVector *Q, Epetra_MultiVector *MQ, Epetra_MultiVector *Y,
                      Epetra_MultiVector *V, double *H, int ldh, opType _optype, double *work1) {

    int myPid = MyComm.MyPID();

    double VAL_TOL = 0.0;
    callLAPACK.LAMCH('E', VAL_TOL);
    VAL_TOL *= 100.0;

    double s = 0.0;
    double maxS;

    int r, c;

    int k = Q->NumVectors();
    int j = V->NumVectors();

    Epetra_Vector w1(View, Q->Map(), work1);

    // Check MQ = M*Q
    if (M) {
        maxS = 0.0;
        for (c = 0; c < k; c ++) {
            M->Apply(*((*Q)(c)), w1);
            w1.Update(1.0, *((*MQ)(c)), -1.0);
            w1.NormInf(&s);
            maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
        }
        if ((myPid == 0) && (maxS > VAL_TOL))
            rWarning("Norm of MQ - M*Q = %.4e", maxS);
    }

    // Check Q'*Q = I or Q'*MQ = I resp. //
    if (M == 0) {
        maxS = 0.0;
        for (r = 0; r < k; r ++) {
            for (c = 0; c < k; c ++) {
                (*Q)(r)->Dot(*((*Q)(c)), &s);
                if (r == c)
                    s -= 1.0;
                maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
            }
        }
        if ((myPid == 0) && (maxS > VAL_TOL))
            rWarning("Norm of Q'*Q - I = %.4e", maxS);
    } else {
        maxS = 0.0;
        for (r = 0; r < k; r ++) {
            for (c = 0; c < k; c ++) {
                (*MQ)(r)->Dot(*((*Q)(c)), &s);
                if (r == c)
                    s -= 1.0;
                maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
            }
        }
        if ((myPid == 0) && (maxS > VAL_TOL))
            rWarning("Norm of MQ'*Q - I = %.4e", maxS);
    }

    // Check V'*V = I or V'*M*V = I resp. //
    if (M == 0) {
        maxS = 0.0;
        for (r = 0; r < j; r ++) {
            for (c = 0; c < j; c ++) {
                (*V)(r)->Dot(*((*V)(c)), &s);
                if (r == c)
                    s -= 1.0;
                maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
            }
        }
        if ((myPid == 0) && (maxS > VAL_TOL))
            rWarning("Norm of V'*V - I = %.4e", maxS);
    } else {
        for (r = 0; r < j; r ++) {
            for (c = 0; c < j; c ++) {
                M->Apply(*((*V)(r)), w1);
                w1.Dot(*((*V)(c)), &s);
                if (r == c)
                    s -= 1.0;
                maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
            }
        }
        if ((myPid == 0) && (maxS > VAL_TOL))
            rWarning("Norm of V'*M*V - I = %.4e", maxS);
    }

    // Check Y = inv(K)*Q or Y = inv(K)*MQ resp. //
    if (Prec) {
        maxS = 0.0;
        if ((M == 0) || (_optype == OP_UNSYM1)) {
            for (c = 0; c < k; c ++) {
                Prec->ApplyInverse(*((*Q)(c)), w1);
                w1.Update(1.0, *((*Y)(c)), -1.0);
                w1.NormInf(&s);
                maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
            }
            if ((myPid == 0) && (maxS > VAL_TOL))
                rWarning("Norm of Y - Prec^{-1}*Q = %.4e", maxS);
        } else {
            for (c = 0; c < k; c ++) {
                Prec->ApplyInverse(*((*MQ)(c)), w1);
                w1.Update(1.0, *((*Y)(c)), -1.0);
                w1.NormInf(&s);
                maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
            }
            if ((myPid == 0) && (maxS > VAL_TOL))
                rWarning("Norm of Y - Prec^{-1}*M*Q = %.4e", maxS);
        }
    }

    // Check H = Q'*Y or H = MQ'*Y resp. //
    if (Prec) {
        maxS = 0.0;
        if (M == 0) {
            for (r = 0; r < k; r ++) {
                for (c = 0; c < k; c ++) {
                    (*Q)(r)->Dot(*((*Y)(c)), &s);
                    s -= H[c*ldh + r];
                    maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
                }
            }
            if ((myPid == 0) && (maxS > VAL_TOL))
                rWarning("Norm of Q'*Y - H = %.4e", maxS);
        } else {
            for (r = 0; r < k; r ++) {
                for (c = 0; c < k; c ++) {
                    (*MQ)(r)->Dot(*((*Y)(c)), &s);
                    s -= H[c*ldh + r];
                    maxS = (fabs(s) > maxS) ? fabs(s) : maxS;
                }
            }
            if ((myPid == 0) && (maxS > VAL_TOL))
                rWarning("Norm of MQ'*Y - H = %.4e", maxS);
        }
    }

}


// PRINT_STATUS - print status line (called for each outer iteration)
void JDBSYM::print_status(int clvl, int it, int k, int j, int kmax,
                          int _blkSize, int actblkSize,
                          double *s, double *resnrm, int *actcorrits) {
    const int max_vals = 5;
    ostringstream buf;

    if (clvl > 0) {
        if (_blkSize == 1) {
            if (it == 1) {
                buf << "  IT   K   J       RES CGIT  RITZVALS(1:" << max_vals << ")";
                rInfo(buf.str().c_str());
                buf.str("");
                int idummy = 27 + ( 13 > max_vals*10 ? 13 : max_vals*10);
                for (int i = 0; i < idummy; i ++)
                    buf << '-';
                rInfo(buf.str().c_str());
                buf.str("");
            }
            buf << right << scientific
                << setw(4) << it << ' '
                << setw(3) << k << ' '
                << setw(3) << j << ' '
                << setw(9) << setprecision(2) << resnrm[0] << ' '
                << setw(4) << actcorrits[0];
            for (int i = 0; i < (j < max_vals ? j : max_vals); i ++)
                buf << ' ' << setw(9) << setprecision(2) << s[i];
        } else {
            buf << " Iteration " << it;
            buf << " - Number of converged eigenvectors " << k << endl;

            if (clvl > 2) {
                for (int i = 0; i < _blkSize; ++i) {
                    resHistory[historyCount*_blkSize+i] = (i < actblkSize) ? resnrm[i] : 0.0;
                    spaceSizeHistory[historyCount] = j;
                    iterPCGHistory[historyCount*_blkSize+i] = (i < actblkSize) ? actcorrits[i] : 0;
                }
                historyCount += 1;
            }

            if (clvl > 1) {
                buf.precision(2);
                buf.setf(ios::scientific, ios::floatfield);
                for (int i = 0; i < actblkSize; ++i) {
                    buf << " Iteration " << it;
                    buf << " - Norm of Residual " << i;
                    buf << " = " << resnrm[i] << endl;
                }
                for (int i = 0; i < actblkSize; ++i) {
                    buf << " Iteration " << it << " - Inner iterations for vector " << i;
                    buf << " = " << actcorrits[i] << endl;
                }
                for (int i = 0; i < actblkSize; ++i) {
                    buf << " Iteration " << it << " - Ritz eigenvalue " << i;
                    buf << " = " << s[i];
                }
            }
        }
    }
    rInfo(buf.str().c_str());
}


/*
 * SORTEIG
 *
 * Default behaviour (strategy == 0):
 *
 *   Sort eigenpairs (S(i),U(:,i)), such that
 *
 *       |S(i) - tau| <= |S(i+1) -tau| for i=1..j-1.
 *
 *     j  : dimension of S
 *     ldu: leading dimension of U
 *   dtemp: double array of length j
 *     idx: int array of length j
 *
 * Alternate behaviour (strategy == 1):
 *
 *   Same as above but put all S(i) < tau to the end. This is used to
 *   avoid computation of zero eigenvalues.
 */


void JDBSYM::sorteig(int j, double S[], double U[], int ldu, double tau,
                     double dtemp[], int idx1[], int idx2[], int strategy) {

    int i;

    /* setup vector to be sorted and index vector */
    switch (strategy) {
    case 0:
        for (i = 0; i < j; i ++)
            dtemp[i] = fabs(S[i] - tau);
        break;
    case 1:
        for (i = 0; i < j; i ++)
            if (S[i] < tau)
                dtemp[i] = Epetra_MaxDouble;
            else
                dtemp[i] = fabs(S[i] - tau);
        break;
    default:
        assert(0);
    }
    for (i = 0; i < j; i ++)
        idx1[i] = i;

    /* sort dtemp in ascending order carrying itemp along */
    quicksort(j, dtemp, idx1);

    /* compute 'inverse' index vector */
    for (i = 0; i < j; i ++)
        idx2[idx1[i]] = i;

    /* sort eigenvalues */
    memcpy(dtemp, S, j * sizeof(double));
    for (i = 0; i < j; i ++)
        S[i] = dtemp[idx1[i]];

    /* sort eigenvectors (in place) */
    for (i = 0; i < j; i ++) {
        if (i != idx1[i]) {
            memcpy(dtemp, U+i*ldu, j*sizeof(double));
            memcpy(U+i*ldu, U+idx1[i]*ldu, j*sizeof(double));
            memcpy(U+idx1[i]*ldu, dtemp, j*sizeof(double));
            idx1[idx2[i]] = idx1[i];
            idx2[idx1[i]] = idx2[i];
        }
    }

}


/*
 * QUICKSORT 
 *
 * Sorts a double array using a non-recursive quicksort algorithm in
 * ascending order carrying along an int array.
 *  
 */


void JDBSYM::quicksort(int n, double arr[], int idx[]) {

    double v, td;
    int i, j, l, r, ti, tos, stack[32];

    l = 0; r = n-1; tos = -1;
    for (;;) {
        while (r > l) {
            v = arr[r]; i = l; j = r-1;
            for (;;) {
                while (arr[i] < v) i ++;
                /* j > l prevents underflow */
                while (arr[j] >= v && j > l) j --;
                if (i >= j) break;
                td = arr[i]; arr[i] = arr[j]; arr[j] = td;
                ti = idx[i]; idx[i] = idx[j]; idx[j] = ti;
            }
            td = arr[i]; arr[i] = arr[r]; arr[r] = td;
            ti = idx[i]; idx[i] = idx[r]; idx[r] = ti;
            if (i-l > r-i) { stack[++tos] = l; stack[++tos] = i-1; l = i+1; } else { stack[++tos] = i+1; stack[++tos] = r; r = i-1; }
            assert(tos < 32);
        }
        if (tos == -1) break;
        r = stack[tos--]; l = stack[tos--];
    }

}

void JDBSYM::set_defaults(Teuchos::ParameterList& params, int kmax) {
    // The tol parameter is set directly before constructing the solver object
    // params.set("tol", 1e-8);
    params.set("jmin", kmax < 5? 5 : kmax+1);
    params.set("jmax", kmax+10);
    params.set("tau", 0.0);
    params.set("eps_tr", 1e-3);
    params.set("tol_decay", 1.5);
    params.set("max_it", 200);
    params.set("max_inner_it", 50);
    params.set("block_size", 1);
    params.set("block_wise", 0);
    params.set("strategy", 0);
    params.set("op_type", "sym");
    params.set("clvl", 1);
}

int JDBSYM::solve(int numEigen, Epetra_MultiVector &Q, double *lambda) {

    int info = 0;

    info = myVerify.inputArguments(numEigen, K, M, Prec, Q, minimumSpaceDimension(numEigen));
    if (info < 0)
        return info;

    if ((params_.get<int>("clvl") > 2) && (MyComm.MyPID() == 0)) {
        resHistory = new double[params_.get<int>("max_it")*params_.get<int>("block_size")];
        spaceSizeHistory = new int[params_.get<int>("max_it")];
        iterPCGHistory = new int[params_.get<int>("max_it")*params_.get<int>("block_size")];
    }

    JDBSYM::opType optype;
    if (params_.get<string>("op_type") == "sym")
        optype = OP_SYM;
    else if (params_.get<string>("op_type") == "unsym1")
        optype = OP_UNSYM1;
    else if (params_.get<string>("op_type") == "unsym2")
        optype = OP_UNSYM2;
    else
        assert(false);

    int knownEV = 0;

    info = jdbsym(params_.get<double>("tau"),
                  params_.get<double>("tol"),
                  numEigen,
                  params_.get<int>("jmax"),
                  params_.get<int>("jmin"),
                  params_.get<int>("block_size"),
                  params_.get<int>("block_wise"),
                  0, // V0
                  0, // V0dim
                  optype,
                  params_.get<double>("eps_tr"),
                  params_.get<double>("tol_decay"),
                  params_.get<int>("strategy"),
                  params_.get<int>("max_it"),
                  params_.get<int>("max_inner_it"),
                  params_.get<int>("clvl"),
                  knownEV,
                  Q.Values(),
                  lambda,
                  outerIter);
    //  info = jdbsym(tau, tolEigenSolve, numEigen, 2*numEigen, numEigen, 0, 0, 0, OP_SYM,
    //                sqrt(tolEigenSolve), 1.5, 0, verbose, knownEV, Q.Values(), lambda, outerIter);

    // Sort the eigenpairs
    timePostProce -= MyWatch.WallTime();
    if ((info == 0) && (knownEV > 0)) {
        mySort.sortScalars_Vectors(knownEV, lambda, Q.Values(), Q.MyLength());
    }
    timePostProce += MyWatch.WallTime();

    return (info == 0) ? knownEV : info;
}


void JDBSYM::initializeCounters() {

    historyCount = 0;
    if (resHistory) {
        delete[] resHistory;
        resHistory = 0;
    }

    maxSpaceSize = 0;
    sumSpaceSize = 0;
    if (spaceSizeHistory) {
        delete[] spaceSizeHistory;
        spaceSizeHistory = 0;
    }

    maxIterPCG = 0;
    sumIterPCG = 0;
    if (iterPCGHistory) {
        delete[] iterPCGHistory;
        iterPCGHistory = 0;
    }

    memRequested = 0.0;
    highMem = 0.0;

    numCorrectionSolve = 0;
    numOrtho = 0;
    outerIter = 0;
    precOp = 0;
    residual = 0;
    stifOp = 0;

    timeBuildG = 0.0;
    timeBuildH = 0.0;
    timeCorrectionRHS = 0.0;
    timeCorrectionSolve = 0.0;
    timeLocalProj = 0.0;
    timeLocalSolve = 0.0;
    timeLocalUpdate = 0.0;
    timeNorm = 0.0;
    timeOrtho = 0.0;
    timeOuterLoop = 0.0;
    timePostProce = 0.0;
    timePrecOp = 0.0;
    timeResidual = 0.0;
    timeRestart = 0.0;
    timeStifOp = 0.0;
}


void JDBSYM::algorithmInfo() const {
    rDebug("Algorithm: Jacobi-Davidson algorithm (Geus' implementation)");
    rDebug("Block Size: %d", params_.get<int>("block_size"));
}

void JDBSYM::historyInfo() const {

    if (resHistory) {
        ostringstream buf;
        int j;
        buf << " Block Size: " << params_.get<int>("block_size") << endl;
        buf << endl;
        buf << " Residuals    Search Space Dim.   Inner Iter. \n";
        buf << endl;
        buf.precision(4);
        buf.setf(ios::scientific, ios::floatfield);
        for (j = 0; j < historyCount; ++j) {
            int ii;
            for (ii = 0; ii < params_.get<int>("block_size"); ++ii) {
                buf << resHistory[params_.get<int>("block_size")*j + ii] << "      ";
                buf.width(4);
                buf << spaceSizeHistory[j] << "          ";
                buf.width(4);
                buf << iterPCGHistory[j] << endl;
            }
        }
        rDebug(buf.str().c_str());
    }

}

void JDBSYM::memoryInfo() const {
    double maxHighMem = 0.0;
    double tmp = highMem;
    MyComm.MaxAll(&tmp, &maxHighMem, 1);

    double maxMemRequested = 0.0;
    tmp = memRequested;
    MyComm.MaxAll(&tmp, &maxMemRequested, 1);

    rDebug("Memory requested per processor by the eigensolver (EST) = %.2f MB", maxMemRequested);
    rDebug("High water mark in eigensolver                    (EST) = %.2f MB", maxHighMem);
}


void JDBSYM::operationInfo() const {
    ostringstream buf;
    buf << "JDBSYM operations:\n";
    buf << "\tTotal number of stiffness matrix operations      = ";
    buf.width(9);
    buf << stifOp << endl;
    buf << "\tTotal number of preconditioner applications      = ";
    buf.width(9);
    buf << precOp << endl;
    buf << "\tTotal number of linear solves                    = ";
    buf.width(9);
    buf << numCorrectionSolve << endl;
    buf.setf(ios::fixed, ios::floatfield);
    buf.precision(2);
    buf.width(9);
    buf << "\t\tAverage number of iter. per solve   : ";
    buf.width(9);
    buf << ((double) sumIterPCG)*params_.get<int>("block_size")/((double) numCorrectionSolve) << endl;
    buf << "\t\tMaximum number of iter. per solve   : ";
    buf.width(9);
    buf << maxIterPCG << endl;
    buf << "\tTotal number of computed eigen-residuals         = ";
    buf.width(9);
    buf << residual << endl;
    buf << "\tTotal number of orthogonalization steps          = ";
    buf.width(9);
    buf << numOrtho << endl;
    buf << "\tTotal number of outer iterations                 = ";
    buf.width(9);
    buf << outerIter << endl;
    buf << "\tMaximum size of search space                     = ";
    buf.width(9);
    buf << maxSpaceSize << endl;
    buf << "\tAverage size of search space                     = ";
    buf.setf(ios::fixed, ios::floatfield);
    buf.precision(2);
    buf.width(9);
    buf << ((double) sumSpaceSize)/((double) outerIter);
    rInfo(buf.str().c_str());
}


void JDBSYM::timeInfo() const {
    ostringstream buf;
    buf << "JDBSYM timings:\n";
    buf.setf(ios::fixed, ios::floatfield);
    buf.precision(2);
    buf << "\tTotal time for outer loop                         = ";
    buf.width(9);
    buf << timeOuterLoop << " s" << endl;
    buf << "\t\tTime for stiffness matrix operations      = ";
    buf.width(9);
    buf << timeStifOp << " s     ";
    buf.width(5);
    buf << 100*timeStifOp/timeOuterLoop << " %\n";
    buf << "\t\tTime for local projection                 = ";
    buf.width(9);
    buf << timeLocalProj << " s     ";
    buf.width(5);
    buf << 100*timeLocalProj/timeOuterLoop << " %\n";
    buf << "\t\tTime for local update                     = ";
    buf.width(9);
    buf << timeLocalUpdate << " s     ";
    buf.width(5);
    buf << 100*timeLocalUpdate/timeOuterLoop << " %\n";
    buf << "\t\tTime for preconditioner applications      = ";
    buf.width(9);
    buf << timePrecOp << " s     ";
    buf.width(5);
    buf << 100*timePrecOp/timeOuterLoop << " %\n";
    buf << "\t\tTime for making the correction RHS        = ";
    buf.width(9);
    buf << timeCorrectionRHS << " s     ";
    buf.width(5);
    buf << 100*timeCorrectionRHS/timeOuterLoop << " %\n";
    buf << "\t\tTime for solving the correction equation  = ";
    buf.width(9);
    buf << timeCorrectionSolve << " s     ";
    buf.width(5);
    buf << 100*timeCorrectionSolve/timeOuterLoop << " %\n";
    buf << "\t\t\tPreconditioner Mult.     : ";
    buf.width(9);
    buf << projP->getTimeApplyInverse() << " s" << endl;
    buf << "\t\t\tLinear Operator Mult.    : ";
    buf.width(9);
    buf << projK->getTimeApply() << " s" << endl;
    buf << "\t\tTime for orthogonalizations               = ";
    buf.width(9);
    buf << timeOrtho << " s     ";
    buf.width(5);
    buf << 100*timeOrtho/timeOuterLoop << " %\n";
    buf << "\t\tTime for restart                          = ";
    buf.width(9);
    buf << timeRestart << " s     ";
    buf.width(5);
    buf << 100*timeRestart/timeOuterLoop << " %\n";
    buf << "\t\tTime for for building Q^T*Q               = ";
    buf.width(9);
    buf << timeBuildG << " s     ";
    buf.width(5);
    buf << 100*timeBuildG/timeOuterLoop << " %\n";
    buf << "\t\tTime for for building Q^T*M*P^{-1}*M*Q    = ";
    buf.width(9);
    buf << timeBuildH << " s     ";
    buf.width(5);
    buf << 100*timeBuildH/timeOuterLoop << " %\n";
    buf << "\t\tTime for residual computations            = ";
    buf.width(9);
    buf << timeResidual << " s     ";
    buf.width(5);
    buf << 100*timeResidual/timeOuterLoop << " %\n";
    buf << "\tTotal time for post-processing                    = ";
    buf.width(9);
    buf << timePostProce << " s";
    rInfo(buf.str().c_str());
}

---