Given A, an m-by-n matrix with m >=n, and B a vector of size m, the routine DGELS (or S/D/C/Z-GELS) solves the least squares problem
minimize || B - A*X ||_2 over all vector X of size n.
(see
http://www.netlib.org/lapack/double/dgels.f.)
Multiple linear regression problems end up with a linear least squares problem so you can use LAPACK routine DGELS to eventually solve your problem.
Note that you are responsible for constructing the matrix A and B.
For example, if you simply want to solve the multiple linear regression problem
Y = a + b1*X1 + b2*X2 + ... + bp*Xp
then you set
A = [ ones, X1, X2, ... Xp ] , where A is a m-by-(p+1) matrix, ones is the vector of m ones, each column j of A represents the values of the variable Xj at the m different points
B = Y, B is a vector of size m with the values of the variable Y at the m different points
DGELS returns you the vector X of size (p+1) with the regression coefficient:
X = [ a, b1, b2, ... bp ] (in a column not in a line as written).