Reporting mistake in the GSL reference, release 2.7

[顾加银]

Dear GSL developers

     I am Jiayin Gu, a GSL user from China. I recently found maybe a mistake in 
the GSL reference, 2.7.  In the description of LU decomposition (Section 14.1). 
It states that

"The j-th column of the matrix P is given by the k-th column of the identity 
matrix, where k = p_j the j-th element of the permutation vector."

However, I tested this LU decomposition routine and found that the statement 
should be replaced by

"The j-th column of the identity matrix is placed as the k-th column of the 
matrix P, where k=p_j the j-th element of the permutation vector."




I tested this on Ubuntu 20.04 with g++ compiler. The GSL library is installed 
with the command "sudo apt install gsl-bin". Executing the command "sudo apt 
install --upgrade gsl-bin" gives  "gsl-bin is already the newest version 
(2.5+dfsg-6build1)." in the print. The attachment is the test code.




I hope the developers can fix this issue.

Best regards




Jiayin Gu
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>


int main(int argc, char *argv[])
{


	double A_data[]={0.18, 0.60, 0.57, 0.96, 0.41, 0.24, 0.99, 0.58, 0.14, 0.30, 0.97, 0.66, 0.51, 0.13, 0.19, 0.85};
	double AA_data[]={0.18, 0.60, 0.57, 0.96, 0.41, 0.24, 0.99, 0.58, 0.14, 0.30, 0.97, 0.66, 0.51, 0.13, 0.19, 0.85};

	gsl_matrix_view A=gsl_matrix_view_array(A_data, 4, 4);
	gsl_matrix_view AA=gsl_matrix_view_array(AA_data, 4, 4);


	int signum;
	gsl_permutation *perm=gsl_permutation_alloc(4);
	gsl_linalg_LU_decomp(&AA.matrix, perm, &signum); // this changes the value in AA array



	gsl_matrix *L=gsl_matrix_alloc(4, 4);
	gsl_matrix_set_zero(L);
	for(int i=0; i<4; i++)
	{
		gsl_matrix_set(L, i, i, 1);
		for(int j=0; j<i; j++)
		{
			double value=gsl_matrix_get(&AA.matrix, i, j);
			gsl_matrix_set(L, i, j, value);
		}
	}


	gsl_matrix *U=gsl_matrix_alloc(4, 4);
	gsl_matrix_set_zero(U);
	for(int i=0; i<4; i++)
	{
		for(int j=i; j<4; j++)
		{
			double value=gsl_matrix_get(&AA.matrix, i, j);
			gsl_matrix_set(U, i, j, value);
		}
	}


	gsl_matrix *P=gsl_matrix_alloc(4, 4);
	gsl_matrix_set_zero(P);
	for(int j=0; j<4; j++)
	{
		int k=gsl_permutation_get(perm, j);
		gsl_matrix_set(P, j, k, 1); // the j-th column of identity matrix is placed in the k-th column of P matrix; there is a mistake in the reference
	}


	gsl_matrix *mat=gsl_matrix_alloc(4, 4);


	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, L, U, 0.0, mat);
	std::cout<<gsl_matrix_get(mat, 3, 0)<<std::endl;


	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, P, &A.matrix, 0.0, mat);
	std::cout<<gsl_matrix_get(mat, 3, 0)<<std::endl;


	gsl_permutation_free(perm);
	gsl_matrix_free(L);
	gsl_matrix_free(U);
	gsl_matrix_free(mat);


	return 0;
}