Re: complex arrays in octfiles

[John W. Eaton]

On 29-Jun-2011, c. wrote:

| Here is an example of how to build A ComplexNDArray:
| 
| #include <octave/oct.h>
| 
| DEFUN_DLD (pippo, args, , "test return statement")
| {
| 
|   dim_vector dims;
|   dims.resize (3);
|   dims(2) = 4;
|   dims(1) = 3;
|   dims(0) = 2;
| 
| 
|   ComplexNDArray a (dims, std::complex<double> (0.0, 0.0));
|   // or 
|   //Array<std::complex<double> > a (dims, std::complex<double> (0.0, 0.0));
|   for (octave_idx_type i = 0; i < dims(0); i++)
|     for (octave_idx_type j = 0; j < dims(1); j++)
|       for (octave_idx_type k = 0; k < dims(2); k++)
|         a (i, j, k) = std::complex<double> (i, j + k);
| 
|   return octave_value (a);
| 
| }

This should work, but using the indexing operator for assignment is
relatively slow because of Octave's copy-on-write semantics.  So if
you want better performance, you can write it this way:

  ComplexNDArray a (dims, Complex (0.0, 0.0));

  Complex *cptr = a.fortran_vec ();

  for (octave_idx_type k = 0; k < dims(2); k++)
    for (octave_idx_type j = 0; j < dims(1); j++)
      for (octave_idx_type i = 0; i < dims(0); i++)
        *cptr++ = Complex (i, j + k);

Note that the data is stored in column-major order, so I reversed the
loop order.

Complex is a typedef for std::complex<double> in the Octave sources.

jwe