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

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，结论

遇到非负整数，比如阶数、数组下标等，尽量用 unsigned long int 类型；