anchor.diff:          subroutine decode_tri_ind(ind, i, j)
   
               ! Decode a combined index.
  +            ! The combined index, p, is given by:
  +            !   p = i(i-1)/2 + j
  +            ! Following the procedure by Rifkin, it is easy to show:
  +            !   i = int(0.5 + sqrt( 2p - 7/4 )
  +            ! and j hence follows.
  +            !
               ! In:
               !    ind: combined index based upon i,j.
               ! Out:
  
positive.code:             integer, intent(in) :: ind
              integer, intent(out) :: i, j
  
              i = int(0.50_dp + sqrt(2*ind-0.750_dp))
              j = ind - (i*(i-1))/2
  
          end subroutine decode_tri_ind
  
negative.code:         subroutine decode_rect_ind(ind, n, i, j)
  
              ! Decode a combined index.
              ! The combined index, p, is given by:
              !   p = p = n(n+1)/2 + (i-n-1)*n + j
              ! Following the procedure by Rifkin, it is easy to show:
              !   i = int[ 1/n (p - n(n+1)/2 - 1 )] + n + 1
              ! and j hence follows.
              !
              ! In:
              !    ind: combined index based upon i,j.
              !    n: max value of j.
  