Re: modifying matrix allocation functions for use with R

[Faheem Mitha <faheem at email dot unc dot edu>]

On Sat, 6 Oct 2001, Brian Gough wrote:

> Faheem Mitha writes:
>  > One option I was considering was rewriting gsl_matrix_alloc (leaving out
>  > the error checking) in terms of R_alloc. However, I can't find the code
>  > for FUNCTION(gsl_block, alloc) (I assume this means gsl_block_alloc)
>  > used in the code for gsl_matrix_alloc. gsl_matrix_alloc itself seems to be
>  > in matrix/init_source.c. Can someone point me to the code for
>  > gsl_block_alloc, or suggest any other options?
>
> Based on my reading of the R documentation, I would suggest using the
> standard gsl_matrix_alloc in this case, and handling any errors that
> occur with the R libray function error().
>
> e.g.,
>          m = gsl_matrix_alloc(...)
>          if (m == 0) { error ("could not allocate matrix"); } ;
>          ....
>          gsl_matrix_free(m);
>
> This would only add a few lines to your code, and will work with the
> standard GSL library so your extension will be more portable.

Thanks, this is a good sugestion. After considering all the alternatives,
I decided to go with something like this.

I think I spotted a typo. In the section LU Decomposition, it is written

 - Function: int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation
          * P, int *SIGNUM)
     This function factorizes the square matrix A into the LU
     decomposition PA = LU.  On output the diagonal and upper
     triangular part of the input matrix A contain the matrix R. The
     lower triangular part of the input matrix (excluding the diagonal)
     contains L.  The diagonal elements of L are unity, and are not
     stored.

Surely, this should be "...the input matrix A contain the matrix U"?

Incidentally, it is a little awkward to compute the inverse of a matrix,
since this involves calling both int gsl_linalg_LU_decomp and int
gsl_linalg_LU_invert. Would it not be useful to have a "wrapper" function
which uses both together?

I have one (slighly idle) question. I have noticed that macros are used
quite liberally in the gsl source code, for example you have a macro
called TYPE which seems to do nothing at all, there is another macro
called PUNCTION(foo, bar) which concatenates foo and bar etc. I'm just
wondering what purpose some of these serve, like the TYPE macro for
example. Do they just make the code easier to write?

                                    Sincerely, Faheem Mitha.