bddcml_wrapper.cc

Go to the documentation of this file.

//! @author Jakub Sistek
//! @date   7/2012

#include "la/bddcml_wrapper.hh"

//------------------------------------------------------------------------------
/** Constructor for the BDDCML solver object
 */
la::BddcmlWrapper::BddcmlWrapper( const unsigned numDofs,
                                  const unsigned numDofsSub,
                                  const MatrixType matrixType, 
                                  const MPI_Comm comm,
                                  int numSubLoc )
  : numDofs_( static_cast<int>( numDofs ) ),
    numDofsSub_( static_cast<int>( numDofsSub ) ),
    matrixType_( matrixType ),
    comm_ (comm),
    numSubLoc_ ( numSubLoc )
{

    // Orient in the communicator
    MPI_Comm_size( comm_, &nProc_ );
    MPI_Comm_rank( comm_, &rank_ );

    // prepare RHS vector
    rhsVec_.resize( numDofsSub_, 0. );

    // prepare arrays to mark boundary conditions
    ifix_.resize( numDofsSub_, 0 );
    fixv_.resize( numDofsSub_, 0. );

    // initialize subdomain solution vector
    sol_.resize( numDofsSub_, 0. );
    normSol_ = 0.;
    numIter_ = 0;
    convergedReason_ = 0;
    condNumber_ = 0.;

    // initialize state of mesh and matrix
    meshLoaded_     = false;
    isMatAssembled_ = false;
}
//------------------------------------------------------------------------------
/** Destructor frees the BDDC solver structures
 */
la::BddcmlWrapper::~BddcmlWrapper()
{
    // clear RHS
    rhsVec_.clear();

    // clear BC
    ifix_.clear();
    fixv_.clear();

    // clear solution
    sol_.clear();

    // clear mesh data
    inet_.clear();            
    nnet_.clear();            
    nndf_.clear();            
    xyz_.clear();             
    isngn_.clear();         
    isvgvn_.clear();         
    isegn_.clear();           

    // clear global2local map
    global2LocalDofMap_.clear();

    // clear triplet
    coo_.clear();

}

//------------------------------------------------------------------------------
/** Routine for loading raw mesh data into the solver - for cases of strange meshes, 
 * where these are not Mesh or Grid objects, user can create own raw description
 * to exploit the flexibility of mesh format underlaying BDDCML.
 */
void la::BddcmlWrapper::loadRawMesh( const int nDim, const int numNodes,
                                     const std::vector<int> & inet, 
                                     const std::vector<int> & nnet, 
                                     const std::vector<int> & nndf, 
                                     const std::vector<int> & isegn, 
                                     const std::vector<int> & isngn, 
                                     const std::vector<int> & isvgvn,
                                     const std::vector<double> & xyz,
                                     const std::vector<double> & element_data,
                                     const int meshDim )
{
    // space dimension
    nDim_ = nDim;

    // set topological dimension of mesh
    if ( meshDim == 0 ) { // if not specified, use space dimension
       meshDim_ = nDim_;
    }
    else {
       meshDim_ = static_cast<int>( meshDim ); // if specified, use it
    }

    // fill internal variables about mesh
    numElemSub_  = isegn.size( );
    numNodesSub_ = isngn.size( );

    // global number of nodes
    numNodes_ = numNodes;

    std::vector<int> procElemStarts( nProc_ + 1 );
    int numElemSubInt = static_cast<int>( numElemSub_ );
    MPI_Allgather( &numElemSubInt,       1, MPI_INT,
                   &(procElemStarts[0]), 1, MPI_INT, comm_ );
    // now procElemStarts_ contains counts of local elements at each process
    // change it to starts
    for ( unsigned i = 1; i < static_cast<unsigned>( nProc_ ); i++ )
        procElemStarts[i] = procElemStarts[i-1] + procElemStarts[i];
    // shift it one back
    for ( unsigned i = nProc_; i > 0; i-- )
            procElemStarts[i] = procElemStarts[i-1];
    // start from zero
    procElemStarts[0] = 0;

    // find global number of elements - assumes nonoverlapping division
    numElem_ = procElemStarts[ nProc_ ];

    // check sizes of arrays
    OLD_ASSERT( std::accumulate( nnet.begin(), nnet.end(), 0 ) == (int)(inet.size()),
                            "array inet size mismatch \n " );
    OLD_ASSERT( std::accumulate( nndf.begin(), nndf.end(), 0 ) == numDofsSub_,
                            "array nndf content mismatch: %d %d \n ",  std::accumulate( nndf.begin(), nndf.end(), 0 ), numDofsSub_ );
    OLD_ASSERT( static_cast<int>(isvgvn.size()) == numDofsSub_,
                            "array isvgvn size mismatch \n " );
    OLD_ASSERT( static_cast<int>(nnet.size()) == numElemSub_,
                            "array nnet size mismatch \n " );
    OLD_ASSERT( static_cast<int>(nndf.size()) == numNodesSub_,
                            "array nndf size mismatch \n " );
    OLD_ASSERT( static_cast<int>(xyz.size()) == numNodesSub_ * nDim_,
                            "array xyz size mismatch: %d %d \n ", static_cast<int>(xyz.size()), numNodesSub_ * nDim_ );

    // Simply copy input data into the private object equivalents
    inet_.resize( inet.size() );
    std::copy( inet.begin(), inet.end(), inet_.begin() );
    nnet_.resize( nnet.size() );
    std::copy( nnet.begin(), nnet.end(), nnet_.begin() );
    nndf_.resize( nndf.size() );
    std::copy( nndf.begin(), nndf.end(), nndf_.begin() );
    isegn_.resize( isegn.size() );
    std::copy( isegn.begin(), isegn.end(), isegn_.begin() );
    isngn_.resize( isngn.size() );
    std::copy( isngn.begin(), isngn.end(), isngn_.begin() );
    isvgvn_.resize( isvgvn.size() );
    std::copy( isvgvn.begin(), isvgvn.end(), isvgvn_.begin() );
    xyz_.resize( xyz.size() );
    std::copy( xyz.begin(), xyz.end(), xyz_.begin() );
    element_data_.resize( element_data.size() );
    std::copy( element_data.begin(), element_data.end(), element_data_.begin() );

    // prepare auxiliary map for renumbering dofs in the global array
    for ( unsigned ind = 0; ind < isvgvn_.size(); ++ind ) {
        global2LocalDofMap_.insert( std::make_pair( static_cast<unsigned>( isvgvn_[ind] ), ind ) );
    }

    // change the state
    meshLoaded_     = true;

    return;
}

//------------------------------------------------------------------------------
/** Routine for preparing MATRIX assembly. It reserves memory in triplet.
 */
void la::BddcmlWrapper::prepareMatAssembly( unsigned numElements, unsigned elMatSize )
{
    unsigned length = numElements * elMatSize;
    coo_.prepareAssembly( length );

    return;
}
//------------------------------------------------------------------------------
/** Insert element stiffness matrix into global system matrix
 */
void la::BddcmlWrapper::insertToMatrix( const SubMat_  & subMat, 
                                        const VecUint_ & rowIndices,
                                        const VecUint_ & colIndices )
{
    // copy dense matrix into triplet
    for ( unsigned i = 0; i < rowIndices.size(); i++) {
        for ( unsigned j = 0; j < colIndices.size(); j++) {

            // insert value
            if (subMat(i,j) != 0.0)
              coo_.insert( rowIndices[i], colIndices[j], subMat(i,j) );
        }
    }

    // set flag
    isMatAssembled_ = false;

    return;
}

//------------------------------------------------------------------------------
/** Routine for triplet finalization.
 */
void la::BddcmlWrapper::finishMatAssembly( )
{
    // do nothing if already called
    if ( isMatAssembled_ ) return;

    // finish assembly of triplet
    coo_.finishAssembly( );

    // set flag
    isMatAssembled_ = true;

    return;
}
//------------------------------------------------------------------------------
/** Insert subvector to the system's RHS vector
 */
void la::BddcmlWrapper::insertToRhs( const SubVec_  & subVec, 
                                     const VecUint_ & dofIndices )
{
    for ( unsigned i = 0; i < dofIndices.size( ); i ++ ) {

        // map global dof index to subdomain dof
        Global2LocalMap_::iterator pos = global2LocalDofMap_.find( dofIndices[i] );
        /*
        OLD_ASSERT( pos != global2LocalDofMap_.end(),
                                "Cannot remap index %d to local indices in right-hand side," 
                                " perhaps missing call to loadMesh() member function. \n ", dofIndices[i] );
        */
        if (pos == global2LocalDofMap_.end()) {
            DebugOut().every_proc() << "RHS, Missing map to local for idx: " << dofIndices[i];
        }
        const unsigned indLoc = pos -> second;

        rhsVec_[ indLoc ] += subVec[i];
    }
    return;
}

//------------------------------------------------------------------------------
/** Apply boundary conditions
 *  Since BDDCML handles constraints, just fill proper BC arrays.
 *  Scalar is not used, it is kept for compatibility with constraint functors.
 */
void la::BddcmlWrapper::applyConstraints( ConstraintVec_ & constraints,
                                          const double factor, FMT_UNUSED const double scalar ) 
{
    // Constraint iterators
    ConstraintVec_::const_iterator cIter = constraints.begin( );
    ConstraintVec_::const_iterator cEnd  = constraints.end( );
    // collect global dof indices and the correpsonding values
    for ( ; cIter != cEnd; ++cIter ) {
        unsigned globalDof =  cIter -> first;
        double   value    = (cIter -> second );

        // map global dof index to subdomain dof
        Global2LocalMap_::iterator pos = global2LocalDofMap_.find( globalDof );
        if ( pos != global2LocalDofMap_.end() ) {
            // I contain the dof
            unsigned indLoc = pos -> second;

            ifix_[ indLoc ] = 1;
            fixv_[ indLoc ] = factor * value;
        }
    }

    return;
}
//------------------------------------------------------------------------------
/** Simply set the dof with index 'index' to zero
 */
void la::BddcmlWrapper::fixDOF( const unsigned index, FMT_UNUSED const double scalar )
{
    std::vector< std::pair<unsigned,double> > thisConstraint;
    thisConstraint.push_back( std::make_pair( index, 0. ) );
    this -> applyConstraints( thisConstraint, 1., scalar );
    thisConstraint.clear();
    return;
}

//------------------------------------------------------------------------------
/** Insert elements on diagonal for weights.
 */
void la::BddcmlWrapper::insertToDiagonalWeights( const int & index,
                                                 const double & value )
{
    // insert value
    diagWeightsCoo_.insert( index, index, value );

    // set flag
    isDiagAssembled_ = false;

    return;
}

//------------------------------------------------------------------------------
/** Routine for diagonal finalization.
 */
void la::BddcmlWrapper::finishDiagAssembly( )
{
    // do nothing if already called
    if ( isDiagAssembled_ ) return;

    // finish assembly of triplet
    diagWeightsCoo_.finishAssembly( );

    // set flag
    isDiagAssembled_ = true;

    return;
}

//------------------------------------------------------------------------------
/** Solve the system by BDDCML
 */
void la::BddcmlWrapper::solveSystem( double tol, int  numLevels, std::vector<int> *  numSubAtLevels, 
                                     int verboseLevel, int  maxIt, int  ndecrMax, bool use_adaptive )
{

    // check that mesh was loaded 
    OLD_ASSERT( meshLoaded_ ,
                            "Subdomain mesh was not loaded, perhaps missing call to loadMesh() member function. \n " );
   
    //============= initialize BDDC
    int numSub;                  //!< number of subdomains on first level   
    // Number of subdomains at first level - if not given, deduce it from MPI communicator
    MPI_Allreduce ( &numSubLoc_, &numSub, 1, MPI_INT, MPI_SUM, comm_ );

    // Number of subdomains at levels
    std::vector<int> numSubLev( numLevels ); //!< numbers of subdomains on levels ( should be monotonically decreasing, e.g. for 3 levels 1024, 32, 1 )
    if ( numSubAtLevels != NULL ) {
        // numbers are given 
        
        OLD_ASSERT( numLevels == static_cast<int>((*numSubAtLevels).size()),
                                "Inconsistent size of numbers of subdomains at levels: %d %d \n",
                                numLevels, static_cast<int>((*numSubAtLevels).size()) );
        OLD_ASSERT( (*numSubAtLevels)[0] == numSub, 
                                "Inconsistent number of subdomains at first level: %d %d \n",
                                (*numSubAtLevels)[0], numSub );

        std::copy( (*numSubAtLevels).begin(), (*numSubAtLevels).end(),
                    numSubLev.begin() );
    }
    else {
        // assume 2-level BDDC
        // generate default vector for two levels
        if ( numLevels == 2 ) {
            numSubLev[0] = numSub;
            numSubLev[1] = 1;
        }
        else if ( numLevels > 2 ) {
            double coarsening = pow( numSub, 1. / double(numLevels - 1) );
            // prescribe number of subdomains on levels so that coarsening is fixed between levels
            numSubLev[0] = numSub;
            for (int i = 1; i <= numLevels-2; i++) {
                int ir = numLevels - i - 1;
                numSubLev[i] = (int) pow( coarsening, ir );
                if ( numSubLev[i]%2 != 0 ) {
                    numSubLev[i] = numSubLev[i] + 1;
                }
            }
            numSubLev[numLevels-1] = 1;
        }
        else {
            OLD_ASSERT( false, "Missing numbers of subdomains at levels. \n" );
        }
    }
    //std::cout << "numSubLev ";
    //std::copy( numSubLev.begin(), numSubLev.end(), std::ostream_iterator<int>( std::cout, " " ) );
    //std::cout << std::endl;

    int lnumSubLev   = numSubLev.size();
    int commInt      = MPI_Comm_c2f( comm_ ); // translate C MPI communicator hanlder to Fortran integer
    int numBase      = 0;                     // indexing starts with 0 for C++
    int justDirectSolve = 0;

    bddcml_init_c( &numLevels, &(numSubLev[0]), &lnumSubLev, &numSubLoc_, &commInt, &verboseLevel, &numBase, &justDirectSolve );


    //============= uploading data for subdomain
    // at the moment, number of subdomains equals number of processors
    int iSub = rank_;

    int linet   = inet_.size();
    int lnnet   = nnet_.size();
    int lnndf   = nndf_.size();
    int lisngn  = isngn_.size();
    int lisvgvn = isvgvn_.size();
    int lisegn  = isegn_.size();

    int lxyz1   = numNodesSub_;
    int lxyz2   = nDim_;

    int lifix   = ifix_.size();
    int lfixv   = fixv_.size();

    int lrhsVec = rhsVec_.size();
    int isRhsCompleteInt = 0; // since only local subassembly of right-hand side was performed, it does not contain 
                              // repeated entries

    std::vector<int>    i_sparse;
    std::vector<int>    j_sparse;
    std::vector<double> a_sparse;

    // specify if matrix is symmetric, if yes, use only upper triangle
    bool onlyUpperTriangle;
    if ( matrixType_ == GENERAL ) {
        onlyUpperTriangle = false;
    }
    else {
        onlyUpperTriangle = true;
    }

    // copy out arrays from triplet
    this -> finishMatAssembly();
    coo_.extractArrays( i_sparse, j_sparse, a_sparse, onlyUpperTriangle );
    coo_.clear();

    // map global row and column indices to local subdomain indices
    int indRow    = -1;
    int indCol    = -1;
    int indRowLoc = -1;
    int indColLoc = -1;
    for ( unsigned inz = 0; inz < a_sparse.size(); inz++ ) {
        if ( i_sparse[inz] != indRow ) {
           indRow = i_sparse[inz];
           Global2LocalMap_::iterator pos = global2LocalDofMap_.find( static_cast<unsigned> ( indRow ) );
           OLD_ASSERT( pos != global2LocalDofMap_.end(),
                                   "Cannot remap index %d to local indices. \n ", indRow );
           indRowLoc = static_cast<int> ( pos -> second );
        }
        if ( j_sparse[inz] != indCol ) {
           indCol = j_sparse[inz];
           Global2LocalMap_::iterator pos = global2LocalDofMap_.find( static_cast<unsigned> ( indCol ) );
           OLD_ASSERT( pos != global2LocalDofMap_.end(),
                                   "Cannot remap index %d to local indices. \n ", indCol );
           indColLoc = static_cast<int> ( pos -> second );
        }

        i_sparse[ inz ] = indRowLoc;
        j_sparse[ inz ] = indColLoc;

        //std::cout << " row: " << i_sparse[ inz ] << " col: " << j_sparse[ inz ] << std::endl;
    }
    // due to remap, it is not guaranteed that the triplet is still correctly ordered
    int isMatAssembledInt = 0;

    int la = a_sparse.size();

    // diagonal weights for BDDC loaded by user
    OLD_ASSERT( diagWeightsCoo_.nnz() > 0,
            "It appears that diagonal weights for BDDC are not loaded. This is currently mandatory. \n " );

    std::vector<int>    i_diag_sparse;
    std::vector<int>    j_diag_sparse;
    std::vector<double> diag_sparse;

    // copy out arrays from triplet
    this -> finishDiagAssembly();
    diagWeightsCoo_.extractArrays( i_diag_sparse, j_diag_sparse, diag_sparse );
    diagWeightsCoo_.clear();

    // map global row indices to local subdomain indices
    indRow    = -1;
    indRowLoc = -1;
    for ( unsigned inz = 0; inz < diag_sparse.size(); inz++ ) {
        indRow = i_diag_sparse[inz];
        Global2LocalMap_::iterator pos = global2LocalDofMap_.find( static_cast<unsigned> ( indRow ) );
        OLD_ASSERT( pos != global2LocalDofMap_.end(),
                                "Cannot remap index %d to local indices. \n ", indRow );
        indRowLoc = static_cast<int> ( pos -> second );

        i_diag_sparse[ inz ] = indRowLoc;
    }

    int lsub_diagonal = numDofsSub_;
    OLD_ASSERT( lsub_diagonal == static_cast<int>(diag_sparse.size()),
            "Array length mismatch: %d %d . \n ", lsub_diagonal, static_cast<int>(diag_sparse.size()) );

    std::vector<double> sub_diagonal( lsub_diagonal, -1. );
    // permute the vector according to subdomain indexing
    for ( int i = 0; i < lsub_diagonal; i++ ){
        indRowLoc    = i_diag_sparse[i];
        double value = diag_sparse[i];
        sub_diagonal[indRowLoc] = value;
    }
    OLD_ASSERT( std::find(sub_diagonal.begin(), sub_diagonal.end(), -1. ) == sub_diagonal.end(),
            "There are missing entries in the diagonal of weights. \n " );

    // remove const attribute
    int numDofsInt    = static_cast<int> ( numDofs_ );
    int numDofsSubInt = static_cast<int> ( numDofsSub_ );

    int matrixTypeInt;
    if      ( matrixType_ == GENERAL )                  matrixTypeInt = 0 ;
    else if ( matrixType_ == SPD )                      matrixTypeInt = 1 ;
    else if ( matrixType_ == SYMMETRICGENERAL )         matrixTypeInt = 2 ;
    else if ( matrixType_ == SPD_VIA_SYMMETRICGENERAL ) matrixTypeInt = 2 ; 
    else    {OLD_ASSERT( false, "Illegal matrixType \n " );}

    // download local solution
    int lsol = sol_.size();

    // user constraints are not used
    std::vector<double> userConstraints(1,0.);
    int lUserConstraints1 = 0;
    int lUserConstraints2 = 0;

    int lelement_data1 = 1;
    int lelement_data2 = element_data_.size();
    
    // if > 0, BDDCML will check if the mesh is not composed of several components and if
    // yes, BDDCML will handle them as independent subdomains for selection of coarse 
    // degrees of freedom
    int find_components_int = 1;

    // upload local data to the solver
    bddcml_upload_subdomain_data_c( &numElem_, &numNodes_, &numDofsInt, &nDim_, &meshDim_, 
                                  &iSub, &numElemSub_, &numNodesSub_, &numDofsSubInt, 
                                  &(inet_[0]),&linet, &(nnet_[0]),&lnnet, &(nndf_[0]),&lnndf, 
                                  &(isngn_[0]),&lisngn, &(isvgvn_[0]),&lisvgvn, &(isegn_[0]),&lisegn, 
                                  &(xyz_[0]),&lxyz1,&lxyz2, 
                                  &(ifix_[0]),&lifix, &(fixv_[0]),&lfixv, 
                                  &(rhsVec_[0]),&lrhsVec, &isRhsCompleteInt,
                                  &(sol_[0]), &lsol,
                                  &matrixTypeInt, &(i_sparse[0]), &(j_sparse[0]), &(a_sparse[0]), &la, &isMatAssembledInt,
                                  &(userConstraints[0]),&lUserConstraints1,&lUserConstraints2,
                                  &(element_data_[0]),&lelement_data1,&lelement_data2,
                                  &(sub_diagonal[0]), &lsub_diagonal,
                                  &find_components_int);
    i_sparse.clear();
    j_sparse.clear();
    a_sparse.clear();

    
    //============= setup of BDDC preconditioner
    
    int use_defaults_int          = 0;
    int parallel_division_int     = 0;
    int use_arithmetic_int        = 1;
    int use_adaptive_int;
    if ( use_adaptive ) {
        use_adaptive_int = 1;
    }
    else {
        use_adaptive_int = 0;
    }
    int use_user_constraints_int  = 0;
    int weights_type              = 3;
    
    // if > 0, BDDCML will use constraints on connectivity at corners
    int use_corner_constraints_int = 1;

    // setup multilevel BDDC preconditioner
    START_TIMER("BDDC preconditioner setup");
    bddcml_setup_preconditioner_c( & matrixTypeInt, & use_defaults_int, 
                                 & parallel_division_int, & use_arithmetic_int, 
                                 & use_corner_constraints_int,
                                 & use_adaptive_int, & use_user_constraints_int,
                                 & weights_type );
    END_TIMER("BDDC preconditioner setup");

    //============= compute the norm on arrays with overlapping entries - apply weights
    double normSquaredLoc = 0.;
    if (nProc_ > 1 ) {
        bddcml_dotprod_subdomain_c( &iSub, &(rhsVec_[0]), &lrhsVec, &(rhsVec_[0]), &lrhsVec, &normSquaredLoc );
    } else {
        for ( double elem: rhsVec_ )  normSquaredLoc += elem*elem;
    }

    // sum over processes
    double normSquared = 0.;
    MPI_Allreduce ( &normSquaredLoc, &normSquared, 1, MPI_DOUBLE, MPI_SUM, comm_ );
    normRhs_ = std::sqrt( normSquared );

    //============= solution of the problem
    int method;
    if      ( matrixType_ == SPD )
        method = 0; // PCG - set for SPD matrix
    else if ( matrixType_ == SPD_VIA_SYMMETRICGENERAL)
        method = 0; // PCG - set for benign matrix
    else
        method = 1; // BICGSTAB - set for any other matrix

    int recycling_int = 0;
    int max_number_of_stored_vectors = 100;

    // call iterative solver
    START_TIMER("BDDC solve");
    bddcml_solve_c( & commInt, & method, & tol, & maxIt, & ndecrMax, 
                  & recycling_int, & max_number_of_stored_vectors,
                  & numIter_, & convergedReason_, & condNumber_);
    END_TIMER("BDDC solve");


    //============= downloading solution from BDDCML solver
    
    // download local solution
    lsol = sol_.size();

    bddcml_download_local_solution_c( &iSub, &(sol_[0]), &lsol );

    //============= compute the norm on arrays with overlapping entries - apply weights
    normSquaredLoc = 0.;
    if (nProc_ > 1 ) {
       bddcml_dotprod_subdomain_c( &iSub, &(sol_[0]), &lsol, &(sol_[0]), &lsol, &normSquaredLoc );
    } else {
        for ( double elem: sol_ )  normSquaredLoc += elem*elem;
    }

    // sum over processes
    normSquared = 0.;
    MPI_Allreduce ( &normSquaredLoc, &normSquared, 1, MPI_DOUBLE, MPI_SUM, comm_ );
    normSol_ = std::sqrt( normSquared );

    //============= erase memory used by BDDCML solver
    bddcml_finalize_c();

    return;
}