function [L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr]  =  ludecp(A,INDEX,iorder,ierr)
%  SUBROUTINE COMPUTES L-U DECOMPOSITION ELEMENTS

%     ierr = 0   good result
%     ierr = -1  first diagonal element zero (failure! but
%                not necessarily singular)
%     ierr = -2  failure singular matrix

%A=G
ierr = 0;
L11 = A(INDEX,1,1);
if L11 == 0.0
  ierr = -1
  return;
end
V1 = 1./L11;
U12 =  V1.*A(INDEX,1,2);
U13 =  V1.*A(INDEX,1,3);
L21 = A(INDEX,2,1);
L22 = A(INDEX,2,2) - L21.*U12;
if L22 == 0.0
  ierr = -2
  return;
end
V2 = 1./L22;
U23 = (A(INDEX,2,3) -L21.*U13).* V2;
L31 = A(INDEX,3,1);
L32 = A(INDEX,3,2) - L31.*U12;
L33 = A(INDEX,3,3) - L31.*U13 -L32.*U23;
if L33 == 0.0
  ierr = -2;
  return;
end
V3 = 1./L33;

if iorder == 0
  L11 = L11'; L21 = L21'; L22 = L22'; L31 = L31'; L32 = L32'; L33 = L33';
  V1 = V1';   V2 = V2';   V3 = V3';
  U12 = U12'; U13 = U13'; U23 = U23';
end