cblas_dgemm - 仅当（beta）是二的幂时才有效

Question

我完全被难住了。我有一个用 c 编写的相当大的递归程序，它调用 cblas_dgemm()。结果由正确运行的程序独立验证。

C = alpha*A*B + beta*C 

在使用随机矩阵和所有可能的参数组合进行重复测试时，只有当 abs(beta) = 2^n (1,2,4,8..) 时，程序才会给出正确答案。任何值都适用于 alpha。任何其他正/负、奇/偶的 beta 值在 10-30% 的时间内给出正确的黑白答案。

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

任何提示或建议将不胜感激。我对潜伏在这个网站上的非常慷慨（而且聪明）的人们感到惊讶。

预先感谢大家，

拉斯

Answer 1

两个完全不相关的错误共同产生了一幅虚幻的画面。这让我在错误的地方寻找问题。

(1) 调用 dgemm 的函数逻辑存在简单错误。如果我没有追寻错误的问题，本来很容易解决。

（2）我的双重比较函数：AlmostEqual2sComplement()的双重版本（http: //www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm）使用了错误大小的整数 - 在某些罕见的情况下导致错误的 TRUE。这是我第一次被这个错误咬到！

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

拉斯

Answer 2

是的，一个完整的例子会很方便。这是我使用 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;    
}