cblas_dgemm - работает ТОЛЬКО если (бета) имеет степень двойки

Question

Я полностью в тупике. У меня есть довольно большая рекурсивная программа, написанная на c, которая вызывает cblas_dgemm (). Результат проверяется независимо программой, которая работает правильно.

C = alpha*A*B + beta*C 

На повторных тестах с использованием случайных матриц и всех возможных комбинаций параметров программа дает правильный ответ ТОЛЬКО если abs (бета) = 2 ^ n (1,2,4,8 ..). Любое значение работает на альфа. Любое другое положительное / отрицательное, нечетное / четное значение для беты дает правильный ответ ч / б в 10-30% случаев.

Я использую Ubuntu 10.04, GCC 4.4.x, я пробовал установить систему blas / cblas / atlas, а также скомпилированный вручную атлас.

Любые советы или предложения будут с благодарностью. Я поражен удивительно щедрыми (и умными) людьми, скрывающимися на этом сайте.

Заранее благодарю всех,

Расс

Answer 1

Две совершенно не связанные ошибки сговорились, чтобы создать иллюзорную картину. Это заставило меня искать проблемы не в том месте.

(1) В логике функции, вызывающей dgemm, произошла простая ошибка. Было бы легко исправить, если бы я не преследовал неправильную проблему.

(2) Моя функция двойного сравнения: двойная версия почтиEqual2sComplement () (http://www.cygnus -software.com /apers / comparingfloats / comparingfloats.htm ) использовала целочисленное значение неправильного размера - в результате неверное ИСТИНА при определенных редких обстоятельствах. Это был первый раз, когда ошибка укусила меня!

Еще раз спасибо за полезное предложение использовать научный метод при попытке отладки программы.

Russ

Answer 2

Да, полный пример будет полезен.Вот старый пример, который я использовал, используя вариант sgemm GSL;должно быть легко исправить на 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;    
}