From: John Nolan <jpnolan_at_american.edu>

Date: Fri, 12 Jun 2009 08:44:27 -0400

9 10 11 12

13 14 15 16

17 18 19 20

}

}

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel

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

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