Eigen Datatypes

[sriharikrishna]

Hi,
I am trying to differentiate a code that uses Eigen::ArrayXd Eigen::ArrayXXd and Eigen::MatrixXd . It appears that CoDiPack supports Eigen::Matrix at the moment. Is there a workaround that you can suggest or is this an extension that must be added to LinearSystemInterfaceTypes.hpp?
Thanks
Krishna

[MaxSagebaum]

Hi, if you use one of the CoDiPack types as the template argument it should work out of the box. The EigenLinearSystem class in CoDiPack is just tested with Eigen::Vector and Eigen::MatrixXd. If it also works with Eigen::ArrayXXd I do not know.

What problems do you have? Did you have a look at https://scicomp.rptu.de/codi/d2/d8e/Example_21_Special_handling_of_linear_system_solvers.html ?

[sriharikrishna]

Thanks @MaxSagebaum. The code uses Eigen data structures as a convenient way to create and manipulate matrices and vectors before passing them off to LAPACK/MKL calls like cblas_dsyr2k. Is this usage of Eigen supported?

[MaxSagebaum]

If you can formulate this with CoDiPack types, than it is supported. So Eigen::Matrix, Eigen::Map or so should not be any problem. Passing the CoDiPack structures through LAPACK/MKL requires the external functions of CoDiPack. If you just want to solve linear systems, then you can implement the codi::LinearSystemInterface interface for this. (See https://scicomp.rptu.de/codi/d2/d8e/Example_21_Special_handling_of_linear_system_solvers.html and have a look at the implementation of codi::EigenLinearSystem.

It seems that cblas_dsyr2k is not a linear system solve. You have to options:
Option 1: Create an external function for each of these methods. If they are not called that often, the conversion overhead from CoDiPack types to regular types is not that problematic.

Option 2: A faster and more thorough approach is along the lines of adjoint_petsc https://github.com/SciCompKL/adjoint-PETSc . Here, the primal data and AD data is split. So you basically have a matrix with your primal data and a matrix with the same layout as the primal for the AD data. This allows to call functions like cblas_dsyr2k on the primal data and afterwards add the adjoint formulation via external functions. (E.g. https://github.com/SciCompKL/adjoint-PETSc/blob/6dfccb3dca2e6e419d9a777134d3764f7a909425/src/adjoint_petsc/ksp.cpp#L183). No conversion overhead is necessary. The downside on this approach is, that you do not have the regular CoDiPack types for the tracking of regular operations. You then have to reprogram the logic, so that a pseudo CoDiPack type is returned (https://github.com/SciCompKL/adjoint-PETSc/blob/6dfccb3dca2e6e419d9a777134d3764f7a909425/src/adjoint_petsc/vec.cpp#L643, https://github.com/SciCompKL/adjoint-PETSc/blob/6dfccb3dca2e6e419d9a777134d3764f7a909425/src/adjoint_petsc/util/vec_iterator_util.hpp#L75, https://github.com/SciCompKL/adjoint-PETSc/blob/6dfccb3dca2e6e419d9a777134d3764f7a909425/include/adjoint_petsc/util/codi_def.hpp#L41 (sorry the generalization made this quite convoluted.))
Here a short example for the wrapper:

double a = ...;
int id = ...;

auto a_codi = codi::ActiveTypeWrapper<codi::RealReverse>(a, id); // <- behaves like a regular CoDiPack type.

Or you can program something along the lines of https://github.com/SciCompKL/adjoint-PETSc/blob/master/include/adjoint_petsc/util/wrapper_array.hpp.

Hope this helps. If you want we can discuss this in video call.

[MaxSagebaum]

Did this solve your issue?

[sriharikrishna]

Working on it. Do you handle Eigen complex types?

[MaxSagebaum]

I was not aware that eigen has a special complex type. What I could find, is that they support std::complex. Do you have a link?

It was already possible to use std::complex with CoDiPack types. Since version 3.0, there is specialization of std::complex for CoDiPack types. This adds complex types to the expression templates of CoDiPack.