Sparse Hessian of a mass-spring system

sparse_hessian_adolc.cpp

// 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;
  }
}