LocalADOLCStiffness: use "hessian2" instead of "hess_vec"

With this commit, we do two things:

• replace "hess_vec" with a simpler call of "hessian" which calls internally "hess_vec".

• use "hessian2" by default instead of "hessian". "hessian2" calls internally "hess_mat" which is faster than "hessian/hess_vec". If the user encounters probems with the new "hessian2" driver, a boolean "useHessian2" can be set to false in the constructor.

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -110,6 +110,7 @@ assembleGradientAndHessian(const LocalView& localView,
 
     localHessian.setSize(nDoubles,nDoubles);
 
+#if 1 // use hess_vec
     std::vector<double> v(nDoubles);
     std::vector<double> w(nDoubles);
 
@@ -131,8 +132,109 @@ assembleGradientAndHessian(const LocalView& localView,
       // Make v the null vector again
       v[i] = 0.;
     }
+#elif 0 // use hessian
+    // allocate raw doubles
+    double* rawHessian[nDoubles];
+    for(size_t i=0; i<nDoubles; i++)
+    {
+      // we need only the lower half of the hessian, thus allocate i+1 doubles only
+      rawHessian[i] = (double*)malloc((i+1)*sizeof(double));
+    }
+
+    // here comes the magic
+    // ADOL-C won't accept pointers with const values, therefore the const_cast
+    hessian(rank,nDoubles,const_cast<double*>(localConfiguration.data()),rawHessian);
+
+    // copy the raw values to the localHessian
+    for (size_t i=0; i<nDoubles; i++)
+    {
+      for (size_t j=0; j<i; j++)
+      {
+        localHessian[i][j] = rawHessian[i][j];
+        localHessian[j][i] = rawHessian[i][j];
+      }
+      localHessian[i][i] = rawHessian[i][i];
+    }
+
+    // don't forget the clean everything up
+    for(size_t i=0; i<nDoubles; i++)
+      free(rawHessian[i]);
+
+#elif 0 // use hess_mat
+    // allocate raw doubles
+    double* rawHessian[nDoubles];
+    double* rawIdentity[nDoubles];
+    for(size_t i=0; i<nDoubles; i++)
+    {
+      rawHessian[i] = (double*)malloc(nDoubles*sizeof(double));
+      rawIdentity[i] = (double*)malloc(nDoubles*sizeof(double));
+      for(size_t j=0; j<nDoubles; j++)
+        rawIdentity[i][j] = i==j ? 1. : 0.;
+    }
 
-    // TODO: implement ADOLC_VECTOR_MODE variant
+    // here comes the magic
+    // ADOL-C won't accept pointers with const values, therefore the const_cast
+    hess_mat(rank,nDoubles,nDoubles,const_cast<double*>(localConfiguration.data()),rawIdentity,rawHessian);
+
+    // copy the raw values to the localHessian
+    for (size_t i=0; i<nDoubles; i++)
+    {
+      for (size_t j=0; j<i; j++)
+      {
+        localHessian[i][j] = rawHessian[i][j];
+        localHessian[j][i] = rawHessian[i][j];
+      }
+      localHessian[i][i] = rawHessian[i][i];
+    }
+
+    // don't forget the clean everything up
+    for(size_t i=0; i<nDoubles; i++)
+    {
+      free(rawHessian[i]);
+      free(rawIdentity[i]);
+    }
+#else // use hess_mat with slices
+    // slice the hessian in col blocks of size "blocksize" and repeat
+    int blocksize = 6;
+    // allocate raw doubles
+    double* rawHessian[nDoubles];
+    double* rawIdentity[nDoubles];
+    for(size_t i=0; i<nDoubles; i++)
+    {
+      rawHessian[i] = (double*)malloc(blocksize*sizeof(double));
+      rawIdentity[i] = (double*)malloc(blocksize*sizeof(double));
+    }
+
+    for( size_t block=0; block*blocksize<nDoubles; block++ )
+    {
+      for(size_t i=0; i<nDoubles; i++)
+      {
+        for(size_t j=0; j<blocksize; j++)
+          rawIdentity[i][j] = i==(block*blocksize+j) ? 1. : 0.;
+      }
+
+      // here comes the magic
+      // ADOL-C won't accept pointers with const values, therefore the const_cast
+      hess_mat(rank,nDoubles,blocksize,const_cast<double*>(localConfiguration.data()),rawIdentity,rawHessian);
+
+      // copy the raw values to the localHessian
+      for (size_t i=0; i<nDoubles; i++)
+      {
+        for (size_t j=0; j<blocksize; j++)
+        {
+          if ( block*blocksize+j < nDoubles )
+            localHessian[i][block*blocksize+j] = rawHessian[i][j];
+        }
+      }
+    }
+
+    // don't forget the clean everything up
+    for(size_t i=0; i<nDoubles; i++)
+    {
+      free(rawHessian[i]);
+      free(rawIdentity[i]);
+    }
+#endif
 }
 
 } // namespace Dune::Elasticity