anchor.diff:  
   ! Solver subroutines
   public :: solve_nnls
  +public :: solve_jacobi
   
   #ifdef USE_LAPACK
   public :: solve_dgels
   public :: solve_dgesv
  +public :: solve_dsyev
   public :: solve_dsysv
   #endif
   
  
positive.code: double precision :: alocal(nrows,ncols), blocal(nrows,1)
  
  ierr = 0
  
  ! Use local arrays for inversion
  alocal = Ain
  
negative.code: 
  #endif
  
  !
  ! !--------------------------------------------------------------------------------------------------!
  !
  ! subroutine jacobi(a,n,d,v,nrot)
  ! !----
  ! ! Numerical Recipes subroutine JACOBI:
  ! !
  ! ! Computes all eigenvalues and eigenvectors of a real symmetric matrix a, which is of size n by n,
  ! ! stored in a physical np by np array. On output, elements of a above the diagonal are destroyed. d
  ! ! returns the eigenvalues of a in its first n elements. v is a matrix with the same logical and
  ! ! physical dimensions as a, whose columns contain, on output, the normalized eigenvectors of a.
  ! ! nrot returns the number of Jacobi rotations that were required.
  ! !----
  ! implicit none
  !
  ! ! Arguments
  ! integer :: n, nrot
  ! double precision :: a(n,n), d(n), v(n,n)
  !
  ! ! Local variables
  ! integer :: i, ip, iq, j
  ! double precision :: c, g, h, s, sm, t, tau, theta, tresh, b(:), z(:)
  !
  !
  ! ! Initialize v as the identity matrix
  ! v = 0.0d0
  ! do ip = 1,n
  !     v(ip,ip) = 1.0d0
  ! enddo
  ! !
  ! !
  ! !
  ! !
  ! !
  !
  ! return
  ! end subroutine
  !
  !
  !--------------------------------------------------------------------------------------------------!
  !--------------------------------------------------------------------------------------------------!
  ! NON-NEGATIVE LEAST SQUARES SUBROUTINES - EDIT AT YOUR OWN RISK                                   !
  