Re: [Rd] Can't get F77_CALL(dgemm) to work [SEC=UNCLASSIFIED]

From: John Nolan <jpnolan_at_american.edu>
Date: Fri, 12 Jun 2009 08:44:27 -0400

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_at_american.edu
202.885.3140 voice
202.885.3155 fax
http://academic2.american.edu/~jpnolan

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



-----r-devel-bounces_at_r-project.org wrote: -----

To: r-devel_at_r-project.org
From: Daniel Elazar <daniel.elazar_at_abs.gov.au> Sent by: r-devel-bounces_at_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_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

        [[alternative HTML version deleted]]



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 12 Jun 2009 - 12:52:12 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 12 Jun 2009 - 15:05:40 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive