2010-09-19 28 views
1

Je suis totalement perplexe. J'ai un assez gros programme récursif écrit en c qui appelle cblas_dgemm(). Le résultat est vérifié indépendamment par un programme qui fonctionne correctement.cblas_dgemm - fonctionne UNIQUEMENT si (beta) est power-of-two

C = alpha*A*B + beta*C 

sur des tests répétés à l'aide de matrices aléatoires et toutes les combinaisons possibles de paramètres du programme donne la réponse correcte seulement si abs (bêta) = 2^n (1,2,4,8 ..). Toute valeur fonctionne pour alpha. Toute autre valeur positive/négative, impaire/paire pour la bêta donne une réponse correcte de 10 à 30% du temps. J'utilise Ubuntu 10.04, GCC 4.4.x, j'ai essayé blas/cblas/atlas installé par le système aussi bien qu'atlas compilé manuellement.

Tous les conseils ou suggestions seraient grandement appréciés. Je suis étonné par les gens merveilleusement généreux (et intelligents) qui se cachent sur ce site.

vous Remerciant tout à l'avance,

Russ

+0

Pouvez-vous fournir le code exact avec lequel vous avez appelé cblas_dgemm? Et aussi, avez-vous essayé d'appeler directement dans la routine fortran DGEMM? Je ne me souviens pas avoir eu des problèmes avec dgemm avant ... –

Répondre

1

Oui, un exemple complet serait à portée de main. Voici un vieil exemple que j'ai traîné en utilisant la variante sgemm de GSL; devrait être facile à corriger à double. S'il vous plaît essayer de voir si cela donne le résultat indiqué dans le manuel GSL:

/* from the gsl info documentation in node 'gsl cblas examples' */ 
/* compile via 'gcc -o $file $file.c -lgslcblas' */ 
/* edd 15 Nov 2003 */ 

#include <stdio.h>  
#include <gsl/gsl_cblas.h> 

int 
main (void)  
{  
    int lda = 3; 
    float A[] = { 0.11, 0.12, 0.13, 
       0.21, 0.22, 0.23 }; 
    int ldb = 2;      
    float B[] = { 1011, 1012, 
       1021, 1022,              
       1031, 1032 }; 
    int ldc = 2; 
    float C[] = { 0.00, 0.00, 
       0.00, 0.00 };  
    /* Compute C = A B */ 
    cblas_sgemm (CblasRowMajor, 
       CblasNoTrans, CblasNoTrans, 2, 2, 3, 
       1.0, A, lda, B, ldb, 0.0, C, ldc); 
    printf ("[ %g, %g\n", C[0], C[1]);   
    printf (" %g, %g ]\n", C[2], C[3]); 

    return 0;  
}   
+0

s/float/double et s/sgemm/dgemm bien sûr. –

+0

Dès que j'ai posté la question, j'ai réalisé que je devais m'assurer qu'un appel direct à cblas_dgemm() fonctionne correctement. Donc, j'ai essayé un exemple comme le vôtre et en effet il n'y a pas de problèmes quand dgemm est appelé directement (en utilisant des directives include/link identiques). Ce qui exclut les problèmes avec les erreurs de liaison blas/cblas etc. Cela rend la situation encore plus bizarre. La fonction appelant gemm sélectionne différents paramètres (transpose/lda/beta etc.). AFAIK, ceux-ci sont cueillis correctement - fonctionne pour beta (0,1,2,4,8 ..). Je vais essayer d'isoler le problème plus loin. Va poster des résultats. Merci. Russ – user151410

2

Deux erreurs sans aucun lien conspiré pour produire une image trompeuse. Cela m'a fait chercher des problèmes au mauvais endroit.

(1) Il y avait une erreur simple dans la logique de la fonction appelée dgemm. Aurait été facilement réparé si je ne cherchais pas le mauvais problème.

(2) Ma double fonction de comparaison: la version double de AlmostEqual2sComplement() (http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm) a utilisé un entier de taille incorrecte - ce qui donne un TRUE incorrect dans de rares circonstances. C'était la première fois que l'erreur me mordait!

Merci encore pour la suggestion utile d'utiliser la méthode scientifique lorsque vous essayez de déboguer un programme.

Russ

+0

Parfois, le simple fait de montrer votre problème à quelqu'un vous permet de repérer la solution. –