basic BLAS library only provides level 1 (vec-vec), 2 (mat-vec) and 3 (mat-mat) operations for dense matrices.
There is an extension to the BLAS lib called sparse BLAS which is available in intel MKL and some other implementations but not sure if it is freely available too,
BiCGSTAB can be implemented quite easily following the guidelines in the book
Templates for the Solution of Linear Systems,
you only need to populate your matrices using COO format (which is the easiest) and call the functions to change it to CRS or CCS and you only need a few levels one and 2 functions calls to finish the iterations. for the convergence tests there are again good suggestions in that book,
If you are looking for sophisticated parallel implementation I suggest you consider using software such as HYPRE or PETSc and use the linear solvers they provide. HYPER has both multigrid and krylov subspace type solvers but I have no experience with PETSc.