I have a code deeply embedded with GNU Scientific Library (GSL) matrix arithmetic, the main computation of this code is solving a large system of linear equations that takes a very long time in serial and with GSL and BLAS functions, is there a way to parallelize this computation or convert it for use in an already parallel library like ScaLAPACK?
If your matrix is sparse, i.e. it contains a lot of zero entries, then you can easily implement many sparse matrix algebra packages without too much trouble. Unfortunately this will require you to store your matrices in sparse format which, to my knowledge, gsl does not do. Once you have your matrix stored in some sparse format, you should be able to handle large systems without too much trouble, even in serial applications.
I suggest using UMFPACK because it requires the least amount of work to implement as it doesn't require you to put your data into their structures.
A note on paralleism: If your code is currently serial, going to a parallel solver is NOT trivial. It is possible that it may be simple to implement a multi-threaded package, but I don't have much experience with threaded programs. Additionally, truly parallel (distributed memory) direct solvers are not all that efficient, since each processor needs its own copy of the full matrix, and it is better to use iterative methods.
A little more detail would be helpful: How long is a long time? Do you need the inverse for some reason, or are you just solving a system of equations?