Loading ...
Sorry, an error occurred while loading the content.

Re: [Boost-Users] ublas: questions from a newcomer

Expand Messages
  • Kresimir Fresl
    ... Atlas bindings work for me in the `real code (but I am biased ;o) I am also using SuperLU bindings for sparse matrices, but they are not in boost-sandbox
    Message 1 of 17 , Mar 25, 2003
    • 0 Attachment
      stuart macgregor wrote:

      > I am interested in ublas for linear algebra. In addition to the simple
      > operations supported by ublas I will need comprehensive access to lapack
      > for the real work.

      > I see in the wiki that some work has been done in this area. Is this
      > numeric/binding stuff working,

      Atlas bindings work for me in the `real' code (but I am biased ;o)
      I am also using SuperLU bindings for sparse matrices, but they are
      not in boost-sandbox yet.

      > and can it be readily extended to the full
      > set of lapack routines?

      Some time ago Toon wrote in boost-ml:

      ``Note that the bindings are not complete but the 'framework' is there
      and it's just a matter of time to complete them.''

      Bindings are mostly a free-time activity. Or, one can say that we
      try to generalize some more or less ad-hoc parts of our `regular' work.
      So `matter of time' probably means `some (considerable) time' :o(

      > Some kind of meta program to generate lapack interfaces from a parameter
      > spec would be nice (ruby/perl?

      Interesting idea.

      > or just a description of what needs to be
      > done) Then I could happily bring up the interfaces I need.

      As I said, bindings are free-time activity, so any help is
      welcomed.

      Regarding (c)lapack part of atlas bindings, there are interfaces
      for all (interesting) functions which atlas provides (gesv, getrf,
      getrs, getri, posv, potrf, potrs, potri; only lauum and trtri are
      missing), see

      http://math-atlas.sourceforge.net/psdoc/lapackqref.ps.gz

      Lapack bindings are Toon's realm, but I hope that my understanding
      is good enough to (try to) explain what needs to be done.

      As an example, I'll use function `getrf' (which is already there).
      As you know, there are four versions of each lapack function
      -- float, double, complex, complex double. AFAICS Toon decided
      to write bindings only for double and complex double versions,
      but other two can be easily added following the same pattern.
      I'll use only dgetrf here.

      Include files are in

      boost/numeric/bindings/lapack

      and the `traits.cpp' file is in

      libs/numeric/bindings/lapack

      (1) `lapack.h' contains declarations of lapack functions:

      As there are several conventions for calling Fortran from C/C++
      (some compilers require trailing underscore in function names,
      some don't, see excelent page

      http://www.aei.mpg.de/~jthorn/c2f.html

      and, in particular, section `Name mangling'), first a macro
      LAPACK_DGETRF is defined, depending on naming convention:
      =============================================
      #if defined (BIND_FORTRAN_LOWERCASE_UNDERSCORE)
      #define LAPACK_DGETRF dgetrf_
      #elif defined(BIND_FORTRAN_LOWERCASE)
      #define LAPACK_DGETRF dgetrf
      #endif
      =============================================
      (macros BIND_FORTRAN_... are defined in
      boost/numeric/bindings/traits/fortran.h).

      Now lapack function can be declared:
      =============================================
      void LAPACK_DGETRF(const int* m, const int* n,
      double* a, const int* lda, int* ipiv, int*
      info) ;
      =============================================
      All parameters are pointers.

      (2) `traits.hpp' contains class template `traits'.

      For each lapack there is a typedef of a function pointer, which
      corresponds to the function declaration in `lapack.h':
      =============================================
      typedef void (*getrf_type)(const int* m, const int* n,
      bind_type* a, const int* lda, int* ipiv, int* info) ;
      =============================================
      but instead of `double*' there is a `bind_type*', which depends
      on the template parameter `T' of the class `traits'. `bind_type' is
      defined using `value_traits' class (which is defined and specialized
      in boost/numeric/bindings/traits/value_traits.hpp).

      Class `traits' also contains a (static) function pointer:
      =============================================
      static getrf_type getrf ;
      =============================================

      (3) in `traits.cpp' the function pointer is initialized with
      (the address of) the lapack function:
      =============================================
      template<> traits<double>::getrf_type
      traits<double>::getrf = LAPACK_DGETRF ;
      =============================================

      `traits.cpp' must be separately compiled and linked with your
      program.

      But, in fact, here we have a problem -- see comment [4] below.

      (4) `lapack.hpp' contains the C++ interfaces (or bindings).

      All bindings (blas, lapack, atlas) are written in terms of traits
      classes which are interfaces between these functions and
      vector/matrix classes. Namely, idea is to make bindings
      independent of the concrete vector/matrix libraries as much
      as possible.

      Lapack functions require the starting address of the matrix,
      number of rows and columns and leading dimension.
      So, matrix_traits class, that is, its specializations provide this
      functions and some typedefs. For example, size1() (i.e. number
      of rows) in matrix_traits specialization for ublas::matrix traits is

      static int size1 (matrix_type& m) { return m.size1(); }

      and in specialization for Array2D from Roldan Pozo's TNT

      static int size1 (matrix_type& m) { return m.dim1(); }

      If you use ublas, you don't need to worry about matrix_traits,
      because they are already written -- you just have to include
      boost/numeric/bindings/traits/ublas_matrix.hpp
      for general matrices and matrix proxies and
      boost/numeric/bindings/traits/ublas_symmetric.hpp
      for symmetric matrices and symmetric adaptors.
      (Traits classes for Roldan Pozo's TNT are also
      there, but I didn't test them thoroughly.)
      Note that bindings don't know and need not know anything
      about traits specializations --
      boost/numeric/bindings/traits/ublas_matrix.hpp etc.
      must be included in your program, not in bindings
      headers.

      To simplify the use of traits classes, Toon recently introduced
      some free accessor function, so instead of writting

      matrix_traits<matrix_type>::size1 (m)

      you can write simply

      matrix_size1 (m)

      Now, let's see how all this works in the interface for getrf
      (explanations are given as comments; this is Toon's code and I'll
      also comment some implementation details -- longer comments,
      not really important for functionality, are given after the code):
      ==============================================
      template < typename matrix_type, typename IpivIterator >
      int getrf (matrix_type& a, IpivIterator begin_ipiv, IpivIterator end_ipiv)
      // instead of 6 parameters of lapack's (sdcz)getrf, we have only
      // 3 parameters: matrix to factor and iterator range for pivot vector
      // (but see also comment [1] below)
      {
      using namespace boost::numeric::bindings::traits ;
      // why? all names below are fully namespace qualified
      typedef typename matrix_type::value_type value_type ;
      // see comment [2] below
      typedef typename
      boost::numeric::bindings::traits::value_traits< value_type
      >::value_type
      bind_type ;
      // as in `traits.hpp' in (2) above

      // int m = boost::numeric::bindings::traits::size1( a ) ;
      // int n = boost::numeric::bindings::traits::size2( a ) ;
      // number of matrix rows and columns;
      // free accessor functions are used, but in new version they are
      // matrix_size1() and matrix_size2(), so this will not compile :o(
      // correct calls are:
      int m = boost::numeric::bindings::traits::matrix_size1( a ) ;
      int n = boost::numeric::bindings::traits::matrix_size2( a ) ;

      value_type *a_ptr =
      boost::numeric::bindings::traits::matrix_storage( a ) ;
      const int lda =
      boost::numeric::bindings::traits::leading_dimension( a ) ;
      // starting address of matrix storage and its leading dimension

      assert( std::distance( begin_ipiv, end_ipiv ) >= std::min( m, n ) );
      // pivot vector must be large enough

      int info = 0; // return value; in short: is matrix singular?
      int& ipiv_ptr = *begin_ipiv ;
      // reference to first element of the pivot vector,
      // &ipiv_ptr is then its starting address;
      // but see also comment [3]

      traits<value_type>::getrf( &m, &n, (bind_type*)a_ptr, &lda,
      &ipiv_ptr, &info );
      // finally, call lapack's function

      return info;
      }
      ==============================================

      Comments:

      [1] I think that this interface [i.e. func (matrix&, begin_iter,
      end_iter)] is
      mixing of two idioms -- why STL's idiom of iterator ranges [begin, end)
      here, and not for other vectors in e.g. blas 1&2 bindings?
      Interface of getrf in atlas bindings (atlas/clapack.hpp) is simply:

      int getrf (MatrA& a, IVec& ipiv);

      [2] It is true that both ublas::matrix<> and TNT::Array2D<> define
      value_type, but (in ideal world ;o) bindings should work for other
      matrix types which can have val_t or something else.
      And that's why we have matrix_traits; so, instead of

      typedef typename matrix_type::value_type value_type ;

      it should be

      typedef typename matrix_traits<matrix_type>::value_type value_type;

      [3] See the interface of getrf in atlas bindings (atlas/clapack.hpp)
      for IMHO more elegant way to do this (yes, I'm probably biased ;o)

      [4] The biggest problem with blas & lapack bindings is that
      although they work with g++ (at least with 3.2, I didn't try later
      versions), they cannot be compiled with Comeau C++ 4.3
      --- they are not standard conforming.

      At first sight, class `traits<>', defined in `lapack/traits.hpp'
      (and, similarly, class `blas<>' in `blas/blaspp.hpp') is very elegant
      solution for the function overloading -- one class for all four
      versions of blas/lapack functions (float, double, complex,
      complex double). But function pointers defined in `traits<>'
      and `blas<>' have "C++" linkage and blas/lapack functions
      have (extern) "C" linkage; therefore, this may, but need not,
      work:

      template<> traits<double>::getrf_type
      traits<double>::getrf = LAPACK_DGETRF ;

      Here is a quote from the recent post by Greg Comeau in
      comp.lang.c++.moderated (`Re: use of extern "C"'):
      ===============================================
      One reason why this has probably been acceptable in the code
      you've written is that your compilers are allowing implicit
      conversions between extern "C" and extern "C++" entities,
      actually pointers to some such entities, since often the
      same internal/underlying conventions (calling, representations, etc)
      are used for each, but from a strict standpoint, the standard
      does not support that.
      ===============================================

      Steve Clamage gave an example in the same thread:
      ===============================================
      The extern "C" linkage declaration says that C
      linkage rules are in effect. If C and C++ functions have different
      calling sequences, the code generation will be different depending on
      the declared linkage. The same goes for pointer-to function, since
      the linkage is part of the pointer type. For example:
      int f(int); // C++ linkage
      extern "C" int g( int(*)(int) );
      ...
      g(f); // error
      Function g is declared to take a pointer to a C-linkage function, and
      passing a pointer to a C++-linkage function is not valid. Not all
      compilers complain about the error, however.
      ===============================================

      Our example is inverted -- traits<double>::getrf_type is a pointer
      to a C++-linkage function and LAPACK_DGETRF is (a pointer)
      to a C-linkage function -- but the same rule obviously applies.

      I am afraid that blas and lapack bindings need complete rewriting
      as in atlas bindings which use less elegant and more conservative
      (but also standard conforming) approach -- simple function overloading:
      ===============================================
      int gesv (CBLAS_ORDER const Order, int const N, int const NRHS,
      float* A, int const lda, int* ipiv, float* B, int const ldb)
      {
      return clapack_sgesv (Order, N, NRHS, A, lda, ipiv, B, ldb);
      }
      int gesv (CBLAS_ORDER const Order, int const N, int const NRHS,
      double* A, int const lda, int* ipiv, double* B, int const
      ldb)
      {
      return clapack_dgesv (Order, N, NRHS, A, lda, ipiv, B, ldb);
      }
      ===============================================

      > Is there any way to generate vectore and matrices which do not own
      the data
      > buffer which they provide access to - 'view' or 'sub' objects in
      other lin
      > alg packages? This would be very valuable to me for interfacing
      existing
      > code, without duplicating storage and doing unnecessary copies. I
      found it
      > difficult to determine from the current documentation.

      Joerg answered this one ...

      > I have one worry: Ublas, and most of the boost modules, seems to me
      not for
      > the faint hearted. They are clearly powerful, but seem to expect a lot
      > from a potential user.

      Well, some boost libraries are certainly frightening, but I don't think
      that ublas is among them -- its interface seems rather intuitive.

      > Are these facilities aimed at library writers rather than users whoes
      skill
      > focus may lie in other domains?

      Some very simple, self-explanatory ublas examples can be found in
      subdirectory

      libs/numeric/ublas/doc/samples

      of the boost distribution.

      > I expect that a user level doc may eventually make ublas/lapack much
      more
      > approachable by a non computer science mathematician or engineer.

      I agree. ublas documentation is good and detailed reference
      (when you know what you are looking for ;o), but a tutorial with
      few well commented examples is probably missing.

      Regarding bindings, I can only hope that I'll find some time in the near
      future to write some documentation (at least for traits classes and atlas
      bindings). In the meantime, there are some examples in

      libs/numeric/bindings/atlas/

      (in the boost-sandbox); I tried to write meaningful comments ;o).
      Begin with ublas_gesv2.cc, then ublas_gesv.cc and
      ublas_posv.cc.

      Also, don't hesitate to ask further questions.

      (BTW, few days ago Toon updated bindings in the sandbox,
      so it now contains the latest version -- that is, the note on
      Wiki page is obsolete.)

      > Forgive me if the above sounds critical and ungracious for free
      software,

      No, not at all ... I had (and have) the same problem with some other
      libraries ;o)

      > but I want to use boost and I am worried that it may be too demanding
      for
      > its benefites to be justified for use by the group of people I work
      with.

      AFAICS, the main problem is the lack of documentation, in particular
      tutorials and introductory examples. Also, as Joerg suggested,
      ``compiler diagnostics clearly are a pain sometimes'', but this is
      a problem with all templated code.

      Regards,

      fres
    • Toon Knapen
      ... I don t believe it s feasible. When implementing the bindings I found it to be very subtle. For instance, lapack functions might have different signatures
      Message 2 of 17 , Apr 1, 2003
      • 0 Attachment
        On Tuesday 25 March 2003 17:26, Kresimir Fresl wrote:
        > Some time ago Toon wrote in boost-ml:
        >
        > ``Note that the bindings are not complete but the 'framework' is there
        > and it's just a matter of time to complete them.''
        >
        > Bindings are mostly a free-time activity. Or, one can say that we
        > try to generalize some more or less ad-hoc parts of our `regular' work.
        > So `matter of time' probably means `some (considerable) time' :o(
        >
        > > Some kind of meta program to generate lapack interfaces from a parameter
        > > spec would be nice (ruby/perl?
        >
        > Interesting idea.
        I don't believe it's feasible. When implementing the bindings I found it to be very subtle.
        For instance, lapack functions might have different signatures depending on the
        value_type.

        > > or just a description of what needs to be
        > > done) Then I could happily bring up the interfaces I need.
        >
        > As I said, bindings are free-time activity, so any help is
        > welcomed.
        All help is welcome. As Kresimir will mention later in this email,
        the current implementation is not standard conforming. So if you care
        to help, we can launch the discussion on how we best implement
        the bindings and start working on it together.

        > Lapack bindings are Toon's realm, but I hope that my understanding
        > is good enough to (try to) explain what needs to be done.
        Thanks for responding to this question while I was enjoying snowboarding.
        You practically write the whole documentation in a single email ;-)

        > [1] I think that this interface [i.e. func (matrix&, begin_iter,
        > end_iter)] is
        > mixing of two idioms -- why STL's idiom of iterator ranges [begin, end)
        > here, and not for other vectors in e.g. blas 1&2 bindings?
        > Interface of getrf in atlas bindings (atlas/clapack.hpp) is simply:
        >
        > int getrf (MatrA& a, IVec& ipiv);
        Agreed !


        > [2] It is true that both ublas::matrix<> and TNT::Array2D<> define
        > value_type, but (in ideal world ;o) bindings should work for other
        > matrix types which can have val_t or something else.
        > And that's why we have matrix_traits; so, instead of
        >
        > typedef typename matrix_type::value_type value_type ;
        >
        > it should be
        >
        > typedef typename matrix_traits<matrix_type>::value_type value_type;
        Agreed !


        > [3] See the interface of getrf in atlas bindings (atlas/clapack.hpp)
        > for IMHO more elegant way to do this (yes, I'm probably biased ;o)
        I admit that it's a bit clumsy while trying to avoid code duplication.

        > > I have one worry: Ublas, and most of the boost modules, seems to me
        > not for
        > > the faint hearted. They are clearly powerful, but seem to expect a lot
        > > from a potential user.
        >
        > Well, some boost libraries are certainly frightening, but I don't think
        > that ublas is among them -- its interface seems rather intuitive.
        Joerg,
        You know what I appreciate even more than the documentation: a small
        comment on all classes and their member function.


        > > but I want to use boost and I am worried that it may be too demanding
        > for
        > > its benefites to be justified for use by the group of people I work
        > with.
        >
        powerfull things might be difficult to understand at first but ... hey they are powerfull.
        You might achieve the same goal writing more code but that takes time.
        You can just as well walk your whole life but take the time to learn how to drive
        a car and you will get there much faster (but you first need to invest the time to
        learn to drive a car). I understand your concern but I'm convinced that investing
        in a well established library is worth it.

        toon
      • Matthias Troyer
        Hi, After some discussions in my group when we had Jeremy Siek visit this weekend we agreed that we would like to see generic C++ lapack interfaces that work
        Message 3 of 17 , Apr 1, 2003
        • 0 Attachment
          Hi,

          After some discussions in my group when we had Jeremy Siek visit this
          weekend we agreed that we would like to see generic C++ lapack
          interfaces that work with ublas and other matrix libraries submitted to
          boost as a library. One could start with the high level driver
          functions only and add more specialized functions later. I thus would
          like to ask what the current status of the ublas lapack bindings is and
          if somebody is working on a generic version for boost, or if my
          research group could help here?

          Now I just saw this e-mail ...

          On Tuesday, April 1, 2003, at 01:03 PM, Toon Knapen wrote:

          > On Tuesday 25 March 2003 17:26, Kresimir Fresl wrote:
          >>
          >>> Some kind of meta program to generate lapack interfaces from a
          >>> parameter
          >>> spec would be nice (ruby/perl?
          >>
          >> Interesting idea.
          > I don't believe it's feasible. When implementing the bindings I found
          > it to be very subtle.
          > For instance, lapack functions might have different signatures
          > depending on the
          > value_type.

          It might still be possible, but only after a number of interfaces have
          been coded manually to see how much could be automated.

          >>> or just a description of what needs to be
          >>> done) Then I could happily bring up the interfaces I need.
          >>
          >> As I said, bindings are free-time activity, so any help is
          >> welcomed.
          > All help is welcome. As Kresimir will mention later in this email,
          > the current implementation is not standard conforming. So if you care
          > to help, we can launch the discussion on how we best implement
          > the bindings and start working on it together.

          We are certainly interested and could help with discussions or some
          coding.

          All the best,

          Matthias
        • Kresimir Fresl
          Hi Matthias, hi all, ... This is the main idea behind the vector and matrix traits classes in the bindings library -- for example, currently blas level 1
          Message 4 of 17 , Apr 1, 2003
          • 0 Attachment
            Hi Matthias, hi all,

            > After some discussions in my group when we had Jeremy Siek visit this
            > weekend we agreed that we would like to see generic C++ lapack
            > interfaces that work with ublas and other matrix libraries submitted to
            > boost as a library.

            This is the main idea behind the vector and matrix traits classes
            in the bindings library -- for example, currently blas level 1 functions
            can be applied to ublas::vector, std::vector, std::valarray, boost::array,
            TNT::Array1D and plain C array; e.g.:
            ================================================
            #include <boost/numeric/bindings/traits/std_vector.hpp>
            #include <boost/numeric/bindings/traits/std_valarray.hpp>
            #include <boost/numeric/bindings/traits/boost_array.hpp>
            #include <boost/numeric/bindings/traits/c_array.hpp>
            #include <boost/numeric/bindings/atlas/cblas1.hpp>

            // ...

            std::vector<double> sv (n);
            std::valarray<double> va (n);
            // ...
            cout << atlas::dot (sv, va) << endl;

            boost::array<double, 10> ba;
            atlas::set (0.1, ba);
            double ca[10];
            atlas::set (1., ca);
            atlas::axpy (0.1, ba, ca);
            ================================================
            (for complete examples see

            libs/numeric/bindings/atlas/other.cc
            libs/numeric/bindings/atlas/ublas_vct.cc

            in boost-sandbox or in `bindings_2003_01_30.zip' in the files
            section of this ml).

            Similarly, (c)lapack functions included in Atlas can be currently
            applied to ublas::matrix, TNT::Array2D and TNT::FortranArray2D.

            Above `currently' means that traits classes are already written
            for above-mentioned `vector' and `matrix' classes, see directory

            boost/numeric/bindings/traits

            Of course, the framework can be extended to other `vector' and
            `matrix' classes by writting appropriate traits classes.

            > One could start with the high level driver
            > functions only and add more specialized functions later. I thus would
            > like to ask what the current status of the ublas lapack bindings is and
            > if somebody is working on a generic version for boost, or if my
            > research group could help here?

            See

            http://groups.yahoo.com/group/ublas-dev/message/593

            But you already know about it ;o) :

            > Now I just saw this e-mail ...

            [...]

            >>>>Some kind of meta program to generate lapack interfaces from a
            >>>>parameter
            >>>>spec would be nice (ruby/perl?
            >>>
            >>>Interesting idea.
            >>
            >>I don't believe it's feasible. When implementing the bindings I found
            >>it to be very subtle.
            >>For instance, lapack functions might have different signatures
            >>depending on the
            >>value_type.
            >
            > It might still be possible, but only after a number of interfaces have
            > been coded manually to see how much could be automated.

            As I said:

            >>>Interesting idea.

            [...]

            >>All help is welcome. As Kresimir will mention later in this email,
            >>the current implementation is not standard conforming. So if you care
            >>to help, we can launch the discussion on how we best implement
            >>the bindings and start working on it together.

            > We are certainly interested and could help with discussions or some
            > coding.

            Goood.

            I think that the most important thing is to unify the interface so that
            individual functions can be independently added, that is, so that one
            can write function interface when he needs it and add it painlessly
            to the library.

            By `interface' I mean e.g. whether `getrf' should be

            template <typename MatrA, typename IpivIterator >
            int getrf (MatrA& a, IpivIterator begin_ipiv, IpivIterator end_ipiv);

            as in `boost/numeric/bindings/lapack/lapack.hpp', or

            template <typename MatrA, typename IVec>
            int getrf (MatrA& a, IVec& ipiv);

            as in `boost/numeric/bindings/atlas/clapack.hpp'. Also, although
            probably less important, should (forwarding) functions with more
            descriptive names be added, e.g. (from `atlas/clapack.hpp'):

            template <typename MatrA, typename IVec>
            inline
            int lu_factor (MatrA& a, IVec& ipiv) { return getrf (a, ipiv); }

            Furthermore, for `higher level' driver functions, should pivot vectors
            and/or work spaces be `hidden'? `atlas/clapack.hpp' provides
            both options:

            template <typename MatrA, typename MatrB, typename IVec>
            int gesv (MatrA& a, IVec& ipiv, MatrB& b);

            template <typename MatrA, typename MatrB>
            int gesv (MatrA& a, MatrB& b);
            // with:
            // int *ipiv = new int[n];


            Regards,

            fres
          • Toon Knapen
            Kresimir, Mathias, all, Already from the experience of implementing the current bindings I learned a lot of things through cooperation with Kresimir and I
            Message 5 of 17 , Apr 4, 2003
            • 0 Attachment
              Kresimir, Mathias, all,

              Already from the experience of implementing the current bindings
              I learned a lot of things through cooperation with Kresimir and I think
              we have knowledge about how it should be implemented
              (as Kresimir recently posted, the current blas and lapack bindings are
              not standard conforming and can't be without being reworked totally ;
              initially we also changed the directory organisation a few times, ...)

              And it would be great to do this job all together. So ideally we should
              revisit all components of the bindings library, discuss it, document it
              and implement it. We should start from the lowest level (which are the
              traits).

              So to put my money where my mouth is, I started a small doc
              in boost-sandbox/libs/numeric/bindings/traits/doc/index.html
              about the traits library.
              (here's a link : http://cvs.sourceforge.net/cgi-bin/viewcvs.cgi/boost-sandbox/boost-sandbox/libs/numeric/bindings/traits/doc/index.html?rev=1.1&content-type=text/vnd.viewcvs-markup )

              I started of with commenting about the fortran.h file. But in the mean time
              I also looked at the link to the webpage provided by Kresimir recently and
              I started thinking. Should we still have our own fortran.h our would we reuse
              the fortran.h provided on this page ?

              I like the system he implements but I fail to see how you can use
              a macro to change the macro-argument (which should be in lowercase)
              to uppercase.

              toon
            • Matthias Troyer
              ... That s a good idea - it also helped me in other projects to discuss interface issues with others and am happy to participate. ... I just checked it out but
              Message 6 of 17 , Apr 4, 2003
              • 0 Attachment
                On Friday, April 4, 2003, at 06:51 PM, Toon Knapen wrote:
                > And it would be great to do this job all together. So ideally we should
                > revisit all components of the bindings library, discuss it, document it
                > and implement it. We should start from the lowest level (which are the
                > traits).

                That's a good idea - it also helped me in other projects to discuss
                interface issues with others and am happy to participate.

                > So to put my money where my mouth is, I started a small doc
                > in boost-sandbox/libs/numeric/bindings/traits/doc/index.html
                > about the traits library.
                > (here's a link :
                > http://cvs.sourceforge.net/cgi-bin/viewcvs.cgi/boost-sandbox/boost-
                > sandbox/libs/numeric/bindings/traits/doc/index.html?rev=1.1&content-
                > type=text/vnd.viewcvs-markup )

                I just checked it out but did not find anything about the fortran.h file
                >
                > I started of with commenting about the fortran.h file. But in the mean
                > time
                > I also looked at the link to the webpage provided by Kresimir recently
                > and
                > I started thinking. Should we still have our own fortran.h our would
                > we reuse
                > the fortran.h provided on this page ?


                Is that the link from http://www.aei.mpg.de/~jthorn/c2f.html? I took a
                look at that file but think that we definitively need our own fortran.h
                and can't just build on this incomplete one.

                We need typedefs for the following Fortran types (where they exist)

                REAL
                REAL*4
                REAL*8
                REAL*16
                COMPLEX
                COMPLEX*8
                COMPLEX*16
                COMPLEX*32
                INTEGER
                INTEGER*4
                INTEGER*8
                LOGICAL

                Next we need typedefs for the single and double precision floating
                point types used in BLAS and LAPACK. Note that BLAS single precision is
                usually 32 bit which on all platforms so far seems to be float in
                C/C++. On the (vector) Crays however single precision is 64 bit and
                thus double. The BLAS double precision type on the Cray is 128 bit and
                long double in C++, while it is 64 bit and double on most other
                platforms.

                Next, I do not understand why he has these macros (again qualified by
                my looking at the right fortran.h file)

                #define FORTRAN_INTEGER_IS_INT TRUE
                #define FORTRAN_LOGICAL_TRUE 1
                #define FORTRAN_LOGICAL_FALSE 0

                The last two can be constants in C++, or do we also want to provide C
                bindings in addition to the C++ ones?


                > I like the system he implements but I fail to see how you can use
                > a macro to change the macro-argument (which should be in lowercase)
                > to uppercase.

                The only way I see is the one taken by Atlas, e.g. to define separate
                macros for the name. You could thus define e.g.

                #ifdef FORTRAN_UPPER_CASE
                #define daxpy_ DAXPY_
                ...
                #endif

                A macro clearly cannot do that.

                What is the next step? Will we first discuss a bit more and then code
                the fortran.h file, shall we discuss more topics before coding, or is
                somebody already working on an implementation?

                Best regards,

                Matthias
              • Toon Knapen
                ... in boost-sandbox/libs/numeric/bindings/traits/doc/index.html ? I doublechecked and it s there ?! ... yep , that s him. ... Where do we need types with
                Message 7 of 17 , Apr 7, 2003
                • 0 Attachment
                  On Friday 04 April 2003 22:16, Matthias Troyer wrote:

                  > I just checked it out but did not find anything about the fortran.h file
                  in boost-sandbox/libs/numeric/bindings/traits/doc/index.html ?
                  I doublechecked and it's there ?!

                  > Is that the link from http://www.aei.mpg.de/~jthorn/c2f.html?
                  yep , that's him.

                  > I took a
                  > look at that file but think that we definitively need our own fortran.h
                  > and can't just build on this incomplete one.
                  >
                  > We need typedefs for the following Fortran types (where they exist)
                  >
                  > REAL
                  > REAL*4
                  > REAL*8
                  > REAL*16
                  > COMPLEX
                  > COMPLEX*8
                  > COMPLEX*16
                  > COMPLEX*32
                  > INTEGER
                  > INTEGER*4
                  > INTEGER*8
                  > LOGICAL

                  Where do we need types with specific size for. BLAS and LAPACK
                  for instance are formulated purely in terms of REAL and COMPLEX,
                  using the 'natural' size.
                  A question I have though is if we'll bind std::complex to the fortran
                  COMPLEX or would users be able to use other complex types in
                  C++ (remember that C++ does not guarantee the two real's in the
                  complex to be consecutive in memory and stored as real() and imag()).

                  > Next we need typedefs for the single and double precision floating
                  > point types used in BLAS and LAPACK. Note that BLAS single precision is
                  > usually 32 bit which on all platforms so far seems to be float in
                  > C/C++. On the (vector) Crays however single precision is 64 bit and
                  > thus double. The BLAS double precision type on the Cray is 128 bit and
                  > long double in C++, while it is 64 bit and double on most other
                  > platforms.
                  OK.

                  >
                  > Next, I do not understand why he has these macros (again qualified by
                  > my looking at the right fortran.h file)
                  >
                  > #define FORTRAN_INTEGER_IS_INT TRUE
                  > #define FORTRAN_LOGICAL_TRUE 1
                  > #define FORTRAN_LOGICAL_FALSE 0
                  >
                  > The last two can be constants in C++, or do we also want to provide C
                  > bindings in addition to the C++ ones?
                  Well, since I propose to always bind via C, we can just as well provide
                  bindings for C. But I think the C++ layer is another layer on top of the
                  C bindings.

                  >
                  > > I like the system he implements but I fail to see how you can use
                  > > a macro to change the macro-argument (which should be in lowercase)
                  > > to uppercase.
                  >
                  > The only way I see is the one taken by Atlas, e.g. to define separate
                  > macros for the name. You could thus define e.g.
                  >
                  > #ifdef FORTRAN_UPPER_CASE
                  > #define daxpy_ DAXPY_
                  > ...
                  > #endif
                  >
                  > A macro clearly cannot do that.
                  Is there really no way to do that? Maybe we should ask the
                  guys from the preprocessor library because it's troublesome to
                  do this ifdef for adding underscores _and_ capitalising.


                  > What is the next step? Will we first discuss a bit more and then code
                  > the fortran.h file, shall we discuss more topics before coding, or is
                  > somebody already working on an implementation?
                  Well, I can start adapting the existing fortran.h one I know the answer
                  to questions above.

                  t
                • Matthias Troyer
                  ... I found oth the index.html and the fortran.h file but did not see fortran.h mentioned in the documentation ... If we just want to bind to blas/lapack/atlas
                  Message 8 of 17 , Apr 7, 2003
                  • 0 Attachment
                    On Monday, April 7, 2003, at 04:11 PM, Toon Knapen wrote:

                    > On Friday 04 April 2003 22:16, Matthias Troyer wrote:
                    >
                    >> I just checked it out but did not find anything about the fortran.h
                    >> file
                    > in boost-sandbox/libs/numeric/bindings/traits/doc/index.html ?
                    > I doublechecked and it's there ?!

                    I found oth the index.html and the fortran.h file but did not see
                    fortran.h mentioned in the documentation

                    >
                    >> I took a
                    >> look at that file but think that we definitively need our own
                    >> fortran.h
                    >> and can't just build on this incomplete one.
                    >>
                    >> We need typedefs for the following Fortran types (where they exist)
                    >>
                    >> REAL
                    >> REAL*4
                    >> REAL*8
                    >> REAL*16
                    >> COMPLEX
                    >> COMPLEX*8
                    >> COMPLEX*16
                    >> COMPLEX*32
                    >> INTEGER
                    >> INTEGER*4
                    >> INTEGER*8
                    >> LOGICAL
                    >
                    > Where do we need types with specific size for. BLAS and LAPACK
                    > for instance are formulated purely in terms of REAL and COMPLEX,
                    > using the 'natural' size.

                    If we just want to bind to blas/lapack/atlas for now then we can of
                    course limit ourselves to two types for the platform specific single
                    and double precision real and complex types. However we could make it
                    more general and provide complete mappings to help designers of
                    interfaces to other libraries.

                    > A question I have though is if we'll bind std::complex to the fortran
                    > COMPLEX or would users be able to use other complex types in
                    > C++ (remember that C++ does not guarantee the two real's in the
                    > complex to be consecutive in memory and stored as real() and imag()).

                    That's a valid point and is being addressed by the standards committee
                    as far as I know. Since all implementations that I know of currently
                    store real and imaginary parts consecutively I would suggest a to e.g.
                    typedef std::complex<fortran_real> fortran_complex, so that we could
                    later (if ever needed on some obscure platforms) provide alternative
                    implementations.

                    >> Next, I do not understand why he has these macros (again qualified by
                    >> my looking at the right fortran.h file)
                    >>
                    >> #define FORTRAN_INTEGER_IS_INT TRUE
                    >> #define FORTRAN_LOGICAL_TRUE 1
                    >> #define FORTRAN_LOGICAL_FALSE 0
                    >>
                    >> The last two can be constants in C++, or do we also want to provide C
                    >> bindings in addition to the C++ ones?
                    > Well, since I propose to always bind via C, we can just as well provide
                    > bindings for C. But I think the C++ layer is another layer on top of
                    > the
                    > C bindings.

                    How would you do the bindings for complex in C? Doesn't this just make
                    things harder? I guess we need separate C and C++ interfaces to the
                    Fortran bindings (maybe in the same source using preprocessor tricks
                    and macros like COMPLEX_TYPE), and then build higher level C++
                    interfaces on top of them as you suggested.


                    >>> I like the system he implements but I fail to see how you can use
                    >>> a macro to change the macro-argument (which should be in lowercase)
                    >>> to uppercase.
                    >>
                    >> The only way I see is the one taken by Atlas, e.g. to define separate
                    >> macros for the name. You could thus define e.g.
                    >>
                    >> #ifdef FORTRAN_UPPER_CASE
                    >> #define daxpy_ DAXPY_
                    >> ...
                    >> #endif
                    >>
                    >> A macro clearly cannot do that.
                    > Is there really no way to do that? Maybe we should ask the
                    > guys from the preprocessor library because it's troublesome to
                    > do this ifdef for adding underscores _and_ capitalising.

                    Underscores can be added by preprocessing tricks. Regarding
                    capitalizing we should indeed make sure there is no to do so by asking
                    the preprocessor library gurus. [I guess something like
                    UPPERCASE(d,a,x,p,y) should be doable].

                    > Well, I can start adapting the existing fortran.h one I know the answer
                    > to questions above.

                    Great, I'm looking forward to it,

                    Matthias
                  • Peter Schmitteckert
                    Salut, ... at least one could do something ugly like #define USE_LAPACK_UNDERSCORE_IN_NAME #define USE_LAPAPCK_CAPITAL #ifdef USE_LAPACK_UNDERSCORE_IN_NAME
                    Message 9 of 17 , Apr 7, 2003
                    • 0 Attachment
                      Salut,


                      > > #define daxpy_ DAXPY_
                      > > ...
                      > > #endif
                      > >
                      > > A macro clearly cannot do that.
                      > Is there really no way to do that? Maybe we should ask the
                      > guys from the preprocessor library because it's troublesome to
                      > do this ifdef for adding underscores _and_ capitalising.

                      at least one could do something ugly like

                      #define USE_LAPACK_UNDERSCORE_IN_NAME
                      #define USE_LAPAPCK_CAPITAL

                      #ifdef USE_LAPACK_UNDERSCORE_IN_NAME
                      #define LAPACK_UNDERSCORE(x) x##_
                      #else
                      #define LAPACK_UNDERSCORE_IN_NAME_STRING(x) x
                      #endif

                      #ifdef USE_LAPAPCK_CAPITAL
                      #define Daxpy DAXPY
                      #else
                      #define Daxpy daxpy
                      #endif


                      LAPACK_UNDERSCORE(Daxpy) (...);
                      LAPACK_UNDERSCORE(dgemm) (...);


                      best wishes,
                      Peter
                    • Peter Schmitteckert
                      Salut, ... sorry, I just realized that the Daxpy in the second last line does not get substituted. Maybe it s easier to write a small C++ programm that
                      Message 10 of 17 , Apr 7, 2003
                      • 0 Attachment
                        Salut,

                        On Mon, 7 Apr 2003, Peter Schmitteckert wrote:


                        > #define USE_LAPACK_UNDERSCORE_IN_NAME
                        > #define USE_LAPAPCK_CAPITAL
                        >
                        > #ifdef USE_LAPACK_UNDERSCORE_IN_NAME
                        > #define LAPACK_UNDERSCORE(x) x##_
                        > #else
                        > #define LAPACK_UNDERSCORE_IN_NAME_STRING(x) x
                        > #endif
                        >
                        > #ifdef USE_LAPAPCK_CAPITAL
                        > #define Daxpy DAXPY
                        > #else
                        > #define Daxpy daxpy
                        > #endif
                        >
                        >
                        > LAPACK_UNDERSCORE(Daxpy) (...);
                        > LAPACK_UNDERSCORE(dgemm) (...);
                        >


                        sorry, I just realized that the Daxpy in the second
                        last line does not get substituted.

                        Maybe it's easier to write a small C++ programm that
                        produces various header files. In real life there aren"t
                        to many version, Uppercase / lowercase. with "_" , "$" or
                        nothing, is all I ever encountered.

                        Best wishes,
                        Peter
                      • Kresimir Fresl
                        Hi Toon, Hi Mathias, hi all, Toon wrote: [...] ... [...] ... [...] I must admit that Fortran is not my strong point, so I ll leave that part to you -- but, of
                        Message 11 of 17 , Apr 7, 2003
                        • 0 Attachment
                          Hi Toon, Hi Mathias, hi all,

                          Toon wrote:
                          [...]
                          > And it would be great to do this job all together. So ideally we should
                          > revisit all components of the bindings library, discuss it, document it
                          > and implement it. We should start from the lowest level (which are the
                          > traits).
                          [...]
                          > I started of with commenting about the fortran.h file.
                          [...]

                          I must admit that Fortran is not my strong point, so I'll leave
                          that part to you -- but, of course, I hope that I'll learn something
                          from yopur discussion.

                          I think that vector and matrix traits classes are `orthogonal'
                          to Fortran issues (after all, they are used in atlas cblas/clapack
                          bindings), so I'll try to write a description and/or explanation
                          of them in next few days.

                          fres
                        • Matthias Troyer
                          ... They are indeed orthogonal, but you need to know the requirements (which you probably due): pointer to first element, stride, leading dimension, size.
                          Message 12 of 17 , Apr 7, 2003
                          • 0 Attachment
                            On Monday, April 7, 2003, at 06:21 PM, Kresimir Fresl wrote:
                            >
                            > I think that vector and matrix traits classes are `orthogonal'
                            > to Fortran issues (after all, they are used in atlas cblas/clapack
                            > bindings), so I'll try to write a description and/or explanation
                            > of them in next few days.

                            They are indeed orthogonal, but you need to know the requirements
                            (which you probably due): pointer to first element, stride, leading
                            dimension, size.

                            Looking forward to your results,

                            Matthias
                          • Matthias Troyer
                            ... This might indeed be the easiest but it is such a pity that we have to use C++ executables to create C++ code instead of using preprocessor hacks or TMP
                            Message 13 of 17 , Apr 7, 2003
                            • 0 Attachment
                              On Monday, April 7, 2003, at 06:17 PM, Peter Schmitteckert wrote:
                              >
                              > Maybe it's easier to write a small C++ programm that
                              > produces various header files. In real life there aren"t
                              > to many version, Uppercase / lowercase. with "_" , "$" or
                              > nothing, is all I ever encountered.

                              This might indeed be the easiest but it is such a pity that we have to
                              use C++ executables to create C++ code instead of using preprocessor
                              hacks or TMP tricks.
                            • Kresimir Fresl
                              ... Well, yes, I know the requirements. Vector and matrix traits are already written, see boost/numeric/bindings/traits/traits.hpp and various
                              Message 14 of 17 , Apr 7, 2003
                              • 0 Attachment
                                Matthias Troyer wrote:

                                > They are indeed orthogonal, but you need to know the requirements
                                > (which you probably due): pointer to first element, stride, leading
                                > dimension, size.

                                Well, yes, I know the requirements. Vector and matrix traits are
                                already written, see

                                boost/numeric/bindings/traits/traits.hpp

                                and various specializations:

                                - for ublas vectors, vector ranges and vector slices in

                                boost/numeric/bindings/traits/ublas_vector.hpp

                                - for std vector and valarray, boost array and C array in

                                boost/numeric/bindings/traits/std_vector.hpp
                                boost/numeric/bindings/traits/std_valarray.hpp
                                boost/numeric/bindings/traits/boost_array.hpp
                                boost/numeric/bindings/traits/c_array.hpp

                                - for ublas matrices in

                                boost/numeric/bindings/traits/ublas_matrix.hpp

                                - for ublas symmetric and hermitian matrices and adaptors in

                                boost/numeric/bindings/traits/ublas_symmetric.hpp
                                boost/numeric/bindings/traits/ublas_hermitian.hpp

                                - for TNT 1D and 2D arrays in

                                boost/numeric/bindings/traits/tnt.hpp

                                I don't think that they can be significantly changed or that
                                the need for this exists -- maybe some name changes.
                                After all, they are already used in existing blas, lapack and
                                atlas bindings.

                                fres
                              • Leonildo Tincani
                                Hi, I m not sure, because I never used it, but GNU AutoGen can convert characters: `string-downcase - lower case a new string It creates a new SCM string
                                Message 15 of 17 , Apr 7, 2003
                                • 0 Attachment
                                  Hi,
                                  I'm not sure, because I never used it, but GNU AutoGen can convert characters:
                                  `string-downcase' - lower case a new string
                                  It creates a new SCM string containing the same text as the original, only all
                                  the upper case letters are changed to lower case.

                                  You could try to verify,

                                  Leo


                                  On Monday April 7 2003 19:15, Matthias Troyer wrote:
                                  > On Monday, April 7, 2003, at 06:17 PM, Peter Schmitteckert wrote:
                                  > > Maybe it's easier to write a small C++ programm that
                                  > > produces various header files. In real life there aren"t
                                  > > to many version, Uppercase / lowercase. with "_" , "$" or
                                  > > nothing, is all I ever encountered.
                                  >
                                  > This might indeed be the easiest but it is such a pity that we have to
                                  > use C++ executables to create C++ code instead of using preprocessor
                                  > hacks or TMP tricks.
                                  >
                                  >
                                  >
                                  > To unsubscribe from this group, send an email to:
                                  > ublas-dev-unsubscribe@yahoogroups.com
                                  >
                                  >
                                  >
                                  > Your use of Yahoo! Groups is subject to http://docs.yahoo.com/info/terms/
                                • jhr.walter@t-online.de
                                  Hi Matthias (hi all), ... I ve been thinking about the lapack bindings for some time now and I believe it s firstly necessary to limit the scope of
                                  Message 16 of 17 , Apr 7, 2003
                                  • 0 Attachment
                                    Hi Matthias (hi all),

                                    you wrote:

                                    > After some discussions in my group when we had Jeremy Siek visit this
                                    > weekend we agreed that we would like to see generic C++ lapack
                                    > interfaces that work with ublas and other matrix libraries submitted to
                                    > boost as a library. One could start with the high level driver
                                    > functions only and add more specialized functions later. I thus would
                                    > like to ask what the current status of the ublas lapack bindings is and
                                    > if somebody is working on a generic version for boost, or if my
                                    > research group could help here?

                                    I've been thinking about the lapack bindings for some time now and I believe
                                    it's firstly necessary to limit the scope of investigation.

                                    The LAPACK users guide distinguishes the following parts:

                                    <cite>
                                    (1) Linear Equations
                                    (2) Linear Least Squares (LLS) Problems
                                    (3) Generalized Linear Least Squares (LSE and GLM) Problems
                                    (4) Standard Eigenvalue and Singular Value Problems
                                    (a) Symmetric Eigenproblems (SEP)
                                    (b) Nonsymmetric Eigenproblems (NEP)
                                    (c) Singular Value Decomposition (SVD)
                                    (5) Generalized Eigenvalue and Singular Value Problems
                                    (a) Generalized Symmetric Definite Eigenproblems (GSEP)
                                    (b) Generalized Nonsymmetric Eigenproblems (GNEP)
                                    (c) Generalized Singular Value Decomposition (GSVD)
                                    </cite>

                                    This structure seems to suggest some splitting into sublibraries or
                                    namespaces. Which parts are most important?

                                    Golub & van Loan reference the following procedures:

                                    Chapter 3:
                                    (1) _TRSV, _TRSM, _TRCON, _TRRFS, _TRTRS, _TRTRI
                                    (1) _GESV, _GECON, _GERFS, _GESVX, _GETRF, _GETRS, _GETRI, _GEEQU
                                    Chapter 4:
                                    (1) _GBSV, _CGBCON, _GBRFS, _GBSVX, _GBTRF, _GBTRS, _GBEQU
                                    (1) _GTSV, _GTCON, _GTRFS, _GTSVX, _GTTRF, _GBTRS
                                    (1) _POSV, _POCON, _PORFS, _POSVX, _POTRF, _POTRS, _POTRI, _POEQU
                                    (1) _PBSV, _PBCON, _PBRFS, _PBSVX, _PBTRF, _PBTRS
                                    (1) _PTSV, _PTCON, _PTRFS, _PTSVX, _PTTRF, _PTTRS
                                    (1) _SYSV, _SYCON, _SYRFS, _SYSVX, _SYTRF, _SYTRS, _SYTRI
                                    (1) _TBCON, _TBRFS, _TBTRS
                                    Chapter 5:
                                    (?) _LARFG, _LARF, _LARFX, _LARFB, _LARFT
                                    (?) _LARTG, _LARGV, _LARTV, _LARSR, CSROT, CROT, CLACGV
                                    (?) _GEQRF, _GEQPF, _ORMQR, _UNMQR, _ORGQR, _UNGQR
                                    (?) _GERQF, _GELQF, _GEQLF, _TZRQF
                                    (4.c) _GESVD, _BDSQR, _GEBRD, _ORGBR, _GBBRD
                                    (2) _GELS, _GELSS, _GELSX, _GEEQU
                                    Chapter 7:
                                    (4 b) _GEBAL, _GEBAK
                                    (4 b) _GEHRD, _ORMHR, _ORGHR, _UNMHR, _UNGHR
                                    (4 b) _HSEQR, _HSEIN
                                    (4 b) _GEES, _GEESX, _GEEV, _GEEVX
                                    (4 b) _TREVC, _TRSNA, _TREXC, _TRSEN, _TRSYL
                                    (5 b) _GGBAL, _GGHRD, _HGEQZ, _TGEVC, _GGBAK
                                    Chapter 8:
                                    (4 a) _SYEV, _SYEVD, _SYEVX
                                    (4 a) _SYTRD, _SBTRD, _SPTRD
                                    (4 a) _STEQR, _STEDC, _STERF, _PTEQR, _STEBZ, _STEIN
                                    (?) _SYGST, _PBSTF, _SBGST
                                    (4 c) _GESVD, _BDSQR, _GEBRD, _ORGBR, _GBBRD
                                    (5) _GGSVP, _TGSJA

                                    My classification currently is an educated guess only ;-) If I count these
                                    function correctly, I still find an astonishing number of around 130 ;-( To
                                    get any clue of how to proceed I believe we have no other chance than to
                                    follow Kresimir's suggestion to invent names for the corresponding families
                                    of overloaded functions or specialized classes. If nobody else is
                                    interested, I'll try to come up with a proposal after some more scanning...

                                    All the best,

                                    Joerg

                                    P.S.: I didn't see any definitive motion for your idea to wrap up the driver
                                    routines first, so this is kind of an alternative ansatz.

                                    P.P.S.: How are we able to test these wrappers?
                                  Your message has been successfully submitted and would be delivered to recipients shortly.