diff --git a/NEWS b/NEWS
index bee4c00..2c79814 100644
--- a/NEWS
+++ b/NEWS
@@ -1,10 +1,23 @@
-Changes in 1.2.0:
-~~~~~~~~~~~~~~~~~
-This branch introduces support for gsl 2.3. An overview:
+Changes in 1.2.0 (by R Bader):
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+This release introduces support for gsl 2.3. An overview:
 
 1. added support for the multifit_nlinear and multilarge_nlinear interfaces, 
    including example programs
-
+2. added fgsl_rstat_rms, fgsl_rstat_quantile_reset to rstat interface
+3. updated the linalg interface with additional Cholesky routines 
+   (fgsl_linalg_cholesky_decomp1, fgsl_linalg_cholesky_rcond, 
+   fgsl_linalg_cholesky_*2, fgsl_linalg_cholesky_scale*, 
+   fgsl_linalg_pcholesky_*, fgsl_linalg_mcholesky_*), with 
+   additional QRPT routines, with the linalg_tri triangular matrix routines, 
+   with the linalg_cod complete orthogonal decomposition routines
+4. added fgsl_permute_matrix and test case
+5. updated the randist interface with the multivariate Gaussian distribution
+6. updated the multifit_linear interface with fgsl_multifit_linear_rank,  
+   fgsl_multifit_linear_gcv* and the fgsl_multifit_*linear_tsvd routines.
+7. Travis CI configuration updated by T. Schoonjans
+8. added fgsl_spblas_dgem[v,m] to sparse matrix routines
+9. minor changes to work around compiler bugs
 
 
 Changes in 1.1.0 (by Tom Schoonjans):
diff --git a/README b/README
index 9a76b65..838fd84 100644
--- a/README
+++ b/README
@@ -124,4 +124,4 @@ Releases:
   * 1.0.0: February 11, 2014 (revision 34)
   --- Migrated to Github ---
   * 1.1.0: February, 2016 
-  
+  * 1.2.0: January, 2017
diff --git a/README_WORKFLOW b/README_WORKFLOW
index 2441509..2f84234 100644
--- a/README_WORKFLOW
+++ b/README_WORKFLOW
@@ -1,6 +1,8 @@
 Some notes on the source management workflow
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+(A) Trivial changes
+
 * generate local clone of repo
   git clone https://github.com/reinh-bader/fgsl.git
   produces a folder fgsl that contains the clone
@@ -11,3 +13,22 @@ Some notes on the source management workflow
   git commit -m "<commit message>" 
 * push updated repository back to master
   git push
+  
+
+(B) Non-trivial changes
+
+* generate local clone of repo
+  git clone https://github.com/reinh-bader/fgsl.git
+* create a new branch and check it out:
+  git checkout -b fgsl_devel_<new_version>
+  (this can be done even if some files have been modified in advance)
+* do some work
+* commit
+* push the new branch with the commit to github and track it
+  git push -u origin fgsl_devel_<new_version>
+* open a pull-request for this branch in the repository, and check the results of 
+  the Travis-CI builds: if one of them fails, fix the underlying cause in the code, 
+  commit and push, thereby triggering a new build
+* if needed, keep making changes to the code, commit them and check the build results.
+* when done, merge the PR in and delete the local and remote working branch.
+
diff --git a/api/linalg.finc b/api/linalg.finc
index ccc99d4..33441a5 100644
--- a/api/linalg.finc
+++ b/api/linalg.finc
@@ -283,6 +283,27 @@
     fgsl_linalg_qrpt_svx = gsl_linalg_qrpt_svx(qr%gsl_matrix, &
          tau%gsl_vector, p%gsl_permutation, x%gsl_vector)
   end function fgsl_linalg_qrpt_svx
+  function fgsl_linalg_qrpt_lssolve (qr, tau, p, b, x, residual)
+    type(fgsl_matrix), intent(in) :: qr
+    type(fgsl_vector), intent(in) :: tau, b
+    type(fgsl_permutation), intent(in) :: p
+    type(fgsl_vector), intent(inout) :: x, residual
+    integer(fgsl_int) :: fgsl_linalg_qrpt_lssolve
+    fgsl_linalg_qrpt_lssolve = gsl_linalg_qrpt_lssolve(qr%gsl_matrix, &
+         tau%gsl_vector, p%gsl_permutation, b%gsl_vector, &
+         x%gsl_vector, residual%gsl_vector)
+  end function fgsl_linalg_qrpt_lssolve
+  function fgsl_linalg_qrpt_lssolve2 (qr, tau, p, b, rank, x, residual)
+    type(fgsl_matrix), intent(in) :: qr
+    type(fgsl_vector), intent(in) :: tau, b
+    integer(fgsl_size_t), intent(in) :: rank
+    type(fgsl_permutation), intent(in) :: p
+    type(fgsl_vector), intent(inout) :: x, residual
+    integer(fgsl_int) :: fgsl_linalg_qrpt_lssolve2
+    fgsl_linalg_qrpt_lssolve2 = gsl_linalg_qrpt_lssolve2(qr%gsl_matrix, &
+         tau%gsl_vector, p%gsl_permutation, b%gsl_vector, rank, &
+         x%gsl_vector, residual%gsl_vector)
+  end function fgsl_linalg_qrpt_lssolve2
   function fgsl_linalg_qrpt_qrsolve (q, r, p, b, x)
     type(fgsl_matrix), intent(in) :: q, r
     type(fgsl_vector), intent(in) :: b
@@ -318,6 +339,67 @@
     fgsl_linalg_qrpt_rsvx = gsl_linalg_qrpt_rsvx(qr%gsl_matrix, &
          p%gsl_permutation, x%gsl_vector)
   end function fgsl_linalg_qrpt_rsvx
+  function fgsl_linalg_qrpt_rank (qr, tol)
+    type(fgsl_matrix), intent(in) :: qr
+    real(fgsl_double), intent(in) :: tol
+    integer(fgsl_size_t) :: fgsl_linalg_qrpt_rank
+    fgsl_linalg_qrpt_rank = gsl_linalg_qrpt_rank(qr%gsl_matrix, tol)
+  end function fgsl_linalg_qrpt_rank
+  function fgsl_linalg_qrpt_rcond (qr, rcond, work)
+    type(fgsl_matrix), intent(in) :: qr
+    real(fgsl_double), intent(inout) :: rcond
+    type(fgsl_vector), intent(inout) :: work
+    integer(fgsl_int) :: fgsl_linalg_qrpt_rcond
+    fgsl_linalg_qrpt_rcond = gsl_linalg_qrpt_rcond(qr%gsl_matrix, rcond, work%gsl_vector)
+  end function fgsl_linalg_qrpt_rcond
+!
+! Complete orthogonal decomposition (COD)
+!
+  integer(fgsl_int) function fgsl_linalg_cod_decomp(a, tau_q, tau_z, p, rank, work)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_vector), intent(inout) :: tau_q, tau_z, work
+    type(fgsl_permutation), intent(inout) :: p
+    integer(fgsl_size_t), intent(inout) :: rank
+    fgsl_linalg_cod_decomp = gsl_linalg_cod_decomp(a%gsl_matrix, tau_q%gsl_vector, &
+          tau_z%gsl_vector, p%gsl_permutation, rank, work%gsl_vector)
+  end function
+  integer(fgsl_int) function fgsl_linalg_cod_decomp_e(a, tau_q, tau_z, p, tol, rank, work)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_vector), intent(inout) :: tau_q, tau_z, work
+    real(fgsl_double), intent(in) :: tol
+    type(fgsl_permutation), intent(inout) :: p
+    integer(fgsl_size_t), intent(inout) :: rank
+    fgsl_linalg_cod_decomp_e = gsl_linalg_cod_decomp_e(a%gsl_matrix, tau_q%gsl_vector, &
+          tau_z%gsl_vector, p%gsl_permutation, tol, rank, work%gsl_vector)
+  end function
+  integer(fgsl_int) function fgsl_linalg_cod_lssolve(qrz, tau_q, tau_z, p, rank, b, x, residual)
+    type(fgsl_matrix), intent(in)  :: qrz
+    type(fgsl_vector), intent(in) :: tau_q, tau_z, b
+    type(fgsl_vector), intent(inout) :: x, residual
+    type(fgsl_permutation), intent(in) :: p
+    integer(fgsl_size_t), intent(in) :: rank
+    fgsl_linalg_cod_lssolve = gsl_linalg_cod_lssolve(qrz%gsl_matrix, tau_q%gsl_vector, &
+          tau_z%gsl_vector, p%gsl_permutation, rank, b%gsl_vector, x%gsl_vector, &
+          residual%gsl_vector)
+  end function
+  integer(fgsl_int) function fgsl_linalg_cod_unpack(qrz, tau_q, tau_z, p, rank, q, r, z)
+    type(fgsl_matrix), intent(in)  :: qrz
+    type(fgsl_vector), intent(in) :: tau_q, tau_z
+    type(fgsl_permutation), intent(in) :: p
+    integer(fgsl_size_t), intent(in) :: rank
+    type(fgsl_matrix), intent(inout)  :: q, r, z
+    fgsl_linalg_cod_unpack = gsl_linalg_cod_unpack(qrz%gsl_matrix, tau_q%gsl_vector, &
+          tau_z%gsl_vector, p%gsl_permutation, rank, q%gsl_matrix, r%gsl_matrix, z%gsl_matrix)
+  end function
+  integer(fgsl_int) function fgsl_linalg_cod_matz(qrz, tau_z, rank, a, work)
+    type(fgsl_matrix), intent(in)  :: qrz
+    type(fgsl_vector), intent(in) ::  tau_z
+    integer(fgsl_size_t), intent(in) :: rank
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_vector), intent(inout) ::  work
+    fgsl_linalg_cod_matz = gsl_linalg_cod_matz(qrz%gsl_matrix, tau_z%gsl_vector, &
+          rank, a%gsl_matrix, work%gsl_vector)
+  end function
 !
 ! SVD
 !
@@ -359,6 +441,11 @@
 !
 ! Cholesky
 !
+  function fgsl_linalg_cholesky_decomp1(a)
+    type(fgsl_matrix), intent(inout)  :: a
+    integer(fgsl_int) :: fgsl_linalg_cholesky_decomp1
+    fgsl_linalg_cholesky_decomp1 = gsl_linalg_cholesky_decomp1(a%gsl_matrix)
+  end function fgsl_linalg_cholesky_decomp1
   function fgsl_linalg_cholesky_decomp(a)
     type(fgsl_matrix), intent(inout)  :: a
     integer(fgsl_int) :: fgsl_linalg_cholesky_decomp
@@ -402,6 +489,30 @@
          gsl_linalg_complex_cholesky_svx(chol%gsl_matrix_complex, &
          x%gsl_vector_complex)
   end function fgsl_linalg_complex_cholesky_svx
+  function fgsl_linalg_cholesky_decomp2(a,s)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_vector), intent(inout) :: s
+    integer(fgsl_int) :: fgsl_linalg_cholesky_decomp2
+    fgsl_linalg_cholesky_decomp2 = gsl_linalg_cholesky_decomp2(a%gsl_matrix, &
+         s%gsl_vector)
+  end function fgsl_linalg_cholesky_decomp2
+  function fgsl_linalg_cholesky_solve2(chol, s, b, x)
+    type(fgsl_matrix), intent(in) :: chol
+    type(fgsl_vector), intent(in) :: s
+    type(fgsl_vector), intent(in) :: b
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_cholesky_solve2
+    fgsl_linalg_cholesky_solve2 = gsl_linalg_cholesky_solve2(chol%gsl_matrix, &
+         s%gsl_vector, b%gsl_vector, x%gsl_vector)
+  end function fgsl_linalg_cholesky_solve2
+  function fgsl_linalg_cholesky_svx2(chol, s, x)
+    type(fgsl_matrix), intent(in) :: chol
+    type(fgsl_vector), intent(in) :: s
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_cholesky_svx2
+    fgsl_linalg_cholesky_svx2 = gsl_linalg_cholesky_svx2(chol%gsl_matrix, &
+         s%gsl_vector, x%gsl_vector)
+  end function fgsl_linalg_cholesky_svx2
   function fgsl_linalg_cholesky_invert (chol)
     integer(fgsl_int) :: fgsl_linalg_cholesky_invert
     type(fgsl_matrix), intent(inout) :: chol
@@ -412,6 +523,143 @@
     type(fgsl_matrix_complex), intent(inout) :: chol
     fgsl_linalg_complex_cholesky_invert = gsl_linalg_complex_cholesky_invert(chol%gsl_matrix_complex)
   end function fgsl_linalg_complex_cholesky_invert
+  function fgsl_linalg_cholesky_scale(a,s)
+    type(fgsl_matrix), intent(in)  :: a
+    type(fgsl_vector), intent(inout) :: s
+    integer(fgsl_int) :: fgsl_linalg_cholesky_scale
+    fgsl_linalg_cholesky_scale = gsl_linalg_cholesky_scale(a%gsl_matrix, &
+         s%gsl_vector)
+  end function fgsl_linalg_cholesky_scale
+  function fgsl_linalg_cholesky_scale_apply(a,s)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_vector), intent(in) :: s
+    integer(fgsl_int) :: fgsl_linalg_cholesky_scale_apply
+    fgsl_linalg_cholesky_scale_apply = gsl_linalg_cholesky_scale_apply(a%gsl_matrix, &
+         s%gsl_vector)
+  end function fgsl_linalg_cholesky_scale_apply
+  function fgsl_linalg_cholesky_rcond (chol, rcond, work)
+    integer(fgsl_int) :: fgsl_linalg_cholesky_rcond
+    type(fgsl_matrix), intent(in) :: chol
+    real(fgsl_double), intent(inout) :: rcond
+    type(fgsl_vector), intent(inout) :: work
+    fgsl_linalg_cholesky_rcond = gsl_linalg_cholesky_rcond(chol%gsl_matrix, rcond, &
+         work%gsl_vector)
+  end function fgsl_linalg_cholesky_rcond
+!
+! pivoted Cholesky
+!
+  function fgsl_linalg_pcholesky_decomp(a,p)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_permutation), intent(inout)  :: p
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_decomp
+    fgsl_linalg_pcholesky_decomp = gsl_linalg_pcholesky_decomp(a%gsl_matrix, p%gsl_permutation)
+  end function fgsl_linalg_pcholesky_decomp
+  function fgsl_linalg_pcholesky_solve(ldlt, p, b, x)
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in)  :: p
+    type(fgsl_vector), intent(in) :: b
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_solve
+    fgsl_linalg_pcholesky_solve = gsl_linalg_pcholesky_solve(ldlt%gsl_matrix, &
+         p%gsl_permutation, b%gsl_vector, x%gsl_vector)
+  end function fgsl_linalg_pcholesky_solve
+  function fgsl_linalg_pcholesky_svx(ldlt, p, x)
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in)  :: p
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_svx
+    fgsl_linalg_pcholesky_svx = gsl_linalg_pcholesky_svx(ldlt%gsl_matrix, &
+         p%gsl_permutation, x%gsl_vector)
+  end function fgsl_linalg_pcholesky_svx
+  function fgsl_linalg_pcholesky_decomp2(a,p,s)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_permutation), intent(inout)  :: p
+    type(fgsl_vector), intent(inout) :: s
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_decomp2
+    fgsl_linalg_pcholesky_decomp2 = gsl_linalg_pcholesky_decomp2(a%gsl_matrix, &
+         p%gsl_permutation, s%gsl_vector)
+  end function fgsl_linalg_pcholesky_decomp2
+  function fgsl_linalg_pcholesky_solve2(ldlt, p, s, b, x)
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in)  :: p
+    type(fgsl_vector), intent(in) :: s
+    type(fgsl_vector), intent(in) :: b
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_solve2
+    fgsl_linalg_pcholesky_solve2 = gsl_linalg_pcholesky_solve2(ldlt%gsl_matrix, &
+         p%gsl_permutation, s%gsl_vector, b%gsl_vector, x%gsl_vector)
+  end function fgsl_linalg_pcholesky_solve2
+  function fgsl_linalg_pcholesky_svx2(ldlt, p, s, x)
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in)  :: p
+    type(fgsl_vector), intent(in) :: s
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_svx2
+    fgsl_linalg_pcholesky_svx2 = gsl_linalg_pcholesky_svx2(ldlt%gsl_matrix, &
+         p%gsl_permutation, s%gsl_vector, x%gsl_vector)
+  end function fgsl_linalg_pcholesky_svx2
+  function fgsl_linalg_pcholesky_invert (ldlt, p, ainv)
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_invert
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in) :: p
+    type(fgsl_matrix), intent(inout) :: ainv
+    fgsl_linalg_pcholesky_invert = gsl_linalg_pcholesky_invert(ldlt%gsl_matrix, &
+         p%gsl_permutation, ainv%gsl_matrix)
+  end function fgsl_linalg_pcholesky_invert
+  function fgsl_linalg_pcholesky_rcond (ldlt, p, rcond, work)
+    integer(fgsl_int) :: fgsl_linalg_pcholesky_rcond
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in) :: p
+    real(fgsl_double), intent(inout) :: rcond
+    type(fgsl_vector), intent(inout) :: work
+    fgsl_linalg_pcholesky_rcond = gsl_linalg_pcholesky_rcond(ldlt%gsl_matrix, &
+         p%gsl_permutation, rcond, work%gsl_vector)
+  end function fgsl_linalg_pcholesky_rcond
+!
+! modified Cholesky
+!
+  function fgsl_linalg_mcholesky_decomp(a,p,e)
+    type(fgsl_matrix), intent(inout)  :: a
+    type(fgsl_permutation), intent(inout)  :: p
+    type(fgsl_vector), intent(inout)  :: e
+    integer(fgsl_int) :: fgsl_linalg_mcholesky_decomp
+    fgsl_linalg_mcholesky_decomp = gsl_linalg_mcholesky_decomp(a%gsl_matrix, &
+         p%gsl_permutation, e%gsl_vector)
+  end function fgsl_linalg_mcholesky_decomp
+  function fgsl_linalg_mcholesky_solve(ldlt, p, b, x)
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in)  :: p
+    type(fgsl_vector), intent(in) :: b
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_mcholesky_solve
+    fgsl_linalg_mcholesky_solve = gsl_linalg_mcholesky_solve(ldlt%gsl_matrix, &
+         p%gsl_permutation, b%gsl_vector, x%gsl_vector)
+  end function fgsl_linalg_mcholesky_solve
+  function fgsl_linalg_mcholesky_svx(ldlt, p, x)
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in)  :: p
+    type(fgsl_vector), intent(inout) :: x
+    integer(fgsl_int) :: fgsl_linalg_mcholesky_svx
+    fgsl_linalg_mcholesky_svx = gsl_linalg_mcholesky_svx(ldlt%gsl_matrix, &
+         p%gsl_permutation, x%gsl_vector)
+  end function fgsl_linalg_mcholesky_svx
+  function fgsl_linalg_mcholesky_invert (ldlt, p, ainv)
+    integer(fgsl_int) :: fgsl_linalg_mcholesky_invert
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in) :: p
+    type(fgsl_matrix), intent(inout) :: ainv
+    fgsl_linalg_mcholesky_invert = gsl_linalg_mcholesky_invert(ldlt%gsl_matrix, &
+         p%gsl_permutation, ainv%gsl_matrix)
+  end function fgsl_linalg_mcholesky_invert
+  function fgsl_linalg_mcholesky_rcond (ldlt, p, rcond, work)
+    integer(fgsl_int) :: fgsl_linalg_mcholesky_rcond
+    type(fgsl_matrix), intent(in) :: ldlt
+    type(fgsl_permutation), intent(in) :: p
+    real(fgsl_double), intent(inout) :: rcond
+    type(fgsl_vector), intent(inout) :: work
+    fgsl_linalg_mcholesky_rcond = gsl_linalg_mcholesky_rcond(ldlt%gsl_matrix, &
+         p%gsl_permutation, rcond, work%gsl_vector)
+  end function fgsl_linalg_mcholesky_rcond
 !
 ! Tridiag
 !
@@ -686,3 +934,37 @@
     real(fgsl_double), intent(in) :: c, s
     call gsl_linalg_givens_gv(v%gsl_vector, i, j, c, s)
   end subroutine fgsl_linalg_givens_gv
+!
+! Triangular
+!
+  integer(fgsl_int) function fgsl_linalg_tri_upper_invert(t)
+    type(fgsl_matrix), intent(inout) :: t
+    fgsl_linalg_tri_upper_invert = gsl_linalg_tri_upper_invert(t%gsl_matrix)
+  end function
+  integer(fgsl_int) function fgsl_linalg_tri_lower_invert(t)
+    type(fgsl_matrix), intent(inout) :: t
+    fgsl_linalg_tri_lower_invert = gsl_linalg_tri_lower_invert(t%gsl_matrix)
+  end function
+  integer(fgsl_int) function fgsl_linalg_tri_upper_unit_invert(t)
+    type(fgsl_matrix), intent(inout) :: t
+    fgsl_linalg_tri_upper_unit_invert = gsl_linalg_tri_upper_unit_invert(t%gsl_matrix)
+  end function
+  integer(fgsl_int) function fgsl_linalg_tri_lower_unit_invert(t)
+    type(fgsl_matrix), intent(inout) :: t
+    fgsl_linalg_tri_lower_unit_invert = gsl_linalg_tri_lower_unit_invert(t%gsl_matrix)
+  end function
+  integer(fgsl_int) function fgsl_linalg_tri_upper_rcond(t, rcond, work)
+    type(fgsl_matrix), intent(inout) :: t
+    real(fgsl_double), intent(inout) :: rcond
+    type(fgsl_vector), intent(inout) :: work
+    fgsl_linalg_tri_upper_rcond = gsl_linalg_tri_upper_rcond(t%gsl_matrix, rcond, &
+         work%gsl_vector)
+  end function
+  integer(fgsl_int) function fgsl_linalg_tri_lower_rcond(t, rcond, work)
+    type(fgsl_matrix), intent(inout) :: t
+    real(fgsl_double), intent(inout) :: rcond
+    type(fgsl_vector), intent(inout) :: work
+    fgsl_linalg_tri_lower_rcond = gsl_linalg_tri_lower_rcond(t%gsl_matrix, rcond, &
+         work%gsl_vector)
+  end function
+  
diff --git a/api/multifit.finc b/api/multifit.finc
index 28d5094..8f4650e 100644
--- a/api/multifit.finc
+++ b/api/multifit.finc
@@ -404,6 +404,19 @@
     fgsl_multifit_linear = gsl_multifit_linear(x%gsl_matrix, y%gsl_vector, &
          c%gsl_vector, cov%gsl_matrix, chisq, work%gsl_multifit_linear_workspace)
   end function fgsl_multifit_linear
+  function fgsl_multifit_linear_tsvd(x, y, tol, c, cov, chisq, rank, work)
+    type(fgsl_matrix), intent(in) :: x
+    type(fgsl_vector), intent(in) :: y
+    real(fgsl_double), intent(in) :: tol
+    type(fgsl_vector), intent(inout) :: c
+    type(fgsl_matrix), intent(inout) :: cov
+    integer(fgsl_size_t), intent(inout) :: rank
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    real(fgsl_double), intent(inout) :: chisq
+    integer(fgsl_int) :: fgsl_multifit_linear_tsvd
+    fgsl_multifit_linear_tsvd = gsl_multifit_linear_tsvd(x%gsl_matrix, y%gsl_vector, tol, &
+         c%gsl_vector, cov%gsl_matrix, chisq, rank, work%gsl_multifit_linear_workspace)
+  end function fgsl_multifit_linear_tsvd
   function fgsl_multifit_linear_svd(x, work)
     type(fgsl_matrix), intent(in) :: x
     type(fgsl_multifit_linear_workspace), intent(inout) :: work
@@ -550,6 +563,46 @@
     fgsl_multifit_linear_lcorner2 = gsl_multifit_linear_lcorner2(&
     reg_param%gsl_vector, eta%gsl_vector, idx)
   end function fgsl_multifit_linear_lcorner2
+  integer(fgsl_int) function fgsl_multifit_linear_gcv_init(y, reg_param, uty, delta0, work)
+    type(fgsl_vector), intent(in) :: y
+    type(fgsl_vector), intent(inout) :: reg_param, uty
+    real(fgsl_double), intent(inout) :: delta0
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    fgsl_multifit_linear_gcv_init = gsl_multifit_linear_gcv_init(y%gsl_vector, &
+          reg_param%gsl_vector, uty%gsl_vector, delta0, work%gsl_multifit_linear_workspace)
+  end function
+  integer(fgsl_int) function fgsl_multifit_linear_gcv_curve(reg_param, uty, delta0, g, work)
+    type(fgsl_vector), intent(in) :: reg_param, uty
+    type(fgsl_vector), intent(inout) :: g
+    real(fgsl_double), intent(in) :: delta0
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    fgsl_multifit_linear_gcv_curve = gsl_multifit_linear_gcv_curve(reg_param%gsl_vector, &
+          uty%gsl_vector, delta0, g%gsl_vector, work%gsl_multifit_linear_workspace)
+  end function
+  integer(fgsl_int) function fgsl_multifit_linear_gcv_min(reg_param, uty, delta0, g, & 
+          lambda, work)
+    type(fgsl_vector), intent(in) :: reg_param, uty, g
+    real(fgsl_double), intent(in) :: delta0
+    real(fgsl_double), intent(inout) :: lambda
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    fgsl_multifit_linear_gcv_min = gsl_multifit_linear_gcv_min(reg_param%gsl_vector, &
+          uty%gsl_vector, delta0, g%gsl_vector, lambda, work%gsl_multifit_linear_workspace)
+  end function
+  real(fgsl_double) function fgsl_multifit_linear_gcv_calc(lambda, uty, delta0, work)
+    type(fgsl_vector), intent(in) :: uty
+    real(fgsl_double), intent(in) :: delta0, lambda
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    fgsl_multifit_linear_gcv_calc = gsl_multifit_linear_gcv_calc(lambda, &
+          uty%gsl_vector, delta0, work%gsl_multifit_linear_workspace)
+  end function
+  integer(fgsl_int) function fgsl_multifit_linear_gcv(y, reg_param, g, lambda, g_lambda, work)
+    type(fgsl_vector), intent(in) :: y
+    type(fgsl_vector), intent(inout) :: reg_param, g
+    real(fgsl_double), intent(inout) :: lambda, g_lambda
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    fgsl_multifit_linear_gcv = gsl_multifit_linear_gcv(y%gsl_vector, reg_param%gsl_vector, &
+          g%gsl_vector, lambda, g_lambda, work%gsl_multifit_linear_workspace)
+  end function
   function fgsl_multifit_linear_lk(p, k, l)
     integer(fgsl_size_t), intent(in) :: p, k
     type(fgsl_matrix), intent(inout) :: l
@@ -608,6 +661,20 @@
     fgsl_multifit_wlinear = gsl_multifit_wlinear(x%gsl_matrix, w%gsl_vector, y%gsl_vector, &
          c%gsl_vector, cov%gsl_matrix, chisq, work%gsl_multifit_linear_workspace)
   end function fgsl_multifit_wlinear
+  function fgsl_multifit_wlinear_tsvd(x, w, y, tol, c, cov, chisq, rank, work)
+    type(fgsl_matrix), intent(in) :: x
+    type(fgsl_vector), intent(in) :: w, y
+    real(fgsl_double), intent(in) :: tol
+    type(fgsl_vector), intent(inout) :: c
+    type(fgsl_matrix), intent(inout) :: cov
+    integer(fgsl_size_t), intent(inout) :: rank
+    type(fgsl_multifit_linear_workspace), intent(inout) :: work
+    real(fgsl_double), intent(inout) :: chisq
+    integer(fgsl_int) :: fgsl_multifit_wlinear_tsvd
+    fgsl_multifit_wlinear_tsvd = gsl_multifit_wlinear_tsvd( &
+         x%gsl_matrix, w%gsl_vector, y%gsl_vector, tol, &
+         c%gsl_vector, cov%gsl_matrix, chisq, rank, work%gsl_multifit_linear_workspace)
+  end function fgsl_multifit_wlinear_tsvd
   function fgsl_multifit_wlinear_svd(x, w, y, tol, rank, c, cov, chisq, work)
     type(fgsl_matrix), intent(in) :: x
     type(fgsl_vector), intent(in) :: w, y
@@ -652,6 +719,11 @@
     fgsl_multifit_linear_residuals = gsl_multifit_linear_residuals( &
          x%gsl_matrix, y%gsl_vector, c%gsl_vector, r%gsl_vector)
   end function fgsl_multifit_linear_residuals
+  integer(fgsl_size_t) function fgsl_multifit_linear_rank(tol, work)
+    real(fgsl_double), intent(in) :: tol
+    type(fgsl_multifit_linear_workspace), intent(in) :: work
+    fgsl_multifit_linear_rank = gsl_multifit_linear_rank(tol, work%gsl_multifit_linear_workspace)
+  end function
   function fgsl_multifit_status(multifit)
     type(fgsl_multifit_linear_workspace), intent(in) :: multifit
     logical :: fgsl_multifit_status
diff --git a/api/nlfit.finc b/api/nlfit.finc
index 89eb251..98cf2dc 100644
--- a/api/nlfit.finc
+++ b/api/nlfit.finc
@@ -10,6 +10,10 @@
     character(kind=fgsl_char, len=*) :: s 
     fgsl_multifit_nlinear_setup%gsl_multifit_nlinear_type = gsl_multifit_nlinear_setup(s(1:1))
   end function
+  type(fgsl_multilarge_nlinear_type) function fgsl_multilarge_nlinear_setup(s)
+    character(kind=fgsl_char, len=*) :: s
+    fgsl_multilarge_nlinear_setup%gsl_multilarge_nlinear_type = gsl_multilarge_nlinear_setup(s(1:1))
+  end function
   function fgsl_multifit_nlinear_alloc(t, params, n, p)
     type(fgsl_multifit_nlinear_type), intent(in) :: t
     type(fgsl_multifit_nlinear_parameters), intent(in) :: params
@@ -305,3 +309,77 @@
       if (present(h_fvv)) p%h_fvv = h_fvv
     end associate
   end subroutine
+  function fgsl_multilarge_nlinear_fdf_init(ndim, p, params, func, dfunc, fvv)
+    procedure(fgsl_nlinear_fdf_func), optional :: func
+    procedure(fgsl_nlinear_fdf_dlfunc), optional :: dfunc
+    procedure(fgsl_nlinear_fdf_fvv), optional :: fvv
+    integer(fgsl_size_t), intent(in) :: ndim, p
+    type(c_ptr), intent(in) :: params
+    type(fgsl_multilarge_nlinear_fdf) :: fgsl_multilarge_nlinear_fdf_init
+!
+    type(c_funptr) :: fp, dfp, fvvp
+    fp = c_funloc(func)
+    if (present(dfunc)) then
+       dfp = c_funloc(dfunc)
+    else
+       dfp = c_null_funptr
+    end if
+    if (present(fvv)) then
+       fvvp = c_funloc(fvv)
+    else
+       fvvp = c_null_funptr
+    end if
+    fgsl_multilarge_nlinear_fdf_init%gsl_multilarge_nlinear_fdf = &
+         fgsl_multilarge_nlinear_fdf_cinit(ndim, p, params, fp, dfp, fvvp)
+  end function fgsl_multilarge_nlinear_fdf_init
+  subroutine fgsl_multilarge_nlinear_fdf_free(fun)
+    type(fgsl_multilarge_nlinear_fdf), intent(inout) :: fun
+    call fgsl_multilarge_nlinear_fdf_cfree(fun%gsl_multilarge_nlinear_fdf)
+  end subroutine fgsl_multilarge_nlinear_fdf_free
+  subroutine fgsl_multilarge_nlinear_fdf_get(fdf, func, dfunc, fvv, n, p, params, &
+         nevalf, nevaldfu, nevaldf2, nevalfvv)
+    type(fgsl_multilarge_nlinear_fdf), intent(in) :: fdf
+    procedure(fgsl_nlinear_fdf_func), pointer, optional :: func
+    procedure(fgsl_nlinear_fdf_dlfunc), pointer, optional :: dfunc
+    procedure(fgsl_nlinear_fdf_fvv), pointer, optional :: fvv
+    integer(fgsl_size_t), intent(out), optional :: n, p, nevalf, nevaldfu, nevaldf2, nevalfvv
+    type(c_ptr), intent(out), optional :: params
+    type(c_funptr) :: pfunc, pdfunc, pfvv
+    integer(fgsl_size_t) :: nl, pl, nevalfl, nevaldful, nevaldf2l, nevalfvvl
+    type(c_ptr) :: paramsl
+    call gsl_multilarge_nlinear_fdf_get(fdf%gsl_multilarge_nlinear_fdf, pfunc, pdfunc, pfvv, &
+          nl, pl, paramsl, nevalfl, nevaldful, nevaldf2l, nevalfvvl)
+    if (present(func)) call c_f_procpointer(pfunc, func)
+    if (present(dfunc)) call c_f_procpointer(pdfunc, dfunc)
+    if (present(fvv)) call c_f_procpointer(pfvv, fvv)
+    if (present(n)) n = nl
+    if (present(p)) p = pl
+    if (present(params)) params = paramsl
+    if (present(nevalf)) nevalf = nevalfl
+    if (present(nevaldfu)) nevaldfu = nevaldful
+    if (present(nevaldf2)) nevaldf2 = nevaldf2l
+    if (present(nevalfvv)) nevalfvv = nevalfvvl
+  end subroutine fgsl_multilarge_nlinear_fdf_get
+  subroutine fgsl_multilarge_nlinear_parameters_set(params, trs, scale, solver, &
+         fdtype, factor_up, factor_down, avmax, h_df, h_fvv, max_iter, tol)
+    type(fgsl_multilarge_nlinear_parameters) :: params
+    type(fgsl_multilarge_nlinear_trs), optional :: trs
+    type(fgsl_multilarge_nlinear_scale), optional :: scale
+    type(fgsl_multilarge_nlinear_solver), optional :: solver
+    integer(fgsl_int), optional :: fdtype
+    integer(fgsl_size_t), optional :: max_iter
+    real(c_double), optional :: factor_up, factor_down, avmax, h_df, h_fvv, tol
+    associate(p => params%gsl_multilarge_nlinear_parameters)
+      if (present(trs)) p%trs = gsl_multilarge_nlinear_get_trs(trs%which)
+      if (present(scale)) p%scale = gsl_multilarge_nlinear_get_scale(scale%which)
+      if (present(solver)) p%solver = gsl_multilarge_nlinear_get_solver(solver%which)
+      if (present(fdtype)) p%fdtype = fdtype
+      if (present(factor_up)) p%factor_up = factor_up
+      if (present(factor_down)) p%factor_down = factor_down
+      if (present(avmax)) p%avmax = avmax
+      if (present(h_df)) p%h_df = h_df
+      if (present(h_fvv)) p%h_fvv = h_fvv
+      if (present(max_iter)) p%max_iter = max_iter
+      if (present(tol)) p%tol = tol
+    end associate
+  end subroutine
diff --git a/api/permutation.finc b/api/permutation.finc
index 4fb50cf..42c90bd 100644
--- a/api/permutation.finc
+++ b/api/permutation.finc
@@ -129,6 +129,13 @@
     fgsl_permute_vector_inverse = &
          gsl_permute_vector_inverse(p%gsl_permutation,v%gsl_vector)
   end function fgsl_permute_vector_inverse
+  function fgsl_permute_matrix(p,a)
+    type(fgsl_permutation), intent(in) :: p
+    type(fgsl_matrix), intent(inout) :: a
+    integer(fgsl_int) :: fgsl_permute_matrix
+    fgsl_permute_matrix = &
+         gsl_permute_matrix(p%gsl_permutation,a%gsl_matrix)
+  end function fgsl_permute_matrix
   function fgsl_permutation_mul(p, pa, pb)
     type(fgsl_permutation), intent(inout) :: p
     type(fgsl_permutation),  intent(in) :: pa, pb
diff --git a/api/rng.finc b/api/rng.finc
index 47f2a9c..7391f2f 100644
--- a/api/rng.finc
+++ b/api/rng.finc
@@ -252,6 +252,47 @@
     fgsl_ran_bivariate_gaussian_pdf = &
          gsl_ran_bivariate_gaussian_pdf(x, y, sigma_x, sigma_y, rho)
   end function fgsl_ran_bivariate_gaussian_pdf
+  function fgsl_ran_multivariate_gaussian(r, mu, l, result)
+    type(fgsl_rng), intent(in) :: r
+    type(fgsl_vector), intent(in) :: mu
+    type(fgsl_vector), intent(inout) :: result
+    type(fgsl_matrix), intent(in) :: l
+    integer(fgsl_int) :: fgsl_ran_multivariate_gaussian
+    fgsl_ran_multivariate_gaussian = gsl_ran_multivariate_gaussian(r%gsl_rng, &
+         mu%gsl_vector, l%gsl_matrix, result%gsl_vector)
+  end function fgsl_ran_multivariate_gaussian
+  function fgsl_ran_multivariate_gaussian_pdf(x, mu, l, result, work)
+    type(fgsl_vector), intent(in) :: x, mu
+    real(fgsl_double), intent(inout) :: result
+    type(fgsl_vector), intent(inout) :: work
+    type(fgsl_matrix), intent(in) :: l
+    integer(fgsl_int) :: fgsl_ran_multivariate_gaussian_pdf
+    fgsl_ran_multivariate_gaussian_pdf = gsl_ran_multivariate_gaussian_pdf(x%gsl_vector, &
+         mu%gsl_vector, l%gsl_matrix, result, work%gsl_vector)
+  end function fgsl_ran_multivariate_gaussian_pdf
+  function fgsl_ran_multivariate_gaussian_log_pdf(x, mu, l, result, work)
+    type(fgsl_vector), intent(in) :: x, mu
+    real(fgsl_double), intent(inout) :: result
+    type(fgsl_vector), intent(inout) :: work
+    type(fgsl_matrix), intent(in) :: l
+    integer(fgsl_int) :: fgsl_ran_multivariate_gaussian_log_pdf
+    fgsl_ran_multivariate_gaussian_log_pdf = gsl_ran_multivariate_gaussian_log_pdf( &
+         x%gsl_vector, mu%gsl_vector, l%gsl_matrix, result, work%gsl_vector)
+  end function fgsl_ran_multivariate_gaussian_log_pdf
+  function fgsl_ran_multivariate_gaussian_mean(x, mu_hat)
+    type(fgsl_matrix), intent(in) :: x
+    type(fgsl_vector), intent(inout) :: mu_hat
+    integer(fgsl_int) :: fgsl_ran_multivariate_gaussian_mean
+    fgsl_ran_multivariate_gaussian_mean = gsl_ran_multivariate_gaussian_mean( &
+         x%gsl_matrix, mu_hat%gsl_vector )
+  end function fgsl_ran_multivariate_gaussian_mean
+  function fgsl_ran_multivariate_gaussian_vcov(x, sigma_hat)
+    type(fgsl_matrix), intent(in) :: x
+    type(fgsl_matrix), intent(inout) :: sigma_hat
+    integer(fgsl_int) :: fgsl_ran_multivariate_gaussian_vcov
+    fgsl_ran_multivariate_gaussian_vcov = gsl_ran_multivariate_gaussian_vcov( &
+         x%gsl_matrix, sigma_hat%gsl_matrix )
+  end function fgsl_ran_multivariate_gaussian_vcov
   function fgsl_ran_exponential(r, mu)
     type(fgsl_rng), intent(in) :: r
     real(fgsl_double), intent(in) :: mu
diff --git a/api/rstat.finc b/api/rstat.finc
index c9fa253..d6a0460 100644
--- a/api/rstat.finc
+++ b/api/rstat.finc
@@ -14,6 +14,10 @@ subroutine fgsl_rstat_quantile_free(w)
   type(fgsl_rstat_quantile_workspace), intent(inout) :: w
   call gsl_rstat_quantile_free(w%gsl_rstat_quantile_workspace)
 end subroutine fgsl_rstat_quantile_free
+integer(fgsl_int) function fgsl_rstat_quantile_reset(w)
+  type(fgsl_rstat_quantile_workspace), intent(inout) :: w
+  fgsl_rstat_quantile_reset = gsl_rstat_quantile_reset(w%gsl_rstat_quantile_workspace)
+end function fgsl_rstat_quantile_reset
 function fgsl_rstat_quantile_add(x, w)
   real(fgsl_double), intent(in) :: x
   type(fgsl_rstat_quantile_workspace), intent(inout) :: w
@@ -59,6 +63,11 @@ function fgsl_rstat_mean(w)
   real(fgsl_double) :: fgsl_rstat_mean
   fgsl_rstat_mean = gsl_rstat_mean(w%gsl_rstat_workspace)
 end function fgsl_rstat_mean
+function fgsl_rstat_rms(w)
+  type(fgsl_rstat_workspace), intent(inout) :: w
+  real(fgsl_double) :: fgsl_rstat_rms
+  fgsl_rstat_rms = gsl_rstat_rms(w%gsl_rstat_workspace)
+end function fgsl_rstat_rms
 function fgsl_rstat_variance(w)
   type(fgsl_rstat_workspace), intent(inout) :: w
   real(fgsl_double) :: fgsl_rstat_variance
diff --git a/api/spmatrix.finc b/api/spmatrix.finc
index 3348bd6..6efc1c5 100644
--- a/api/spmatrix.finc
+++ b/api/spmatrix.finc
@@ -16,6 +16,11 @@ function fgsl_spmatrix_alloc_nzmax(n1, n2, nzmax, flags)
   fgsl_spmatrix_alloc_nzmax%gsl_spmatrix = &
   gsl_spmatrix_alloc_nzmax(n1, n2, nzmax, flags)
 end function fgsl_spmatrix_alloc_nzmax
+subroutine fgsl_spmatrix_size(m, n1, n2)
+  type(fgsl_spmatrix), intent(in) :: m
+  integer(fgsl_size_t), intent(inout) :: n1, n2
+  call gsl_spmatrix_size(m%gsl_spmatrix, n1, n2)
+end subroutine
 subroutine fgsl_spmatrix_free(m)
   type(fgsl_spmatrix), intent(in) :: m
   call gsl_spmatrix_free(m%gsl_spmatrix)
@@ -116,3 +121,19 @@ function fgsl_spmatrix_transpose_memcpy(dest, src)
   integer(fgsl_int) :: fgsl_spmatrix_transpose_memcpy
   fgsl_spmatrix_transpose_memcpy = gsl_spmatrix_transpose_memcpy(dest%gsl_spmatrix, src%gsl_spmatrix)
 end function fgsl_spmatrix_transpose_memcpy
+integer(fgsl_int) function fgsl_spblas_dgemv(transa, alpha, a, x, beta, y)
+  integer(fgsl_int), intent(in) :: transa
+  real(fgsl_double), intent(in) :: alpha, beta
+  type(fgsl_spmatrix), intent(in) :: a
+  type(fgsl_vector), intent(in) :: x
+  type(fgsl_vector), intent(inout) :: y
+  fgsl_spblas_dgemv = gsl_spblas_dgemv(transa, alpha, a%gsl_spmatrix, x%gsl_vector, &
+       beta, y%gsl_vector)
+end function fgsl_spblas_dgemv
+integer(fgsl_int) function fgsl_spblas_dgemm(alpha, a, b, c)
+  real(fgsl_double), intent(in) :: alpha
+  type(fgsl_spmatrix), intent(in) :: a, b
+  type(fgsl_spmatrix), intent(inout) :: c
+  fgsl_spblas_dgemm = gsl_spblas_dgemm(alpha, a%gsl_spmatrix, b%gsl_spmatrix, &
+       c%gsl_spmatrix)
+end function fgsl_spblas_dgemm
diff --git a/configure.ac b/configure.ac
index 53c161a..5e54fa9 100644
--- a/configure.ac
+++ b/configure.ac
@@ -22,7 +22,7 @@ AC_USE_SYSTEM_EXTENSIONS
 m4_ifdef([AM_PROG_AR],[AM_PROG_AR])
 
 
-LT_PREREQ([2.3.0])
+LT_PREREQ([2.1.0])
 LT_INIT([win32-dll])
 
 AC_PROG_CC
@@ -32,7 +32,7 @@ fi
 #initialize pkg-config
 PKG_PROG_PKG_CONFIG
 
-PKG_CHECK_MODULES([gsl],gsl >= 2.3)
+PKG_CHECK_MODULES([gsl], [gsl >= 2.3])
 
 local_gsl_version=`$PKG_CONFIG --modversion gsl`
 AC_SUBST(GSL_VERSION,$local_gsl_version)
diff --git a/doc/examples/Makefile.am b/doc/examples/Makefile.am
index ac0a6d1..03e61a6 100644
--- a/doc/examples/Makefile.am
+++ b/doc/examples/Makefile.am
@@ -11,7 +11,7 @@ check_PROGRAMS = cdf cheb combination diff dwt \
 		 fitting fitting3 fitting2 histogram \
 		 histogram2d ieee ieeeround integration \
 		 interp interpp intro matrix min monte \
-                 multifit nlfit nlfit2 \
+                 multifit nlfit nlfit2 nlfit4 \
 		 ntuplew ntupler \
 		 linalglu eigen eigen_nonsymm fft \
 		 ode-initval permutation permshuffle \
@@ -69,6 +69,8 @@ nlfit_SOURCES = nlfit.f90
 
 nlfit2_SOURCES = nlfit2.f90
 
+nlfit4_SOURCES = nlfit4.f90
+
 ntupler_SOURCES = ntupler.f90
 
 ntuplew_SOURCES = ntuplew.f90
@@ -128,7 +130,7 @@ dist_fgsl_example_DATA = $(data_files) cdf.f90 cheb.f90 combination.F90 diff.f90
 			fitting.f90 fitting3.f90 fitting2.f90 histogram.F90 \
 			histogram2d.f90 ieee.f90 ieeeround.f90 integration.f90 \
 			interp.f90 interpp.f90 intro.f90 matrix.f90 min.f90 monte.f90 \
-			nlfit.f90 ntuplew.f90 ntupler.f90 \
+			nlfit.f90 nlfit2.f90 nlfit4.f90 ntuplew.f90 ntupler.f90 \
 			linalglu.f90 eigen.f90 eigen_nonsymm.f90 fft.f90 \
 			ode-initval.F90 permutation.f90 permshuffle.f90 \
 			polyroots.f90 qrng.f90 randpoisson.f90 randwalk.f90 \
diff --git a/doc/examples/fitting2.f90 b/doc/examples/fitting2.f90
index 0797354..ac5f078 100644
--- a/doc/examples/fitting2.f90
+++ b/doc/examples/fitting2.f90
@@ -2,7 +2,7 @@ program fitting2
   use fgsl
   implicit none
   integer(fgsl_int) :: status
-  integer(fgsl_size_t) :: i, n
+  integer(fgsl_size_t) :: i, n, rk
   real(fgsl_double) :: xi, yi, ei, chisq
   type(fgsl_matrix) :: x, cov
   type(fgsl_vector) :: y, w, c
@@ -42,14 +42,16 @@ program fitting2
   close(20)
   work = fgsl_multifit_linear_alloc(n, 3_fgsl_size_t)
   status = fgsl_multifit_wlinear(x, w, y, c, cov, chisq, work)
+  rk = fgsl_multifit_linear_rank(1.0d-4, work)
   call fgsl_multifit_linear_free(work)
-  write(6, '(''# best fit:  Y = '',F12.5,'' + '',F12.5, &
+  write(*, '(''# best fit:  Y = '',F12.5,'' + '',F12.5, &
        & ''*X + '',F12.5,''*X^2'')') c_v
-  write(6, '(''# covariance matrix: '')') 
-  write(6, '(3(F12.5,1X))')  cov_m(1:3,1)
-  write(6, '(3(F12.5,1X))')  cov_m(1:3,2)
-  write(6, '(3(F12.5,1X))')  cov_m(1:3,3)
-  write(6, '(''# chisq = '',F12.5)') chisq
+  write(*, '(''# covariance matrix: '')') 
+  write(*, '(3(F12.5,1X))')  cov_m(1:3,1)
+  write(*, '(3(F12.5,1X))')  cov_m(1:3,2)
+  write(*, '(3(F12.5,1X))')  cov_m(1:3,3)
+  write(*, '(''# chisq = '',F12.5)') chisq
+  write(*, '(''# rank = '',I0)') rk 
 
   call fgsl_vector_free(y)
   call fgsl_vector_free(w)
diff --git a/doc/examples/histogram.F90 b/doc/examples/histogram.F90
index 44a97f9..0c2d10b 100644
--- a/doc/examples/histogram.F90
+++ b/doc/examples/histogram.F90
@@ -4,12 +4,12 @@ program histogram
   real(fgsl_double) :: a, b, x
   integer(fgsl_size_t) :: n
   integer(fgsl_int) :: status
-  integer(8) :: i
+  integer :: i
   type(fgsl_histogram) :: h
   type(fgsl_file) :: stdout
   namelist / histpar / a, b, n
 !
-  write(6, *) 'Reading histogram parameters from file histogram.dat ... '
+  write(*, *) 'Reading histogram parameters from file histogram.dat ... '
   open(20, file=HISTOGRAM_DAT, form='formatted', status='old')
   read(20, nml=histpar) 
   close(20)
@@ -22,7 +22,7 @@ program histogram
      status = fgsl_histogram_increment (h, x)
   end do
   close(20)
-  write(6, *) 'Printing histogram:'
+  write(*, *) 'Printing histogram:'
   status = fgsl_histogram_fprintf(stdout, h,'%12.5f','%7.0f')
   call fgsl_histogram_free(h)
 end program histogram
diff --git a/doc/examples/nlfit4.f90 b/doc/examples/nlfit4.f90
new file mode 100644
index 0000000..5507436
--- /dev/null
+++ b/doc/examples/nlfit4.f90
@@ -0,0 +1,265 @@
+module timer
+! Timer implementation based on system_clock intrinsic
+! Note that this facility is not thread-safe
+  implicit none
+  private
+  public :: dwalltime
+  integer, parameter :: ik = selected_int_kind(6)
+  logical, save :: first = .true.
+  integer(ik), save ::  count_rate, count_max
+  double precision, save :: conversion = 0.0d0
+contains
+  double precision function dwalltime()
+    integer(ik) :: count
+    if (first) then
+       first = .false.
+       call system_clock(count, count_rate, count_max)
+       conversion = 1.0d0 / dble(count_rate)
+    else
+       call system_clock(count)
+    end if
+    dwalltime = count * conversion
+  end function dwalltime
+end module timer
+module mod_nlfit4 
+  use fgsl
+  use, intrinsic :: iso_c_binding
+  implicit none
+  integer(fgsl_size_t), parameter :: max_iter = 200
+  real(fgsl_double), parameter :: xtol = 1.0e-8_fgsl_double, &
+       gtol = 1.0e-8_fgsl_double, ftol = 1.0e-8_fgsl_double
+! 
+! parameters for functions
+  type :: model_params
+     real(fgsl_double) :: alpha
+     type(fgsl_spmatrix) :: J
+  end type model_params
+contains
+!
+! penalty function
+  integer(c_int) function penalty_f(x, params, f) BIND(C)
+    type(c_ptr), value :: x, params, f
+
+    type(model_params), pointer :: par
+    real(fgsl_double) :: sqrt_alpha, sx
+    type(fgsl_vector) :: xx, ff
+    real(fgsl_double), pointer :: xxp(:), ffp(:)
+    integer(fgsl_size_t) :: i
+    integer(fgsl_int) :: status
+
+    call c_f_pointer(params, par)
+    sqrt_alpha = sqrt(par%alpha)
+    call fgsl_obj_c_ptr(xx, x)
+    status = fgsl_vector_align(xxp, xx)
+    call fgsl_obj_c_ptr(ff, f)
+    status = fgsl_vector_align(ffp, ff)
+    sx = 0.0_fgsl_double
+    do i = 1, size(xxp)
+       ffp(i) = sqrt_alpha * (xxp(i) - 1.0_fgsl_double)
+       sx = sx + xxp(i)**2
+    end do
+    ffp(size(xxp)+1) = sx - 0.25_fgsl_double
+    penalty_f = fgsl_success
+  end function penalty_f
+  integer(c_int) function penalty_df(transj, x, u, params, &
+       v, jtj) BIND(C)
+    integer(c_int), value :: transj
+    type(c_ptr), value :: x, u, params, v, jtj
+    
+    type(model_params), pointer :: par
+    type(fgsl_vector) :: xx, uu, vv
+    real(fgsl_double), pointer :: xxp(:), uup(:)
+    type(fgsl_matrix) :: jj
+    real(fgsl_double), pointer :: jtjp(:,:)
+    integer(fgsl_size_t) :: i, j
+    integer(fgsl_int) :: status
+
+    call c_f_pointer(params, par)
+    call fgsl_obj_c_ptr(xx, x)
+    status = fgsl_vector_align(xxp, xx)
+    call fgsl_obj_c_ptr(uu, u)
+    status = fgsl_vector_align(uup, uu)
+!
+!   store 2*x in last row of J
+    do i = 1, size(xxp)
+       status = fgsl_spmatrix_set(par%J, size(xxp, KIND=fgsl_size_t), i-1, &
+            2.0_fgsl_double * xxp(i))
+    end do
+!
+!   compute v = op(J) * u
+    if (c_associated(v)) then
+       call fgsl_obj_c_ptr(vv, v)
+       call fgsl_spmatrix_size(par%J, i, j)
+       
+       status = fgsl_spblas_dgemv(transj, 1.0_fgsl_double, &
+            par%J, uu, 0.0_fgsl_double, vv)
+    end if
+
+!
+!   compute JTJ = alpha*I_p + 4 * x * x^T
+!   note that "Lower" requires i <= j for GSL
+    if (c_associated(jtj)) then
+       call fgsl_obj_c_ptr(jj, jtj)
+       status = fgsl_matrix_align(jtjp, jj)
+       do j = 1, size(jtjp, 2)
+          jtjp(j, j) = 4.0_fgsl_double * xxp(j)**2 + par%alpha
+          do i = 1, j-1
+             jtjp(i, j) = 4.0_fgsl_double * xxp(i) * xxp(j)
+          end do
+       end do
+    end if
+    penalty_df = fgsl_success
+  end function penalty_df
+  integer(fgsl_int) function penalty_fvv(x, v, params, fvv) BIND(C)
+    type(c_ptr), value :: x, v, params, fvv
+    type(fgsl_vector) :: xx, f_v, f_fvv
+    real(fgsl_double), pointer :: p_xx(:), p_v(:), p_fvv(:)
+    integer(fgsl_int) :: status
+
+    call fgsl_obj_c_ptr(f_v, v)
+    status = fgsl_vector_align(p_v, f_v)
+    call fgsl_obj_c_ptr(f_fvv, fvv)
+    status = fgsl_vector_align(p_fvv, f_fvv)
+    call fgsl_obj_c_ptr(xx, x)
+    status = fgsl_vector_align(p_xx, xx)
+   
+    p_fvv = 0.0_fgsl_double
+    p_fvv( size(p_xx) + 1 ) = 2.0_fgsl_double * sum(p_v * p_v)
+    
+    penalty_fvv = fgsl_success
+  end function penalty_fvv
+
+  subroutine solve_system(x0, fdf, params)
+    use timer
+    type(fgsl_vector) :: x0
+    type(fgsl_multilarge_nlinear_fdf) :: fdf
+    type(fgsl_multilarge_nlinear_parameters) :: params
+
+    type(fgsl_multilarge_nlinear_type) :: t
+    type(fgsl_multilarge_nlinear_workspace) :: work
+    integer(fgsl_size_t) :: n, p
+    type(fgsl_vector) :: f, x
+    real(fgsl_double), pointer :: p_x(:), p_f(:)
+    integer(fgsl_int) :: info, status
+    integer(fgsl_size_t) :: nevalf, nevaldfu, nevaldf2, nevalfvv
+    real(fgsl_double) :: chisq0, chisq, rcond, xsq, ti
+
+
+    t = fgsl_multilarge_nlinear_type('trust')
+    call fgsl_multilarge_nlinear_fdf_get(fdf, N=n, P=p)
+    work = fgsl_multilarge_nlinear_alloc(t, params, n, p)
+    f = fgsl_multilarge_nlinear_residual(work)
+    x = fgsl_multilarge_nlinear_position(work)
+    status = fgsl_vector_align(p_x, x)
+    status = fgsl_vector_align(p_f, f)
+
+   ! write(*, *) 'SIzes n, p: ', n, p
+   !write(*, *) 'SIzes f, x: ', size(p_f), size(p_x)
+
+    ti = dwalltime()
+!
+! initialize solver
+    status = fgsl_multilarge_nlinear_init(x0, fdf, work)
+!
+! store initial cost
+    chisq0 = dot_product(p_f,p_f)
+!
+! iterate to convergence
+    status = fgsl_multilarge_nlinear_driver(max_iter, xtol, gtol, ftol, &
+         CALLBACK_PARAMS=c_null_ptr, INFO=info, W=work)
+
+    ti = dwalltime() - ti
+!
+! store final cost
+    chisq = dot_product(p_f,p_f)
+!
+! store final x^2
+    xsq = dot_product(p_x,p_x)
+!
+! store cond(J(x))
+    status = fgsl_multilarge_nlinear_rcond(rcond, work)
+
+!
+! print summary
+    call fgsl_multilarge_nlinear_fdf_get(fdf, NEVALF=nevalf, &
+         NEVALDFU=nevaldfu, NEVALDF2=nevaldf2, NEVALFVV=nevalfvv)
+
+    write(*, fmt='(A25,I4,I4,I4,I4,I4,1P,4(E11.4),1X,0P,F6.2)') &
+         fgsl_multilarge_nlinear_trs_name(work), &
+         fgsl_multilarge_nlinear_niter(work), &
+         nevalf, nevaldfu, nevaldf2, nevalfvv, &
+         chisq0, chisq, 1.0_fgsl_double/rcond, xsq, ti
+
+    call fgsl_multilarge_nlinear_free(work)
+  end subroutine solve_system
+end module mod_nlfit4
+program nlfit4
+  use mod_nlfit4
+  implicit none
+! Problem size reduced to decrease run time
+  integer(fgsl_size_t), parameter :: p = 500, n = p+1
+  integer(fgsl_int) :: status
+  integer(fgsl_size_t) :: i
+
+  type(fgsl_vector) :: f, x
+  real(fgsl_double), target :: tf(n), tx(p)
+  type(fgsl_multilarge_nlinear_fdf) :: fdf
+  type(fgsl_multilarge_nlinear_parameters) :: fdf_params
+  type(model_params), target :: params
+
+  fdf_params = fgsl_multilarge_nlinear_default_parameters()
+
+  f = fgsl_vector_init(type = 1.0_fgsl_double)
+  x = fgsl_vector_init(type = 1.0_fgsl_double)
+  status = fgsl_vector_align(tf, n, f, n, 0_fgsl_size_t, 1_fgsl_size_t)
+  status = fgsl_vector_align(tx, p, x, p, 0_fgsl_size_t, 1_fgsl_size_t)
+!
+! model params: sparse Jacobian matrix with 2*p non-zero elements
+! in triplet format
+  params%alpha = 1.0e-5_fgsl_double
+  params%J = fgsl_spmatrix_alloc_nzmax(n, p, 2*p, FGSL_SPMATRIX_TRIPLET)
+!
+! define function to be minimized
+  fdf = fgsl_multilarge_nlinear_fdf_init(n, p, c_loc(params), &
+       penalty_f, penalty_df, penalty_fvv)
+
+  do i=1, p
+!    starting point
+     tx(i) = real(i, fgsl_double)
+!    store sqrt(alpha)*I_p in upper p-by-p block of J
+     status = fgsl_spmatrix_set( params%J, i-1, i-1, sqrt(params%alpha) )
+  end do
+
+  write(*, '(''Method                NITER   NFEV  NJUEV NJTJEV NAEV&
+       &Init Cost Final Cost cond(J) Final x^2 Time (s)'')')
+
+  call fgsl_multilarge_nlinear_parameters_set(fdf_params, &
+       SCALE=fgsl_multilarge_nlinear_scale_levenberg, &
+       TRS=fgsl_multilarge_nlinear_trs_lm)
+  call solve_system(x, fdf, fdf_params)
+
+  call fgsl_multilarge_nlinear_parameters_set(fdf_params, &
+       TRS=fgsl_multilarge_nlinear_trs_lmaccel)
+  call solve_system(x, fdf, fdf_params)
+
+  call fgsl_multilarge_nlinear_parameters_set(fdf_params, &
+       TRS=fgsl_multilarge_nlinear_trs_dogleg)
+  call solve_system(x, fdf, fdf_params)
+
+  call fgsl_multilarge_nlinear_parameters_set(fdf_params, &
+       TRS=fgsl_multilarge_nlinear_trs_ddogleg)
+  call solve_system(x, fdf, fdf_params)
+
+  call fgsl_multilarge_nlinear_parameters_set(fdf_params, &
+       TRS=fgsl_multilarge_nlinear_trs_subspace2d)
+  call solve_system(x, fdf, fdf_params)
+
+  call fgsl_multilarge_nlinear_parameters_set(fdf_params, &
+       TRS=fgsl_multilarge_nlinear_trs_cgst)
+  call solve_system(x, fdf, fdf_params)
+
+  call fgsl_vector_free(x)
+  call fgsl_vector_free(f)
+  call fgsl_multilarge_nlinear_fdf_free(fdf)
+  call fgsl_spmatrix_free(params%J)
+end program nlfit4
diff --git a/doc/examples/rstat.f90 b/doc/examples/rstat.f90
index ea5f824..8fc50a9 100644
--- a/doc/examples/rstat.f90
+++ b/doc/examples/rstat.f90
@@ -3,7 +3,7 @@ subroutine rstat1()
   use, intrinsic :: iso_fortran_env
   implicit none
   real(fgsl_double) :: data(5) = [17.2, 18.1, 16.5, 18.3, 12.6]
-  real(fgsl_double) :: mean, variance, largest, smallest
+  real(fgsl_double) :: mean, variance, largest, smallest, rms
   type(fgsl_rstat_workspace) :: rstat_p
   integer(fgsl_size_t) :: i
   integer(fgsl_int) :: status
@@ -18,12 +18,13 @@ subroutine rstat1()
   variance = fgsl_rstat_variance(rstat_p)
   largest  = fgsl_rstat_max(rstat_p)
   smallest = fgsl_rstat_min(rstat_p)
+  rms      = fgsl_rstat_rms(rstat_p)
 
-  write(output_unit, '(A,5G10.4)') 'The dataset is ', (data(i), i=1,5)
-  write(output_unit, '(A,G10.4)') 'The sample mean is ', mean
-  write(output_unit, '(A,G10.4)') 'The estimated variance is ', variance
-  write(output_unit, '(A,G10.4)') 'The largest value is ', largest
-  write(output_unit, '(A,G10.4)') 'The smallest value is ', smallest
+  write(output_unit, '(A,5G11.4)') 'The dataset is ', (data(i), i=1,5)
+  write(output_unit, '(A,G11.4)') 'The sample mean is ', mean
+  write(output_unit, '(A,G11.4)') 'The estimated variance is ', variance
+  write(output_unit, '(A,G11.4)') 'The root mean square is ', rms
+  write(output_unit, '(A,G11.4)') 'The largest value is ', largest
 
   status = fgsl_rstat_reset(rstat_p)
   call fgsl_rstat_free(rstat_p)
@@ -69,19 +70,17 @@ subroutine rstat2
   val_p75 = fgsl_rstat_quantile_get(work_75)
 
   write(output_unit, '(A,5G12.4,A)') 'The dataset is ', (data(i), i=1,5), ' ...'
-  write(output_unit, '(A,G10.4,A,G10.4,A,G10.4)') '0.25 quartile: exact = ', &
+  write(output_unit, '(A,G11.4,A,G11.4,A,G11.4)') '0.25 quartile: exact = ', &
     exact_p25, ', estimated = ', val_p25, ', error = ', (val_p25 - exact_p25) / exact_p25
-  write(output_unit, '(A,G10.4,A,G10.4,A,G10.4)') '0.50 quartile: exact = ', &
+  write(output_unit, '(A,G11.4,A,G11.4,A,G11.4)') '0.50 quartile: exact = ', &
     exact_p50, ', estimated = ', val_p50, ', error = ', (val_p50 - exact_p50) / exact_p50
-  write(output_unit, '(A,G10.4,A,G10.4,A,G10.4)') '0.75 quartile: exact = ', &
+  write(output_unit, '(A,G11.4,A,G11.4,A,G11.4)') '0.75 quartile: exact = ', &
     exact_p75, ', estimated = ', val_p75, ', error = ', (val_p75 - exact_p75) / exact_p75
 
-  !  printf ("0.25 quartile: exact = %.5f, estimated = %.5f, error = %.6e\n",
-  !          exact_p25, val_p25, (val_p25 - exact_p25) / exact_p25);
-  !  printf ("0.50 quartile: exact = %.5f, estimated = %.5f, error = %.6e\n",
-  !          exact_p50, val_p50, (val_p50 - exact_p50) / exact_p50);
-  !  printf ("0.75 quartile: exact = %.5f, estimated = %.5f, error = %.6e\n",
-  !          exact_p75, val_p75, (val_p75 - exact_p75) / exact_p75);
+
+  status = fgsl_rstat_quantile_reset(work_25)
+  val_p25 = fgsl_rstat_quantile_get(work_25)
+  write(output_unit, '(A,G11.4)') 'estimate after reset: ', val_p25
 
   call fgsl_rstat_quantile_free(work_25)
   call fgsl_rstat_quantile_free(work_50)
diff --git a/fgsl.F90 b/fgsl.F90
index 23bea2c..82456a7 100644
--- a/fgsl.F90
+++ b/fgsl.F90
@@ -304,7 +304,7 @@ module fgsl
        fgsl_permutation_valid, fgsl_permutation_reverse, fgsl_permutation_inverse, &
        fgsl_permutation_next, fgsl_permutation_prev, &
        fgsl_permute, fgsl_permute_inverse, fgsl_permute_vector, &
-       fgsl_permute_vector_inverse, fgsl_permutation_mul, &
+       fgsl_permute_vector_inverse, fgsl_permute_matrix, fgsl_permutation_mul, &
        fgsl_permutation_linear_to_canonical, fgsl_permutation_canonical_to_linear, &
        fgsl_permutation_inversions, fgsl_permutation_linear_cycles, &
        fgsl_permutation_canonical_cycles, fgsl_permutation_fwrite, &
@@ -347,13 +347,28 @@ module fgsl
        fgsl_linalg_qrpt_decomp2, fgsl_linalg_qrpt_solve, fgsl_linalg_qrpt_svx, &
        fgsl_linalg_qrpt_qrsolve, fgsl_linalg_qrpt_update, &
        fgsl_linalg_qrpt_rsolve, fgsl_linalg_qrpt_rsvx, &
+       fgsl_linalg_qrpt_lssolve, fgsl_linalg_qrpt_lssolve2, &
+       fgsl_linalg_qrpt_rank, fgsl_linalg_qrpt_rcond, &
+       fgsl_linalg_cod_decomp, fgsl_linalg_cod_decomp_e, fgsl_linalg_cod_lssolve, &
+       fgsl_linalg_cod_unpack, fgsl_linalg_cod_matz, &
        fgsl_linalg_sv_decomp, fgsl_linalg_sv_decomp_mod, &
        fgsl_linalg_sv_decomp_jacobi, fgsl_linalg_sv_solve, &
-       fgsl_linalg_cholesky_decomp, fgsl_linalg_complex_cholesky_decomp, &
+       fgsl_linalg_cholesky_decomp1, fgsl_linalg_complex_cholesky_decomp, &
+       fgsl_linalg_cholesky_decomp, &
        fgsl_linalg_cholesky_solve, fgsl_linalg_complex_cholesky_solve, &
        fgsl_linalg_cholesky_svx, fgsl_linalg_complex_cholesky_svx, &
-       fgsl_linalg_cholesky_invert, &
+       fgsl_linalg_cholesky_invert, fgsl_linalg_cholesky_rcond, &
        fgsl_linalg_complex_cholesky_invert, &
+       fgsl_linalg_cholesky_decomp2, fgsl_linalg_cholesky_solve2, &
+       fgsl_linalg_cholesky_svx2, &
+       fgsl_linalg_cholesky_scale, fgsl_linalg_cholesky_scale_apply, &
+       fgsl_linalg_pcholesky_decomp, fgsl_linalg_pcholesky_solve, &
+       fgsl_linalg_pcholesky_svx, fgsl_linalg_pcholesky_decomp2, &
+       fgsl_linalg_pcholesky_solve2, fgsl_linalg_pcholesky_svx2, &
+       fgsl_linalg_pcholesky_invert, fgsl_linalg_pcholesky_rcond, &
+       fgsl_linalg_mcholesky_decomp, fgsl_linalg_mcholesky_solve, &
+       fgsl_linalg_mcholesky_svx, fgsl_linalg_mcholesky_rcond, &
+       fgsl_linalg_mcholesky_invert, &
        fgsl_linalg_symmtd_decomp, fgsl_linalg_symmtd_unpack, &
        fgsl_linalg_symmtd_unpack_t, fgsl_linalg_hermtd_decomp, &
        fgsl_linalg_hermtd_unpack, fgsl_linalg_hermtd_unpack_t, &
@@ -370,7 +385,10 @@ module fgsl
        fgsl_linalg_hh_solve, fgsl_linalg_hh_svx, fgsl_linalg_solve_tridiag, &
        fgsl_linalg_solve_symm_tridiag, fgsl_linalg_solve_cyc_tridiag, &
        fgsl_linalg_solve_symm_cyc_tridiag, fgsl_linalg_balance_matrix, &
-       fgsl_linalg_qr_matq, fgsl_linalg_givens, fgsl_linalg_givens_gv
+       fgsl_linalg_qr_matq, fgsl_linalg_givens, fgsl_linalg_givens_gv, &
+       fgsl_linalg_tri_upper_invert, fgsl_linalg_tri_lower_invert, &
+       fgsl_linalg_tri_upper_unit_invert, fgsl_linalg_tri_lower_unit_invert, &
+       fgsl_linalg_tri_upper_rcond, fgsl_linalg_tri_lower_rcond
   public :: fgsl_linalg_sv_leverage
 ! eigensystems
   public :: fgsl_eigen_symm_alloc, fgsl_eigen_symm_free, fgsl_eigen_symm, &
@@ -438,7 +456,12 @@ module fgsl
        fgsl_cdf_gaussian_pinv, fgsl_cdf_gaussian_qinv, fgsl_cdf_ugaussian_p, &
        fgsl_cdf_ugaussian_q, fgsl_cdf_ugaussian_pinv, fgsl_cdf_ugaussian_qinv, &
        fgsl_ran_gaussian_tail, fgsl_ran_gaussian_tail_pdf, fgsl_ran_ugaussian_tail, &
-       fgsl_ran_ugaussian_tail_pdf, fgsl_ran_bivariate_gaussian, fgsl_ran_exponential, &
+       fgsl_ran_ugaussian_tail_pdf, fgsl_ran_bivariate_gaussian, &
+       fgsl_ran_bivariate_gaussian_pdf,  &
+       fgsl_ran_multivariate_gaussian, fgsl_ran_multivariate_gaussian_pdf, &
+       fgsl_ran_multivariate_gaussian_log_pdf, &
+       fgsl_ran_multivariate_gaussian_mean, fgsl_ran_multivariate_gaussian_vcov, &
+       fgsl_ran_exponential, &
        fgsl_ran_exponential_pdf, fgsl_cdf_exponential_p, fgsl_cdf_exponential_q, &
        fgsl_cdf_exponential_pinv, fgsl_cdf_exponential_qinv, fgsl_ran_laplace, &
        fgsl_ran_laplace_pdf, fgsl_cdf_laplace_p, fgsl_cdf_laplace_q, &
@@ -621,17 +644,21 @@ module fgsl
        fgsl_fit_wmul, fgsl_fit_mul_est
 
   public :: fgsl_multifit_linear_alloc, fgsl_multifit_linear_free, &
-       fgsl_multifit_linear, fgsl_multifit_linear_svd, &
-       fgsl_multifit_wlinear, fgsl_multifit_wlinear_svd, &
+       fgsl_multifit_linear, fgsl_multifit_linear_tsvd, fgsl_multifit_linear_svd, &
+       fgsl_multifit_wlinear,fgsl_multifit_wlinear_tsvd, fgsl_multifit_wlinear_svd, &
        fgsl_multifit_wlinear_usvd, fgsl_multifit_linear_bsvd, &
        fgsl_multifit_linear_est, fgsl_multifit_linear_solve, &
-       fgsl_multifit_linear_residuals, fgsl_multifit_linear_applyw, &
+       fgsl_multifit_linear_residuals, fgsl_multifit_linear_rank, &
+       fgsl_multifit_linear_applyw, &
        fgsl_multifit_linear_stdform1, fgsl_multifit_linear_wstdform1, &
        fgsl_multifit_linear_l_decomp, fgsl_multifit_linear_stdform2, &
        fgsl_multifit_linear_wstdform2, fgsl_multifit_linear_genform1, &
        fgsl_multifit_linear_genform2, fgsl_multifit_linear_wgenform2, &
        fgsl_multifit_linear_lreg, fgsl_multifit_linear_lcurve, &
        fgsl_multifit_linear_lcorner, fgsl_multifit_linear_lcorner2, &
+       fgsl_multifit_linear_gcv_init, fgsl_multifit_linear_gcv_curve, &
+       fgsl_multifit_linear_gcv_min, fgsl_multifit_linear_gcv_calc, &
+       fgsl_multifit_linear_gcv, &
        fgsl_multifit_linear_lk, fgsl_multifit_linear_lsobolev, &
        fgsl_multifit_linear_rcond, fgsl_multifit_robust_maxiter, &
        fgsl_multifit_robust_weights, fgsl_multifit_robust_residuals
@@ -679,8 +706,9 @@ module fgsl
        fgsl_multifit_nlinear_driver, fgsl_multilarge_nlinear_driver, &
        fgsl_multifit_nlinear_covar, fgsl_multilarge_nlinear_covar, &
        fgsl_multifit_nlinear_fdf_init, fgsl_multifit_nlinear_fdf_free, &
-       fgsl_multifit_nlinear_fdf_get, fgsl_multifit_nlinear_parameters_set
-       
+       fgsl_multifit_nlinear_fdf_get, fgsl_multifit_nlinear_parameters_set, &
+       fgsl_multilarge_nlinear_fdf_init, fgsl_multilarge_nlinear_fdf_free, &
+       fgsl_multilarge_nlinear_fdf_get, fgsl_multilarge_nlinear_parameters_set
 !
 ! large linear least squares systems
   public :: fgsl_multilarge_linear_alloc, fgsl_multilarge_linear_free, &
@@ -717,13 +745,14 @@ module fgsl
   public :: fgsl_bspline_knots_greville
 
 ! sparse matrices
-  public :: fgsl_spmatrix_alloc, fgsl_spmatrix_alloc_nzmax, &
+  public :: fgsl_spmatrix_alloc, fgsl_spmatrix_alloc_nzmax, fgsl_spmatrix_size, &
        fgsl_spmatrix_free, fgsl_spmatrix_realloc, fgsl_spmatrix_set_zero, &
        fgsl_spmatrix_nnz, fgsl_spmatrix_compare_idx, fgsl_spmatrix_memcpy, &
        fgsl_spmatrix_get, fgsl_spmatrix_set, fgsl_spmatrix_compcol, &
        fgsl_spmatrix_cumsum, fgsl_spmatrix_scale, fgsl_spmatrix_minmax, &
        fgsl_spmatrix_add, fgsl_spmatrix_d2sp, fgsl_spmatrix_sp2d, &
-       fgsl_spmatrix_equal, fgsl_spmatrix_transpose_memcpy
+       fgsl_spmatrix_equal, fgsl_spmatrix_transpose_memcpy, &
+       fgsl_spblas_dgemv, fgsl_spblas_dgemm
 
 ! sparse matrix linear algebra
   public :: fgsl_splinalg_itersolve_alloc, fgsl_splinalg_itersolve_free, &
@@ -733,9 +762,10 @@ module fgsl
 ! running statistics
   public :: fgsl_rstat_quantile_alloc, fgsl_rstat_quantile_free, &
        fgsl_rstat_quantile_add, fgsl_rstat_quantile_get, &
+       fgsl_rstat_quantile_reset, &
        fgsl_rstat_alloc, fgsl_rstat_free, fgsl_rstat_n, &
        fgsl_rstat_add, fgsl_rstat_min, fgsl_rstat_max, &
-       fgsl_rstat_mean, fgsl_rstat_variance, fgsl_rstat_sd, &
+       fgsl_rstat_mean, fgsl_rstat_rms, fgsl_rstat_variance, fgsl_rstat_sd, &
        fgsl_rstat_sd_mean, fgsl_rstat_median, fgsl_rstat_skew, &
        fgsl_rstat_kurtosis, fgsl_rstat_reset
 
@@ -1824,7 +1854,7 @@ module fgsl
      type(c_ptr) :: trs, scale, solver
      integer(c_int) :: fdtype
      real(c_double) :: factor_up, factor_down, avmax, h_df, h_fvv
-     integer(c_size_t) :: maxiter    
+     integer(c_size_t) :: max_iter    
      real(c_double) :: tol
   end type
   type, public :: fgsl_multilarge_nlinear_parameters
@@ -1853,6 +1883,12 @@ integer(c_int) function fgsl_nlinear_fdf_dfunc(x, params, df) bind(c)
       import :: c_ptr, c_int
       type(c_ptr), value :: x, params, df
     end function 
+    integer(c_int) function fgsl_nlinear_fdf_dlfunc(t, x, u, params, v, jtj) bind(c)
+      import :: c_ptr, c_int
+!     assuming int is the correct integer for enum type
+      integer(c_int), value :: t
+      type(c_ptr), value :: x, u, params, v, jtj
+    end function 
     integer(c_int) function fgsl_nlinear_fdf_fvv(x, v, params, vv) bind(c)
       import :: c_ptr, c_int
       type(c_ptr), value :: x, v, params, vv
@@ -1869,6 +1905,16 @@ integer(c_int) function fgsl_nlinear_fdf_fvv(x, v, params, vv) bind(c)
           fgsl_multifit_nlinear_trs_dogleg = fgsl_multifit_nlinear_trs(3), &
           fgsl_multifit_nlinear_trs_ddogleg = fgsl_multifit_nlinear_trs(4), &
           fgsl_multifit_nlinear_trs_subspace2d = fgsl_multifit_nlinear_trs(5)
+  type, private :: fgsl_multilarge_nlinear_trs
+    integer(c_int) :: which = 0
+  end type
+  type(fgsl_multilarge_nlinear_trs), public, parameter :: &
+          fgsl_multilarge_nlinear_trs_lm = fgsl_multilarge_nlinear_trs(1), &
+          fgsl_multilarge_nlinear_trs_lmaccel = fgsl_multilarge_nlinear_trs(2), &
+          fgsl_multilarge_nlinear_trs_dogleg = fgsl_multilarge_nlinear_trs(3), &
+          fgsl_multilarge_nlinear_trs_ddogleg = fgsl_multilarge_nlinear_trs(4), &
+          fgsl_multilarge_nlinear_trs_subspace2d = fgsl_multilarge_nlinear_trs(5), &
+          fgsl_multilarge_nlinear_trs_cgst = fgsl_multilarge_nlinear_trs(6)
 ! 
 ! scaling matrix strategies
   type, private :: fgsl_multifit_nlinear_scale
@@ -1878,6 +1924,13 @@ integer(c_int) function fgsl_nlinear_fdf_fvv(x, v, params, vv) bind(c)
           fgsl_multifit_nlinear_scale_levenberg = fgsl_multifit_nlinear_scale(1), &
           fgsl_multifit_nlinear_scale_marquardt = fgsl_multifit_nlinear_scale(2), &
           fgsl_multifit_nlinear_scale_more = fgsl_multifit_nlinear_scale(3)
+  type, private :: fgsl_multilarge_nlinear_scale
+    integer(c_int) :: which = 0
+  end type
+  type(fgsl_multilarge_nlinear_scale), public, parameter :: &
+          fgsl_multilarge_nlinear_scale_levenberg = fgsl_multilarge_nlinear_scale(1), &
+          fgsl_multilarge_nlinear_scale_marquardt = fgsl_multilarge_nlinear_scale(2), &
+          fgsl_multilarge_nlinear_scale_more = fgsl_multilarge_nlinear_scale(3)
 !
 ! linear solvers
   type, private :: fgsl_multifit_nlinear_solver
@@ -1889,6 +1942,11 @@ integer(c_int) function fgsl_nlinear_fdf_fvv(x, v, params, vv) bind(c)
           fgsl_multifit_nlinear_solver_svd = fgsl_multifit_nlinear_solver(3)
   integer(fgsl_int), parameter, public :: FGSL_MULTIFIT_NLINEAR_FWDIFF = 0, & 
                                   FGSL_MULTIFIT_NLINEAR_CTRDIFF = 1
+  type, private :: fgsl_multilarge_nlinear_solver
+    integer(c_int) :: which = 0
+  end type
+  type(fgsl_multilarge_nlinear_solver), public, parameter :: &
+          fgsl_multilarge_nlinear_solver_cholesky = fgsl_multilarge_nlinear_solver(1)
 !
 ! nonlinear fitting legacy interface
   type, public :: fgsl_multifit_function
@@ -1934,8 +1992,8 @@ integer(c_int) function fgsl_nlinear_fdf_fvv(x, v, params, vv) bind(c)
 !
 ! Types: sparse matrices
 !
-  integer(fgsl_int), public, parameter :: fgsl_spmatrix_triplet = 0
-  integer(fgsl_int), public, parameter :: fgsl_spmatrix_ccs = 1
+  integer(fgsl_size_t), public, parameter :: fgsl_spmatrix_triplet = 0
+  integer(fgsl_size_t), public, parameter :: fgsl_spmatrix_ccs = 1
   type, public :: fgsl_spmatrix
     private
     type(c_ptr) :: gsl_spmatrix = c_null_ptr
diff --git a/fgsl_utils.c b/fgsl_utils.c
index fe79c92..9286b75 100644
--- a/fgsl_utils.c
+++ b/fgsl_utils.c
@@ -34,6 +34,7 @@
 #include <gsl/gsl_multifit.h>
 #include <gsl/gsl_multifit_nlinear.h>
 #include <gsl/gsl_multilarge.h>
+#include <gsl/gsl_multilarge_nlinear.h>
 #include <gsl/gsl_interp2d.h>
 #include <gsl/gsl_spmatrix.h>
 #include <gsl/gsl_splinalg.h>
@@ -972,7 +973,45 @@ void gsl_multifit_nlinear_fdf_get(gsl_multifit_nlinear_fdf *fdf,
 gsl_multifit_nlinear_type *gsl_multifit_nlinear_setup(char *s) {
     return gsl_multifit_nlinear_trust;
 }
+gsl_multilarge_nlinear_type *gsl_multilarge_nlinear_setup(char *s) {
+    return gsl_multilarge_nlinear_trust;
+}
+
+gsl_multilarge_nlinear_fdf *fgsl_multilarge_nlinear_fdf_cinit(
+       size_t ndim, size_t p, void *params, 
+       int (*f)(const gsl_vector *x, void *params, gsl_vector *f), 
+       int (*df)(const gsl_vector *x, void *params, gsl_matrix *df), 
+       int (*fvv)(const gsl_vector *x, const gsl_vector *v, void *params, gsl_vector *vv) 
+       ) {
+    gsl_multilarge_nlinear_fdf *result;
+    result = (gsl_multilarge_nlinear_fdf *) malloc(sizeof(gsl_multilarge_nlinear_fdf));
+    result->f = f;
+    result->df = df;
+    result->fvv = fvv;
+    result->n = ndim;
+    result->p = p;
+    result->params = params;
+    return result;
+}
+
+void fgsl_multilarge_nlinear_fdf_cfree(gsl_multilarge_nlinear_fdf *fun) {
+    free(fun);
+}
 
+void gsl_multilarge_nlinear_fdf_get(gsl_multilarge_nlinear_fdf *fdf, 
+       int (**fp)(const gsl_vector *x, void *params, gsl_vector *f),
+       int (**dfp)(CBLAS_TRANSPOSE_t TransJ, const gsl_vector *x, const gsl_vector *u, 
+           void *params, gsl_vector *v, gsl_matrix *jtj),
+       int (**fvvp)(const gsl_vector *x, const gsl_vector *v, void *params, gsl_vector *vv),
+       size_t *n, size_t *p, void **params, size_t *nevalf, size_t *nevaldfu, 
+       size_t *nevaldf2, size_t *nevalfvv) {
+
+       *fp = fdf->f; *dfp = fdf->df; *fvvp = fdf->fvv;
+       *n = fdf->n; *p = fdf->p; 
+       *params = fdf->params;
+       *nevalf = fdf->nevalf; *nevaldfu = fdf->nevaldfu; 
+       *nevaldf2 = fdf->nevaldf2; *nevalfvv = fdf->nevalfvv;
+}
 
 const gsl_multiroot_fsolver_type *fgsl_aux_multiroot_fsolver_alloc(int i) {
     const gsl_multiroot_fsolver_type *res;
@@ -994,7 +1033,66 @@ const gsl_multiroot_fsolver_type *fgsl_aux_multiroot_fsolver_alloc(int i) {
 	    break;
     }
     return res;
+}
 
+const gsl_multilarge_nlinear_trs *gsl_multilarge_nlinear_get_trs(int i) {
+    const gsl_multilarge_nlinear_trs *res;
+    switch(i) {
+	case 1:
+	    res = gsl_multilarge_nlinear_trs_lm;
+	    break;
+	case 2:
+	    res = gsl_multilarge_nlinear_trs_lmaccel;
+	    break;
+	case 3:
+	    res = gsl_multilarge_nlinear_trs_dogleg;
+	    break;
+	case 4:
+	    res = gsl_multilarge_nlinear_trs_ddogleg;
+	    break;
+	case 5:
+	    res = gsl_multilarge_nlinear_trs_subspace2D;
+	    break;
+	case 6:
+	    res = gsl_multilarge_nlinear_trs_cgst;
+	    break;
+	default:
+	    res = NULL;
+	    break;
+    }
+    return res;
+}
+
+const gsl_multilarge_nlinear_scale *gsl_multilarge_nlinear_get_scale(int i) {
+    const gsl_multilarge_nlinear_scale *res;
+    switch(i) {
+	case 1:
+	    res = gsl_multilarge_nlinear_scale_levenberg;
+	    break;
+	case 2:
+	    res = gsl_multilarge_nlinear_scale_marquardt;
+	    break;
+	case 3:
+	    res = gsl_multilarge_nlinear_scale_more;
+	    break;
+	default:
+	    res = NULL;
+	    break;
+    }
+    return res;
+}
+
+const gsl_multilarge_nlinear_solver *gsl_multilarge_nlinear_get_solver(int i) {
+    const gsl_multilarge_nlinear_solver *res;
+    switch(i) {
+	case 1:
+	    res = gsl_multilarge_nlinear_solver_cholesky;
+	    break;
+	default:
+	    res = NULL;
+	    break;
+    }
+    return res;
 }
 
 const gsl_multiroot_fdfsolver_type *fgsl_aux_multiroot_fdfsolver_alloc(int i) {
@@ -1017,7 +1115,6 @@ const gsl_multiroot_fdfsolver_type *fgsl_aux_multiroot_fdfsolver_alloc(int i) {
 	    break;
     }
     return res;
-
 }
 
 const gsl_multilarge_linear_type *fgsl_aux_multilarge_linear_alloc(int i) {
@@ -1191,6 +1288,11 @@ const gsl_splinalg_itersolve_type *fgsl_aux_splinalg_itersolve_alloc(int i) {
   return res;
 }
 
+void gsl_spmatrix_size(gsl_spmatrix *m, size_t *n1, size_t *n2) {
+  *n1 = m->size1;
+  *n2 = m->size2;
+}
+
 
 gsl_ntuple_select_fn *fgsl_ntuple_select_fn_cinit(int (*func)(void *data, void *params), void *params) {
     gsl_ntuple_select_fn *result;
diff --git a/interface/generics.finc b/interface/generics.finc
index 524db6a..3054751 100644
--- a/interface/generics.finc
+++ b/interface/generics.finc
@@ -68,6 +68,9 @@
   interface fgsl_multifit_nlinear_type
      module procedure fgsl_multifit_nlinear_setup
   end interface
+  interface fgsl_multilarge_nlinear_type
+     module procedure fgsl_multilarge_nlinear_setup
+  end interface
   interface fgsl_sizeof
      module procedure fgsl_sizeof_double
      module procedure fgsl_sizeof_float
diff --git a/interface/linalg.finc b/interface/linalg.finc
index 4eba13a..f379597 100644
--- a/interface/linalg.finc
+++ b/interface/linalg.finc
@@ -215,6 +215,19 @@
     type(c_ptr), value :: qr, tau, p, x
     integer(c_int) :: gsl_linalg_qrpt_svx
   end function gsl_linalg_qrpt_svx
+  function gsl_linalg_qrpt_lssolve (qr, tau, p, b, x, r) &
+       bind(c, name='gsl_linalg_QRPT_lssolve')
+    import :: c_ptr, c_int
+    type(c_ptr), value :: qr, tau, p, b, x, r
+    integer(c_int) :: gsl_linalg_qrpt_lssolve
+  end function gsl_linalg_qrpt_lssolve
+  function gsl_linalg_qrpt_lssolve2 (qr, tau, p, b, rank, x, r) &
+       bind(c, name='gsl_linalg_QRPT_lssolve2')
+    import :: c_ptr, c_int, c_size_t
+    type(c_ptr), value :: qr, tau, p, b, x, r
+    integer(c_size_t), value :: rank
+    integer(c_int) :: gsl_linalg_qrpt_lssolve2
+  end function gsl_linalg_qrpt_lssolve2
   function gsl_linalg_qrpt_qrsolve (q, r, p, b, x) &
        bind(c, name='gsl_linalg_QRPT_QRsolve')
     import :: c_ptr, c_int
@@ -239,6 +252,59 @@
     type(c_ptr), value :: qr, p, x
     integer(c_int) :: gsl_linalg_qrpt_rsvx
   end function gsl_linalg_qrpt_rsvx
+  function gsl_linalg_qrpt_rank (qr, tol) &
+       bind(c, name='gsl_linalg_QRPT_rank')
+    import :: c_ptr, c_size_t, c_double
+    type(c_ptr), value :: qr
+    real(c_double), value :: tol
+    integer(c_size_t) :: gsl_linalg_qrpt_rank
+  end function gsl_linalg_qrpt_rank
+  function gsl_linalg_qrpt_rcond (qr, rcond, wk) &
+       bind(c, name='gsl_linalg_QRPT_rcond')
+    import :: c_ptr, c_int, c_double
+    type(c_ptr), value :: qr, wk
+    real(c_double) :: rcond
+    integer(c_int) :: gsl_linalg_qrpt_rcond
+  end function gsl_linalg_qrpt_rcond
+!
+! Complete orthogonal decomposition (COD)
+!
+  function gsl_linalg_cod_decomp (a, tau_q, tau_z, p, rank, work) &
+       bind(c, name='gsl_linalg_COD_decomp')
+    import :: c_ptr, c_int, c_size_t
+    type(c_ptr), value :: a, tau_q, tau_z, p, work
+    integer(c_size_t) :: rank
+    integer(c_int) :: gsl_linalg_cod_decomp
+  end function gsl_linalg_cod_decomp
+  function gsl_linalg_cod_decomp_e (a, tau_q, tau_z, p, tol, rank, work) &
+       bind(c, name='gsl_linalg_COD_decomp_e')
+    import :: c_ptr, c_int, c_size_t, c_double
+    type(c_ptr), value :: a, tau_q, tau_z, p, work
+    real(c_double), value :: tol
+    integer(c_size_t) :: rank
+    integer(c_int) :: gsl_linalg_cod_decomp_e
+  end function gsl_linalg_cod_decomp_e
+  function gsl_linalg_cod_lssolve (qrz, tau_q, tau_z, p, rank, b, x, residual) &
+       bind(c, name='gsl_linalg_COD_lssolve')
+    import :: c_ptr, c_int, c_size_t
+    type(c_ptr), value :: qrz, tau_q, tau_z, p, b, x, residual
+    integer(c_size_t), value :: rank
+    integer(c_int) :: gsl_linalg_cod_lssolve
+  end function gsl_linalg_cod_lssolve
+  function gsl_linalg_cod_unpack (qrz, tau_q, tau_z, p, rank, q, r, z) &
+       bind(c, name='gsl_linalg_COD_unpack')
+    import :: c_ptr, c_int, c_size_t
+    type(c_ptr), value :: qrz, tau_q, tau_z, p, q, r, z
+    integer(c_size_t), value :: rank
+    integer(c_int) :: gsl_linalg_cod_unpack
+  end function gsl_linalg_cod_unpack
+  function gsl_linalg_cod_matz (qrz, tau_z, rank, a, work) &
+       bind(c, name='gsl_linalg_COD_matZ')
+    import :: c_ptr, c_int, c_size_t
+    type(c_ptr), value :: qrz, tau_z, a, work
+    integer(c_size_t), value :: rank
+    integer(c_int) :: gsl_linalg_cod_matz
+  end function gsl_linalg_cod_matz
 !
 ! SVD
 !
@@ -279,6 +345,12 @@
 !
 ! Cholesky
 !
+  function gsl_linalg_cholesky_decomp1 (a) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a
+    integer(c_int) :: gsl_linalg_cholesky_decomp1
+  end function gsl_linalg_cholesky_decomp1
   function gsl_linalg_cholesky_decomp (a) &
        bind(c)
     import :: c_ptr, c_int
@@ -325,6 +397,127 @@
     integer(c_int) :: gsl_linalg_complex_cholesky_invert
     type(c_ptr), value :: chol
   end function gsl_linalg_complex_cholesky_invert
+  function gsl_linalg_cholesky_decomp2 (a,s) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a,s
+    integer(c_int) :: gsl_linalg_cholesky_decomp2
+  end function gsl_linalg_cholesky_decomp2
+  function gsl_linalg_cholesky_solve2 (chol, s, b, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: chol, s, b, x
+    integer(c_int) :: gsl_linalg_cholesky_solve2
+  end function gsl_linalg_cholesky_solve2
+  function gsl_linalg_cholesky_svx2 (chol, s, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: chol, s, x
+    integer(c_int) :: gsl_linalg_cholesky_svx2
+  end function gsl_linalg_cholesky_svx2
+  function gsl_linalg_cholesky_scale (a,s) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a,s
+    integer(c_int) :: gsl_linalg_cholesky_scale
+  end function gsl_linalg_cholesky_scale
+  function gsl_linalg_cholesky_scale_apply (a,s) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a,s
+    integer(c_int) :: gsl_linalg_cholesky_scale_apply
+  end function gsl_linalg_cholesky_scale_apply
+  function gsl_linalg_cholesky_rcond (chol, rcond, work) &
+       bind(c)
+    import :: c_ptr, c_int, c_double
+    type(c_ptr), value :: chol, work
+    real(c_double), intent(inout) :: rcond
+    integer(c_int) :: gsl_linalg_cholesky_rcond
+  end function gsl_linalg_cholesky_rcond
+!
+! pivoted Cholesky
+!
+  function gsl_linalg_pcholesky_decomp (a, p) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a, p
+    integer(c_int) :: gsl_linalg_pcholesky_decomp
+  end function gsl_linalg_pcholesky_decomp
+  function gsl_linalg_pcholesky_solve (ldlt, p, b, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: ldlt, p, b, x
+    integer(c_int) :: gsl_linalg_pcholesky_solve
+  end function gsl_linalg_pcholesky_solve
+  function gsl_linalg_pcholesky_svx (ldlt, p, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: ldlt, p, x
+    integer(c_int) :: gsl_linalg_pcholesky_svx
+  end function gsl_linalg_pcholesky_svx
+  function gsl_linalg_pcholesky_invert (ldlt, p, ainv) bind(c)
+    import :: c_ptr, c_int
+    integer(c_int) :: gsl_linalg_pcholesky_invert
+    type(c_ptr), value :: ldlt, p, ainv
+  end function gsl_linalg_pcholesky_invert
+  function gsl_linalg_pcholesky_decomp2 (a, p,s) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a, p,s
+    integer(c_int) :: gsl_linalg_pcholesky_decomp2
+  end function gsl_linalg_pcholesky_decomp2
+  function gsl_linalg_pcholesky_solve2 (ldlt, p, s, b, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: ldlt, p, s, b, x
+    integer(c_int) :: gsl_linalg_pcholesky_solve2
+  end function gsl_linalg_pcholesky_solve2
+  function gsl_linalg_pcholesky_svx2 (ldlt, p, s, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: ldlt, p, s, x
+    integer(c_int) :: gsl_linalg_pcholesky_svx2
+  end function gsl_linalg_pcholesky_svx2
+  function gsl_linalg_pcholesky_rcond (ldlt, p, rcond, work) &
+       bind(c)
+    import :: c_ptr, c_int, c_double
+    type(c_ptr), value :: ldlt, p, work
+    real(c_double), intent(inout) :: rcond
+    integer(c_int) :: gsl_linalg_pcholesky_rcond
+  end function gsl_linalg_pcholesky_rcond
+!
+! modified Cholesky
+!
+  function gsl_linalg_mcholesky_decomp (a, p, e) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: a, p, e
+    integer(c_int) :: gsl_linalg_mcholesky_decomp
+  end function gsl_linalg_mcholesky_decomp
+  function gsl_linalg_mcholesky_solve (ldlt, p, b, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: ldlt, p, b, x
+    integer(c_int) :: gsl_linalg_mcholesky_solve
+  end function gsl_linalg_mcholesky_solve
+  function gsl_linalg_mcholesky_svx (ldlt, p, x) &
+       bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: ldlt, p, x
+    integer(c_int) :: gsl_linalg_mcholesky_svx
+  end function gsl_linalg_mcholesky_svx
+  function gsl_linalg_mcholesky_invert (ldlt, p, ainv) bind(c)
+    import :: c_ptr, c_int
+    integer(c_int) :: gsl_linalg_mcholesky_invert
+    type(c_ptr), value :: ldlt, p, ainv
+  end function gsl_linalg_mcholesky_invert
+  function gsl_linalg_mcholesky_rcond (ldlt, p, rcond, work) &
+       bind(c)
+    import :: c_ptr, c_int, c_double
+    type(c_ptr), value :: ldlt, p, work
+    real(c_double), intent(inout) :: rcond
+    integer(c_int) :: gsl_linalg_mcholesky_rcond
+  end function gsl_linalg_mcholesky_rcond
 !
 ! Tridiag
 !
@@ -530,3 +723,32 @@
     integer(c_size_t), value :: i, j
     real(c_double), value :: c, s
   end subroutine gsl_linalg_givens_gv
+!
+! Triangular matrices
+! 
+  integer(c_int) function gsl_linalg_tri_upper_invert(t) bind(c)
+    import :: c_int, c_ptr
+    type(c_ptr), value :: t
+  end function
+  integer(c_int) function gsl_linalg_tri_lower_invert(t) bind(c)
+    import :: c_int, c_ptr
+    type(c_ptr), value :: t
+  end function
+  integer(c_int) function gsl_linalg_tri_upper_unit_invert(t) bind(c)
+    import :: c_int, c_ptr
+    type(c_ptr), value :: t
+  end function
+  integer(c_int) function gsl_linalg_tri_lower_unit_invert(t) bind(c)
+    import :: c_int, c_ptr
+    type(c_ptr), value :: t
+  end function
+  integer(c_int) function gsl_linalg_tri_upper_rcond(t, rcond, work) bind(c)
+    import :: c_int, c_double, c_ptr
+    type(c_ptr), value :: t, work
+    real(c_double) :: rcond
+  end function
+  integer(c_int) function gsl_linalg_tri_lower_rcond(t, rcond, work) bind(c)
+    import :: c_int, c_double, c_ptr
+    type(c_ptr), value :: t, work
+    real(c_double) :: rcond
+  end function
diff --git a/interface/multifit.finc b/interface/multifit.finc
index b6a8176..3172f5a 100644
--- a/interface/multifit.finc
+++ b/interface/multifit.finc
@@ -252,6 +252,14 @@
     real(c_double) :: chisq
     integer(c_int) :: gsl_multifit_linear
   end function gsl_multifit_linear
+  function gsl_multifit_linear_tsvd(x, y, tol, c, cov, chisq, rank, work) bind(c)
+    import :: c_ptr, c_int, c_double, c_size_t
+    type(c_ptr), value :: x, y, c, cov, work
+    real(c_double), value :: tol
+    real(c_double) :: chisq
+    integer(c_size_t) :: rank
+    integer(c_int) :: gsl_multifit_linear_tsvd
+  end function gsl_multifit_linear_tsvd
   function gsl_multifit_linear_svd(x, work) bind(c)
     import :: c_ptr, c_int, c_double, c_size_t
     type(c_ptr), value :: x, work
@@ -339,6 +347,33 @@
     integer(c_size_t) :: idx
     integer(c_int) :: gsl_multifit_linear_lcorner2
   end function gsl_multifit_linear_lcorner2
+  integer(c_int) function gsl_multifit_linear_gcv_init(y, reg_param, uty, delta0, work) bind(c)
+    import :: c_int, c_ptr, c_double
+    type(c_ptr), value :: y, reg_param, uty, work
+    real(c_double) :: delta0
+  end function
+  integer(c_int) function gsl_multifit_linear_gcv_curve(reg_param, uty, delta0, g, work) bind(c)
+    import :: c_int, c_ptr, c_double
+    type(c_ptr), value :: reg_param, uty, g, work
+    real(c_double), value :: delta0
+  end function
+  integer(c_int) function gsl_multifit_linear_gcv_min(reg_param, uty, delta0, g, &
+          lambda, work) bind(c)
+    import :: c_int, c_ptr, c_double
+    type(c_ptr), value :: reg_param, uty, g, work
+    real(c_double), value :: delta0
+    real(c_double) :: lambda
+  end function
+  real(c_double) function gsl_multifit_linear_gcv_calc(lambda, uty, delta0, work) bind(c)
+    import :: c_ptr, c_double
+    type(c_ptr), value :: uty, work
+    real(c_double), value :: delta0, lambda
+  end function
+  integer(c_int) function gsl_multifit_linear_gcv(y, reg_param, g, lambda, g_lambda, work) bind(c)
+    import :: c_int, c_ptr, c_double
+    type(c_ptr), value :: y, reg_param, g, work
+    real(c_double) :: lambda, g_lambda
+  end function
   function gsl_multifit_linear_lk(p, k, L) &
     bind(c, name='gsl_multifit_linear_Lk')
     import :: c_size_t, c_int, c_ptr
@@ -383,6 +418,14 @@
     real(c_double) :: chisq
     integer(c_int) :: gsl_multifit_wlinear
   end function gsl_multifit_wlinear
+  function gsl_multifit_wlinear_tsvd(x, w, y, tol, c, cov, chisq, rank, work) bind(c)
+    import :: c_ptr, c_int, c_double, c_size_t
+    type(c_ptr), value :: x, w, y, c, cov, work
+    real(c_double), value :: tol
+    real(c_double) :: chisq
+    integer(c_size_t) :: rank
+    integer(c_int) :: gsl_multifit_wlinear_tsvd
+  end function gsl_multifit_wlinear_tsvd
   function gsl_multifit_wlinear_svd(x, w, y, tol, rank, c, cov, chisq, work) bind(c)
     import :: c_ptr, c_int, c_double, c_size_t
     type(c_ptr), value :: x, w, y, c, cov, work
@@ -412,6 +455,11 @@
     type(c_ptr), value :: r
     integer(c_int) :: gsl_multifit_linear_residuals
   end function gsl_multifit_linear_residuals
+  integer(c_size_t) function gsl_multifit_linear_rank(tol, work) bind(c)
+    import :: c_ptr, c_size_t, c_double
+    type(c_ptr), value :: work
+    real(c_double), value :: tol
+  end function
   function gsl_multifit_fdfridge_alloc(T, n, p) bind(c)
     import :: c_ptr, c_size_t
     type(c_ptr), value :: T
diff --git a/interface/nlfit.finc b/interface/nlfit.finc
index ee4aa5e..71e372b 100644
--- a/interface/nlfit.finc
+++ b/interface/nlfit.finc
@@ -6,6 +6,10 @@
     import
     character(c_char) :: s
  end function
+ type(c_ptr) function gsl_multilarge_nlinear_setup(s) BIND(C)
+    import
+    character(c_char) :: s
+ end function
  function gsl_multifit_nlinear_alloc(t, params, n, p) BIND(C)
     import
     type(c_ptr), value :: t
@@ -124,20 +128,20 @@
     type(c_ptr), value :: w
  end function
  integer(c_int) function gsl_multifit_nlinear_test(xtol, gtol, ftol, info, w) BIND(C)
-    import
+    import :: c_ptr, c_double, c_int
     real(c_double), value :: xtol, gtol, ftol
     integer(c_int), intent(inout) :: info
     type(c_ptr), value :: w
  end function
  integer(c_int) function gsl_multilarge_nlinear_test(xtol, gtol, ftol, info, w) BIND(C)
-    import
+    import :: c_ptr, c_double, c_int
     real(c_double), value :: xtol, gtol, ftol
     integer(c_int), intent(inout) :: info
     type(c_ptr), value :: w
  end function
  integer(c_int) function gsl_multifit_nlinear_driver(maxiter, xtol, gtol, ftol, &
                        callback, callback_params, info, w) BIND(C)
-    import
+    import :: c_ptr, c_double, c_int, c_size_t, c_funptr
     integer(c_size_t), value :: maxiter
     real(c_double), value :: xtol, gtol, ftol
     type(c_funptr), value :: callback
@@ -146,7 +150,7 @@
  end function
  integer(c_int) function gsl_multilarge_nlinear_driver(maxiter, xtol, gtol, ftol, &
                        callback, callback_params, info, w) BIND(C)
-    import
+    import :: c_ptr, c_double, c_int, c_size_t, c_funptr
     integer(c_size_t), value :: maxiter
     real(c_double), value :: xtol, gtol, ftol
     type(c_funptr), value :: callback
@@ -192,3 +196,33 @@
     import
     integer(c_int), value :: which
  end function
+ type(c_ptr) function fgsl_multilarge_nlinear_fdf_cinit(ndim, p, params, fp, dfp, fvvp) BIND(C)
+    import
+    integer(c_size_t), value :: ndim, p
+    type(c_ptr), value :: params
+    type(c_funptr), value :: fp, dfp, fvvp
+ end function
+ subroutine fgsl_multilarge_nlinear_fdf_cfree(fdf) BIND(C)
+    import
+    type(c_ptr), value :: fdf
+ end subroutine
+ subroutine gsl_multilarge_nlinear_fdf_get(fdf, fp, dfp, fvvp, &
+          n, p, params, nevalf, nevaldfu, nevaldf2, nevalfvv) BIND(C)
+    import
+    type(c_ptr), value :: fdf
+    type(c_funptr) :: fp, dfp, fvvp
+    integer(c_size_t) :: n, p, nevalf, nevaldfu, nevaldf2, nevalfvv
+    type(c_ptr) :: params
+ end subroutine
+ type(c_ptr) function gsl_multilarge_nlinear_get_trs(which) BIND(C)
+    import
+    integer(c_int), value :: which
+ end function
+ type(c_ptr) function gsl_multilarge_nlinear_get_scale(which) BIND(C)
+    import
+    integer(c_int), value :: which
+ end function
+ type(c_ptr) function gsl_multilarge_nlinear_get_solver(which) BIND(C)
+    import
+    integer(c_int), value :: which
+ end function
diff --git a/interface/ode.finc b/interface/ode.finc
index eed7b58..9f58e48 100644
--- a/interface/ode.finc
+++ b/interface/ode.finc
@@ -137,8 +137,8 @@
      function gsl_odeiv2_evolve_apply_fixed_step(e, con, step, dydt, t, h0, y) bind(c)
        import
        type(c_ptr), value :: e, con, step, dydt,  y
-       real(c_double) :: t, h0
-       value :: h0
+       real(c_double) :: t
+       real(c_double), value :: h0
        integer(c_int) :: gsl_odeiv2_evolve_apply_fixed_step
      end function gsl_odeiv2_evolve_apply_fixed_step
      function gsl_odeiv2_evolve_reset(s) bind(c)
diff --git a/interface/permutation.finc b/interface/permutation.finc
index 00092bb..f6aa651 100644
--- a/interface/permutation.finc
+++ b/interface/permutation.finc
@@ -111,6 +111,11 @@
     type(c_ptr), value :: p,v
     integer(c_int) :: gsl_permute_vector_inverse
   end function gsl_permute_vector_inverse
+  function gsl_permute_matrix(p,v) bind(c)
+    import :: c_ptr, c_int
+    type(c_ptr), value :: p,v
+    integer(c_int) :: gsl_permute_matrix
+  end function gsl_permute_matrix
   function gsl_permutation_mul(p, pa, pb) bind(c)
     import
     type(c_ptr), value :: p
diff --git a/interface/rng.finc b/interface/rng.finc
index 2bab6ca..e11e928 100644
--- a/interface/rng.finc
+++ b/interface/rng.finc
@@ -244,6 +244,32 @@
     real(c_double), value :: x, y, sigma_x, sigma_y, rho
     real(c_double) :: gsl_ran_bivariate_gaussian_pdf
   end function gsl_ran_bivariate_gaussian_pdf
+  integer(c_int) function gsl_ran_multivariate_gaussian(r, mu, l, result) bind(c)
+    import
+    type(c_ptr), value :: r, mu, l, result
+  end function gsl_ran_multivariate_gaussian
+  integer(c_int) function gsl_ran_multivariate_gaussian_pdf( &
+       x, mu, l, result, work) bind(c)
+    import
+    type(c_ptr), value :: x, mu, l, work
+    real(fgsl_double) ::  result
+  end function gsl_ran_multivariate_gaussian_pdf
+  integer(c_int) function gsl_ran_multivariate_gaussian_log_pdf( &
+       x, mu, l, result, work) bind(c)
+    import
+    type(c_ptr), value :: x, mu, l, work
+    real(fgsl_double) ::  result
+  end function gsl_ran_multivariate_gaussian_log_pdf
+  integer(c_int) function gsl_ran_multivariate_gaussian_mean( &
+       x, mu_hat) bind(c)
+    import
+    type(c_ptr), value :: x, mu_hat
+  end function gsl_ran_multivariate_gaussian_mean
+  integer(c_int) function gsl_ran_multivariate_gaussian_vcov( &
+       x, sigma_hat) bind(c)
+    import
+    type(c_ptr), value :: x, sigma_hat
+  end function gsl_ran_multivariate_gaussian_vcov
   function gsl_ran_exponential(r, mu) bind(c)
     import
     type(c_ptr), value :: r
diff --git a/interface/rstat.finc b/interface/rstat.finc
index 01ed5a9..61a7e21 100644
--- a/interface/rstat.finc
+++ b/interface/rstat.finc
@@ -7,6 +7,11 @@ function gsl_rstat_quantile_alloc(p) bind(c)
   real(c_double), value :: p
   type(c_ptr) :: gsl_rstat_quantile_alloc
 end function gsl_rstat_quantile_alloc
+function gsl_rstat_quantile_reset(w) bind(c)
+  import :: c_ptr, c_int
+  type(c_ptr), value  :: w
+  integer(c_int) :: gsl_rstat_quantile_reset
+end function gsl_rstat_quantile_reset
 subroutine gsl_rstat_quantile_free(w) bind(c)
   import :: c_ptr
   type(c_ptr), value :: w
@@ -56,6 +61,11 @@ function gsl_rstat_mean(w) bind(c)
   type(c_ptr), value :: w
   real(c_double) :: gsl_rstat_mean
 end function gsl_rstat_mean
+function gsl_rstat_rms(w) bind(c)
+  import :: c_ptr, c_double
+  type(c_ptr), value :: w
+  real(c_double) :: gsl_rstat_rms
+end function gsl_rstat_rms
 function gsl_rstat_variance(w) bind(c)
   import :: c_ptr, c_double
   type(c_ptr), value :: w
diff --git a/interface/spmatrix.finc b/interface/spmatrix.finc
index 4c4f925..4241b3d 100644
--- a/interface/spmatrix.finc
+++ b/interface/spmatrix.finc
@@ -12,6 +12,11 @@ function gsl_spmatrix_alloc_nzmax(n1, n2, nzmax, flags) bind(c)
   integer(c_size_t), value :: n1, n2, nzmax, flags
   type(c_ptr) :: gsl_spmatrix_alloc_nzmax
 end function gsl_spmatrix_alloc_nzmax
+subroutine gsl_spmatrix_size(m, n1, n2) bind(c)
+  import :: c_ptr, c_size_t
+  type(c_ptr), value :: m
+  integer(c_size_t) :: n1, n2
+end subroutine gsl_spmatrix_size
 subroutine gsl_spmatrix_free(m) bind(c)
   import :: c_ptr
   type(c_ptr), value :: m
@@ -57,8 +62,8 @@ function gsl_spmatrix_set(m, i, j, x) bind(c)
 end function gsl_spmatrix_set
 function gsl_spmatrix_compcol(T) bind(c)
   import :: c_ptr
-  type(c_ptr) :: T, gsl_spmatrix_compcol
-  value :: T
+  type(c_ptr) :: gsl_spmatrix_compcol
+  type(c_ptr), value :: T
 end function gsl_spmatrix_compcol
 subroutine gsl_spmatrix_cumsum(n, c) bind(c)
   import :: c_size_t, c_ptr
@@ -102,3 +107,16 @@ function gsl_spmatrix_transpose_memcpy(dest, src) bind(c)
   type(c_ptr), value :: dest, src
   integer(c_int) :: gsl_spmatrix_transpose_memcpy
 end function gsl_spmatrix_transpose_memcpy
+function gsl_spblas_dgemv(transa, alpha, a, x, beta, y) bind(c)
+  import :: c_ptr, c_int, c_double
+  integer(c_int), value :: transa
+  type(c_ptr), value :: a, x, y
+  real(c_double), value :: alpha, beta
+  integer(c_int) :: gsl_spblas_dgemv
+end function gsl_spblas_dgemv
+function gsl_spblas_dgemm(alpha, a, b, c) bind(c)
+  import :: c_ptr, c_int, c_double
+  type(c_ptr), value :: a, b, c
+  real(c_double), value :: alpha
+  integer(c_int) :: gsl_spblas_dgemm
+end function gsl_spblas_dgemm
diff --git a/tests/linalg.F90 b/tests/linalg.F90
index b56688c..4d585b0 100644
--- a/tests/linalg.F90
+++ b/tests/linalg.F90
@@ -6,20 +6,26 @@ program linalg
   complex(fgsl_double), parameter :: ai = (0.0d0, 1.0d0), ui=(1.0d0, 0.0d0)
   type(fgsl_matrix) :: a, a_orig, inv, q, r
   type(fgsl_matrix_complex) :: ac, ac_orig
-  type(fgsl_vector) :: b, x, res, tau
+  type(fgsl_vector) :: b, x, res, tau, sd, work
   type(fgsl_vector_complex) :: bc, xc, resc
   real(fgsl_double), target :: af(3, 3), af_orig(3, 3), bf(3), xf(3), resf(3), &
-       invf(3, 3), tauf(3), qf(3, 3), rf(3, 3), mf(3, 3)
+       invf(3, 3), tauf(3), qf(3, 3), rf(3, 3), mf(3, 3), sdf(3), workf(9)
   complex(fgsl_double) :: acf(3, 3), xcf(3), bcf(3), rescf(3), acf_orig(3, 3)
-  real(fgsl_double) :: det, lndet
+  real(fgsl_double) :: det, lndet, rcond
   type(fgsl_permutation) :: p
   integer(fgsl_int) :: status, signum, ip(3), sgn
-  integer(fgsl_size_t) :: i
+  integer(fgsl_size_t) :: i, rank
 !
 ! Test linear algebra
 ! remember that matrices are transposed vs. usual Fortran convention
 !
   call unit_init(200)
+  sd = fgsl_vector_init(1.0_fgsl_double)
+  work = fgsl_vector_init(1.0_fgsl_double)
+  status = fgsl_vector_align(sdf, 3_fgsl_size_t, sd, 3_fgsl_size_t, &
+       0_fgsl_size_t, 1_fgsl_size_t)
+  status = fgsl_vector_align(workf, 9_fgsl_size_t, work, 9_fgsl_size_t, &
+       0_fgsl_size_t, 1_fgsl_size_t)
 !
 ! LU - real
 !
@@ -248,6 +254,12 @@ program linalg
   call unit_assert_equal('fgsl_linalg_qrpt_svx:status',fgsl_success,status)
   call unit_assert_equal_within('fgsl_linalg_qrpt_svx:x',&
        (/1.0d0/3.0d0, 2.0d0/3.0d0, 5.0d0/3.0d0 /),xf,eps10)
+  status = fgsl_linalg_qrpt_lssolve(a, tau, p, b, x, sd)
+  call unit_assert_equal('fgsl_linalg_qrpt_lssolve:status',fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_qrpt_lssolve:x',&
+       (/1.0d0/3.0d0, 2.0d0/3.0d0, 5.0d0/3.0d0 /),xf,eps10)
+  call unit_assert_equal_within('fgsl_linalg_qrpt_lssolve:residual',&
+       (/0.0d0, 0.0d0, 0.0d0 /),sdf,eps10)
   status = fgsl_linalg_qrpt_rsolve(a, p, b, x)
   call unit_assert_equal('fgsl_linalg_qrpt_rsolve:status',fgsl_success,status)
   call unit_assert_equal_within('fgsl_linalg_qrpt_rsolve:x',&
@@ -278,6 +290,34 @@ program linalg
   af(3, 1:3) = af_orig(3, 1:3) + invf(3, 1:3)
   call unit_assert_equal_within('fgsl_linalg_qrpt_update:mf',&
        af,mf,eps10)
+
+  af = reshape((/ 1.0d0, 1.0d0, 0.0d0, 1.0d0, 0.0d0, 1.0d0, &
+       2.0d0, 1.0d0, 1.0d0 /), (/3, 3/))
+  status = fgsl_linalg_qrpt_decomp(a, tau, p, signum, res)
+  call unit_assert_equal('fgsl_linalg_qrpt_decomp(singular):status',fgsl_success,status)
+  rank = fgsl_linalg_qrpt_rank(a, tol=-1.0d0)
+  call unit_assert_equal('fgsl_linalg_qrpt_rank',2_fgsl_size_t,rank)
+  status = fgsl_linalg_qrpt_lssolve2(a, tau, p, b, rank, x, sd)
+  call unit_assert_equal('fgsl_linalg_qrpt_lssolve2:status',fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_qrpt_lssolve2:x',&
+       (/2.0d0, -1.0d0, 0.0d0 /),xf,eps10)
+  call unit_assert_equal_within('fgsl_linalg_qrpt_lssolve2:residual',&
+       (/0.0d0, 0.0d0, 0.0d0 /),sdf,eps10)
+!
+! COD
+!
+  af = af_orig
+  status = fgsl_linalg_cod_decomp(a, tau, sd, p, rank, x)
+  call unit_assert_equal('fgsl_linalg_cod_decomp:status',fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_cod_decomp:tau_q',&
+       (/ 1.44721359549996d0, 1.74535599249993d0, 0.00000000000000d0 /),tauf,eps10)
+  call unit_assert_equal_within('fgsl_linalg_cod_decomp:tau_z',&
+       (/ -3.525431591703201d-16, -3.525431591703199d-16, 3.525431591703201d-16 /),sdf,eps10)
+  status = fgsl_linalg_cod_lssolve(a, tau, sd, p, rank, b, x, res)
+  call unit_assert_equal_within('fgsl_linalg_cod_decomp:x',&
+       (/1.0d0/3.0d0, 2.0d0/3.0d0, 5.0d0/3.0d0 /),xf,eps10)
+  call unit_assert_equal_within('fgsl_linalg_cod_decomp:res',&
+       (/0.0d0, 0.0d0, 0.0d0 /),resf,eps10)
 !
 ! SVD
 !
@@ -315,21 +355,115 @@ program linalg
   af_orig = reshape((/2.0d0, 0.0d0, 0.0d0, 0.0d0, 2.0d0, 0.0d0, &
        1.0d0, 1.0d0, 2.0d0 /), (/3, 3/))
   af = af_orig
-  status = fgsl_linalg_cholesky_decomp(a)
-  call unit_assert_equal('fgsl_linalg_cholesky_decomp:status', &
+  status = fgsl_linalg_cholesky_decomp1(a)
+  call unit_assert_equal('fgsl_linalg_cholesky_decomp1:status', &
        fgsl_success,status)
   status = fgsl_linalg_cholesky_solve (a, b, x)
   call unit_assert_equal('fgsl_linalg_cholesky_solve:status',&
        fgsl_success,status)
-  call unit_assert_equal_within('fgsl_linalg_sv_solve:x',&
+  call unit_assert_equal_within('fgsl_linalg_cholesky_solve:x',&
        (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
   xf = bf
   status = fgsl_linalg_cholesky_svx (a, x)
   call unit_assert_equal('fgsl_linalg_cholesky_svx:status',&
        fgsl_success,status)
-  call unit_assert_equal_within('fgsl_linalg_sv_svx:x',&
+  call unit_assert_equal_within('fgsl_linalg_cholesky_svx:x',&
+       (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
+  status = fgsl_linalg_cholesky_rcond (a, rcond, work)
+  call unit_assert_equal('fgsl_linalg_cholesky_rcond:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_cholesky_rcond:rcond',&
+       .25d0,rcond,eps10)
+
+  af = af_orig
+  status = fgsl_linalg_cholesky_decomp2(a, sd)
+  call unit_assert_equal('fgsl_linalg_cholesky_decomp2:status', &
+       fgsl_success,status)
+  status = fgsl_linalg_cholesky_solve2 (a, sd, b, x)
+  call unit_assert_equal('fgsl_linalg_cholesky_solve2:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_cholesky_solve2:x',&
+       (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
+  xf = bf
+  status = fgsl_linalg_cholesky_svx2 (a, sd, x)
+  call unit_assert_equal('fgsl_linalg_cholesky_svx2:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_cholesky_svx2:x',&
+       (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
+!
+!  complex Cholesky
+!
+  acf = af_orig * ui
+  bcf = bf * ui
+  status = fgsl_linalg_complex_cholesky_decomp(ac)
+  call unit_assert_equal('fgsl_linalg_complex_cholesky_decomp:status', &
+       fgsl_success,status)
+  status = fgsl_linalg_complex_cholesky_solve (ac, bc, xc)
+  call unit_assert_equal('fgsl_linalg_complex_cholesky_solve:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_complex_cholesky_solve:xc',&
+       (/-.25d0, .25d0, 1.5d0 /)*ui,xcf,eps10)
+  xcf = bf * ui
+  status = fgsl_linalg_complex_cholesky_svx (ac, xc)
+  call unit_assert_equal('fgsl_linalg_complex_cholesky_svx:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_complex_cholesky_svx:x',&
+       (/-.25d0, .25d0, 1.5d0 /)*ui,xcf,eps10)
+!
+!  pivoted Cholesky
+!
+  af = af_orig
+  status = fgsl_linalg_pcholesky_decomp(a, p)
+  call unit_assert_equal('fgsl_linalg_pcholesky_decomp:status', &
+       fgsl_success,status)
+  status = fgsl_linalg_pcholesky_solve (a, p, b, x)
+  call unit_assert_equal('fgsl_linalg_pcholesky_solve:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_pcholesky_solve:x',&
        (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
-! FIXME complex Cholesky
+  xf = bf
+  status = fgsl_linalg_pcholesky_svx (a, p, x)
+  call unit_assert_equal('fgsl_linalg_pcholesky_svx:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_pcholesky_svx:x',&
+       (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
+  af = af_orig
+  status = fgsl_linalg_pcholesky_decomp2(a, p, sd)
+  call unit_assert_equal('fgsl_linalg_pcholesky_decomp2:status', &
+       fgsl_success,status)
+  status = fgsl_linalg_pcholesky_solve2 (a, p, sd, b, x)
+  call unit_assert_equal('fgsl_linalg_pcholesky_solve2:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_pcholesky_solve2:x',&
+       (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
+  xf = bf
+  status = fgsl_linalg_pcholesky_svx2 (a, p, sd, x)
+  call unit_assert_equal('fgsl_linalg_pcholesky_svx2:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_pcholesky_svx2:x',&
+       (/-.25d0, .25d0, 1.5d0 /),xf,eps10)
+!
+! modified Cholesky
+!
+  af_orig = reshape((/0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, &
+       1.0d0, 1.0d0, 2.0d0 /), (/3, 3/))
+  af = af_orig
+  status = fgsl_linalg_mcholesky_decomp(a, p, sd)
+  call unit_assert_equal('fgsl_linalg_mcholesky_decomp:status', &
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_mcholesky_decomp:work',&
+       (/1.d0, 2.d0, 0.d0 /),sdf,eps10)
+  status = fgsl_linalg_mcholesky_solve (a, p, b, x)
+  call unit_assert_equal('fgsl_linalg_mcholesky_solve:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_mcholesky_solve:x',&
+       (/-1.d0, 0.d0, 2.d0 /),xf,eps10)
+  xf = bf
+  status = fgsl_linalg_mcholesky_svx (a, p, x)
+  call unit_assert_equal('fgsl_linalg_mcholesky_svx:status',&
+       fgsl_success,status)
+  call unit_assert_equal_within('fgsl_linalg_mcholesky_svx:x',&
+       (/-1.d0, 0.d0, 2.d0 /),xf,eps10)
 
 !
 ! cleanup
@@ -344,6 +478,8 @@ program linalg
   call fgsl_matrix_free(ac_orig)
   call fgsl_vector_free(b)
   call fgsl_vector_free(bc)
+  call fgsl_vector_free(sd)
+  call fgsl_vector_free(work)
   call fgsl_vector_free(x)
   call fgsl_vector_free(xc)
   call fgsl_vector_free(res)
diff --git a/tests/mod_unit.f90 b/tests/mod_unit.f90
index b195841..3e7b7ae 100644
--- a/tests/mod_unit.f90
+++ b/tests/mod_unit.f90
@@ -5,6 +5,7 @@ module mod_unit
   public :: unit_debug, unit_init, unit_finalize, unit_assert_equal, &
        unit_assert_equal_within, unit_assert_true
   integer, parameter :: ik = selected_int_kind(6)
+  integer, parameter :: lk = selected_int_kind(14)
   integer, parameter :: dk = selected_real_kind(15)
 !FIXME - PGI does not accept  logical, protected :: unit_debug = .false.
   logical :: unit_debug = .false.
@@ -40,12 +41,12 @@ subroutine unit_init(num_tests)
   end subroutine unit_init
   subroutine unit_assert_equal_integer(description, ival_target, ival_check)
     character(len=*), intent(in) :: description
-    integer(kind=4), intent(in) :: ival_target, ival_check
+    integer(kind=ik), intent(in) :: ival_target, ival_check
     call unit_toggle(description, ival_target == ival_check)
   end subroutine unit_assert_equal_integer
   subroutine unit_assert_equal_long(description, ival_target, ival_check)
     character(len=*), intent(in) :: description
-    integer(kind=8), intent(in) :: ival_target, ival_check
+    integer(kind=lk), intent(in) :: ival_target, ival_check
     call unit_toggle(description, ival_target == ival_check)
   end subroutine unit_assert_equal_long
   subroutine unit_assert_equal_integer_array(description, &
diff --git a/tests/permutation.F90 b/tests/permutation.F90
index 4ca9ee9..5e36e61 100644
--- a/tests/permutation.F90
+++ b/tests/permutation.F90
@@ -8,11 +8,13 @@ program permutation
   type(fgsl_multiset) :: m1, m2
   type(fgsl_error_handler_t) :: std
   type(fgsl_vector) :: vec
+  type(fgsl_matrix) :: a
   type(fgsl_file) :: pfile
   integer(fgsl_int) :: status
   integer(fgsl_size_t) :: index
   integer(fgsl_size_t), pointer, dimension(:) :: p_data
   real(fgsl_double) :: da(8)
+  real(fgsl_double) :: af(4,4)
   integer(fgsl_long) :: ida(8)
 !
 ! Test permutations, combinations and multisets
@@ -102,6 +104,21 @@ program permutation
   call unit_assert_equal_within('fgsl_permute_inverse',&
        (/4.0d0,3.0d0,2.0d0,1.0d0/),da(1:4),eps10)
   call fgsl_vector_free(vec)
+  a = fgsl_matrix_init(1.0_fgsl_double)
+  status = fgsl_matrix_align(af, 4_fgsl_size_t, 4_fgsl_size_t, &
+       4_fgsl_size_t,a)
+  af = reshape( [ 1.0d0, 0.0d0, 0.0d0, 0.0d0, &
+                  0.0d0, 1.0d0, 0.0d0, 0.0d0, &
+                  0.0d0, 0.0d0, 1.0d0, 0.0d0, &
+                  0.0d0, 0.0d0, 0.0d0, 1.0d0 ], shape = [ 4, 4 ] )
+  status = fgsl_permute_matrix(p1, a)
+  call unit_assert_equal_within('fgsl_permute_matrix',&
+       reshape( [ 0.0d0, 0.0d0, 0.0d0, 1.0d0, &
+                  0.0d0, 0.0d0, 1.0d0, 0.0d0, &
+                  0.0d0, 1.0d0, 0.0d0, 0.0d0, &
+                  1.0d0, 0.0d0, 0.0d0, 0.0d0 ], shape = [ 4, 4 ] ),af,eps10)
+  call fgsl_matrix_free(a)
+
 ! multiply
   status = fgsl_permutation_mul(p3,p1,p2)
   call unit_assert_equal('fgsl_permute_mul:status',fgsl_success,status)
diff --git a/tests/rng.f90 b/tests/rng.f90
index 1284694..057acb1 100644
--- a/tests/rng.f90
+++ b/tests/rng.f90
@@ -5,8 +5,12 @@ program rng
   real(fgsl_double), parameter :: eps10 = 1.0d-10
   character(kind=fgsl_char, len=fgsl_strmax) :: name
   integer(fgsl_long) :: rn
-  real(fgsl_double) :: rd, rd_copy, re, rid, rie, qd(1), qd_copy(1)
+  real(fgsl_double) :: rd, rd_copy, re, rid, rie, qd(1), qd_copy(1), &
+       rho, sigma, x, y
   type(fgsl_rng) :: r, copy, clone
+  type(fgsl_vector) :: xx, mu, result, work
+  type(fgsl_matrix) :: l
+  real(fgsl_double) :: xxf(2), muf(2), resultf(2), lf(2,2)
   type(fgsl_qrng) :: q, qcopy, qclone
   type(fgsl_rng_type) :: t
   type(fgsl_file) :: rfile
@@ -118,6 +122,53 @@ program rng
 !  rd = fgsl_ran_ugaussian_ratio_method(r)
 !  call unit_assert_equal_within('fgsl_ran_ugaussian_ratio_method',&
 !       0.7630954425575337d0,rd,eps10)
+
+! FIXME: tests for some functions are missing here
+
+  sigma = 1.0d0
+  rho = 0.4d0
+  call fgsl_ran_bivariate_gaussian(r, sigma, 2.0d0*sigma, rho, x, y)
+  call unit_assert_equal_within('fgsl_ran_bivariate_gaussian',&
+       [ 1.10738714772252d0, -1.47314485089665d0 ],[ x, y ],eps10)
+  x = 0.1d0
+  y = 0.2d0
+  rd = fgsl_ran_bivariate_gaussian_pdf(x, y, sigma, 2.0d0*sigma, rho)
+  call unit_assert_equal_within('fgsl_ran_bivariate_gaussian',&
+        8.620816273107094d-2,rd,eps10)
+
+  xx = fgsl_vector_init(1.0_fgsl_double)
+  mu = fgsl_vector_init(1.0_fgsl_double)
+  result = fgsl_vector_init(1.0_fgsl_double)
+  work = fgsl_vector_init(1.0_fgsl_double)
+  l = fgsl_matrix_init(1.0_fgsl_double)
+  status = fgsl_vector_align(xxf, 2_fgsl_size_t, xx, 2_fgsl_size_t, &
+       0_fgsl_size_t, 1_fgsl_size_t)
+  status = fgsl_vector_align(muf, 2_fgsl_size_t, mu, 2_fgsl_size_t, &
+       0_fgsl_size_t, 1_fgsl_size_t)
+  status = fgsl_vector_align(resultf, 2_fgsl_size_t, result, 2_fgsl_size_t, &
+       0_fgsl_size_t, 1_fgsl_size_t)
+  status = fgsl_vector_align(resultf, 2_fgsl_size_t, work, 2_fgsl_size_t, &
+       0_fgsl_size_t, 1_fgsl_size_t)
+! Note result and work are aliased against resultf
+  status = fgsl_matrix_align(lf, 2_fgsl_size_t, 2_fgsl_size_t, &
+       2_fgsl_size_t,l)
+  muf = [ 0.2d0 , 0.8d0 ]
+  lf = reshape( [ 2.0d0, 0.0d0, 0.0d0, 1.0d0 ], [ 2, 2 ] )
+  status = fgsl_ran_multivariate_gaussian(r, mu, l, result)
+  call unit_assert_equal('fgsl_ran_multivariate_gaussian:status',fgsl_success,status)
+  call unit_assert_equal_within('fgsl_ran_multivariate_gaussian:result',&
+       [ 0.381341008927370d0, 0.705879273922799d0 ],resultf,eps10)
+  xxf = [ 0.7, -1.2 ]
+  status = fgsl_ran_multivariate_gaussian_pdf(xx, mu, l, rd, work)
+  call unit_assert_equal('fgsl_ran_multivariate_gaussian_pdf:status',fgsl_success,status)
+  call unit_assert_equal_within('fgsl_ran_multivariate_gaussian_pdf:result',&
+        1.043829169309120d-2,rd,eps10)
+  status = fgsl_ran_multivariate_gaussian_log_pdf(xx, mu, l, rd, work)
+  call unit_assert_equal('fgsl_ran_multivariate_gaussian_log_pdf:status',fgsl_success,status)
+  call unit_assert_equal_within('fgsl_ran_multivariate_gaussian_log_pdf:result',&
+        -4.56227434084661d0,rd,eps10)
+! FIXME tests for gaussian_mean and gaussian_vcov missing
+
   rd = fgsl_cdf_gaussian_p(2.0d0,1.0d0)
   call unit_assert_equal_within('fgsl_cdf_gaussian_p',&
        0.9772498680518208d0,rd,eps10)
@@ -147,6 +198,11 @@ program rng
   call fgsl_rng_free(r)
   call fgsl_rng_free(copy)
   call fgsl_rng_free(clone)
+  call fgsl_vector_free(xx)
+  call fgsl_vector_free(result)
+  call fgsl_vector_free(work)
+  call fgsl_vector_free(mu)
+  call fgsl_matrix_free(l)
   call fgsl_qrng_free(q)
   call fgsl_qrng_free(qcopy)
   call fgsl_qrng_free(qclone)