Re: type discrepancy in rational_fraction operator *

[Konstantinos Poulios]

Hi,

Not sure how well maintained the QD support is. I have never compiled GetFEM with qd support. You should better build it without QD unless you have some specific reason for using high precision arithmetics.

Best regards

Kostas

On Thu, Apr 22, 2021 at 6:53 PM edgar <edgarlux@cryptolab.net> wrote:

Hello, how do I get QD to work with GetFEM?

The compiler says that we are trying to multiply a `rational_fraction'
(vector of type `qd_real') with a `double',
┌────
│ src/bgeot_geometric_trans.cc:894:54: error: ambiguous overload for
‘operator*’ (operand types are ‘bgeot::rational_fraction<qd_real>’ and
‘double’)
│   894 |         trans[0] = (read_base_poly(3, "1-x-y-z") + Q)*0.25;
│       |                    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~
│       |                                                  |    |
│       |                                                  |    double
│       |                                                 
bgeot::rational_fraction<qd_real>
└────

but that is not going to happen, because, `qd_real' is not
`double'. Correct? (the * operator is implicitly defined for the same
type that the vector of `rational_fraction').
┌────
│ template<typename T> class rational_fraction : public std::vector<T> {
│
│     // ---- 8< snip (removed content) ----
│
│     /// Multiply P with the scalar a.
│     rational_fraction operator *(const T &e) const
│     { rational_fraction R = *this; R *= e; return R; }
│
│    // ---- 8< snip (removed content) ----
└────
Listing 1: Excerpts from `src/getfem/bgeot_poly.h'

If I cast the 0.25 as qd_real
┌────
│ trans[0] = (read_base_poly(3, "1-x-y-z") + Q)*qd_real(0.25);
└────

the return type is wrong (the function needs `double').
┌────
│ src/bgeot_poly_composite.cc: In member function ‘bgeot::scalar_type
bgeot::polynomial_composite::eval(const base_node&, bgeot::size_type)
const’:
│ /src/bgeot_poly_composite.cc:126:59: error: cannot convert ‘qd_real’
to ‘bgeot::scalar_type’ {aka ‘double’} in return
│   126 |       if (!local_coordinate) return
poly_of_subelt(l).eval(p.begin());
│       |                                     
~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~
│       |                                                           |
│       |                                                           
qd_real
│ src/bgeot_poly_composite.cc:129:36: error: cannot convert ‘qd_real’ to
‘bgeot::scalar_type’ {aka ‘double’} in return
│   129 |       return poly_of_subelt(l).eval(p1.begin());
│       |              ~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~
│       |                                    |
│       |                                    qd_real
└────

Previously, (naïvely) trying to overload `*' with a double in
`bgeot_poly.h'
┌────
│ rational_fraction operator *(const double &e) const
│ { rational_fraction R = *this; R *= e; return R; }
└────

yields another error, as one can imagine.

I build out of source now (easier to debug the compilation process):
┌────
│ mkdir build
│ cd build
│ ../configure --enable-silent-rules --srcdir=../ --prefix=/usr
--sysconfdir=/etc --bindir=/usr/bin --datadir=/usr/share
--docdir=/usr/share/doc/getfem --with-pic --enable-march --enable-shared
--enable-static=no --enable-python --enable-blas-interface
--enable-blas64-support --enable-metis --enable-paralevel
--enable-multithread-blas --enable-par-mumps --enable-openmp --enable-qd
--with-qd-include-dir=/usr/include/qd/ --enable-qhull --disable-superlu
│ make all
└────

It seems that QD is found:
┌────
│ Libraries Used:
│ ---------------
│
│ - QD library found. High precision (quad-double precision) polynomials
│   and integration methods are enabled.
│ - Qhull found. Using the Qhull library for delaunay triangulations.
│ - Mumps found. A direct solver for large sparse linear systems.
│ - Lapack library found: -llapack
│ - BLAS library found. Link options: -lblas
│   You can give the location of your prefered blas library with either
│   the --with-blas=<lib> option, or the BLAS_LIBS environment variable
│   for example: ./configure BLAS_LIBS="-L/usr/lib/sse2/atlas/ -lblas"
│   Use blas with 64 bits integers or 32/64 bits compatibility mode
└────

I have this in my system
━━━━━━━━━━━━━━━━━━━━━━━━
  qhull         2020.2-4
  qd            2.3.22-3
  mumps          5.3.5-1
  openmp        11.1.0-1
  openmpi        4.0.5-3
  openblas      0.3.13-2
  scalapack      2.1.0-1
  lapack         3.9.0-3
  superlu        5.2.2-2
  metis      5.1.0.p10-1
  python         3.9.2-1
  gcc           10.2.0-6
  make             4.3-3
━━━━━━━━━━━━━━━━━━━━━━━━

Thank you in advance for your help.