1

我完全被难住了。我有一个用 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。

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

先谢谢大家了

拉斯

4

2 回答 2

2

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

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

(2) 我的双重比较功能：AlmostEqual2sComplement() 的双重版本（http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm）使用了不正确大小的整数 - 在某些极少数情况下导致不正确的 TRUE . 这是错误第一次咬我！

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

拉斯

于 2010-09-20T01:00:43.477 回答
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;    
}
于 2010-09-19T22:33:59.697 回答