Daniel,

R apparently uses Fortran order AND storage method when
storing a matrix.  For an (n x m) matrix, Fortran allocates
a single block of nm doubles and stores them in the order
A(1,1),A(2,1),A(3,1),...,A(n,1),A(1,2),A(2,2),...,A(n,m).

In contrast, C allocates a vector of n pointers, each pointing
to a row of the matrix, and then a vector of doubles of length
m for each row. This uses more storage space: nm doubles
and n pointers. Depending on how the matrix is defined, there
is no guarantee that the rows are consecutive memory.  IF the
rows are in consecutive memory and you pass the address of
the first element to Fortran, it will "see" the transpose
of A, not A.

You can force them to be in consecutive memory by allocating
the whole block at once.   It is not clear from the Writing
R extensions manual how Calloc allocates.  The function
dmatrix that is in the free C routines nrutil.c from the
Numerical Recipes books deliberately allocates a block.

Another way to deal with this is to write short routines
to explicitly copy one form of matrix to the other.  This
takes more memory and is a bit slower, but safer and works
in all cases. If anyone wants such code, let me know.

John


...........................................................................

John P. Nolan
Math/Stat Department
227 Gray Hall
American University
4400 Massachusetts Avenue, NW
Washington, DC 20016-8050

jpnolan@american.edu
202.885.3140 voice
202.885.3155 fax
http://academic2.american.edu/~jpnolan
...........................................................................



-----r-devel-bounces@r-project.org wrote: -----


To: r-devel@r-project.org
From: Daniel Elazar <daniel.elazar@abs.gov.au>
Sent by: r-devel-bounces@r-project.org
Date: 06/12/2009 03:03AM
Subject: [Rd] Can't get F77_CALL(dgemm) to work [SEC=UNCLASSIFIED]


Hi

I am new to writing C code and am trying to write an R extension in C.  I
have hit a wall with F77_CALL(dgemm) in that it produces wrong results.
The code below is a simplified example that multiplies the matrices Ab and
Bm to give Cm.  The results below show clearly that Cm is wrong.

 Am=
 1     2    3    4
 5     6    7    8
 9   10  11 12
13 14  15  16
17 18  19  20


 Bm=
 1  1  1
 1  1  1
 1  1  1
 1  1  1


 Cm=
34 38 42
46 50 34
38 42 46
50 34 38
42 46 50


I have searched the internet and suspect it has something to do with column
major matrix format for Fortran being inconsistent with the row major
format for C but I'm not sure how to fix this in C.  One suggestion I came
across (http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=915)  is to
use cblas_dgemm in which the option 'CblasColMajor' can be specified.
However, I would have thought that F77_CALL(dgemm) should work as it has
been used in some R packages.  I'm also not sure that cblas would work from
R.

I tried inputting matrices into dgemm as 2 dimensional arrays and as one
dimensional arrays.  However the results for Cm were the same in both
cases.

Any help would be much appreciated.

Cheers
Daniel



C Code

#include <R.h>
#include <stdio.h>
#include <R_ext/Lapack.h>
#include <R_ext/Applic.h>
#include "math.h"
#define BLAS_H
#include "MPQL.h"
#include <R_ext/PrtUtil.h>

void MPQL(int *iterations) {

      double **Am;
      double  *Am_vec;
      double **Bm;
      double  *Bm_vec;
      double **Cm;
      double  *Cm_vec;
      double  one          =  1.0;
      double  zero         =  0.0;
      int     j, r, c;
      int   rows_A =5;
      int   cols_A =4;
      int   rows_B =4;
      int   cols_B =3;
      int   rows_C =5;
      int   cols_C =3;



Am            = Calloc(rows_A,double *);
Am_vec        = Calloc(rows_A*cols_A,double);
Bm            = Calloc(rows_B,double *);
Bm_vec        = Calloc(rows_B*cols_B,double);
Cm            = Calloc(rows_C,double *);
Cm_vec        = Calloc(rows_C*cols_C,double);

for (j = 0; j < rows_A; j++)  Am[j]  = Am_vec + j * cols_A;
for (j = 0; j < rows_B; j++)  Bm[j]  = Bm_vec + j * cols_B;
for (j = 0; j < rows_C; j++)  Cm[j]  = Cm_vec + j * cols_C;


for (r = 0; r < rows_A; r++)
   for (c = 0; c < cols_A; c++)
       {
       Am[r][c] = r*(cols_A) + c + 1.0;
       };

for (r = 0; r < rows_B; r++)
   for (c = 0; c < cols_B; c++)
       Bm[r][c] = 1.0;

Rprintf("\n\n Am= \n");
for (r = 0; r < rows_A; r++)
   for (c = 0; c < cols_A; c++)
        if (c==(cols_A - 1))  Rprintf("%2.0f \n" ,(double) Am[r][c]);
else
               Rprintf("%2.0f ",(double) Am[r][c]);

Rprintf("\n\n Am_vec= \n");  for (r = 0; r < (rows_A * cols_A); r++)
Rprintf("%2.0f  " ,(double) Am_vec[r]);


Rprintf("\n\n Bm=\n");
for (r = 0; r < rows_B; r++)
   for (c = 0; c < cols_B; c++)
        if (c==(cols_B-1))  Rprintf("%2.0f \n" ,(double) Bm[r][c]);   else
               Rprintf("%2.0f ",(double) Bm[r][c]);

Rprintf("\n\n Bm_vec= \n");  for (r = 0; r < (rows_B * cols_B); r++)
Rprintf("%2.0f  " ,(double) Bm_vec[r]);


/*  Inputting matrices as 2 dimensional arrays gives same results as one
dimensional form below
F77_CALL(dgemm)("N", "N", &rows_A, &cols_B, &cols_A, &one, *Am, &rows_A,
*Bm, &rows_B, &zero, *Cm, &rows_C);
*/

F77_CALL(dgemm)("N", "N", &rows_A, &cols_B, &cols_A, &one, Am_vec, &rows_A,
Bm_vec, &rows_B, &zero, Cm_vec, &rows_C);


Rprintf("\n\n Cm=\n");
for (r = 0; r < rows_C; r++)
   for (c = 0; c < cols_C; c++)
        if (c==(cols_C-1))  Rprintf("%2.0f \n" ,(double) Cm[r][c]);   else
               Rprintf("%2.0f ",(double) Cm[r][c]);

Rprintf("\n\n Cm_vec= \n");  for (r = 0; r < (rows_C * cols_C); r++)
Rprintf("%2.0f  " ,(double) Cm_vec[r]);


Free(Cm_vec);
Free(Cm);
Free(Bm_vec);
Free(Bm);
Free(Am_vec);
Free(Am);
}


Header File

/*  C:\Data\RPackages\MPQL\src\MPQL.h   */

#include <Rmath.h>

/* #include <R_ext/RS.h> */

void MPQL      (int *iterations);



R File compiled with C code above

MPQL <- function(iterations){

Result <-  .C("MPQL",
            as.integer(iterations),
                DUP           = TRUE,
                PACKAGE       = "MPQL"
            )
}


R code to call the compiled C dll


rm(list = ls(all = TRUE))
gc()

# Load R packages.
library(reshape)
library(MPQL)

SAE.MPQL <- function(its) {
MPQL.object <- MPQL(its)
}

BOTT <- SAE.MPQL(its=2)

------------------------------------------------------------------------------------------------

Free publications and statistics available on www.abs.gov.au


    [[alternative HTML version deleted]]

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
	[[alternative HTML version deleted]]

