Using LU Decomposition in C to find the matrix determinant

LUDecomposition.c

//Full walkthrough: https://leetarxiv.substack.com/p/sinkhorn-knopp-algorithm
void PrintMatrix(int rows, int cols, double *matrix)
{
	for(int i = 0; i < rows; i++)
	{
		for(int j = 0; j < cols; j++)
		{
			double element = matrix[INDEX(i,j, cols)];
			printf("%.3f,",element);
		}
		printf("\n");	
	}
	printf("\n");	
}
double LUDecomposition(int n, double *inputMatrix, double *L, double *U, int *P)
{
	double det = 1.0;
	double *A = (double *)malloc(n * n * sizeof(double));
	memcpy(A, inputMatrix, n * n * sizeof(double));

	//Initialize L and U 
	for(int i = 0; i < n; i++){for(int j = 0; j < n; j++){L[i * n + j] = (i == j) ? 1.0 : 0.0;U[i * n + j] = 0.0;}}
	//Initialize P to identity
	for(int i = 0; i < n; i++){P[i] = i;}
    
	for(int i = 0; i < n; i++) 
	{
		//Partial pivoting
		double max = fabs(A[i * n + i]);
		int pivot = i;
		for(int k = i + 1; k < n; k++){double val = fabs(A[k * n + i]);if(val > max){max = val;pivot = k;}}

		if(max < 1e-12)
		{
			printf("WARNING: Singular or near-singular matrix at step %d\n", i);
			printf("Max pivot value: %e\n", max);
			free(A);
			return 0.0;
		}

		if(pivot != i) 
		{
			//Swap rows in A
			for(int j = 0; j < n; j++){double temp = A[i * n + j];A[i * n + j] = A[pivot * n + j];A[pivot * n + j] = temp;}
			//Swap rows in L
			for(int j = 0; j < i; j++){double temp = L[i * n + j];L[i * n + j] = L[pivot * n + j];L[pivot * n + j] = temp;}
			int temp = P[i];
			P[i] = P[pivot];
			P[pivot] = temp;
			det = -det;
		}

		//Compute U[i][j] with NaN check
		for(int j = i; j < n; j++) 
		{
			double sum = 0.0;
			for(int k = 0; k < i; k++){sum += L[i * n + k] * U[k * n + j];}
			U[i * n + j] = A[i * n + j] - sum;

			//CHECK FOR NaN/Inf
			if(!isfinite(U[i * n + j])) 
			{
				printf("ERROR: Non-finite U[%d][%d] = %f\n", i, j, U[i * n + j]);
				printf("A[%d][%d] = %f, sum = %f\n", i, j, A[i * n + j], sum);
				free(A);
				return NAN;
			}
		}

		// Compute L[j][i] with robustness check
		double divisor = U[i * n + i];
		if (fabs(divisor) < 1e-12)
		{
			printf("\nWARNING: Very small divisor at step %d: %e\n", i, divisor);
			free(A);
			return 0.0;
		}

		for(int j = i + 1; j < n; j++) 
		{
			double sum = 0.0;
			for(int k = 0; k < i; k++){sum += L[j * n + k] * U[k * n + i];}
			L[j * n + i] = (A[j * n + i] - sum) / divisor;

			//CHECK FOR NaN/Inf
			if(!isfinite(L[j * n + i]))
			{
				printf("\nERROR: Non-finite L[%d][%d] = %f\n", j, i, L[j * n + i]);
				printf("A[%d][%d] = %f, sum = %f, divisor = %f\n", j, i, A[j * n + i], sum, divisor);
				free(A);
				return NAN;
			}
		}
		det *= U[i * n + i];

		//CHECK DETERMINANT
		if(!isfinite(det))
		{
			printf("\nERROR: Non-finite determinant at step %d: %f\n", i, det);
			printf("U[%d][%d] = %f\n", i, i, U[i * n + i]);
			free(A);
			return det;
		}
	}

	free(A);
	return det;
}


void NaiveMatmul(int rowA, int colA, int colB,double *A,  double *B, double *C)
{
	for(int i = 0; i < rowA; i++)
	{
		for(int j = 0; j < colB; j++)
		{
			double sum = 0.0;
			for(int k = 0; k < colA; k++)
			{
				sum += A[INDEX(i, k, colA)] * B[INDEX(k, j, colB)];
			}
			C[INDEX(i, j, colB)] = sum;
		}
	}
}

void FindPAMatrix(int n, int *P, double *PA, double *A)
{
    for (int i = 0; i < n; i++)
    {
        int srcRow = P[i];  // Row in original A that became row i in the LU step
        for (int j = 0; j < n; j++)
        {
            PA[i * n + j] = A[srcRow * n + j];
        }
    }
}