Preconditioner and solver setup to solve a case using boost::numeric::ublas/CPU/openMP?

[tario]

Hello,

there's something incorrect with the setup of the preconditioner and solver I think. I keep getting compile errors. I tried to follow this example (http://viennacl.sourceforge.net/viennacl-examples-iterative.html) but I can't get it right

Note: This is my base implementation to get the interface right, it's supposed to compute the case on the CPU using OpenMP. GPU computing will be the next step.

...
//
// Solve case using ViennaCL 
//      
    uint n = psi.size();

    boost::numeric::ublas::compressed_matrix<scalar> ublas_matrix(n,n); // create matrix
    ldu2vcl(matrix, ublas_matrix); // load external matrix into ublas_matrix

    boost::numeric::ublas::vector<scalar> vcl_psi(n); // x
    boost::numeric::ublas::vector<scalar> vcl_source(n); // rhs

    // Set solver tolerance and maximum iterations
    viennacl::linalg::cg_tag cg(tolerance_, maxIter_); 

    // Set preconditioner
    typedef viennacl::linalg::ilut_precond< 
                        boost::numeric::ublas::compressed_matrix<scalar> > ublas_ilut_t;

    ublas_ilut_t ublas_ilut(ublas_matrix, viennacl::linalg::ilut_tag());  

    // Solve the linear system
    vcl_psi  = solve(ublas_matrix,   //using ublas objects on CPU
                     vcl_source,
                     viennacl::linalg::cg_tag(), 
                     ublas_ilut_t); // error?
...

What's wrong here?

T.

[Karl Rupp]

It's hard to tell without the error message. Can you please provide the first couple of compiler error lines? Did you include the relevant headers (cf. example)?

Also note that when using ublas, there is no benefit from OpenMP, as ublas doesn't use OpenMP for its operations.

[tario]

Then I'm mislead.

Which configuration offers OpenMP support on the CPU? Is there an example?

[tario]

Here are the relevant compiler error lines, others are warnings about old-style-type-cast:

...
vclPCG.C: In member function ‘virtual Foam::solverPerformance Foam::vclPCG::solve(Foam::scalarField&, const scalarField&, Foam::direction) const’:
vclPCG.C:167:31: error: invalid use of non-static member function
ldu2vcl(matrix, vcl_matrix); // load external matrix into vcl_matrix
^
vclPCG.C:177:64: error: no matching function for call to ‘viennacl::linalg::ilut_precond<viennacl::compressed_matrix<double> >::ilut_precond()’
viennacl::compressed_matrix<ScalarType> > vcl_ilut_t;
^
In file included from vclPCG.C:32:0:
viennacl/linalg/detail/ilu/ilut.hpp:415:3: note: candidate: viennacl::linalg::ilut_precond<viennacl::compressed_matrix<T, AlignmentV=""> >::ilut_precond(const MatrixType&, const viennacl::linalg::ilut_tag&) [with NumericT = double; unsigned int AlignmentV = 1u; viennacl::linalg::ilut_precond<viennacl::compressed_matrix<T, AlignmentV=""> >::MatrixType = viennacl::compressed_matrix<double>]
ilut_precond(MatrixType const & mat, ilut_tag const & tag)
^
viennacl/linalg/detail/ilu/ilut.hpp:415:3: note: candidate expects 2 arguments, 0 provided
viennacl/linalg/detail/ilu/ilut.hpp:410:7: note: candidate: viennacl::linalg::ilut_precond<viennacl::compressed_matrix<double> >::ilut_precond(const viennacl::linalg::ilut_precond<viennacl::compressed_matrix<double> >&)
class ilut_precond< viennacl::compressed_matrix<NumericT, AlignmentV=""> >
^
viennacl/linalg/detail/ilu/ilut.hpp:410:7: note: candidate expects 1 argument, 0 provided
viennacl/linalg/detail/ilu/ilut.hpp:410:7: note: candidate: viennacl::linalg::ilut_precond<viennacl::compressed_matrix<double> >::ilut_precond(viennacl::linalg::ilut_precond<viennacl::compressed_matrix<double> >&&)
viennacl/linalg/detail/ilu/ilut.hpp:410:7: note: candidate expects 1 argument, 0 provided
vclPCG.C:179:16: error: expected ‘;’ before ‘vcl_ilut’
vcl_ilut_t vcl_ilut(vcl_matrix, ilut_conf);
^
vclPCG.C:179:47: warning: statement has no effect [-Wunused-value]
vcl_ilut_t vcl_ilut(vcl_matrix, ilut_conf);
^
vclPCG.C:185:23: error: ‘vcl_ilut’ was not declared in this scope
vcl_ilut);
^
vclPCG.C:187:55: warning: use of old-style cast [-Wold-style-cast]
solverPerf.finalResidual() = (ScalarType)cg.error();
^
vclPCG.C:188:53: warning: use of old-style cast [-Wold-style-cast]
solverPerf.nIterations() = (ScalarType)cg.iters();
^
/opt/openfoam5/wmake/rules/General/transform:25: die Regel für Ziel „Make/linux64GccDPInt32Opt/vclPCG.o“ scheiterte
make: *** [Make/linux64GccDPInt32Opt/vclPCG.o] Fehler 1

The main program code is:

#define VIENNACL_WITH_OPENMP // enables execution on multiple CPU cores via OpenMP backend

#include "vclPCG.H"

// ViennaCL
#include "viennacl/vector.hpp"
#include "viennacl/linalg/detail/ilu/ilut.hpp"
#include "viennacl/linalg/cg.hpp"
#include "viennacl/compressed_matrix.hpp"

typedef double    ScalarType;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(vclPCG, 0);

    lduMatrix::solver::addsymMatrixConstructorToTable<vclPCG>
        addvclPCGSymMatrixConstructorToTable_;
}

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::vclPCG::vclPCG
(
    const word& fieldName,
    const lduMatrix& matrix,
    const FieldField<Field, scalar>& interfaceBouCoeffs,
    const FieldField<Field, scalar>& interfaceIntCoeffs,
    const lduInterfaceFieldPtrsList& interfaces,
    const dictionary& solverControls
)
:
    lduMatrix::solver
    (
        fieldName,
        matrix,
        interfaceBouCoeffs,
        interfaceIntCoeffs,
        interfaces,
        solverControls
    )
{}

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

//
// Convert matrix from LDU format to original format, storing only nonzero values
//

void ldu2org
( 
    const Foam::lduMatrix & matrix, 
    uint* rows, // row and column indexes must be unsigned int for OpenMP and CUDA and cl_uint for OpenCL
    uint* cols, // column                                        
    ScalarType* vals // values                                   
) 
{
    int k = 0; // Index to fill lower triagonal part, diagonal part and upper triagonal part successively into one array

    //  
    // Fill original format matrix arrays with nonzero values
    //

    // Fill original format matrix with lower triangle
    for (int i=0; i < matrix.lower().size() ; i++) 
    {
        rows[k] = matrix.lduAddr().upperAddr()[i];
        cols[k] = matrix.lduAddr().lowerAddr()[i];
        vals[k] = matrix.lower()[i];
    k++;
    } 

    // Fill original format matrix with diagonal
    for (int i=0; i < matrix.diag().size(); i++) 
    {
        rows[k] = matrix.lduAddr().upperAddr()[i];
        cols[k] = matrix.lduAddr().lowerAddr()[i];
        vals[k] = matrix.diag()[i];
    k++;
    } 

    // Fill original format matrix with upper triangle
    for (int i=0; i < matrix.upper().size() ; i++) 
    {
        rows[k] = matrix.lduAddr().lowerAddr()[i];
        cols[k] = matrix.lduAddr().upperAddr()[i];
        vals[k] = matrix.upper()[i];
    k++;
    } 
}

void ldu2vcl
(
    const Foam::lduMatrix & matrix, 
    viennacl::compressed_matrix<ScalarType> & vcl_matrix
)
{
    uint nnz = matrix.lower().size() + matrix.upper().size() + matrix.diag().size(); // number of nonzero values (nnz)

    // Allocate memory for original sparse matrix 
    uint * rows = (uint *)calloc(nnz, sizeof(uint)); // rows
    uint * cols = (uint *)calloc(nnz, sizeof(uint)); // columns
    ScalarType * vals = (ScalarType *)calloc(nnz, sizeof(ScalarType)); // values

    ldu2org(matrix,rows,cols,vals); // row, column, value is the default order when adding values to an vcl_matrix

    // Fill the compressed vcl_matrix
    for (uint i=0; i < nnz; i++)
    {
    vcl_matrix(rows[i],cols[i]) = vals[i];
    }

    // Free and release memory of for original sparse matrix
    free(rows); free(cols); free(vals); //colloc()  
}

Foam::solverPerformance Foam::vclPCG::solve
(
    scalarField& psi,
    const scalarField& source,
    const direction cmpt
) const
{

    // --- Setup class containing solver performance data
    solverPerformance solverPerf
    (
        lduMatrix::preconditioner::getName(controlDict_) + typeName,
        fieldName_
    );

//
// Solve case using ViennaCL 
//
    uint n = psi.size();

    viennacl::compressed_matrix<ScalarType> vcl_matrix; // create matrix
    ldu2vcl(matrix, vcl_matrix); // load external matrix into vcl_matrix

    viennacl::vector<ScalarType> vcl_psi(n); // x
    viennacl::vector<ScalarType> vcl_source(n); // rhs

    // Set solver tolerance and maximum iterations
    viennacl::linalg::cg_tag cg(tolerance_, maxIter_); 

    // Set preconditioner
    viennacl::linalg::ilut_precond<
                    viennacl::compressed_matrix<ScalarType> >  vcl_ilut_t;

    vcl_ilut_t vcl_ilut(vcl_matrix, ilut_conf);

    // Solve the linear system
    vcl_psi  = solve(vcl_matrix,     //using viennacl objects on GPU
                      vcl_source,
                      viennacl::linalg::cg_tag(),
                      vcl_ilut);

    solverPerf.finalResidual() = (ScalarType)cg.error();
    solverPerf.nIterations() = (ScalarType)cg.iters();

    return solverPerf;
}

[Karl Rupp]

Thanks! So in line 177 you are missing a typedef:

typedef viennacl::...    vcl_ilut_t;

The first error in line 167 might be a problem with function visibility, but I'm not entirely sure.

[tario]

Thank you, that helped and at the moment there's one error left:

clPCG_solve.C:47:37: error: ‘ilut_conf’ was not declared in this scope
 vcl_ilut_t vcl_ilut(vcl_matrix, ilut_conf);

 What's the correct syntax for that?.

[Karl Rupp]

You need to declare the configuration tag. That is, either replace ilut_conf with viennacl::linalg::ilut_tag() or declare

viennacl::linalg::ilut_tag ilut_conf;

before line 47.