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 jpno...@american.edu 202.885.3140 voice 202.885.3155 fax http://academic2.american.edu/~jpnolan ........................................................................... -----r-devel-boun...@r-project.org wrote: ----- To: r-devel@r-project.org From: Daniel Elazar <daniel.ela...@abs.gov.au> Sent by: r-devel-boun...@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]] ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
