fetisolver.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  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.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
 
#include "sm/FETISolver/feticommunicator.h"
#include "mathfem.h"
#include "fetisolver.h"
#include "skyline.h"
#include "verbose.h"
#include "dofmanager.h"
#include "engngm.h"
#include "domain.h"
#include "unknownnumberingscheme.h"
#include "classfactory.h"
 
namespace oofem {
REGISTER_SparseLinSolver(FETISolver, ST_Feti);
 
FETISolver :: FETISolver(Domain *d, EngngModel *m) : SparseLinearSystemNM(d, m),
    ni(20),
    err(1.e-6),
    pcbuff(CBT_static),
    processCommunicator(& pcbuff, 0),
    energyNorm_comput_flag(0)
{
}
 
 
FETISolver :: ~FETISolver()
{
    if ( engngModel->giveRank() == 0 ) {
        delete commBuff;
        delete masterCommunicator;
    }
}
 
 
int
FETISolver :: estimateMaxPackSize(IntArray &map, DataStream &buff, int &packUnpackType)
{
    int rank = domain->giveEngngModel()->giveRank();
    int count = 0;
    EModelDefaultEquationNumbering dn;
 
    if ( rank == 0 ) {
        // master comm maps contain boundary dof managers
        for ( int m: map ) {
            count += masterCommunicator->giveDofManager( m )->giveNumberOfDofs();
        }
    } else {
        for ( int m: map ) {
            for ( Dof *dof: *domain->giveDofManager( m ) ) {
                if ( dof->giveEquationNumber(dn) != 0 ) {
                    count++;
                }
            }
        }
    }
 
    return ( buff.givePackSizeOfDouble(1) * count * FETISOLVER_MAX_RBM );
}
 
 
void
FETISolver :: initializeFrom(InputRecord &ir)
{
    IR_GIVE_FIELD(ir, ni, _IFT_FETISolver_maxiter);
    IR_GIVE_FIELD(ir, err, _IFT_FETISolver_maxerr);
    IR_GIVE_FIELD(ir, limit, _IFT_FETISolver_limit);
    energyNorm_comput_flag = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, energyNorm_comput_flag, _IFT_FETISolver_energynormflag);
 
    if ( fabs(limit) < 1.e-20 ) {
        limit = 1.e-20;
    }
 
    if ( err < 1.e-20 ) {
        err = 1.e-20;
    }
 
    if ( err > 0.1 ) {
        err = 0.1;
    }
}
 
void FETISolver :: setUpCommunicationMaps()
{
    int nnodes = domain->giveNumberOfDofManagers();
    int boundaryDofManNum = 0;
    int indx = 1, neq;
    StaticCommunicationBuffer commBuff(MPI_COMM_WORLD);
    IntArray commMap;
    EModelDefaultEquationNumbering dn;
 
    // determine the total number of boundary dof managers
    for ( int i = 1; i <= nnodes; i++ ) {
        if ( domain->giveDofManager(i)->giveParallelMode() == DofManager_shared ) {
            boundaryDofManNum++;
        }
    }
 
    if ( domain->giveEngngModel()->giveRank() != 0 ) {
        commBuff.resize( commBuff.givePackSizeOfInt(1) );
        commBuff.write(boundaryDofManNum);
 
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Sending data to partition 0 (send %d)\n",
                        this->giveEngngModel()->giveRank(),
                        "FETISolver :: setUpCommunicationMaps : send number of boundary dofMans", boundaryDofManNum);
#endif
        commBuff.iSend(0, FETICommunicator :: NumberOfBoundaryDofManagersMsg);
 
        MPI_Barrier(MPI_COMM_WORLD);
    }
 
#ifdef __VERBOSE_PARALLEL
    OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Barrier reached\n",
                    this->giveEngngModel()->giveRank(), "FETISolver :: setUpCommunicationMaps");
#endif
 
 
 
 
    commMap.resize(boundaryDofManNum);
 
    commBuff.resize(commBuff.givePackSizeOfInt(1) * 2 * boundaryDofManNum);
    commBuff.init();
    // determine total number of DOFs per each boundary node
    // and build communication map simultaneously
    indx = 1;
    IntArray locNum;
    if ( domain->giveEngngModel()->giveRank() != 0 ) {
        for ( int i = 1; i <= nnodes; i++ ) {
            if ( domain->giveDofManager(i)->giveParallelMode() == DofManager_shared ) {
                // remember comm map entry
                commMap.at(indx++) = i;
                // determine number of DOFs
                domain->giveDofManager(i)->giveCompleteLocationArray(locNum, dn);
                neq = 0;
                for ( int j = 1; j <= locNum.giveSize(); j++ ) {
                    if ( locNum.at(j) ) {
                        neq++;
                    }
                }
 
                commBuff.write( domain->giveDofManager(i)->giveGlobalNumber() );
                commBuff.write(neq);
            }
        }
    } else {
        for ( int i = 1; i <= nnodes; i++ ) {
            if ( domain->giveDofManager(i)->giveParallelMode() == DofManager_shared ) {
                // remember comm map entry
                commMap.at(indx++) = i;
            }
        }
    }
 
    if ( domain->giveEngngModel()->giveRank() != 0 ) {
        processCommunicator.setToSendArry(this, commMap, 0);
        processCommunicator.setToRecvArry(this, commMap, 0);
        commBuff.iSend(0, FETICommunicator :: BoundaryDofManagersRecMsg);
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Barrier reached\n",
                        this->giveEngngModel()->giveRank(), "FETISolver :: setUpCommunicationMaps");
#endif
        MPI_Barrier(MPI_COMM_WORLD);
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Barrier reached\n",
                        this->giveEngngModel()->giveRank(), "FETISolver :: setUpCommunicationMaps");
#endif
        MPI_Barrier(MPI_COMM_WORLD);
    } else {
        // remember comm map for direct mapping at master
        this->masterCommMap = commMap;
    }
}
 
 
 
 
 
/*******************************************************************/
/*******************************************************************/
/*******************************************************************/
/*******************************************************************/
 
void
FETISolver :: projection(FloatArray &v, FloatMatrix &l, FloatMatrix &l1)
/*
 * funkce provadi projekci v modifikovane metode sdruzenych gradientu
 *
 * vstupy
 * v - vektor pred projekci
 * l - matice L(m,n)
 * l1 - inverzni matice k matici L^T L, rozmery (n,n)
 *
 * vystup
 * v - vektor po projekci
 *
 * 7.12.1998
 */
{
    FloatArray help1, help2;
 
    help1.beTProductOf(l, v);
    help2.beProductOf(l1, help1);
    help1.beProductOf(l, help2);
 
    v.subtract(help1);
}
 
 
 
 
int
FETISolver :: packRBM(ProcessCommunicator &processComm)
{
    int result = 1;
    const IntArray &toSendMap = processComm.giveToSendMap();
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
    IntArray locationArray;
    EModelDefaultEquationNumbering dn;
 
    for ( int inode : toSendMap ) {
        domain->giveDofManager( inode )->giveCompleteLocationArray(locationArray, dn);
        int ndofs = locationArray.giveSize();
        for ( int j = 1; j <= ndofs; j++ ) {
            int eqNum;
            if ( ( eqNum = locationArray.at(j) ) ) {
                for ( int ir = 1; ir <= nse; ir++ ) {
                    result &= send_buff.write( rbm.at(eqNum, ir) );
                }
 
                if ( nse == 0 ) {
                    result &= send_buff.write(0.0);
                }
            }
        }
    }
 
    return result;
}
 
 
int
FETISolver :: masterUnpackRBM(ProcessCommunicator &processComm)
{
    int result = 1;
    const IntArray &toRecvMap = processComm.giveToRecvMap();
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
    IntArray locationArray;
 
    int receivedRank = processComm.giveRank();
    if ( receivedRank != 0 ) {
        //  for (irbm = 1; irbm <= nsem.at(receivedRank+1); irbm++) {
        for ( int to : toRecvMap ) {
            //
            // loop over all dofs
            for ( int idof = 1; idof <= masterCommunicator->giveDofManager(to)->giveNumberOfDofs(); idof++ ) {
                for ( int irbm = 1; irbm <= nsem.at(receivedRank + 1); irbm++ ) {
                    // unpack contribution
                    double value;
                    result &= recv_buff.read(value);
                    if ( masterCommunicator->giveDofManager(to)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                              // corresponding DOFs in boundary node
 
                        // localize to corresponding places
                        int nshared = masterCommunicator->giveDofManager(to)->giveNumberOfSharedPartitions();
                        for ( int j = 1; j <= nshared; j++ ) {
                            int part = masterCommunicator->giveDofManager(to)->giveSharedPartition(j);
                            if ( part == processComm.giveRank() ) {
                                continue;
                            }
 
                            int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(part, idof);
                            l.at(eqNum, rbmAddr.at(receivedRank + 1) + irbm - 1) = value;
                        }
                    } else { // no reference partition
                        int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(receivedRank, idof);
                        l.at(eqNum, rbmAddr.at(receivedRank + 1) + irbm - 1) = ( -1.0 ) * value;
                    }
                }
            }
        }
    }
 
    return result;
}
 
int
FETISolver :: masterMapRBM()
{
    int receivedRank = 0, result;
    IntArray locationArray;
    int size = masterCommMap.giveSize();
    // master will map its own values directly
    int locpos;
    EModelDefaultEquationNumbering dn;
 
    for ( int irbm = 1; irbm <= nsem.at(1); irbm++ ) {
        for ( int i = 1; i <= size; i++ ) {
            int to = masterCommunicator->giveMasterCommMapPtr()->at(i);
            // use receive map, send map is empty to prevent master to send
            // itself any data. Note, however, that send and receive maps are same.
            int from = masterCommMap.at(i);
 
            domain->giveDofManager(from)->giveCompleteLocationArray(locationArray, dn);
            locpos = 1;
            //
            //   ndofs = locationArray.giveSize(); // including supported
            //   for (j=1; j<=ndofs; j++) {
            //   if (eqNum = locationArray.at(j)) {
 
            //
            // loop over all dofs
            for ( int idof = 1; idof <= masterCommunicator->giveDofManager(to)->giveNumberOfDofs(); idof++, locpos++ ) {
                // extract source value
                while ( locationArray.at(locpos) == 0 ) {
                    locpos++;
                    // test if position of nonzero dof is allowed
                    if ( locpos > locationArray.giveSize() ) {
                        OOFEM_ERROR("Consistency dof error");
                    }
                }
 
                double value = rbm.at(locationArray.at(locpos), irbm);
                if ( masterCommunicator->giveDofManager(to)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                          // corresponding DOFs in boundary node
 
                    // localize to corresponding places
                    int nshared = masterCommunicator->giveDofManager(to)->giveNumberOfSharedPartitions();
                    for ( int j = 1; j <= nshared; j++ ) {
                        int part = masterCommunicator->giveDofManager(to)->giveSharedPartition(j);
                        if ( part == 0 ) {
                            continue;
                        }
 
                        int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(part, idof);
                        l.at(eqNum, rbmAddr.at(receivedRank + 1) + irbm - 1) = value;
                    }
                } else { // no reference partition
                    int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(receivedRank, idof);
                    l.at(eqNum, rbmAddr.at(receivedRank + 1) + irbm - 1) = ( -1.0 ) * value;
                }
            }
        }
    }
 
    result = 1;
    return result;
}
 
 
int
FETISolver :: packQQProducts(ProcessCommunicator &processComm)
{
    int result = 1;
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
    IntArray locationArray;
 
    for ( int i = 1; i <= nse; i++ ) {
        result &= send_buff.write( qq.at(i) );
    }
 
    return result;
}
 
 
int
FETISolver :: masterUnpackQQProduct(ProcessCommunicator &processComm)
{
    int result = 1;
    //const IntArray const &toRecvMap = processComm.giveToRecvMap();
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
    IntArray locationArray;
    //double value;
 
    int receivedRank = processComm.giveRank();
 
    if ( receivedRank != 0 ) {
        for ( int i = 1; i <= nsem.at(receivedRank + 1); i++ ) {
            result &= recv_buff.read( q.at(rbmAddr.at(receivedRank + 1) + i - 1) );
        }
    }
 
    return result;
}
 
 
int
FETISolver :: masterMapQQProduct()
{
    int result = 0;
 
    for ( int i = 1; i <= nsem.at(1); i++ ) {
        q.at(rbmAddr.at(1) + i - 1) = qq.at(i);
        result = 1;
    }
 
    return result;
}
 
 
int
FETISolver :: packSolution(ProcessCommunicator &processComm)
{
    // master
    int result = 1;
    const IntArray &toSendMap = processComm.giveToSendMap();
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
    IntArray locationArray;
 
    int rank = processComm.giveRank();
    for ( int from : toSendMap ) {
        int ndofs = masterCommunicator->giveDofManager(from)->giveNumberOfDofs();
        if ( rank == masterCommunicator->giveDofManager(from)->giveReferencePratition() ) {
            // summ corresponding values (multipliers)
            int nshared = masterCommunicator->giveDofManager(from)->giveNumberOfSharedPartitions();
            for ( int k = 1; k <= ndofs; k++ ) {
                double val = 0.0;
                for ( int j = 1; j <= nshared; j++ ) {
                    int part = masterCommunicator->giveDofManager(from)->giveSharedPartition(j);
                    if ( part == processComm.giveRank() ) {
                        continue;
                    }
 
                    int eqNum = masterCommunicator->giveDofManager(from)->giveCodeNumber(part, k);
                    val += w.at(eqNum);
                }
 
                result &= send_buff.write(val);
            }
        } else {
            masterCommunicator->giveDofManager(from)->giveCompleteLocationArray(rank, locationArray);
            for ( int j = 1; j <= ndofs; j++ ) {
                int eqNum;
                if ( ( eqNum = locationArray.at(j) ) ) {
                    result &= send_buff.write( ( -1.0 ) * w.at(eqNum) );
                }
            }
        }
    }
 
    return result;
}
 
 
int
FETISolver :: unpackSolution(ProcessCommunicator &processComm)
{
    // slaves unpack their slotion contributions
    int result = 1;
    const IntArray &toRecvMap = processComm.giveToRecvMap();
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
    IntArray locationArray;
    EModelDefaultEquationNumbering dn;
 
    for ( int to : toRecvMap ) {
        domain->giveDofManager(to)->giveCompleteLocationArray(locationArray, dn);
        int ndofs = locationArray.giveSize();
        for ( int j = 1; j <= ndofs; j++ ) {
            int eqNum;
            if ( ( eqNum = locationArray.at(j) ) ) {
                double value;
                result &= recv_buff.read(value);
                dd.at(eqNum) = value;
 
#ifdef __VERBOSE_PARALLEL
                OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Unpacking solution value %f at %d\n",
                                this->giveEngngModel()->giveRank(), "FETISolver :: solveYourselfAt", value, eqNum);
#endif
            }
        }
    }
 
    return result;
}
 
int
FETISolver :: masterMapSolution()
{
    int receivedRank = 0, result;
    IntArray locationArray;
    int size = masterCommMap.giveSize();
    EModelDefaultEquationNumbering dn;
 
    for ( int i = 1; i <= size; i++ ) {
        int from = masterCommunicator->giveMasterCommMapPtr()->at(i);
        // use receive map, send map is empty to prevent master to send
        // itself any data. Note, however, that send and receive maps are same.
        int to = masterCommMap.at(i);
 
        domain->giveDofManager(to)->giveCompleteLocationArray(locationArray, dn);
        int locpos = 1;
        //
        // loop over all dofs
        for ( int idof = 1; idof <= masterCommunicator->giveDofManager(from)->giveNumberOfDofs(); idof++, locpos++ ) {
            // extract source value
            while ( locationArray.at(locpos) == 0 ) {
                locpos++;
                // test if position of nonzero dof is allowed
                if ( locpos > locationArray.giveSize() ) {
                    OOFEM_ERROR("Consistency dof error");
                }
            }
 
            double value = 0.0;
            if ( masterCommunicator->giveDofManager(from)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                        // corresponding DOFs in boundary node
 
                // localize to corresponding places
                int nshared = masterCommunicator->giveDofManager(from)->giveNumberOfSharedPartitions();
                for ( int j = 1; j <= nshared; j++ ) {
                    int part = masterCommunicator->giveDofManager(from)->giveSharedPartition(j);
                    if ( part == 0 ) {
                        continue;
                    }
 
                    int eqNum = masterCommunicator->giveDofManager(from)->giveCodeNumber(part, idof);
                    value += w.at(eqNum);
                }
            } else { // no reference partition
                int eqNum = masterCommunicator->giveDofManager(from)->giveCodeNumber(receivedRank, idof);
                value = ( -1.0 ) * w.at(eqNum);
            }
 
            // unpack to locaL CORRESPONDING EQUATION
            dd.at( locationArray.at(locpos) ) = value;
        }
    }
 
    result = 1;
    return result;
}
 
 
 
int
FETISolver :: packResiduals(ProcessCommunicator &processComm)
{
    // slaves
 
    int result = 1;
    const IntArray &toSendMap = processComm.giveToSendMap();
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
    IntArray locationArray;
    EModelDefaultEquationNumbering dn;
 
    for ( int inode : toSendMap ) {
        domain->giveDofManager( inode )->giveCompleteLocationArray(locationArray, dn);
        int ndofs = locationArray.giveSize();
        for ( int j = 1; j <= ndofs; j++ ) {
            int eqNum;
            if ( ( eqNum = locationArray.at(j) ) ) {
                result &= send_buff.write( pp.at(eqNum) );
            }
        }
    }
 
    return result;
}
 
 
int
FETISolver :: unpackResiduals(ProcessCommunicator &processComm)
{
    // master
 
    int result = 1;
    const IntArray &toRecvMap = processComm.giveToRecvMap();
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
    IntArray locationArray;
 
    int receivedRank = processComm.giveRank();
 
    if ( receivedRank != 0 ) {
        for ( int to : toRecvMap ) {
            //
            // loop over all dofs
            for ( int idof = 1; idof <= masterCommunicator->giveDofManager(to)->giveNumberOfDofs(); idof++ ) {
                // unpack contribution
                double value;
                result &= recv_buff.read(value);
                if ( masterCommunicator->giveDofManager(to)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                          // corresponding DOFs in boundary node
 
                    // localize to corresponding places
                    int nshared = masterCommunicator->giveDofManager(to)->giveNumberOfSharedPartitions();
                    for ( int j = 1; j <= nshared; j++ ) {
                        int part = masterCommunicator->giveDofManager(to)->giveSharedPartition(j);
                        if ( part == processComm.giveRank() ) {
                            continue;
                        }
 
                        int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(part, idof);
                        g.at(eqNum) += value;
                    }
                } else { // no reference partition
                    int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(receivedRank, idof);
                    g.at(eqNum) += ( -1.0 ) * value;
                }
            }
        }
    }
 
    return result;
}
 
int
FETISolver :: masterMapResiduals()
{ // master will map its own values directly
    int receivedRank = 0, result;
    IntArray locationArray;
    int size = masterCommMap.giveSize();
    EModelDefaultEquationNumbering dn;
 
    for ( int i = 1; i <= size; i++ ) {
        int to = masterCommunicator->giveMasterCommMapPtr()->at(i);
        // use receive map, send map is empty to prevent master to send
        // itself any data. Note, however, that send and receive maps are same.
        int from = masterCommMap.at(i);
 
        domain->giveDofManager(from)->giveCompleteLocationArray(locationArray, dn);
        int locpos = 1;
        //
        //   ndofs = locationArray.giveSize(); // including supported
        //   for (j=1; j<=ndofs; j++) {
        //   if (eqNum = locationArray.at(j)) {
 
        //
        // loop over all dofs
        for ( int idof = 1; idof <= masterCommunicator->giveDofManager(to)->giveNumberOfDofs(); idof++, locpos++ ) {
            // extract source value
            while ( locationArray.at(locpos) == 0 ) {
                locpos++;
                // test if position of nonzero dof is allowed
                if ( locpos > locationArray.giveSize() ) {
                    OOFEM_ERROR("Consistency dof error");
                }
            }
 
            double value = pp.at( locationArray.at(locpos) );
            if ( masterCommunicator->giveDofManager(to)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                      // corresponding DOFs in boundary node
 
                // localize to corresponding places
                int nshared = masterCommunicator->giveDofManager(to)->giveNumberOfSharedPartitions();
                for ( int j = 1; j <= nshared; j++ ) {
                    int part = masterCommunicator->giveDofManager(to)->giveSharedPartition(j);
                    if ( part == 0 ) {
                        continue;
                    }
 
                    int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(part, idof);
                    g.at(eqNum) += value;
                }
            } else { // no reference partition
                int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(receivedRank, idof);
                g.at(eqNum) += ( -1.0 ) * value;
            }
        }
    }
 
    result = 1;
    return result;
}
 
 
 
 
int
FETISolver :: packDirectionVector(ProcessCommunicator &processComm)
{
    // master
    int result = 1;
    const IntArray &toSendMap = processComm.giveToSendMap();
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
    IntArray locationArray;
 
    int rank = processComm.giveRank();
    for ( int from : toSendMap ) {
        int ndofs = masterCommunicator->giveDofManager(from)->giveNumberOfDofs();
        if ( rank == masterCommunicator->giveDofManager(from)->giveReferencePratition() ) {
            // summ corresponding values (multipliers)
            int nshared = masterCommunicator->giveDofManager(from)->giveNumberOfSharedPartitions();
            for ( int k = 1; k <= ndofs; k++ ) {
                double val = 0.0;
                for ( int j = 1; j <= nshared; j++ ) {
                    int part = masterCommunicator->giveDofManager(from)->giveSharedPartition(j);
                    if ( part == processComm.giveRank() ) {
                        continue;
                    }
 
                    int eqNum = masterCommunicator->giveDofManager(from)->giveCodeNumber(part, k);
                    val += d.at(eqNum);
                }
 
                result &= send_buff.write(val);
            }
        } else {
            masterCommunicator->giveDofManager(from)->giveCompleteLocationArray(rank, locationArray);
            for ( int j = 1; j <= ndofs; j++ ) {
                int eqNum;
                if ( ( eqNum = locationArray.at(j) ) ) {
                    result &= send_buff.write( ( -1.0 ) * d.at(eqNum) );
                }
            }
        }
    }
 
    return result;
}
 
 
int
FETISolver :: unpackDirectionVector(ProcessCommunicator &processComm)
{
    // slaves unpack their slotion contributions
    int result = 1;
    const IntArray &toRecvMap = processComm.giveToRecvMap();
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
    IntArray locationArray;
    // int receivedRank = domainComm.giveRank();
    EModelDefaultEquationNumbering dn;
 
    // if (receivedRank != 0) {
    for ( int inode : toRecvMap ) {
        domain->giveDofManager( inode )->giveCompleteLocationArray(locationArray, dn);
        int ndofs = locationArray.giveSize();
        for ( int j = 1; j <= ndofs; j++ ) {
            int eqNum;
            if ( ( eqNum = locationArray.at(j) ) ) {
                result &= recv_buff.read( dd.at(eqNum) );
            }
        }
    }
 
    // }
    return result;
}
 
int
FETISolver :: masterMapDirectionVector()
{
    int receivedRank = 0, result;
    IntArray locationArray;
    int size = masterCommMap.giveSize();
    EModelDefaultEquationNumbering dn;
 
    for ( int i = 1; i <= size; i++ ) {
        int from = masterCommunicator->giveMasterCommMapPtr()->at(i);
        // use receive map, send map is empty to prevent master to send
        // itself any data. Note, however, that send and receive maps are same.
        int to = masterCommMap.at(i);
 
        domain->giveDofManager(to)->giveCompleteLocationArray(locationArray, dn);
        int locpos = 1;
        //
        // loop over all dofs
        for ( int idof = 1; idof <= masterCommunicator->giveDofManager(from)->giveNumberOfDofs(); idof++, locpos++ ) {
            // extract source value
            while ( locationArray.at(locpos) == 0 ) {
                locpos++;
                // test if position of nonzero dof is allowed
                if ( locpos > locationArray.giveSize() ) {
                    OOFEM_ERROR("Consistency dof error");
                }
            }
 
            double value = 0.0;
            if ( masterCommunicator->giveDofManager(from)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                        // corresponding DOFs in boundary node
 
                // localize to corresponding places
                int nshared = masterCommunicator->giveDofManager(from)->giveNumberOfSharedPartitions();
                for ( int j = 1; j <= nshared; j++ ) {
                    int part = masterCommunicator->giveDofManager(from)->giveSharedPartition(j);
                    if ( part == 0 ) {
                        continue;
                    }
 
                    int eqNum = masterCommunicator->giveDofManager(from)->giveCodeNumber(part, idof);
                    value += d.at(eqNum);
                }
            } else { // no reference partition
                int eqNum = masterCommunicator->giveDofManager(from)->giveCodeNumber(receivedRank, idof);
                value = ( -1.0 ) * d.at(eqNum);
            }
 
            // unpack to locaL CORRESPONDING EQUATION
            dd.at( locationArray.at(locpos) ) = value;
        }
    }
 
    result = 1;
    return result;
}
 
 
int
FETISolver :: packPPVector(ProcessCommunicator &processComm)
{
    // slaves
 
    int result = 1;
    const IntArray &toSendMap = processComm.giveToSendMap();
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
    IntArray locationArray;
    EModelDefaultEquationNumbering dn;
 
    for ( int inode : toSendMap ) {
        domain->giveDofManager( inode )->giveCompleteLocationArray(locationArray, dn);
        int ndofs = locationArray.giveSize();
        for ( int j = 1; j <= ndofs; j++ ) {
            int eqNum;
            if ( ( eqNum = locationArray.at(j) ) ) {
                result &= send_buff.write( pp.at(eqNum) );
            }
        }
    }
 
    return result;
}
 
 
int
FETISolver :: unpackPPVector(ProcessCommunicator &processComm)
{
    // master
 
    int result = 1;
    const IntArray &toRecvMap = processComm.giveToRecvMap();
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
    IntArray locationArray;
 
    int receivedRank = processComm.giveRank();
 
    if ( receivedRank != 0 ) {
        for ( int to : toRecvMap ) {
            //
            // loop over all dofs
            for ( int idof = 1; idof <= masterCommunicator->giveDofManager(to)->giveNumberOfDofs(); idof++ ) {
                // unpack contribution
                double value;
                result &= recv_buff.read(value);
                if ( masterCommunicator->giveDofManager(to)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                          // corresponding DOFs in boundary node
 
                    // localize to corresponding places
                    int nshared = masterCommunicator->giveDofManager(to)->giveNumberOfSharedPartitions();
                    for ( int j = 1; j <= nshared; j++ ) {
                        int part = masterCommunicator->giveDofManager(to)->giveSharedPartition(j);
                        if ( part == processComm.giveRank() ) {
                            continue;
                        }
 
                        int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(part, idof);
                        p.at(eqNum) += value;
                    }
                } else { // no reference partition
                    int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(receivedRank, idof);
                    p.at(eqNum) += ( -1.0 ) * value;
                }
            }
        }
    }
 
    return result;
}
 
int
FETISolver :: masterMapPPVector()
{ // master will map its own values directly
    int receivedRank = 0, result;
    IntArray locationArray;
    int size = masterCommMap.giveSize();
    EModelDefaultEquationNumbering dn;
 
    for ( int i = 1; i <= size; i++ ) {
        int to = masterCommunicator->giveMasterCommMapPtr()->at(i);
        // use receive map, send map is empty to prevent master to send
        // itself any data. Note, however, that send and receive maps are same.
        int from = masterCommMap.at(i);
 
        domain->giveDofManager(from)->giveCompleteLocationArray(locationArray, dn);
        int locpos = 1;
        //
        //   ndofs = locationArray.giveSize(); // including supported
        //   for (j=1; j<=ndofs; j++) {
        //   if (eqNum = locationArray.at(j)) {
 
        //
        // loop over all dofs
        for ( int idof = 1; idof <= masterCommunicator->giveDofManager(to)->giveNumberOfDofs(); idof++, locpos++ ) {
            // extract source value
            while ( locationArray.at(locpos) == 0 ) {
                locpos++;
                // test if position of nonzero dof is allowed
                if ( locpos > locationArray.giveSize() ) {
                    OOFEM_ERROR("Consistency dof error");
                }
            }
 
            double value = pp.at( locationArray.at(locpos) );
            if ( masterCommunicator->giveDofManager(to)->giveReferencePratition() == receivedRank ) { // contribution from reference partition localizes to all
                                                                                                      // corresponding DOFs in boundary node
 
                // localize to corresponding places
                int nshared = masterCommunicator->giveDofManager(to)->giveNumberOfSharedPartitions();
                for ( int j = 1; j <= nshared; j++ ) {
                    int part = masterCommunicator->giveDofManager(to)->giveSharedPartition(j);
                    if ( part == 0 ) {
                        continue;
                    }
 
                    int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(part, idof);
                    p.at(eqNum) += value;
                }
            } else { // no reference partition
                int eqNum = masterCommunicator->giveDofManager(to)->giveCodeNumber(receivedRank, idof);
                p.at(eqNum) += ( -1.0 ) * value;
            }
        }
    }
 
    result = 1;
    return result;
}
 
 
int
FETISolver :: packGammas(ProcessCommunicator &processComm)
{
    // master
 
    int result = 1;
    int rank = processComm.giveRank();
    CommunicationBuffer &send_buff = processComm.giveProcessCommunicatorBuff()->giveSendBuff();
 
    for ( int irbm = 1; irbm <= nsem.at(rank + 1); irbm++ ) {
        result &= send_buff.write( gamma.at(rbmAddr.at(rank + 1) + irbm - 1) );
    }
 
    return result;
}
 
 
int
FETISolver :: unpackGammas(ProcessCommunicator &processComm)
{
    // slaves
    int result = 1;
    CommunicationBuffer &recv_buff = processComm.giveProcessCommunicatorBuff()->giveRecvBuff();
 
    localGammas.resize(nse);
    for ( int irbm = 1; irbm <= nse; irbm++ ) {
        result &= recv_buff.read( localGammas.at(irbm) );
    }
 
    return result;
}
 
int
FETISolver :: masterMapGammas()
{
    // slaves
    int result = 1;
 
    localGammas.resize(nse);
    for ( int irbm = 1; irbm <= nse; irbm++ ) {
        localGammas.at(irbm) = gamma.at(rbmAddr.at(1) + irbm - 1);
    }
 
    return result;
}
 
 
ConvergedReason
FETISolver :: solve(SparseMtrx &A, FloatArray &partitionLoad, FloatArray &partitionSolution)
{
    int tnse = 0, rank = domain->giveEngngModel()->giveRank();
    int source, tag;
    int masterLoopStatus;
    double nom = 0.0, denom, alpha, beta, energyNorm = 0.0;
    FloatMatrix l1;
    StaticCommunicationBuffer commBuff(MPI_COMM_WORLD);
    Skyline *partitionStiffness = dynamic_cast< Skyline * >(&A);
    if ( !partitionStiffness ) {
        OOFEM_ERROR("unsuported sparse matrix type");
    }
 
 
    if ( ( partitionSolution.giveSize() ) != partitionLoad.giveSize() ) {
        OOFEM_ERROR("size mismatch");
    }
 
    int neq = partitionStiffness->giveNumberOfRows();
    int size = domain->giveEngngModel()->giveNumberOfProcesses();
    nsem.resize(size);
 
    // initialize communication maps
    this->setUpCommunicationMaps();
    if ( rank == 0 ) {
        this->commBuff = new CommunicatorBuff(size);
        masterCommunicator = new FETICommunicator(engngModel, this->commBuff, engngModel->giveRank(), size);
        masterCommunicator->setUpCommunicationMaps(engngModel);
    }
 
    MPI_Barrier(MPI_COMM_WORLD);
 
#ifdef __VERBOSE_PARALLEL
    OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Barrier reached\n",
                    this->giveEngngModel()->giveRank(), "FETISolver :: solveYourselfAt");
#endif
 
    /*********************************/
    /*  rozklad matice               */
    /*  vypocet rigid body modes     */
    /*********************************/
    /*  pole indexu zavislych rovnic  */
    se.resize(6);
    se.zero();
    /*  pole posunu jako tuheho telesa  */
    // rbm.resize (neq, 6); rbm.zero();
 
#ifdef __VERBOSE_PARALLEL
    OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: rbmodes computation startup, lneq is %d\n",
                    rank, "FETISolver :: solveYourselfAt", neq);
#endif
 
    partitionStiffness->rbmodes(rbm, nse, se, limit, 3);
 
#ifdef __VERBOSE_PARALLEL
    OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Number of Rigid body modes %3d\n",
                    rank, "FETISolver :: solveYourselfAt", nse);
 
    //fprintf(stdout, "\n[process rank %3d]: %-30s: RBM dump",
    //    rank,"FETISolver :: solveYourselfAt");
    //rbm.printYourself();
 
#endif
 
    /*****************************/
    /*  vypocet soucinu R^T . f  */
    /*****************************/
 
    if ( nse != 0 ) {
        qq.resize(nse);
        qq.zero();
 
        qq.beTProductOf(rbm, partitionLoad);
        qq.negated();
    }
 
    /***************************************************************************/
    /*  zjisteni rozmeru matice L, vektoru q, poctu neznamych na podoblastech  */
    /***************************************************************************/
 
    //  commBuff.init();
    if ( rank == 0 ) {
        commBuff.resize( commBuff.givePackSizeOfInt(1) );
        //
        // receive data
        //
        for ( int i = 1; i < size; i++ ) {
            commBuff.iRecv(MPI_ANY_SOURCE, FETISolver :: NumberOfRBMMsg);
            while ( !commBuff.testCompletion(source, tag) ) {
                ;
            }
 
            // unpack data
#ifdef __VERBOSE_PARALLEL
            OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Received data from partition %3d\n",
                            rank, "FETICommunicator :: setUpCommunicationMaps : received number of partition rbm", source);
#endif
            commBuff.init();
            commBuff.read( nsem.at(source + 1) );
            tnse += nsem.at(source + 1);
        }
 
        nsem.at(1) = nse;
        tnse += nse;
 
        OOFEM_LOG_INFO("Number of RBM per partion\npart. rbm\n-------------------------------\n");
        for ( int i = 1; i <= size; i++ ) {
            OOFEM_LOG_INFO( "%-4d %8d\n", i - 1, nsem.at(i) );
        }
 
 
        // initialize partition start indexes into rbm array
        rbmAddr.resize(size);
        rbmAddr.zero();
        rbmAddr.at(1) = 1;
        for ( int i = 2; i <= size; i++ ) {
            rbmAddr.at(i) = rbmAddr.at(i - 1) + nsem.at(i - 1);
        }
    } else { // slave code
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Sending number of Rigid body modes %3d\n",
                        rank, "FETISolver :: solveYourselfAt", nse);
#endif
 
        commBuff.resize( commBuff.givePackSizeOfInt(1) );
        commBuff.write(nse);
        commBuff.iSend(0, FETISolver :: NumberOfRBMMsg);
    }
 
    MPI_Barrier(MPI_COMM_WORLD);
 
    /*****************************/
    /*  sestaveni cele matice L  */
    /*****************************/
 
    if ( rank == 0 ) {
        // assemble l matrix containing rbm localized from partition contibutions
        if ( tnse ) {
            l.resize(masterCommunicator->giveNumberOfDomainEquations(), tnse);
            l.zero();
        }
 
        masterCommunicator->initReceive(FETISolver :: RBMMessage);
        masterCommunicator->unpackAllData(this, & FETISolver :: masterUnpackRBM);
        this->masterMapRBM();
 
        if ( tnse ) {
            l.negated();
        }
    } else {
        // pack RBMs
        processCommunicator.packData(this, & FETISolver :: packRBM);
        processCommunicator.initSend(FETISolver :: RBMMessage);
        MPI_Barrier(MPI_COMM_WORLD);
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: RBMMessage Barrier finished",
                        rank, "FETISolver :: solveYourselfAt");
#endif
    }
 
    if ( rank == 0 ) {
        // receive products of rbm * LoadVector
        if ( tnse ) {
            q.resize(tnse);
            q.zero();
        }
 
        masterCommunicator->initReceive(FETISolver :: QQMessage);
        masterCommunicator->unpackAllData(this, & FETISolver :: masterUnpackQQProduct);
        this->masterMapQQProduct();
    } else {
        // pack QQ products
        processCommunicator.packData(this, & FETISolver :: packQQProducts);
        processCommunicator.initSend(FETISolver :: QQMessage);
        MPI_Barrier(MPI_COMM_WORLD);
    }
 
    /*******************************************/
    /*******************************************/
    /*  sestaveni matice L^T.L a jeji inverze  */
    /*******************************************/
    /*******************************************/
 
    if ( ( rank == 0 ) && ( tnse != 0 ) ) {
        FloatMatrix l2;
        // compute l1, the inverse of (l^Tl)
        l2.beTProductOf(l, l);
        l1.beInverseOf(l2);
    }
 
    /****************************************************************************/
    /****************************************************************************/
    /*  vypocet pocatecni aproximace neznamych do modifikovane met. sdr. grad.  */
    /****************************************************************************/
    /****************************************************************************/
 
    /****************************/
    /*  alokace nekterych poli  */
    /****************************/
    /*  vektor smeru na jedne podoblasti  */
    dd.resize(neq);
    /*  pomocny vektor na jedne podoblasti  */
    pp.resize(neq);
 
 
    if ( rank == 0 ) {
        /*  vektor smeru  */
        d.resize(masterCommunicator->giveNumberOfDomainEquations() + 1);
        /*  pomocny vektor  */
        p.resize( masterCommunicator->giveNumberOfDomainEquations() );
        /*  vektor gradientu  */
        g.resize( masterCommunicator->giveNumberOfDomainEquations() );
        /*  vektor neznamych  */
        w.resize( masterCommunicator->giveNumberOfDomainEquations() );
 
        /*  soucin (L^T.L)^{-1}.q  */
        if ( tnse ) {
            FloatArray help(tnse);
            help.beProductOf(l1, q);
            w.beProductOf(l, help);
        } else {
            w.zero();
        }
 
        // send aproximation of solution to slaves
        masterCommunicator->packAllData(this, & FETISolver :: packSolution);
        masterCommunicator->initSend(FETISolver :: SolutionMessage);
        this->masterMapSolution();
    } else {
        // receive send aproximation of solution
        processCommunicator.initReceive(FETISolver :: SolutionMessage);
        while ( !processCommunicator.receiveCompleted() ) {
            ;
        }
 
        processCommunicator.unpackData(this, & FETISolver :: unpackSolution);
    }
 
    MPI_Barrier(MPI_COMM_WORLD);
 
#ifdef __VERBOSE_PARALLEL
    OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Solution Approx Barrier finished\n",
                    rank, "FETISolver :: solveYourselfAt");
#endif
 
    dd.subtract(partitionLoad);
    partitionStiffness->ldl_feti_sky(pp, dd, nse, limit, se);
 
    if ( rank == 0 ) {
        // receive contributions (to resudals) from slaves
 
        masterCommunicator->initReceive(FETISolver :: ResidualMessage);
        masterCommunicator->unpackAllData(this, & FETISolver :: unpackResiduals);
        this->masterMapResiduals();
    } else {
        // send contributions (to resudals) to master
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Residual contribution packing initiated\n",
                        rank, "FETISolver :: solveYourselfAt");
#endif
        processCommunicator.packData(this, & FETISolver :: packResiduals);
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Residual contribution send initiated\n",
                        rank, "FETISolver :: solveYourselfAt");
#endif
        processCommunicator.initSend(FETISolver :: ResidualMessage);
 
        MPI_Barrier(MPI_COMM_WORLD);
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Residual contribution Barrier finished\n",
                        rank, "FETISolver :: solveYourselfAt");
#endif
    }
 
 
    /*  v tuto chvili je vypocten nulty vektor gradientu  */
 
    /***********************************/
    /*  vypocet nulteho vektoru smeru  */
    /***********************************/
    if ( rank == 0 ) {
        if ( tnse ) {
            this->projection(g, l, l1);
        }
 
        d = g;
        d.negated();
 
        /*  vypocet nulteho citatele  */
        nom = g.computeSquaredNorm();
 
        //#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("\nIteration process\n");
        if ( energyNorm_comput_flag ) {
            OOFEM_LOG_DEBUG("iteration      gradient vector norm      energy norm\n=====================================================================\n");
        } else {
            OOFEM_LOG_DEBUG("iteration      gradient vector norm\n================================================\n");
        }
 
        //#endif
    }
 
    /*  cyklus vyjadrujici pocet iteraci  */
    for ( int i = 0; i < ni; i++ ) {
        dd.zero();
        if ( rank == 0 ) {
            /***************************************/
            /*  vypocet pomocneho vektoru K.d = p  */
            /***************************************/
 
            // send direction vector to slaves
            masterCommunicator->packAllData(this, & FETISolver :: packDirectionVector);
            masterCommunicator->initSend(FETISolver :: DirectionVectorMessage);
            this->masterMapDirectionVector();
        } else {
            // receive direction vector on slaves
            processCommunicator.initReceive(FETISolver :: DirectionVectorMessage);
            while ( !processCommunicator.receiveCompleted() ) {
                ;
            }
 
            processCommunicator.unpackData(this, & FETISolver :: unpackDirectionVector);
        }
 
        MPI_Barrier(MPI_COMM_WORLD);
 
 
        pp.zero();
        partitionStiffness->ldl_feti_sky(pp, dd, nse, limit, se);
 
        if ( rank == 0 ) {
            p.zero();
 
            // receive contributions (to resudals) from slaves
 
            masterCommunicator->initReceive(FETISolver :: PPVectorMessage);
            masterCommunicator->unpackAllData(this, & FETISolver :: unpackPPVector);
            this->masterMapPPVector();
        } else {
            // pack pp on all partitions
            processCommunicator.packData(this, & FETISolver :: packPPVector);
            processCommunicator.initSend(FETISolver :: PPVectorMessage);
            MPI_Barrier(MPI_COMM_WORLD);
        }
 
        if ( rank == 0 ) {
            /*  calculation of the denominator of the fraction defining Alpha */
            denom = d.dotProduct(p);
 
            if ( fabs(denom) < FETISOLVER_ZERONUM ) {
                OOFEM_LOG_RELEVANT("FETISolver::solve :  v modifikovane metode sdruzenych gradientu je nulovy jmenovatel u soucinitele alpha\n");
                commBuff.resize( commBuff.givePackSizeOfInt(1) );
                commBuff.init();
                commBuff.write(FETISolverIterationBreak);
                commBuff.bcast(0);
                break;
            }
 
            /*******************/
            /*  vypocet alpha  */
            /*******************/
            alpha = nom / denom;
 
            /**************************************************************/
            /*  vypocet noveho gradientu g a nove aproximace neznamych x  */
            /**************************************************************/
            for ( int j = 1; j <= masterCommunicator->giveNumberOfDomainEquations(); j++ ) {
                w.at(j) += alpha * d.at(j);
                g.at(j) += alpha * p.at(j);
            }
 
 
 
            if ( tnse ) {
                this->projection(g, l, l1);
            }
 
            denom = nom;
            if ( fabs(denom) < FETISOLVER_ZERONUM ) {
                OOFEM_LOG_RELEVANT("FETISolver::solve : v modifikovane metode sdruzenych gradientu je nulovy jmenovatel u soucinitele beta\n");
                commBuff.resize( commBuff.givePackSizeOfInt(1) );
                commBuff.init();
                commBuff.write(FETISolverIterationBreak);
                commBuff.bcast(0);
                break;
            }
 
            nom = g.computeSquaredNorm();
 
            if ( nom < err ) {
                commBuff.resize( commBuff.givePackSizeOfInt(1) );
                commBuff.init();
                commBuff.write(FETISolverIterationBreak);
                commBuff.bcast(0);
                break;
            }
 
 
            /******************/
            /*  vypocet beta  */
            /******************/
            beta = nom / denom;
 
            /****************************/
            /*  vypocet noveho smeru d  */
            /****************************/
            for ( int j = 1; j <= masterCommunicator->giveNumberOfDomainEquations(); j++ ) {
                d.at(j) = beta * d.at(j) - g.at(j);
            }
        }
 
 
        if ( rank == 0 ) {
            commBuff.resize( commBuff.givePackSizeOfInt(1) );
            commBuff.init();
            commBuff.write(FETISolverIterationContinue);
            commBuff.bcast(0);
        } else {
            commBuff.resize( commBuff.givePackSizeOfInt(1) );
            commBuff.init();
            commBuff.bcast(0);
            commBuff.read(masterLoopStatus);
            if ( masterLoopStatus == FETISolverIterationBreak ) {
#ifdef __VERBOSE_PARALLEL
                OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Received Loop break signal from master\n",
                                rank, "FETISolver :: solveYourselfAt");
#endif
 
                break;
            }
        }
 
 
        // computation of energy norm
        if ( energyNorm_comput_flag ) {
            dd.zero();
            if ( rank == 0 ) {
                // send aproximation of solution to slaves
                masterCommunicator->packAllData(this, & FETISolver :: packSolution);
                masterCommunicator->initSend(FETISolver :: SolutionMessage);
                this->masterMapSolution();
            } else {
                // receive send aproximation of solution
                processCommunicator.initReceive(FETISolver :: SolutionMessage);
                while ( !processCommunicator.receiveCompleted() ) {
                    ;
                }
 
                processCommunicator.unpackData(this, & FETISolver :: unpackSolution);
            }
 
            MPI_Barrier(MPI_COMM_WORLD);
 
            pp.zero();
            partitionStiffness->ldl_feti_sky(pp, dd, nse, limit, se);
 
            if ( rank == 0 ) {
                p.zero();
 
                // receive contributions (to resudals) from slaves
 
                masterCommunicator->initReceive(FETISolver :: PPVectorMessage);
                masterCommunicator->unpackAllData(this, & FETISolver :: unpackPPVector);
                this->masterMapPPVector();
            } else {
                // pack pp on all partitions
                processCommunicator.packData(this, & FETISolver :: packPPVector);
                processCommunicator.initSend(FETISolver :: PPVectorMessage);
                MPI_Barrier(MPI_COMM_WORLD);
            }
 
            if ( rank == 0 ) {
                energyNorm = p.dotProduct(w);
            }
        }
 
        // end of energy norm computation
 
        //#ifdef __VERBOSE_PARALLEL
        if ( rank == 0 ) {
            if ( energyNorm_comput_flag ) {
                OOFEM_LOG_DEBUG("%-9d%15e %15e\n", i, nom, energyNorm);
            } else {
                OOFEM_LOG_DEBUG("%-9d%15e\n", i, nom);
            }
        }
#ifdef __VERBOSE_PARALLEL        
        if ( rank == 0 ) {
          OOFEM_LOG_DEBUG("Konec metody sdruzenych gradientu, nite %d, err %e\n", i, nom);
        }
#endif
 
    } // end iterative loop
 
 
    /*  prelom (break) */
 
    /****************************************************/
    /*  konec modifikovane metody sdruzenych gradientu  */
    /****************************************************/
 
    if ( rank == 0 ) {
        OOFEM_LOG_INFO("End of iteration, reached norm %15e\n", nom);
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("\nVysledne Lagrangeovy multiplikatory\n");
        for ( int i = 1; i <= masterCommunicator->giveNumberOfDomainEquations(); i++ ) {
            OOFEM_LOG_DEBUG( "lambda %4d          %f\n", i, w.at(i) );
        }
 
#endif
    }
 
    dd.zero();
    if ( rank == 0 ) {
        // send solution to slaves
        masterCommunicator->packAllData(this, & FETISolver :: packSolution);
        masterCommunicator->initSend(FETISolver :: SolutionMessage);
        this->masterMapSolution();
    } else {
        // receive solution
        processCommunicator.initReceive(FETISolver :: SolutionMessage);
        while ( !processCommunicator.receiveCompleted() ) {
            ;
        }
 
        processCommunicator.unpackData(this, & FETISolver :: unpackSolution);
    }
 
    MPI_Barrier(MPI_COMM_WORLD);
 
 
#ifdef __VERBOSE_PARALLEL
    //fprintf(stdout, "\n[process rank %3d]: %-30s: Partition lagrange multipliers dump:\n",
    //    rank,"FETISolver :: solveYourselfAt");
    //dd.printYourself();
#endif
 
    partitionLoad.subtract(dd);
#ifdef __VERBOSE_PARALLEL
    // fprintf(stdout, "\n[process rank %3d]: %-30s: Partition load dump:\n",
    //    rank,"FETISolver :: solveYourselfAt");
    //partitionLoad.printYourself();
#endif
    partitionStiffness->ldl_feti_sky(partitionSolution, partitionLoad, nse, limit, se);
    pp = partitionSolution;
 
#ifdef __VERBOSE_PARALLEL
    //fprintf(stdout, "\n[process rank %3d]: %-30s: Partition solution dump:\n",
    //    rank,"FETISolver :: solveYourselfAt");
    //pp.printYourself();
#endif
 
    if ( rank == 0 ) {
        g.zero();
        // receive contributions (to resudals) from slaves
 
        masterCommunicator->initReceive(FETISolver :: ResidualMessage);
        masterCommunicator->unpackAllData(this, & FETISolver :: unpackResiduals);
 
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Received residuals from slaves\n",
                        rank, "FETISolver :: solveYourselfAt");
#endif
 
        this->masterMapResiduals();
    } else {
        // send contributions (to resudals) to master
 
#ifdef __VERBOSE_PARALLEL
        OOFEM_LOG_DEBUG("[process rank %3d]: %-30s: Sending residuals to master\n",
                        rank, "FETISolver :: solveYourselfAt");
#endif
 
        processCommunicator.packData(this, & FETISolver :: packResiduals);
        processCommunicator.initSend(FETISolver :: ResidualMessage);
        MPI_Barrier(MPI_COMM_WORLD);
    }
 
#ifdef __VERBOSE_PARALLEL
    if ( rank == 0 ) {
        //fprintf(stdout, "\n[process rank %3d]: %-30s: Partition residuals dump:\n",
        //    rank,"FETISolver :: solveYourselfAt");
        //g.printYourself();
    }
 
#endif
 
    if ( rank == 0 ) {
        if ( tnse ) {
            FloatArray help1;
            help1.beTProductOf(l, g);
            // coefficient of linear combinations
            gamma.beProductOf(l1, help1);
        } else {
            gamma.zero();
        }
 
        // send corresponding gammas to slaves
        masterCommunicator->packAllData(this, & FETISolver :: packGammas);
        masterCommunicator->initSend(FETISolver :: GammasMessage);
        this->masterMapGammas();
    } else {
        // receive gammas
        processCommunicator.initReceive(FETISolver :: GammasMessage);
        while ( !processCommunicator.receiveCompleted() ) {
            ;
        }
 
        processCommunicator.unpackData(this, & FETISolver :: unpackGammas);
    }
 
    MPI_Barrier(MPI_COMM_WORLD);
 
    if ( nse != 0 ) {
        FloatArray help;
        help.beProductOf(rbm, localGammas);
 
        partitionSolution.add(help);
    }
 
    return CR_CONVERGED;
}
} // end namespace oofem