From: Daniel Elazar <daniel.elazar_at_abs.gov.au>

Date: Fri, 12 Jun 2009 17:03:05 +1000

9 10 11 12

13 14 15 16

17 18 19 20

}

}

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

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

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 <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.
*