alecjacobson · April 23, 2024 16:34 · Apr 23, 2024
diff --git a/sparse_hessian_adolc.cpp b/sparse_hessian_adolc.cpp
@@ -0,0 +1,116 @@
+// Adapted from ADOL-C/examples/additional_examples/sparse/sparse_hessian.cpp 
+
+#include <math.h>
+#include <cstdlib>
+#include <cstdio>   
+
+#include <adolc/adolc.h>
+#include <adolc/adolc_sparse.h>
+
+#define tag 1
+
+void print_symmetric(const char* name, int m, int n, double** M) { 
+    int i,j;
+
+    printf("%s =[\n",name);
+    for(i=0; i<m ;i++) { 
+        for(j=0;j<n ;j++)    
+            printf(" %g ", M[i][j]);
+        printf("\n");        
+    }
+    printf("]; %s = %s + tril(%s,-1)';\n",name,name,name);
+}
+
+template <typename T>
+T spring_energy(int n, int dim, T * x)
+{
+  T f = 0;
+  for(int i = 0;i<n;i++)
+  {
+    int j = (i+1)%n;
+    T sqr_len = 0;
+    for(int d = 0;d<dim;d++)
+    {
+      T dijd = x[dim*i + d]-x[dim*j + d];
+      sqr_len += dijd*dijd;
+    }
+    T len = sqrt(sqr_len); 
+
+    T eij = (len-T(1.0));
+    f += eij*eij;
+  }
+  return f;
+}
+
+template <typename T>
+void circle_of_springs(int n, int dim, T * x)
+{
+  for(int i = 0;i<n;i++)
+  {
+    T th = T(i)/T(n) * T(2.0*M_PI);
+    x[i*dim+0] = cos(th);
+    x[i*dim+1] = sin(th);
+  }
+}
+
+void print_positions(int n,int dim,double * x, std::string name = "x")
+{
+  printf("%s = [\n",name.c_str());
+  for(int i = 0;i<n;i++)
+  {
+    printf("%g %g\n",x[dim*i+0],x[dim*i+1]);
+  }
+  printf("];\n");
+}
+
+int main()
+{
+  int n = 5;
+  int dim = 2;
+  double f;
+  double * x = new double[dim*n]; 
+
+  circle_of_springs(n,dim,x);
+  print_positions(n,dim,x,"x");
+  f = spring_energy(n,dim,x);
+  printf("    f = %g;\n",f);
+
+  adouble fad;
+  adouble * xad = new adouble[dim*n];
+  trace_on(tag);
+    for(int i=0;i<n*dim;i++){ xad[i] <<= x[i];}
+    fad = spring_energy(n,dim,xad);
+    fad >>= f;
+  trace_off();
+  printf("fad = %g;\n",f);
+
+  // Full Hessian
+  double **H;
+  H = myalloc2(n*dim,n*dim);
+  hessian(tag,n*dim,x,H);
+  print_symmetric("H_full",n*dim,n*dim,H);
+
+  // Sparse Hessian
+  {
+    unsigned int    *rind  = NULL;
+    unsigned int    *cind  = NULL;
+    double *values = NULL;
+    int nnz;
+    int options[2];
+
+    options[0] = 0;          /*                       safe mode (default) */ 
+    options[1] = 0;          /*               indirect recovery (default) */ 
+    //options[1] = 1;          /*                           direct recovery */ 
+
+    sparse_hess(tag, n*dim, 0, x, &nnz, &rind, &cind, &values, options);
+
+    printf("H = [\n");   
+    for (int i=0;i<nnz;i++)
+        printf("%d %d %g\n",rind[i]+1,cind[i]+1,values[i]);
+    printf("];H = sparse(H(:,1),H(:,2),H(:,3),%d,%d); H = H + triu(H,1)';\n",n*dim,n*dim);
+
+    free(rind); rind = NULL;
+    free(cind); cind = NULL;
+    free(values); values = NULL;
+  }
+}