레벨 3의 복잡한 BLAS 함수가 잘못된 값 오류를 던졌습니다.

Question

복잡한 행렬 행렬 곱셈을위한 Fortran 함수를 호출하는 함수를 작성하고 있습니다. 복잡한 곱셈을 위해 CGEMM_ 및 ZGEMM_ 함수를 호출합니다. 모든 xGEMM_ 함수는 본질적으로 동일하므로 SGEMM_에서 CGEMM__ 및 ZGEMM_으로 코드를 복사했습니다. 유일한 변경은 각각의 데이터 유형입니다. SGEMM_ 및 DGEMM_ 함수는 정상적으로 작동하지만 CGEMM_은 오류를 발생시킵니다. 모든 입력은 동일합니다.

** On entry to CGEMM parameter number 13 had an illegal value 

및 zgemm_ 정말 무슨 일이 일어나고 있는지 전혀 모르는

** On entry to ZGEMM parameter number 1 had an illegal value 

가 발생합니다. liblapack 패키지에 어떤 종류의 버그가 있습니까? liblapack-dev 패키지를 사용하고 있습니다. 내 큰 코드의 작은 버전을 만들었고 CGEMM_에서 여전히 동일한 오류가 발생합니다.

32 비트 시스템을 실행 중이고 문제인지 궁금합니다.

코드 : 다음과 같이

#include<iostream> 
using namespace std; 
#include<stdlib.h> 
#include<string.h> 
#include<complex> 

typedef complex<float> c_float; 
extern "C" 
{c_float cgemm_(char*,char*,int*,int*,int*,c_float*, c_float[0],int*,c_float[0],int*,c_float*,c_float[0],int*);//Single Complex Matrix Multiplication 
} 

c_float** allocate(int rows, int columns) 
{ 
    c_float** data; 

    // Allocate Space 
    data = new c_float*[columns]; //Allocate memory for using multidimensional arrays in column major format. 
    data[0] = new c_float[rows*columns]; 
    for (int i=0; i<columns; i++) 
    { 
     data[i] = data[0] + i*rows; 
    } 

    // Randomize input 
    for (int i=0; i<columns; i++) 
    {for (int j=0; j<rows; j++) 
     { 
      data[j][i] =complex<double>(drand48()*10 +1,drand48()*10 +1); //Randomly generated matrix with values in the range [1 11) 
      } 
    } 
    return(data); 
} 

// Destructor 
void dest(c_float** data) 
{ 
    delete [] data[0]; 
    delete [] data; 
} 

// Multiplication 
void mult(int rowsA,int columnsA, int rowsB,int columnsB) 
{ 
    c_float **matA,**matB,**matC; 
    char transA, transB; 
    int m,n,k,LDA,LDB,LDC; 
    c_float *A,*B,*C; 
    c_float alpha(1.0,0.0); 
    c_float beta(0.0,0.0); 

    matA = allocate(rowsA,columnsA); 
    matB = allocate(rowsB,columnsB); 
    matC = allocate(rowsA,columnsB); 

    transA = 'N'; 
    transB = 'N'; 
    A = matA[0]; 
    B = matB[0]; 
    m = rowsA; 
    n = columnsB; 
    C = matC[0]; 
    k = columnsA; 
    LDA = m; 
    LDB = k; 
    LDC = m; 
    cout<<"Matrix A"<<endl; 
    for (int i=0; i<rowsA; i++) 
    {for (int j=0; j<columnsA; j++) 
     { 
      cout<<matA[i][j]; 
      cout<<" "; 
     }cout<<endl; 
    } 
    cout<<"Matrix B"<<endl; 
    for (int i=0; i<rowsB; i++) 
    {for (int j=0; j<columnsB; j++) 
     { 
      cout<<matB[i][j]; 
      cout<<" "; 
     }cout<<endl; 
    } 

    cgemm_(&transA,&transB,&m,&n,&k,&alpha,A,&LDA,B,&LDB,&beta,C,&LDC); 
    cout<<"Matrix A*B"<<endl; 
    for (int i=0; i<rowsA; i++) 
    {for (int j=0; j<columnsB; j++) 
     { 
      cout<<matC[i][j]; 
      cout<<""; 
     } 
     cout<<endl; 
    } 
    dest(matA); 
    dest(matB); 
    dest(matC); 
} 



main() 
{ 
    mult (2,2,2,2); 
} 

출력 및 Valgrind의 보고서는 다음과 같습니다

----------------------------------------- 
Compilation using g++ -g -o matrix Matrix_multiplication.cpp -lblas -llapack -lgfortran 
./matrix gives 
Matrix A 
(1.00985,1) (1.91331,4.64602) 
(2.76643,1.41631) (5.87217,1.92298) 
Matrix B 
(5.54433,6.2675) (6.6806,10.3173) 
(9.31292,3.33178) (1.50832,6.56094) 
** On entry to CGEMM parameter number 1 had an illegal value 

Valgrind output looks like 

==4710== Memcheck, a memory error detector 
==4710== Copyright (C) 2002-2013, and GNU GPL'd, by Julian Seward et al. 
==4710== Using Valgrind-3.10.0.SVN and LibVEX; rerun with -h for copyright info 
==4710== Command: ./o 
==4710== Parent PID: 3337 
==4710== 
==4710== Conditional jump or move depends on uninitialised value(s) 
==4710== at 0x46E5096: lsame_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x46DD683: cgemm_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x8048C7E: mult(int, int, int, int) (Matrix_multiplication.cpp:83) 
==4710== by 0x8048D70: main (Matrix_multiplication.cpp:102) 
==4710== Uninitialised value was created by a stack allocation 
==4710== at 0x8048A18: mult(int, int, int, int) (Matrix_multiplication.cpp:43) 
==4710== 
==4710== Conditional jump or move depends on uninitialised value(s) 
==4710== at 0x46DD686: cgemm_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x8048C7E: mult(int, int, int, int) (Matrix_multiplication.cpp:83) 
==4710== by 0x8048D70: main (Matrix_multiplication.cpp:102) 
==4710== Uninitialised value was created by a stack allocation 
==4710== at 0x8048A18: mult(int, int, int, int) (Matrix_multiplication.cpp:43) 
==4710== 
==4710== Conditional jump or move depends on uninitialised value(s) 
==4710== at 0x46E5096: lsame_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x46DD7B1: cgemm_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x8048C7E: mult(int, int, int, int) (Matrix_multiplication.cpp:83) 
==4710== by 0x8048D70: main (Matrix_multiplication.cpp:102) 
==4710== Uninitialised value was created by a stack allocation 
==4710== at 0x8048A18: mult(int, int, int, int) (Matrix_multiplication.cpp:43) 
==4710== 
==4710== Conditional jump or move depends on uninitialised value(s) 
==4710== at 0x46DD7B4: cgemm_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x8048C7E: mult(int, int, int, int) (Matrix_multiplication.cpp:83) 
==4710== by 0x8048D70: main (Matrix_multiplication.cpp:102) 
==4710== Uninitialised value was created by a stack allocation 
==4710== at 0x8048A18: mult(int, int, int, int) (Matrix_multiplication.cpp:43) 
==4710== 
==4710== Conditional jump or move depends on uninitialised value(s) 
==4710== at 0x46E5096: lsame_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x46DD859: cgemm_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x8048C7E: mult(int, int, int, int) (Matrix_multiplication.cpp:83) 
==4710== by 0x8048D70: main (Matrix_multiplication.cpp:102) 
==4710== Uninitialised value was created by a stack allocation 
==4710== at 0x8048A18: mult(int, int, int, int) (Matrix_multiplication.cpp:43) 
==4710== 
==4710== Conditional jump or move depends on uninitialised value(s) 
==4710== at 0x46DD85C: cgemm_ (in /usr/lib/atlas-base/atlas/libblas.so.3.0) 
==4710== by 0x8048C7E: mult(int, int, int, int) (Matrix_multiplication.cpp:83) 
==4710== by 0x8048D70: main (Matrix_multiplication.cpp:102) 
==4710== Uninitialised value was created by a stack allocation 
==4710== at 0x8048A18: mult(int, int, int, int) (Matrix_multiplication.cpp:43) 
==4710== 
==4710== 
==4710== HEAP SUMMARY: 
==4710==  in use at exit: 120 bytes in 6 blocks 
==4710== total heap usage: 43 allocs, 37 frees, 13,897 bytes allocated 
==4710== 
==4710== LEAK SUMMARY: 
==4710== definitely lost: 0 bytes in 0 blocks 
==4710== indirectly lost: 0 bytes in 0 blocks 
==4710==  possibly lost: 0 bytes in 0 blocks 
==4710== still reachable: 120 bytes in 6 blocks 
==4710==   suppressed: 0 bytes in 0 blocks 
==4710== Rerun with --leak-check=full to see details of leaked memory 
==4710== 
==4710== For counts of detected and suppressed errors, rerun with: -v 
==4710== ERROR SUMMARY: 6 errors from 6 contexts (suppressed: 0 from 0) 

편집 : 질문을 실행할 수있는 코드로 수정되었습니다. 문제는 동일하게 남아 있으며 질문의 성격은 변하지 않았습니다.

Answer 1

xgemm_ 호출에 전달되는 또는 transB 문자열의 길이가 표시되지 않습니다.

Character Fortran의 더미에는 '숨김'길이 인수가 수반됩니다. 예를 들어 GCC 4.9.0에서 사용되는 규칙은 여기에서 더 자세히 설명됩니다.

https://gcc.gnu.org/onlinedocs/gcc-4.9.0/gfortran/Argument-passing-conventions.html.

인수 목록에서 숨겨진 인수의 위치는 플랫폼에 따라 다릅니다. 리눅스에서는 모든 명시적인 인수 뒤에 위치합니다.

나는 나의 문자열의 길이를 나타내는, s의 (포트란에 숨겨진) 세 번째 인수로 1를 전달하고 리눅스에서 실행되는 s.f90

Subroutine s(c, r) 
    Character (*), Intent (In) :: c 
    Real, Intent (In) :: r 
    Print '(3A,I0,A)', 'c = "', c, '", (len(c)=', len(c), ')' 
    Print *, 'r = ', r 
End Subroutine 

및 main.c

#include <string.h> 
int main(void) 
{ 
    char c[1+1]; 
    float r=4.2; 
    strcpy(c,"A"); 
    s_(c,&r,1); 
} 

을 고려하십시오.

컴파일 및 gfortran 실행은 나에게

> gfortran -g main.c s.f90 && ./a.out 
c = "A", (len(c)=1) 
r = 4.19999981 

그래서 아마 당신의 xgemm_ 통화 ...,&LDC,1,1);을해야을 제공합니다?

Answer 2

Fortran에서 문자 변수의 길이에 대한 대답은 본질적으로 정확하지만 여기서의 문제는 아닙니다. blas 라이브러리 내의 함수의 고정 길이 문자 변수는 함수 인수에서 길이를 읽지 않을 것입니다.함수에 대해 이것을 검사했고 -O0에서도 길이가 컴파일 타임 상수였습니다.

특정 문제의 원인은 c_float cgemm_(...입니다. 여기서 cgemm_은 c_float을 반환한다는 것을 컴파일러에 알립니다. 일반적으로 반환 값은 레지스터에 저장되지만 너무 큰 경우 스택으로 이동할 수도 있습니다. 귀하의 경우, 32 비트 시스템에서 이것은 8 바이트 c_float의 경우 인 것 같습니다. 함수를 void cgemm_ (반드시 있어야 함) 또는 int cgemm_ (레지스터를 사용할 것임)으로 정의하면이 문제가 해결됩니다.

take-home 메시지는 "하지 마라."이것은 hackish 한 호출 방법이며 다른 플랫폼/컴파일러를 다룰 때 두통을 일으킬 것입니다. blas 작업을 위해 cblas 인터페이스 나 C++ 라이브러리를 사용하는 것이 훨씬 낫습니다.