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

From: Daniel Elazar <daniel.elazar_at_abs.gov.au>
Date: Fri, 12 Jun 2009 17:03:05 +1000

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 Received on Fri 12 Jun 2009 - 08:30:44 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 - 14:35:29 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