2010-09-19 34 views
1

我完全被难住了。我有一个用c编写的相当大的递归程序,它调用cblas_dgemm()。结果通过正常工作的程序独立验证。cblas_dgemm - 只有当(β)是2的幂时才有效

C = alpha*A*B + beta*C 

在使用随机矩阵反复试验和的参数的所有可能的组合程序给出正确答案仅当ABS(测试版)= 2^N(1,2,4,8 ...)。任何值都适用于alpha。对于测试版,任何其他正面/负面,奇数/偶数值都会在10-30%的时间内给出正确的答案。

我使用的是Ubuntu 10.04,GCC 4.4.x,我试过系统安装了blas/cblas/atlas以及手动编译的地图集。

任何提示或建议将不胜感激。我惊讶于潜藏在这个地方的美妙慷慨(聪明)的人们。

感谢你的所有提前,

拉斯

+0

你能提供你称为cblas_dgemm的确切代码吗?而且,您是否尝试直接打入Fortran常规DGEMM?我不记得在dgemm之前... – 2010-09-19 22:28:01

回答

1

是的,一个完整的例子是得心应手。以下是我使用GSL的sgemm变体悬挂的一个旧例子;应该很容易修复到double。请尝试,看看这给在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和s/sgemm/dgemm当然有问题。 – 2010-09-19 22:36:44

+0

只要我发布的问题,我意识到我需要确保直接调用cblas_dgemm()正常工作。所以,我尝试了一个像你一样的例子,当直接调用dgemm(使用相同的include/link指令)时,确实没有问题。它排除了blas/cblas链接错误等问题。它使情况更加奇怪。调用gemm的函数选择各种参数(转置/ lda/beta等)。 AFAIK,这些被正确选择 - 适用于beta(0,1,2,4,8 ..)。我会尝试进一步隔离问题。将发布结果。谢谢。拉斯 – user151410 2010-09-19 23:28:33

2

两个完全不相干的错误密谋产生虚幻的画面。这让我在错误的地方寻找问题。

(1)函数调用dgemm的逻辑有一个简单的错误。如果我不追逐错误的问题,本来会很容易解决的。 (2)我的双重比较函数:AlmostEqual2sComplement()(http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm)的双重版本使用了不正确大小的整数 - 在某些罕见情况下导致TRUE不正确。这是我第一次发现错误!

再次感谢在尝试调试程序时使用科学方法的有用建议。

拉斯

+0

有时只是向某人展示您的问题会让您发现解决方案。 – 2010-09-21 11:02:17

相关问题