Daniel Elazar
2009-Jun-12 07:03 UTC
[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
Cm34 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]]
John Nolan
2009-Jun-12 12:44 UTC
[Rd] Can't get F77_CALL(dgemm) to work [SEC=UNCLASSIFIED]
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@american.edu
202.885.3140 voice
202.885.3155 fax
http://academic2.american.edu/~jpnolan
...........................................................................
-----r-devel-bounces@r-project.org wrote: -----
To: r-devel@r-project.org
From: Daniel Elazar <daniel.elazar@abs.gov.au>
Sent by: r-devel-bounces@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
Cm34 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]]