There are a lot of very good quality fortran routines around. It is silly to re-invent the wheel. Consider the case of solving Ax=b (for x) by least squares.

First go to Netlib (www.nelib.org). Then find linpack . In the list - search for "squares" and you'll find dqrsl.f - select the option to download that and it's dependencies (including BLAS - if you had optimized BLAS, for your machine, you could substitute those but you probably don't). This will create and archive, send it to your machine, and will unpack to a handful of fortran routines. The blurb mentions that this takes the factorization already performed by dqrdc - so go back to linpack main page and find that and download it's dependencies.

Next, compile all fortran rountines using
g77 -c name_of_routine.f
for each of the routines.

Now pack these object files into a library (don't have to but it is neater):
ar r test.lib file1.o file2.o ...... (for as many files as you compiled).
Now download my example file test_ls.c and compile (after READING) with
gcc -o test_ls.exe test_ls.c test.lib
Run and be amazed.

You can also use Visual C for the final compile and link:
cl test_ls.c test.lib

So endeth the lesson.