Instantiation of Vector< Sacado::Fad::DFad<double> >

[Doug]

Hello again,

I am trying to instantiate a Vector with an AD type such as Vector< Sacado::Fad::DFad<double> > by changing

for (SCALAR : REAL_AND_COMPLEX_SCALARS) to for (SCALAR : ALL_SCALAR_TYPES)

in the instantiation file

dealii/source/lac/la_parallel_vector.inst.in

However, it seems 3 issues come up.

1. I can fix this one, so skip ahead if you want.

A bunch of static/dynamic casting of ADvar(unsigned int) through real_type(partitioner->local_size()), but Trilinos doesn't have the unsigned int. Simply need to cast the unsigned int value to a long before casting. For example real_type((long)partitioner->local_size()).

2. I can also fix.

dealii/include/deal.II/base/exceptions.h

Requires the definition of 

   dealii::numbers::is_finite(number)

where I can provide a definition in dealii/include/deal.II/base/numbers.h for the AD types.

3. I don't know how to fix

std::complex<double>(number) needs to be defined. 

Now, this translates to std::complex<double>::complex(Sacado::Rad::ADvar<Sacado::Fad::DFad<double> >&), which of course doesn't exist.

I understand why it's casting into a complex to ensure that we can use the exception for all scalar arguments. But this is limiting the behaviour such that I can't use AD types.

Any suggestions about how to go with this?

This is somewhat linked in the grand scheme of things with

https://groups.google.com/forum/#!searchin/dealii/automatic$20differentiation|sort:date/dealii/9YohjQr1aro/QdtzHHoWAwAJ

where the goal will be to automatically differentiate the entire Jacobian and the sensitivities of the residual with respect to the grid points.

Best regards,

Doug

[Wolfgang Bangerth]

> I am trying to instantiate a Vector with an AD type such as Vector<
> Sacado::Fad::DFad<double> > by changing
>
> for (SCALAR : REAL_AND_COMPLEX_SCALARS) to for (SCALAR : ALL_SCALAR_TYPES)
>
> in the instantiation file
>
> dealii/source/lac/la_parallel_vector.inst.in

This might work, but the scheme is actually supposed to work differently: In
your user code, just #include <lac/la_parallel_vector.templates.h> and
everything will be instantiated in your user code.

> However, it seems 3 issues come up.
>
> 1. I can fix this one, so skip ahead if you want.
>
> A bunch of static/dynamic casting of ADvar(unsigned int)
> through real_type(partitioner->local_size()), but Trilinos doesn't have the

> *unsigned* int. Simply need to cast the unsigned int value to a long before

> casting. For example real_type((long)partitioner->local_size()).

This seems like we are abusing the real_type for something it wasn't intended
to do. Can you open a github issue with the error message you get with the
unmodified code?

> 2. I can also fix.
>
> dealii/include/deal.II/base/exceptions.h
>
> Requires the definition of
>    dealii::numbers::is_finite(number)
>
> where I can provide a definition in dealii/include/deal.II/base/numbers.h for
> the AD types.

I suspect that that would be useful in its own right. Can you open an issue or
pull request for this as well?



> *3. I don't know how to fix*
> *
> *

> std::complex<double>(number) needs to be defined.
>
> Now, this translates
> to std::complex<double>::complex(Sacado::Rad::ADvar<Sacado::Fad::DFad<double>
> >&), which of course doesn't exist.
>
> I understand why it's casting into a complex to ensure that we can use the
> exception for all scalar arguments. But this is limiting the behaviour such
> that I can't use AD types.

Where exactly is this conversion necessary?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Doug]

On Monday, September 9, 2019 at 10:13:32 AM UTC-4, Wolfgang Bangerth wrote:

1.

This seems like we are abusing the real_type for something it wasn't intended
to do. Can you open a github issue with the error message you get with the
unmodified code? 

 

2.

I suspect that that would be useful in its own right. Can you open an issue or
pull request for this as well?

 

Should those be two different issues?

 

> *3. I don't know how to fix* 

Where exactly is this conversion necessary? 

This occurs with AssertIsFinite in include/base/exceptions.h. The condition to be checked uses the include/base/numbers.h function is_finite() function. However, the exception thrown (ExcNumberNotFinite) uses std::complex<double>(number)) to generate the signature. 

[Wolfgang Bangerth]

> 1.
>
> This seems like we are abusing the real_type for something it wasn't intended
> to do. Can you open a github issue with the error message you get with the
> unmodified code?
>
> 2.
>
> I suspect that that would be useful in its own right. Can you open an
> issue or
> pull request for this as well?
>
>
> Should those be two different issues?

Yes, separate issues for separate issues :-)

> > *3. I don't know how to fix*
>
>
> Where exactly is this conversion necessary?
>
>
> This occurs with AssertIsFinite in include/base/exceptions.h. The condition to
> be checked uses the include/base/numbers.h function is_finite() function.
> However, the exception thrown (ExcNumberNotFinite)
> uses std::complex<double>(number)) to generate the signature.

What I meant is: Can you show the compiler error message that illustrates
where the assertion is located, what the template arguments are, how it came
that we called that function with these template arguments, etc?

[Doug]

What I meant is: Can you show the compiler error message that illustrates
where the assertion is located, what the template arguments are, how it came
that we called that function with these template arguments, etc?

/home/ddong/Libraries/dealii/include/deal.II/base/numbers.h:583:3: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘const std::complex<long double>&’

In file included from /home/ddong/Libraries/dealii/include/deal.II/base/cuda.h:21:0,

                 from /home/ddong/Libraries/dealii/include/deal.II/base/memory_space.h:22,

                 from /home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.h:21,

                 from /home/ddong/Libraries/dealii/source/lac/la_parallel_vector.cc:16:

/home/ddong/Libraries/dealii/include/deal.II/base/exceptions.h:1671:42: error: no matching function for call to ‘std::complex<double>::complex(const Sacado::Fad::DFad<double>&)’

          dealii::ExcNumberNotFinite(std::complex<double>(number)))

It is basically part of the same call to AssertIsFinite to generate the Exception. Note that it could hit this assertion much more directly since /home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.templates.h itself calls AssertIsFinite. e.g. line 1450

In the end, it is basically trying to convert the Sacado type into a complex number, which is undefined, whenever it tries to perform some vector operations.

Doug

[Wolfgang Bangerth]

On 9/17/19 6:58 PM, Doug wrote:
>
> /home/ddong/Libraries/dealii/include/deal.II/base/numbers.h:583:3:
> note:   no known conversion for argument 1 from ‘const
> Sacado::Fad::DFad<double>’ to ‘const std::complex<long double>&’
> In file included from
> /home/ddong/Libraries/dealii/include/deal.II/base/cuda.h:21:0,
>                  from
> /home/ddong/Libraries/dealii/include/deal.II/base/memory_space.h:22,
>                  from
> /home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.h:21,
>                  from
> /home/ddong/Libraries/dealii/source/lac/la_parallel_vector.cc:16:
> /home/ddong/Libraries/dealii/include/deal.II/base/exceptions.h:1671:42:
> error: no matching function for call to
> ‘std::complex<double>::complex(const Sacado::Fad::DFad<double>&)’
>           dealii::ExcNumberNotFinite(std::complex<double>(number)))
>
> It is basically part of the same call to AssertIsFinite to generate the
> Exception. Note that it could hit this assertion much more directly
> since /home/ddong/Libraries/dealii/include/deal.II/lac/la_parallel_vector.templates.h
> itself calls AssertIsFinite. e.g. line 1450

But I don't know how we get there. What does the rest of the error
message look like? The place you show looks like this:

template <typename Number, typename MemorySpaceType>
void
Vector<Number, MemorySpaceType>::add(const Number a)
{
AssertIsFinite(a);

so 'a' is of type Sacado::Fad::DFad<double>. It then needs to call

numbers::is_finite (Sacado::Fad::DFad<double>), which doesn't exist.
It probably just tries to go through all of the overloads of
numbers::is_finite() and wants to see whether it can convert the
argument Sacado::Fad::DFad<double> to the argument type of these
overloads. The error message you show then would just explain why this
one possibility (namely, numbers::is_finite(std::complex<long double>&))
does not work. But that doesn't mean that this is the right overload
anyway -- I suspect that your compiler produces similar error messages
above or below the one you show for all of the other overloads, right?

I *think* that the solution is to simply provide an overload for
numbers::is_finite (const Sacado::Fad::DFad<double> &x)
Can you try this? You could declare it in your own .cc file before you
#include <deal.II/lac/la_parallel_vector.templates.h>

[Doug]

so 'a' is of type Sacado::Fad::DFad<double>. It then needs to call

   numbers::is_finite (Sacado::Fad::DFad<double>), which doesn't exist.
It probably just tries to go through all of the overloads of
numbers::is_finite() and wants to see whether it can convert the
argument Sacado::Fad::DFad<double> to the argument type of these
overloads. The error message you show then would just explain why this
one possibility (namely, numbers::is_finite(std::complex<long double>&))
does not work. But that doesn't mean that this is the right overload
anyway -- I suspect that your compiler produces similar error messages
above or below the one you show for all of the other overloads, right?

True, the error message gets long pretty quickly, but the undefined is_finite was another issue. Even if is_finite exists, the complex constructor is still an issue.

 

I *think* that the solution is to simply provide an overload for
   numbers::is_finite (const Sacado::Fad::DFad<double> &x)
Can you try this? You could declare it in your own .cc file before you
#include <deal.II/lac/la_parallel_vector.templates.h>

 Attached is a tiny .cc file that simply instantiates the vector. Copy pasted here here convenience.

#include <Sacado.hpp>

namespace dealii{

namespace numbers{

    bool is_finite(const Sacado::Fad::DFad<double> &x) {

    ¦   (void) x;

    ¦   return true;

    }   

}}

#include <deal.II/lac/vector.h>

#include <deal.II/lac/la_parallel_vector.templates.h>

int main (int /*argc*/, char * /*argv*/[]){

    using namespace dealii;

    dealii::LinearAlgebra::distributed::Vector<double> vector_double;

    using ADtype = Sacado::Fad::DFad<double>;

    dealii::LinearAlgebra::distributed::Vector<ADtype> vector_ad;

    return 0;

}

I also provide an function for numbers::is_finite (const Sacado::Fad::DFad<double> &x) to avoid the first set of errors. However, I still get the error below

In file included from /home/ddong/Libraries/dealii/install/include/deal.II/base/aligned_vector.h:22:0,

                 from /home/ddong/Libraries/dealii/install/include/deal.II/lac/vector.h:22,

                 from /home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:13:

/home/ddong/Libraries/dealii/install/include/deal.II/lac/la_parallel_vector.templates.h: In instantiation of ‘dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>& dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>::operator*=(Number) [with Number = Sacado::Fad::DFad<double>; MemorySpace = dealii::MemorySpace::Host]’:

/home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:26:1:   required from here

/home/ddong/Libraries/dealii/install/include/deal.II/lac/la_parallel_vector.templates.h:1651:7: error: no matching function for call to ‘std::complex<double>::complex(const Sacado::Fad::DFad<double>&)’

       AssertIsFinite(factor);

In file included from /usr/include/trilinos/Teuchos_ConfigDefs.hpp:94:0,

                 from /usr/include/trilinos/Teuchos_PromotionTraits.hpp:45,

                 from /usr/include/trilinos/Sacado_Fad_Exp_GeneralFadTraits.hpp:139,

                 from /usr/include/trilinos/Sacado.hpp:52,

                 from /home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:1:

/usr/include/c++/7/complex:1512:3: note: candidate: constexpr std::complex<double>::complex(const std::complex<long double>&)

   complex<double>::complex(const complex<long double>& __z)

   ^~~~~~~~~~~~~~~

/usr/include/c++/7/complex:1512:3: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘const std::complex<long double>&’

/usr/include/c++/7/complex:1219:26: note: candidate: constexpr std::complex<double>::complex(const std::complex<float>&)

       _GLIBCXX_CONSTEXPR complex(const complex<float>& __z)

                          ^~~~~~~

/usr/include/c++/7/complex:1219:26: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘const std::complex<float>&’

/usr/include/c++/7/complex:1209:26: note: candidate: constexpr std::complex<double>::complex(double, double)

       _GLIBCXX_CONSTEXPR complex(double __r = 0.0, double __i = 0.0)

                          ^~~~~~~

/usr/include/c++/7/complex:1209:26: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘double’

/usr/include/c++/7/complex:1207:26: note: candidate: constexpr std::complex<double>::complex(std::complex<double>::_ComplexT)

       _GLIBCXX_CONSTEXPR complex(_ComplexT __z) : _M_value(__z) { }

                          ^~~~~~~

/usr/include/c++/7/complex:1207:26: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘std::complex<double>::_ComplexT {aka __complex__ double}’

/usr/include/c++/7/complex:1202:12: note: candidate: constexpr std::complex<double>::complex(const std::complex<double>&)

     struct complex<double>

            ^~~~~~~~~~~~~~~

/usr/include/c++/7/complex:1202:12: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘const std::complex<double>&’

/usr/include/c++/7/complex:1202:12: note: candidate: constexpr std::complex<double>::complex(std::complex<double>&&)

/usr/include/c++/7/complex:1202:12: note:   no known conversion for argument 1 from ‘const Sacado::Fad::DFad<double>’ to ‘std::complex<double>&&’

In file included from /home/ddong/Libraries/dealii/install/include/deal.II/base/aligned_vector.h:22:0,

                 from /home/ddong/Libraries/dealii/install/include/deal.II/lac/vector.h:22,

                 from /home/ddong/Codes/PHiLiP/src/instantiate_vector_ad.cpp:13:

So, even if the is_finite function exists, the exception to be thrown still requires the complex constructor std::complex<double>::complex(const Sacado::Fad::DFad<double>&) to exist.

[Jean-Paul Pelteret]

H Doug and Wolfgang,

At the time when we introduced the AD framework I took a stab at quickly adding support for to our Vector class. To summarise, I hit the same issue and couldn’t find an elegant solution past to get past it (so hopefully this is the only outstanding problem to add support to Vector & friends). We have all of the tooling in place to convert AD numbers to something printable, but the issue is that we create circular inclusions in our headers when trying to use them, since the exceptions.h is a very fundamental header and number.h and almost all other headers require it.

My suggestion would be to not attack the problem in exceptions.h, but rather in vector.templates.h and la_parallel_vector.templates.h itself. So each call to AssertIsFinite(), for example over here, should have sent into it a number that’s pre-converted to a std::complex<double>. You should be able to do this by invoking the conversion tool that will extract the value data for all AD numbers (so for free you’d include support for all Sacado and ADOL-C types):

AssertIsFinite(internal::NumberType<std::complex<double>>::value(s));

I’m pretty sure that will work. It’s a bit annoying, but I really think that might be the best way forward. If we go with this approach then I would suggest that we also add a note to the documentation of AssertIsFinite() to highlight this solution.

I hope this helps.

Best,

Jean-Paul

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/7c3d83a7-154c-4bd2-a736-6c5c83f6f829%40googlegroups.com.
<instantiate_vector_ad.cpp>