Manipulation de matrices en C via R

Question

Je veux écrire quelques fonctions de manipulation de matrices en C, puis les passer à R, où la matrice sera fournie par R, et obtenir les résultats des manipulations. J'ai une fonction de test comme indiqué ci-dessous (s'il vous plaît, ne vous souciez pas de ce qu'elle fait, dans mes fonctions originales, j'aurai besoin de choisir un élément aléatoire de chaque ligne, et de faire quelques calculs sur eux, et de retourner un tableau composé de ces éléments choisis au hasard de chaque ligne, en d'autres termes, je dois avoir deux boucles for pour passer par tous les éléments de la matrice).

void multMat(double **A, int *r, int *c, double *s)
{
    int i, j;

    for (i = 0; i < *r; ++i)
    {

       for (j = 0; j < *c; ++j) 
       {
           if (j == 5)
               s[i] = A[i][j] * A[i][0];
       }
    }
}

Je l'ai compilé avec R CMD SHLIB multMat.c et il produit multMat.so pour moi. Ensuite, du côté R, j'ai quelque chose comme ça :

dyn.load("multMat.so")

multMat <- function(A)
{
  .C("multMat", A=as.double(A), r=as.integer(nrow(A)), c=as.integer(ncol(A)), s=as.double(nrow(A)))
}

Ensuite, j'ai créé une matrice de test comme dans R Studio et j'ai appelé cette fonction :

A <- matrix(1:100, 10, 10)
multMat(A)

Le problème est que lorsque j'exécute cette fonction, R Studio se plante. Je suppose qu'il y a un problème avec la façon dont la fonction C est définie. Avez-vous une idée ?

Answer 1

Où vous avez tort

Vous définissez A como double ** ,

void multMat(double **A, int *r, int *c, double *s)

tout en passant un double * :

.C("multMat", A=as.double(A), r=as.integer(nrow(A)), c=as.integer(ncol(A)), s=as.double(nrow(A)))

Vous devez réécrire votre fonction C en utilisant un tableau à une dimension. Définissez votre fonction comme :

void multMat(double *A, int *r, int *c, double *s)

et remplacer A[i][j] par A[j * r + i] . (Si je ne vous ai pas mal compris, r est la dimension principale).

Problème de performance :

Pour l'instant, c'est i-j boucle avec j-loop comme boucle interne, de sorte que vous balayez une ligne de la matrice dans la boucle la plus interne. Ceci est pas facile à mettre en cache . Vous devriez échanger la boucle pour obtenir j-i boucle.

Je pense que vous êtes en fait conscient du problème du cache. En C, les matrices sont stockées dans l'ordre des lignes majeures, donc i-j est optimale ; mais dans R, les matrices sont stockées dans l'ordre des colonnes majeures, donc j-i est optimale.

Peut-être que le style de stockage différent des matrices causera des problèmes dans votre code existant. Vous pouvez y réfléchir à deux fois ici. Votre code C original suppose qu'il prend une matrice qui est stockée par ligne, alors que si vous initialisez votre matrice dans R et que vous l'introduisez dans C, elle est stockée par colonne. Il est possible que certaines modifications soient nécessaires. Si cela entraîne une modification trop importante de votre code, vous pouvez essayer de passer la transposition de votre matrice R à C.

Vous devez également utiliser davantage de variables locales/automatiques dans votre code C, au lieu d'utiliser des pointeurs . Par exemple, remplacez

for (i = 0; i < *r; ++i)

par

int r_local = *r;
for (i = 0; i < r_local; ++i)

Vous obtenez des bottes performantes en Réutilisation des registres du CPU y réduction de l'enseignement (pas besoin de déréférencer le point à chaque itération).

R type de données C/FORTRAN supporté

Il n'y a donc aucun moyen de passer directement la matrice depuis R, puis d'utiliser la convention occasionnelle pour la matrice, c'est-à-dire A[i][j] dans le côté C ?

Non, il n'y en a pas. R ne supporte pas double ** . Les matrices dans R sont stockées en format colonne-majeur dans un long vecteur, et appartiennent à double * type.

R storage mode  C type        FORTRAN type
logical         int *         INTEGER
integer         int *         INTEGER
double          double *      DOUBLE PRECISION
complex         Rcomplex *    DOUBLE COMPLEX
character       char **       CHARACTER*255
raw unsigned    char *        none

Ver section 5.2 de Writing R extensions .