anchor.diff:  ! Subroutines in module:
   !     - load_array            Load a sub-array into a main array
   !     - load_constraints      Load x constraints into A and b
  +!     - solve_nnls            Implement non-negative least squares solver
  +!     - solve_jacobi          Implement Numerical Recipes symmetric eigensystem solver
   !     - solve_dgels           Implement LAPACK general least-squares solver
  -!     - solve_nnls            Implement non-negative least squares solver
   !     - solve_dgesv           Implement LAPACK general matrix solver
  +!     - solve_dsyev           Implement LAPACK symmetric matrix eigensystem solver
  +!     - solve_dsysv           Implement LAPACK symmetric matrix solver
   !     - solve_dgssv           Implement SuperLU sparse matrix solver
   !     - nnls                  Non-negative least squares solver
   !     - g1                    Subroutine used by nnls (orthogonal rotation matrix)
  
positive.code: 
  ! Solver subroutines
  public :: solve_nnls
  
  #ifdef USE_LAPACK
  public :: solve_dgels
  public :: solve_dgesv
  public :: solve_dsysv
  #endif
  
  
negative.code: double precision :: alocal(nrows,ncols), blocal(nrows,1)
  
  ierr = 0
  
  if (ncols.gt.nrows) then
      write(stderr,*) 'solve_dgels: ncols is greater than nrows'
      write(stderr,*) 'Only works for overdetermined systems (nrows>=ncols)'
      ierr = 1
      return
  endif
  
  ! Use local arrays for inversion
  alocal = Ain
  