20                                 {
   23 
   24  boost::numeric::ublas::vector<double> input(3);
   25  input.clear();
   27  t2(0, 0) = 1;
   28  t2(0, 1) = 0.5;
   29  t2(1, 1) = 1;
   30 
   32  one(0, 0) = one(1, 1) = 1;
   33  one(1, 0) = one(0, 1) = 1;
   34 
   35  for (int ii = 0; ii != 2; ii++) {
   36    for (int jj = 0; jj != 2; jj++) {
   37      std::cout << t2(ii, jj) << " ";
   38    }
   39    std::cout << std::endl;
   40  }
   41  std::cout << std::endl;
   42 
   44 
   45  int tag = 1;
   46  int keep = 1;
   47 
   48  trace_on(tag, keep);
   50  
   51  a_t2(
I, 
J) <<= t2(
I, 
J);
 
   52  
   54  
   56  trace_off();
   57 
   58  for (int ii = 0; ii != 2; ii++) {
   59    for (int jj = 0; jj != 2; jj++) {
   60      std::cout << a_t2(ii, jj) << " ";
   61    }
   62    std::cout << std::endl;
   63  }
   64  std::cout << std::endl;
   65 
   66  std::cout << a_r << 
" " << 
r << std::endl << std::endl;
 
   67 
   68  boost::numeric::ublas::vector<double> output(3);
   70                                                 &output[2]);
   71 
   72  ::gradient(tag, 3, &input[0], &output[0]);
   73 
   74  std::cout << input << std::endl;
   75  std::cout << output << std::endl;
   76 
   77  for (int ii = 0; ii != 2; ii++) {
   78    for (int jj = 0; jj != 2; jj++) {
   79      std::cout << jac_t2(ii, jj) << " ";
   80      if (jac_t2(ii, jj) != 2) {
   81        std::cerr << "Wrong value should be 2\n" << std::endl;
   82        exit(-1);
   83      }
   84    }
   85    std::cout << std::endl;
   86  }
   87  std::cout << std::endl;
   88 
   89  boost::numeric::ublas::matrix<double> Hessian(3, 3);
   90  Hessian.clear();
   92  for (int nn = 0; nn != 3; nn++) {
   93    H[nn] = &Hessian(nn, 0);
 
   94  }
   95  ::hessian(tag, 3, &input[0], 
H);
 
   96  std::cout << Hessian << std::endl;
   97 
   99      &Hessian(0, 0), &Hessian(0, 1), &Hessian(0, 2), &Hessian(1, 0),
  100      &Hessian(1, 1), &Hessian(1, 2), &Hessian(2, 0), &Hessian(2, 1),
  101      &Hessian(2, 2));
  102  for (int ii = 0; ii != 2; ii++) {
  103    for (int jj = 0; jj != 2; jj++) {
  104      for (int II = 0; II != 2; II++) {
  105        for (int JJ = 0; JJ != 2; JJ++) {
  106          std::cout << "(" << ii << "," << jj << "," << II << "," << JJ << ") "
  107                    << t4(ii, jj, II, JJ) << std::endl;
  108        }
  109      }
  110    }
  111  }
  112  std::cout << std::endl;
  113 
  114  if (t4(0, 0, 0, 0) != 2 || t4(1, 1, 1, 1) != 2 || t4(1, 0, 1, 0) != 4 ||
  115      t4(0, 1, 0, 1) != 4 || t4(0, 1, 1, 0) != 4 || t4(1, 0, 0, 1) != 4) {
  116    std::cerr << "Wrong value\n" << std::endl;
  117    exit(-1);
  118  }
  119 
  120  return 0;
  121}
FTensor::Index< 'J', DIM1 > J
constexpr IntegrationType I