Don't allow non-differentiable non-const array parameters in forward mode

All parameters w.r.t. which we don't differentiate still require adjoints. Since the user doesn't provide it, we have to construct it ourselves. However, if the parameter is a pointer array (e.g., `double*`), we don't know its size and, therefore, cannot initialize the adjoint.
If the array is const, it's guaranteed that it doesn't require an adjoint. If it's a const array (e.g., `double[3]`), we can initialize it. In other cases, we have to produce an error. We've been doing the same in the reverse mode for a while.

Fixes #1035.

Author: PetroZarytskyi

benchmark/BenchmarkedFunctions.h

@@ -41,7 +41,7 @@ inline double product(double p[], int n) {
 }
 
 ///\returns the weighted sum of the elements in \p
-inline double weightedSum(double p[], double w[], int n) {
+inline double weightedSum(const double p[], const double w[], int n) {
   double sum = 0;
   for (int i = 0; i < n; i++)
     sum += p[i] * w[i];

demos/Arrays.cpp

@@ -14,7 +14,7 @@
 // Necessary for clad to work include
 #include "clad/Differentiator/Differentiator.h"
 
-double weighted_avg(double* arr, const double* weights) {
+double weighted_avg(const double* arr, const double* weights) {
   return (arr[0] * weights[0] + arr[1] * weights[1] + arr[2] * weights[2]) / 3;
 }
 

lib/Differentiator/BaseForwardModeVisitor.cpp

@@ -317,16 +317,34 @@ void BaseForwardModeVisitor::GenerateSeeds(const clang::FunctionDecl* dFD) {
     // relation among them, thus it is safe (correct) to use the corresponding
     // non-reference type for creating the derivatives.
     QualType dParamType = param->getType().getNonReferenceType();
-    // We do not create derived variable for array/pointer parameters.
-    if (!utils::IsDifferentiableType(dParamType) ||
-        utils::isArrayOrPointerType(dParamType))
-      continue;
     Expr* dParam = nullptr;
-    if (dParamType->isRealType()) {
+    if (!utils::IsDifferentiableType(dParamType))
+      continue;
+    // If the parameter type decayed const array type, we can still initialize
+    // it. Therefore, we don't have to produce the error.
+    if (const auto* DT = dyn_cast<DecayedType>(dParamType))
+      dParamType = DT->getOriginalType();
+    if (dParamType->isConstantArrayType()) {
+      if (param == m_IndependentVar)
+        continue;
+      dParam = getZeroInit(dParamType);
+    } else if (dParamType->isRealType()) {
       // If param is independent variable, its derivative is 1, otherwise 0.
       int dValue = (param == m_IndependentVar);
       dParam =
           ConstantFolder::synthesizeLiteral(m_Context.IntTy, m_Context, dValue);
+    } else if (utils::isArrayOrPointerType(dParamType)) {
+      // We cannot initialize a pointer array adjoint ourselves.
+      // Produce an error.
+      if (param != m_IndependentVar &&
+          !utils::GetValueType(dParamType).isConstQualified()) {
+        // FIXME: Use diagnostics style as in #1596
+        diag(DiagnosticsEngine::Error, param->getLocation(),
+             "dependent non-const pointer and array parameters "
+             "are not supported; differentiate w.r.t. '%0' or mark it const",
+             {param->getNameAsString()});
+      }
+      continue;
     }
     // For each function arg, create a variable _d_arg to store derivatives
     // of potential reassignments, e.g.:

test/Arrays/ArrayErrorsForwardMode.C

@@ -14,10 +14,24 @@ double addDoubleArr(Double *arr) { // expected-error {{attempted differentiation
   return arr[0].n + arr[1].n + arr[2].n + arr[3].n;
 }
 
+double nonConstNonDiffArr(double x, double *y) { // expected-error {{dependent non-const pointer and array parameters are not supported; differentiate w.r.t. 'y' or mark it const}}
+    y[0] = 3; // expected-warning {{derivative of an assignment attempts to assign to unassignable expr, assignment ignored}}
+    return x * y[0];
+}
+
+double nonConstNonDiffConstArr(double x, double y[3]) {
+    y[0] = 3;
+    return x * y[0];
+}
+
 int main() {
   clad::differentiate(addArr, "arr[1:2]"); // expected-error {{Forward mode differentiation w.r.t. several parameters at once is not supported, call 'clad::differentiate' for each parameter separately}}
 
   clad::differentiate(addArr, "arr[2:1]"); // expected-error {{Range specified in 'arr[2:1]' is in incorrect format}}
 
   clad::differentiate(addDoubleArr, "arr[1]");
+
+  clad::differentiate(nonConstNonDiffArr, "x");
+
+  clad::differentiate(nonConstNonDiffConstArr, "x");
 }

test/Arrays/ArrayInputsForwardMode.C

@@ -41,7 +41,7 @@ double addArr(const double *arr, int n) {
 //CHECK-NEXT:       return _d_ret;
 //CHECK-NEXT:   }
 
-double numMultIndex(double* arr, size_t n, double x) {
+double numMultIndex(const double* arr, size_t n, double x) {
   // compute x * i, where arr[i] = x
   // if x is not present in arr, return 0
   bool flag = false;
@@ -56,7 +56,7 @@ double numMultIndex(double* arr, size_t n, double x) {
   return flag ? idx * x : 0;
 }
 
-// CHECK:   double numMultIndex_darg2(double *arr, size_t n, double x) {
+// CHECK:   double numMultIndex_darg2(const double *arr, size_t n, double x) {
 // CHECK-NEXT:     size_t _d_n = 0;
 // CHECK-NEXT:     double _d_x = 1;
 // CHECK-NEXT:     bool _d_flag = 0;

test/FirstDerivative/CallArguments.C

@@ -132,24 +132,24 @@ float f_const_args_func_8(const float x, float y) {
 // CHECK-NEXT: return _t0.pushforward + _t1.pushforward - _d_y;
 // CHECK-NEXT: }
 
-float f_literal_helper(float x, char ch, float* p, float* q) {
+float f_literal_helper(float x, char ch, const float* p, float* q) {
   if (ch == 'a')
     return x * x;
   return -x * x;
 }
 
-// CHECK: clad::ValueAndPushforward<float, float> f_literal_helper_pushforward(float x, char ch, float *p, float *q, float _d_x, char _d_ch, float *_d_p, float *_d_q) {
+// CHECK: clad::ValueAndPushforward<float, float> f_literal_helper_pushforward(float x, char ch, const float *p, float *q, float _d_x, char _d_ch, const float *_d_p, float *_d_q) {
 // CHECK-NEXT:     if (ch == 'a')
 // CHECK-NEXT:         return {x * x, _d_x * x + x * _d_x};
 // CHECK-NEXT:     return {-x * x, -_d_x * x + -x * _d_x};
 // CHECK-NEXT: }
 
-float f_literal_args_func(float x, float y, float *z) {
+float f_literal_args_func(float x, float y, const float *z) {
   printf("hello world ");
   return x * f_literal_helper(x, 'a', z, nullptr);
 } // x ^ 3
 
-// CHECK: float f_literal_args_func_darg0(float x, float y, float *z) {
+// CHECK: float f_literal_args_func_darg0(float x, float y, const float *z) {
 // CHECK-NEXT: float _d_x = 1;
 // CHECK-NEXT: float _d_y = 0;
 // CHECK-NEXT: printf("hello world ");