Re: NDArray oct file gives error with 4 dimension

[Kai Torben Ohlhus]

On 2/7/20 6:52 AM, niconeuman wrote:
> Greetings,
> I'm writing an oct-file for performing an operation between a 2D NDArray and
> a 4D NDArray. Apparently trying to do
> 
> retMat(mu,nu) += D(ka,la)*gmnkl(mu,nu,ka,la);
> 
> does not compile. On the other hand, doing
> 
> retMat(mu,nu) += D(ka,ka)*gmnkl(mu,nu,ka);
> 
> compiles with no problem. Can anyone help me? 
> Best regards,
> Nicolas 
> 
> Below is the code for contractgD.cc,
> 
> 
> #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 ();
> 
>   const int lengthmu = gmnkl.dims()(0); //have to use it like this, not
> NDArray.dims(i);
>   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(1,1,1,1) << "\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(mu,nu,ka,la);
>         }
>       }
>     }
>   }
> }
> 
>   octave_value_list retval (nargout);
>   retval(0) = octave_value(retMat);
>   return retval;
> }
> 

Unfortunately, the ease of ND-array handling of the Octave language does
not translate 1:1 to it's C++ API.  You get the "comfort" of ND-indices
by calling

   gmnkl(dim_vector (1,1,1,1).as_array ())

instead of

   gmnkl(1,1,1,1)

Same for

   retMat(mu,nu) += D(ka,la)*
      gmnkl(dim_vector (mu,nu,ka,la).as_array ());

Full program attached.  Tested in Octave 5.2.0 with

   mkoctfile contractgD.cc
   D = eye (2);
   gmnkl = eye (2, 2, 2, 2)
   contractgD (D, gmnkl)

As further suggestion, you should check for dimension at the very
beginning of your program to avoid bad memory access crashes, caused by
improper input arrays.

Read more at


https://octave.org/doxygen/5.2/d1/db1/classdim__vector.html#a0189767745d163089c5145b36adf24d4

https://octave.org/doxygen/5.2/d3/dd5/classNDArray.html#a803cb3dcd25b55bcc1cd5333f6373cb1

HTH,
Kai

contractgD.cc
Description: Text Data