sgetrf M N is 103040 时报错，这是个bug么 lapack and Openblas the same，修复备忘

号外：

$ clang-format   -style="{BasedOnStyle: llvm, IndentWidth: 4}"   -i  hello.cpp

$ clang-format   -style="{BasedOnStyle: llvm, IndentWidth: 4}"   -i  hello.cpp

IndentWidth:4不错，默认2太下了

1,现象

M=N=103040时，调用 sgetrf_ 时，无论是 LAPACK 还是 OpenBLAS，都出错：

openblas：

lapack：

2, 复现代码

出现问题的应该是由于M和N相对数字太大，乘积超出32bit整数的表达范围，而接收此参数的类型其实为 unsigned long int，导致误传非常大的值后造成越界。

如下是已经修复的代码，可以正常运行了：

extern "C" void sgetrf_(int* M, int* N, float* A, int *lda, int* piv, int* info);#include <stdlib.h>
#include <stdio.h>
#include <cmath>
#include <iostream>#define ORDER (10304)
//#define ORDER (51520)void print_matrix(int M, int N, float* A, int lda)
{for(unsigned long int i=0; i<M; i++){for(unsigned long int j=0; j<N; j++){printf(" %7.4f", A[i + j*lda]);}printf("\n");}
}void init_matrix(int M, int N, float* A, unsigned long int lda, int seed)
{srand(seed);for(unsigned long int i=0; i<M; i++){for(unsigned long int j=0; j<N; j++){A[i + j*lda] =((float) rand())/RAND_MAX;}}
}int main()
{float* A = NULL;int M = ORDER;int N = M;int lda = M;unsigned long int MM = M;unsigned long int NN = N;unsigned long int ldaa = lda;unsigned long int min_MN = std::min(M, N);int *piv = NULL;int *info = NULL;printf("lda * N * sizeof(float) bytes = %ld\n", (ldaa * NN * sizeof(float)));printf("lda * N * sizeof(float) bytes = %f GB\n", (ldaa * NN * sizeof(float))/1024.0/1024.0/1024.0);A = (float*)malloc(ldaa * NN * sizeof(float));if(A==NULL){printf("failed malloc()\n");}piv = (int*)malloc(min_MN*sizeof(int));info = (int*)malloc(1*sizeof(int));init_matrix(M, N, A, lda, 2024);//printf("A =\n");	print_matrix(7, 7, A, lda);printf("A[%ld] = %7.3f\n", MM -1 + (NN-1)*ldaa, A[MM -1 + (NN-1)*ldaa]);sgetrf_(&M, &N, A, &lda, piv, info);  printf("LU=\n");	print_matrix(7, 7, A, lda);free(A);free(piv);free(info);return 0;
}

3，结论

遇到非负整数，比如阶数、数组下标等，尽量用 signed long int 类型，代替使用 int 类型，现在的数据量太大了，下标/ 数据量 动辄超过 512M/ 2GB 等；

openblas中的 lapack_int 数据类型：

                             #define lapack_int        int64_t

                                                                       ^

                             typedef __int64_t         int64_t;

                                                 ^

typedef signed long int    __int64_t;

详情如下：

 third-party/openblas/local/include/lapack.h

#ifndef lapack_int
#if defined(LAPACK_ILP64)
#define lapack_int        int64_t
#else
#define lapack_int        int32_t
#endif
#endif

third-party/openblas/OpenBLAS-0.3.27/Makefile.system
LAPACK_CFLAGS +=  -DLAPACK_ILP64

# 1 "/usr/include/x86_64-linux-gnu/bits/stdint-intn.h" 1 3 4
# 24 "/usr/include/x86_64-linux-gnu/bits/stdint-intn.h" 3 4
typedef __int8_t int8_t;
typedef __int16_t int16_t;
typedef __int32_t int32_t;
typedef __int64_t int64_t;

typedef unsigned char __u_char;
typedef unsigned short int __u_short;
typedef unsigned int __u_int;
typedef unsigned long int __u_long;typedef signed char __int8_t;
typedef unsigned char __uint8_t;
typedef signed short int __int16_t;
typedef unsigned short int __uint16_t;
typedef signed int __int32_t;
typedef unsigned int __uint32_t;typedef signed long int __int64_t;
typedef unsigned long int __uint64_t;