Re: single/double precision in C++

[Jaroslav Hajek]

On Tue, Nov 10, 2009 at 11:48 AM, Peter L. <address@hidden> wrote:
>> Templates could help you here. Octave 3.3.50+ provides a templatized function
>> octave_value_extract that extracts a specified matrix type from an 
>> octave_value.
>> If you want to target the 3.2.x series as well (which you probably
>> do), you can steal these from the development sources and define them
>> conditionally.
>>
>> Using octave_value_extract, you could write something like:
>> template <class MT1, class MT2>
>> static octave_value do_comp_dgt_fac (const octave_value_list& args,
>>                                                        const int a, const 
>> int M,
>>                                                        void
>> (*lib_func) (/* TODO */))
>> {
>>   const MT1 f = octave_value_extract<MT1> args(0);
>>   const MT2 gf = octave_value_extract<MT1> args(0);
>>
>>
>>        const int L = f.rows();
>>        const int W = f.columns();
>>        const int R = gf.rows()*gf.columns()/L;
>>
>>        const int N = L/a;
>>
>>        MT2 cout(M,N*W*R);
>>
>>    lib_func (...);
>>
>>   return cout;
>> }
>>
>> and then call it for the 4 cases
>> return do_comp_dgt_fac<Matrix, ComplexMatrix> (args, a, m, dgt_fac_r);
>> etc
>
> Thanks for the suggestion, this is just the sort of stuff I don't know
> about. This made me read up on trait classes.
>
>> The C prototypes pose another problem as these use different complex
>> number type than Octave (Octave uses std::complex and assumes the
>> layout conforms to Fortran and C99).
>>
>> A quick solution is to let these types be implied by the template, but
>> that is quite unsafe, as it will silently crash if you misspell a
>> function's name. A better solution is to provide a trait class
>> converting Octave's types to LTFAT's. A possibly yet better solution
>> is to alter ltfat.h in such a manner that it will use a user-defined
>> complex type if one is provided via a #define:
>>
>> #include <oct.h>
>> #define LTFAT_USER_COMPLEX Complex
>> #define LTFAT_USER_SCOMPLEX FloatComplex
>> #include "ltfat.h" // will now use Complex and FloatComplex for the 
>> prototypes
>>
>> this can be useful in general when interfacing with other software.
> This is actually what is already being done:
>
>   typedef fftw_complex  ltfat_complex;
>
> so the complex type is bit-compatible with Octave's.
>

Yes, that is the assumption (it's not actually guaranteed by C++, but
Octave uses the same assumption), but being compatible is not enough -
C++ still won't allow you to implicitly convert fftw_complex to
std::complex.

I checked out ltfat and hacked together a patch (not actually tested)
demonstrating my meaning:

Index: comp_dgt_fac.cc
===================================================================
--- comp_dgt_fac.cc     (revision 879)
+++ comp_dgt_fac.cc     (working copy)
@@ -2,57 +2,64 @@
 #include "config.h"
 #include "ltfat.h"

+template <class MT1, class MT2>
+static octave_value
+do_comp_dgt_fac (const octave_value_list& args,
+                 const int a, const int M,
+                 void (*lib_func) (typename MT1::element_type *,
+                                   typename MT2::element_type *,
+                                   const int, const int,
+                                   const int, const int,
+                                   const int,
+                                   typename MT2::element_type *))
+{
+ typedef typename MT1::element_type T1;
+ typedef typename MT2::element_type T2;
+
+ const MT1 f = octave_value_extract<MT1> args(0);
+ const MT2 gf = octave_value_extract<MT1> args(1);
+
+ const int L = f.rows();
+ const int W = f.columns();
+ const int R = gf.rows()*gf.columns()/L;
+
+ const int N = L/a;
+
+ MT2 cout(M,N*W*R);
+
+ // const_casts are necessary because LTFAT doesn't specify const on
input data.
+ lib_func (const_cast<T1 *>(f.fortran_vec ()), const_cast<T2
*>(gf.fortran_vec ()),
+           L, W, R, a, M, cout.fortran_vec ());
+
+ return cout;
+}
+
+template <class MT1, class MT2>
 DEFUN_DLD (comp_dgt_fac, args, ,
   "This function calls the C-library\n\
   c=comp_dgt_fac(f,gf,a,M);\n\
   Yeah.")
 {
+  octave_value_list retval;

   const bool f_is_complex  = args(0).is_complex_type();

   const int a = args(2).int_value();
   const int M = args(3).int_value();

-  if (f_is_complex)
-  {
-
-     const ComplexMatrix f = args(0).complex_matrix_value();
-     const ComplexMatrix gf = args(1).complex_matrix_value();
-
-     const int L = f.rows();
-     const int W = f.columns();
-     const int R = gf.rows()*gf.columns()/L;
-
-     const int N = L/a;
-
-     ComplexMatrix cout(M,N*W*R);
-
-     dgt_fac((ltfat_complex*)f.data(),(ltfat_complex*)gf.data(),
-     L, W, R, a, M,(ltfat_complex*)cout.data());
-
-     return octave_value (cout);
-
-
-  }
+  if (args(0).is_single_type ())
+    {
+      if (f_is_complex)
+        do_comp_dgt_fac (args, a, M, sdgt_fac);
+      else
+        do_comp_dgt_fac (args, a, M, sdgt_fac_r);
+    }
   else
-  {
+    {
+      if (f_is_complex)
+        do_comp_dgt_fac (args, a, M, dgt_fac);
+      else
+        do_comp_dgt_fac (args, a, M, dgt_fac_r);
+    }

-     const Matrix f = args(0).matrix_value();
-     const ComplexMatrix gf = args(1).complex_matrix_value();
-
-     const int L = f.rows();
-     const int W = f.columns();
-     const int R = gf.rows()*gf.columns()/L;
-
-     const int N = L/a;
-
-     ComplexMatrix cout(M,N*W*R);
-
-     dgt_fac_r((double*)f.data(),(ltfat_complex*)gf.data(),
-     L, W, R, a, M,(ltfat_complex*)cout.data());
-
-     return octave_value (cout);
-
-  }
-
 }
Index: config.h
===================================================================
--- config.h    (revision 879)
+++ config.h    (working copy)
@@ -5,13 +5,15 @@
 #define CONFIG_H 1

 #include "fftw3.h"
+
+#define LTFAT_USER_COMPLEX
+typedef Complex  ltfat_complex;
+typedef FloatComplex ltfat_scomplex;
+
 #include "ltfat.h"

 #define FFTW_OPTITYPE FFTW_ESTIMATE

-typedef fftw_complex  ltfat_complex;
-typedef fftwf_complex ltfat_scomplex;
-
 void* ltfat_malloc (size_t n)
 {
   return fftw_malloc(n);


in this manner, you write the bulk of the interfacing code just once,
in a template, then instantiate the template 4 times.
The trick in ltfat.h is used to map the correct Octave complex types
into the prototypes.
Note that there are still const_casts for the first two arguments,
because the prototypes in LTFAT do not specify const input pointers
(assuming the input data is unchanged). You can get rid of them (and
get even better prototype checking) by
augmenting the prototypes in ltfat_typeindependent.h by CONST
(#defined to const or nothing) specifiers where a const pointer is
only needed. Another option is to not specify const for f and gf, but
that way you'll get a useless copy of the data (if you query a
writable pointer for Octave's matrix types, the data is internally
unshared).

hth

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz