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
