Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[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. 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>::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

> On 18 Sep 2019, at 03:40, Doug  wrote:
> 
> so 'a' is of type Sacado::Fad::DFad. It then needs to call 
> 
>numbers::is_finite (Sacado::Fad::DFad), 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 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&)) 
> 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 ) 
> Can you try this? You could declare it in your own .cc file before you 
> #include  
> 
>  Attached is a tiny .cc file that simply instantiates the vector. Copy pasted 
> here here convenience.
> 
> #include 
> namespace dealii{
> namespace numbers{
> bool is_finite(const Sacado::Fad::DFad ) {
> ¦   (void) x;
> ¦   return true;
> }   
> }}
> #include 
> #include 
> int main (int /*argc*/, char * /*argv*/[]){
> using namespace dealii;
> dealii::LinearAlgebra::distributed::Vector vector_double;
> using ADtype = Sacado::Fad::DFad;
> dealii::LinearAlgebra::distributed::Vector vector_ad;
> return 0;
> }
> 
> I also provide an function for numbers::is_finite (const 
> Sacado::Fad::DFad ) 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 MemorySpace>& dealii::LinearAlgebra::distributed::Vector MemorySpace>::operator*=(Number) [with Number = Sacado::Fad::DFad; 
> 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::complex(const 
> Sacado::Fad::DFad&)’
>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::complex(const std::complex&)
>complex::complex(const complex& __z)
>^~~
> /usr/include/c++/7/complex:1512:3: note:   no known conversion for argument 1 
> from ‘const Sacado::Fad::DFad’ to ‘const std::complex&’
> 

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Doug]

>
> so 'a' is of type Sacado::Fad::DFad. It then needs to call 
>
>numbers::is_finite (Sacado::Fad::DFad), 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 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&)) 
> 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 ) 
> Can you try this? You could declare it in your own .cc file before you 
> #include  
>

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

#include 
namespace dealii{
namespace numbers{
bool is_finite(const Sacado::Fad::DFad ) {
¦   (void) x;
¦   return true;
}   
}}
#include 
#include 
int main (int /*argc*/, char * /*argv*/[]){
using namespace dealii;
dealii::LinearAlgebra::distributed::Vector vector_double;
using ADtype = Sacado::Fad::DFad;
dealii::LinearAlgebra::distributed::Vector vector_ad;
return 0;
}

I also provide an function for numbers::is_finite (const 
Sacado::Fad::DFad ) 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& dealii::LinearAlgebra::distributed::Vector::operator*=(Number) [with Number = Sacado::Fad::DFad; 
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::complex(const Sacado::Fad::DFad&)’
   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::complex(const std::complex&)
   complex::complex(const complex& __z)
   ^~~
/usr/include/c++/7/complex:1512:3: note:   no known conversion for argument 
1 from ‘const Sacado::Fad::DFad’ to ‘const std::complex&’
/usr/include/c++/7/complex:1219:26: note: candidate: constexpr 
std::complex::complex(const std::complex&)
   _GLIBCXX_CONSTEXPR complex(const complex& __z)
  ^~~
/usr/include/c++/7/complex:1219:26: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to ‘const 
std::complex&’
/usr/include/c++/7/complex:1209:26: note: candidate: constexpr 
std::complex::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’ to ‘double’
/usr/include/c++/7/complex:1207:26: note: candidate: constexpr 
std::complex::complex(std::complex::_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’ to 
‘std::complex::_ComplexT {aka __complex__ double}’
/usr/include/c++/7/complex:1202:12: note: candidate: constexpr 
std::complex::complex(const std::complex&)
 struct complex
^~~
/usr/include/c++/7/complex:1202:12: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to ‘const 
std::complex&’
/usr/include/c++/7/complex:1202:12: note: candidate: constexpr 
std::complex::complex(std::complex&&)
/usr/include/c++/7/complex:1202:12: note:   no known conversion for 
argument 1 from ‘const Sacado::Fad::DFad’ to 
‘std::complex&&’
In file included from 
/home/ddong/Libraries/dealii/install/include/deal.II/base/aligned_vector.h:22:0,
 from 

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[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’ to ‘const std::complex&’
> 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::complex(const Sacado::Fad::DFad&)’
>            dealii::ExcNumberNotFinite(std::complex(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 
 void
 Vector::add(const Number a)
 {
   AssertIsFinite(a);

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

   numbers::is_finite (Sacado::Fad::DFad), 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 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&)) 
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 )
Can you try this? You could declare it in your own .cc file before you 
#include 

Best
  W.

-- 

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

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/846b94a1-9c0c-babc-daa3-11653974365b%40colostate.edu.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[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’ 
to ‘const std::complex&’
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::complex(const Sacado::Fad::DFad&)’
  dealii::ExcNumberNotFinite(std::complex(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

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5fae8c77-687b-4127-b6b3-e62a05b24949%40googlegroups.com.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[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(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?

Best
  W.

-- 

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

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a3b717a8-437e-72bf-22b0-a231795f4f57%40colostate.edu.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[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(number)) to generate the signature. 

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/6231eb16-d704-4862-9647-82dd0f40c49c%40googlegroups.com.

Re: [deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Wolfgang Bangerth]

> I am trying to instantiate a Vector with an AD type such as Vector< 
> Sacado::Fad::DFad > 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  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(number) needs to be defined.
> 
> Now, this translates 
> to std::complex::complex(Sacado::Rad::ADvar 
>  >&), 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: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/efa1f1fd-0609-0dcc-9e1a-476e813d17ee%40colostate.edu.

[deal.II] Instantiation of Vector< Sacado::Fad::DFad >

[Doug]

Hello again,

I am trying to instantiate a Vector with an AD type such as Vector< 
Sacado::Fad::DFad > 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(number) needs to be defined. 

Now, this translates 
to std::complex::complex(Sacado::Rad::ADvar 
>&), 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

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/48459ebb-b016-4757-be9a-6a4766ad1621%40googlegroups.com.