// GenTensors - A simple tensor library
// Copyright (C) 2012-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

// check gentensor ops
#include "genTensors.h"

int main(int argc,char**argv)
{
  
  const int tab1[3]={0,1,2};
  const int tab2[3]={2,0,1};
  SwapIndex<3,3> sw1(tab1);
  SwapIndex<3,3> sw2(tab2);
//  for (int i=0;i<SwapIndex<3,3>::Nele;++i) std::cout << i << " " << sw1[i] <<  " " << sw2[i] << std::endl; 
  const int tabndx[SwapIndex<3,3>::Nele] = {0,3,6,9,12,15,18,21,24,1,4,7,10,13,16,19,22,25,2,5,8,11,14,17,20,23,26};
  for (int i=0;i<SwapIndex<3,3>::Nele;++i) if (sw2[i]!=tabndx[i]) return 1;

  double mat0=7.0;
  const double mat1[]={1.0,5.0,3.0};

  const double mat2[]=
  {
    1.0,0.0,-1.0,
    0.0,2.0,0.0,
    -1.0,0.0,3.0
  };

  const double mat3[]=
  {
    1.0,0.0,-1.0,
    0.0,2.0,0.0,
    -1.0,0.0,3.0,

    0.0,0.0,0.0,
    0.0,2.0,0.0,
    0.0,0.0,0.0,

    3.0,0.0,-10.0,
    0.0,4.0,0.0,
    -10.0,0.0,5.0
  };

  const double mat4[]=
  {
    1.0,0.0,-1.0,
    0.0,2.0,0.0,
    -1.0,0.0,3.0,

    0.0,0.0,0.0,
    0.0,2.0,0.0,
    0.0,0.0,0.0,

    3.0,0.0,-10.0,
    0.0,4.0,0.0,
    -10.0,0.0,5.0,


    1.0,0.0,-1.0,
    0.0,2.0,0.0,
    -1.0,0.0,3.0,

    0.0,0.0,0.0,
    0.0,2.0,0.0,
    0.0,0.0,0.0,

    3.0,0.0,-10.0,
    0.0,4.0,0.0,
    -10.0,0.0,5.0,


    1.0,0.0,-1.0,
    0.0,2.0,0.0,
    -1.0,0.0,3.0,

    0.0,0.0,0.0,
    0.0,2.0,0.0,
    0.0,0.0,0.0,

    3.0,0.0,-10.0,
    0.0,4.0,0.0,
    -10.0,0.0,5.0
  };


  genTensor<0,double> t0(mat0);
  genTensor<1,double> t1(mat1);
  genTensor<2,double> t2(mat2);
  genTensor<3,double> t3(mat3);
  genTensor<4,double> t4(mat4);
  int tab[4]={3,2,1,0};
  SwapIndex<4> sw(tab);

  genTensor<0,double> t00;t00.ContractedProduct(t0,t0);
  genTensor<0,double> t11;t11.ContractedProduct(t1,t1);
  genTensor<0,double> t22;t22.ContractedProduct(t2,t2);
  genTensor<0,double> t33;t33.ContractedProduct(t3,t3.swap(),false);
  genTensor<0,double> t44;t44.ContractedProduct(t4,t4.swap(sw).swap(sw));

  genTensor<1,double> t10;t10.ContractedProduct(t1,t0);
  genTensor<1,double> t01;t01.ContractedProduct(t0,t1);
  genTensor<1,double> t12;t12.ContractedProduct(t1,t2);
  genTensor<1,double> t21;t21.ContractedProduct(t2,t1);
  genTensor<1,double> t23;t23.ContractedProduct(t2,t3);
  genTensor<1,double> t32;t32.ContractedProduct(t3,t2);
  genTensor<1,double> t34;t34.ContractedProduct(t3,t4);
  genTensor<1,double> t43;t43.ContractedProduct(t4,t3);

  genTensor<2,double> t20;t20.ContractedProduct(t2,t0);
  genTensor<2,double> t02;t02.ContractedProduct(t0,t2);
  genTensor<2,double> t31;t31.ContractedProduct(t3,t1);
  genTensor<2,double> t13;t13.ContractedProduct(t1,t3);
  genTensor<2,double> t42;t42.ContractedProduct(t4,t2);
  genTensor<2,double> t24;t24.ContractedProduct(t2,t4);

  genTensor<3,double> t30;t30.ContractedProduct(t3,t0);
  genTensor<3,double> t03;t03.ContractedProduct(t0,t3);
  genTensor<3,double> t41;t41.ContractedProduct(t4,t1);
  genTensor<3,double> t14;t14.ContractedProduct(t1,t4);

  genTensor<4,double> t40;t40.ContractedProduct(t4,t0);
  genTensor<4,double> t04;t04.ContractedProduct(t0,t4);

  std::cout << std::endl ;
  std::cout << std::endl ;
  std::cout.precision(17);
  std::cout << std::showpos;
  std::cout << std::scientific;

  std::cout <<  norm(t0) << " " <<  norm(t1) << " " <<  norm(t2) << " " <<  norm(t3) << " " <<  norm(t4) << std::endl;
  std::cout << norm(t00) << " " << norm(t11) << " " << norm(t22) << " " << norm(t33) << " " << norm(t44) << std::endl;
  std::cout << norm(t10) << " " << norm(t01) << " " << norm(t21) << " " << norm(t12) << " " << norm(t32) << " " << norm(t23) << " " << norm(t43) <<  " " << norm(t34) << std::endl;
  std::cout << norm(t20) << " " << norm(t02) << " " << norm(t31) << " " << norm(t13) << " " << norm(t42) << " " << norm(t24) << std::endl;
  std::cout << norm(t30) << " " << norm(t03) << " " << norm(t41) << " " << norm(t14) << std::endl;
  std::cout << norm(t40) << " " << norm(t04)  << std::endl;

  if (fabs(norm(t0) - 7.00000000000000000e+00)> 1e-12) return 1;
  if (fabs(norm(t1) - 5.91607978309961613e+00)> 1e-12) return 1;
  if (fabs(norm(t2) - 4.00000000000000000e+00)> 1e-12) return 1;
  if (fabs(norm(t3) - 1.64316767251549827e+01)> 1e-12) return 1;
  if (fabs(norm(t4) - 2.84604989415154144e+01)> 1e-12) return 1;

  if (fabs(norm(t00) - 49. ) >1e-12) return 1;
  if (fabs(norm(t11) - 35. ) >1e-12) return 1;
  if (fabs(norm(t22) - 16. ) >1e-12) return 1;
  if (fabs(norm(t33) - 65. ) >1e-12) return 1;
  if (fabs(norm(t44) - 1.  ) >1e-12) return 1;

  if (fabs(norm(t10) - 4.14125584816973102e+01)> 1e-12) return 1;
  if (fabs(norm(t01) - 4.14125584816973102e+01)> 1e-12) return 1;
  if (fabs(norm(t21) - 1.29614813968157208e+01)> 1e-12) return 1;
  if (fabs(norm(t12) - 1.29614813968157208e+01)> 1e-12) return 1;
  if (fabs(norm(t32) - 4.88671668914824622e+01)> 1e-12) return 1;
  if (fabs(norm(t23) - 3.76563407675254638e+01)> 1e-12) return 1;
  if (fabs(norm(t43) - 1.12583302491977022e+02)> 1e-12) return 1;
  if (fabs(norm(t34) - 9.13290753265355306e+01)> 1e-12) return 1;

  if (fabs(norm(t20) - 2.80000000000000000e+01)> 1e-12) return 1;
  if (fabs(norm(t02) - 2.80000000000000000e+01)> 1e-12) return 1;
  if (fabs(norm(t31) - 3.77094152699296075e+01)> 1e-12) return 1;
  if (fabs(norm(t13) - 5.40555270069583145e+01)> 1e-12) return 1;
  if (fabs(norm(t42) - 8.46404158779953093e+01)> 1e-12) return 1;
  if (fabs(norm(t24) - 3.28633534503099654e+01)> 1e-12) return 1;

  if (fabs(norm(t30) - 1.15021737076084889e+02)> 1e-12) return 1;
  if (fabs(norm(t03) - 1.15021737076084889e+02)> 1e-12) return 1;
  if (fabs(norm(t41) - 6.53146231712317302e+01)> 1e-12) return 1;
  if (fabs(norm(t14) - 1.47885090526394862e+02)> 1e-12) return 1;

  if (fabs(norm(t40) - 1.99223492590607890e+02)> 1e-12) return 1;
  if (fabs(norm(t04) - 1.99223492590607890e+02)> 1e-12) return 1;

  return 0; 
}