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:         subroutine decode_rect_ind(ind, n, i, j)
  
              ! Decode a combined index.
              ! In:
              !    ind: combined index based upon i,j.
              !    n: max value of j.
  
negative.code:             integer, intent(in) :: ind, n
              integer, intent(out) :: i, j
  
              i = int( (ind - (n*(n+1))/2 - 1.0_dp)/n ) + n + 1
              j = ind - (i*(i-1))/2
  
          end subroutine decode_rect_ind
  