#include <octave/oct.h>

DEFUN_DLD (contractgD, args, nargout,
           "This function takes as arguments D and gmnkl, and calculates the"
           "matrix Jtemp which results from contracting gmnkl with D along"
           "dimensions 3 and 4")
{

  NDArray D = args(0).array_value ();
  NDArray gmnkl = args(1).array_value ();

  //have to use it like this, not NDArray.dims(i);
  const int lengthmu = gmnkl.dims()(0);
  const int lengthnu = gmnkl.dims()(1);
  const int lengthka = gmnkl.dims()(2);
  const int lengthla = gmnkl.dims()(3);

  octave_stdout << lengthmu << "\n";
  octave_stdout << lengthnu << "\n";
  octave_stdout << lengthka << "\n";
  octave_stdout << lengthla << "\n";
  octave_stdout << "gmnkl.ndims() = "  << "\n";
  octave_stdout << gmnkl.ndims() << "\n";
  octave_stdout << "gmnkl(1,1,1) = "  << "\n";
  octave_stdout << gmnkl(1,1,1) << "\n";
  octave_stdout << "gmnkl(1,1,1,1) = "  << "\n";
  octave_stdout << gmnkl(dim_vector (1,1,1,1).as_array ()) << "\n";
  dim_vector dvmn (lengthmu,lengthnu);

  NDArray retMat (dvmn);

  if (gmnkl.ndims() == 4){
    for(octave_idx_type mu = 0; mu < lengthmu; mu++){
      for(octave_idx_type nu = 0; nu < lengthnu; nu++){
        for(octave_idx_type ka = 0; ka < lengthka; ka++){
          for(octave_idx_type la = 0; la < lengthla; la++){
            retMat(mu,nu) += D(ka,la)*gmnkl(dim_vector (mu,nu,ka,la).as_array ());
          }
        }
      }
    }
  }

  octave_value_list retval (nargout);
  retval(0) = octave_value(retMat);
  return retval;
}