sparse-matrix routines are now moved to ../sparse-matrix/ subdirectories

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@1044 f88db872-0e4f-0410-b76b-b9085cfa78c5

3 files changed, 1 insertions, 1061 deletions

diff --git a/src/elliptic/README b/src/elliptic/README
index 09fa662..64f371d 100644
--- a/src/elliptic/README
+++ b/src/elliptic/README
@@ -31,11 +31,6 @@ lapack_wrapper.F77
 
 ilucg.h
 	Header file defining C/C++ prototypes for ILUCG routines
-ilucg.f77
-	ILUCG (Incomplete LU-decomposition / conjugate gradiant)
-	subroutine.  I don't know the original author of this routine;
-	it was provided to me by Tom Nicol, UBC Computing Center,
-	in c.1985.
 ilucg_wrapper.F77
 	Wrapper routines around the ILUCG routines, to avoid problems
 	with C <--> Fortran passing of Fortran logical return results.
diff --git a/src/elliptic/ilucg.f77 b/src/elliptic/ilucg.f77
deleted file mode 100644
index 68f99f5..0000000
--- a/src/elliptic/ilucg.f77
+++ /dev/null
@@ -1,1055 +0,0 @@
-      LOGICAL FUNCTION SILUCG(N,IA,JA,A,B,X,ITEMP,RTEMP,EPS,ITER,IE)
-      IMPLICIT REAL (A-H,O-Z)
-C
-C     INCOMPLETE LU DECOMPOSITION-CONJUGATE GRADIENT
-C     -          --               -         -
-C WHERE:
-C     |N| IS THE NUMBER OF EQUATIONS.  IF N < 0, ITEMP AND
-C       RTEMP CONTAIN THE ILU FROM A PREVIOUS CALL AND
-C       B AND X ARE THE NEW RHS AND INITIAL GUESS.
-C     IA IS AN INTEGER ARRAY DIMENSIONED |N|+1.  IA(I) IS THE
-C       INDEX INTO ARRAYS JA AND A OF THE FIRST NON-ZERO
-C       ELEMENT IN ROW I.  LET MAX=IA(|N|+1)-IA(1).
-C     JA IS AN INTEGER ARRAY DIMENSIONED MAX.  JA(K) GIVES
-C       THE COLUMN NUMBER OF A(K).
-C     A IS A REAL ARRAY DIMENSIONED MAX.  IT CONTAINS THE
-C       NONZERO ELEMENTS OF THE MATRIX STORED BY ROW.
-C     B CONTAINS THE RHS VECTOR.
-C     X IS A REAL ARRAY DIMENSIONED |N|.  ON ENTRY, IT CONTAINS
-C       AN INITIAL ESTIMATE; ON EXIT, THE SOLUTION.
-C     ITEMP IS AN INTEGER SCRATCH ARRAY DIMENSIONED 3*(|N|+MAX)+2.
-C     RTEMP IS AN REAL SCRATCH ARRAY DIMENSIONED 4*|N|+MAX.
-C     EPS IS THE CONVERGENCE CRITERIA.  IT SPECIFIES THE RELATIVE
-C       ERROR ALLOWED IN THE SOLUTION.  TO BE PRECISE, CONVERGENCE
-C       IS DEEMED TO HAVE OCCURED WHEN THE INFINITY-NORM OF THE
-C       CHANGE IN THE SOLUTION IN ONE ITERATION IS .LE. EPS * THE
-C       INFINITY-NORM OF THE CURRENT SOLUTION.  HOWEVER, IF EPS
-C       .LT. 0.0, IT IS INTERNALLY SCALED BY THE MACHINE PRECISION,
-C       SO THAT, FOR EXAMPLE, EPS = -256.0 WILL ALLOW THE LAST TWO
-C       HEXADECIMAL DIGITS OF THE SOLUTION TO BE IN ERROR.
-C     ITER GIVES THE REQUESTED NUMBER OF ITERATIONS, OR IS 0 FOR "NO LIMIT".
-C     IE IS AN INTEGER VARIABLE.  FOR NORMAL RETURN, IT GIVES
-C       THE NUMBER OF ITERATIONS DONE (NEGATIVE IS NO CONVERGENCE).
-C       ON ERROR RETURN, THE ROW IN ERROR.
-C
-C THE FUNCTION RETURNS .FALSE. IF ALL'S WELL, .TRUE. ON AN ERROR RETURN.
-C
-C (MODIFIED TO RETURN LOGICAL VALUE, J. THORNBURG, 13/MAY/85.)
-C (MODIFIED TO ADD CONVERGENCE CRITERIA, J. THORNBURG, 17/MAY/85.)
-C (MODIFIED CONVERGENCE CRITERIA, J. THORNBURG, 4/AUG/89.)
-C (MODIFIED TO AVOID SKIP RETURNS, J. THORNBURG, 10/SEP/89.)
-C
-C REFERENCE:
-C     D.S.KERSHAW,"THE INCOMPLETE CHOLESKY-CONJUGATE GRADIENT
-C       METHOD FOR INTERATIVE SOLUTION OF LINEAR EQUATIONS",
-C       J.COMPUT.PHYS. JAN 1978 PP 43-65
-C
-      DIMENSION IA(*),JA(*),A(*),B(*),X(*),ITEMP(*),RTEMP(*)
-      LOGICAL SLU0
-      NP=IABS(N)
-      IE=0
-      IF (NP.EQ.0) GO TO 20
-C CALCULATE INDICES FOR BREAKING UP TEMPORARY ARRAYS.
-      N1=NP+1
-      MAX=IA(N1)-IA(1)
-      ILU=1
-      JLU=ILU+N1
-      ID=JLU+MAX
-      IC=ID+NP
-      JC=IC+N1
-      JCI=JC+MAX
-      IR=1
-      IP=IR+NP
-      IS1=IP+NP
-      IS2=IS1+NP
-      IALU=IS2+NP
-      IF (N.LT.0) GO TO 10
-C DO INCOMPLETE LU DECOMPOSITION
-      IF (SLU0(NP,IA,JA,A,ITEMP(IC),ITEMP(JC),ITEMP(JCI),RTEMP(IALU),
-     *    ITEMP(ILU),ITEMP(JLU),ITEMP(ID),RTEMP(IR),IE)) GOTO 20
-C AND DO CONJUGATE GRADIENT ITERATIONS
-C CALL MODIFIED TO ADD ADJUSTABLE CONVERGENCE CRITERIA EPS
-C - J. THORNBURG, 17/MAY/85.
-10    CALL SNCG0(NP,IA,JA,A,B,X,ITEMP(ILU),ITEMP(JLU),ITEMP(ID),
-     *  RTEMP(IALU),RTEMP(IR),RTEMP(IP),RTEMP(IS1),RTEMP(IS2),
-     *  EPS,ITER,IE)
-      SILUCG = .FALSE.
-      RETURN
-20    SILUCG = .TRUE.
-      RETURN
-      END
-      LOGICAL FUNCTION SLU0(N,IA,JA,A,IC,JC,JCI,ALU,ILU,JLU,ID,V,IE)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),IC(*),JC(*),JCI(*),
-     *  ALU(*),ILU(*),JLU(*),ID(N),V(N)
-      LOGICAL NODIAG
-      COMMON /ICBD00/ ICBAD
-C     INCOMPLETE LU DECOMPOSITION
-C WHERE:
-C     N,IA,JA, AND A ARE DESCRIBED IN FUNCTION SILUCG
-C     IC IS AN INTEGER ARRAY DIMENSIONED N+1, IC(J) GIVES THE
-C       INDEX OF THE FIRST NONZERO ELEMENT IN COLMN J IN
-C       ARRAY JC.
-C     JC IS AN INTEGER ARRAY WITH THE SAME DIMENSION AS A.
-C       JC(K) GIVES THE ROW NUMBER OF THE K'TH ELEMENT IN
-C       THE COLUMN STRUCTURE.
-C     JCI IS AN INTEGER ARRAY WITH THE SAME DIMENSION AS A.
-C       JCI(K) GIVES THE INDEX INTO ARRAY A OF THE K'TH ELEMENT
-C       OF THE COLUMN STRUCTURE.
-C     ALU HAS THE SAME DIMENSION AS A.  ON EXIT, IT WILL
-C       CONTAIN THE INCOMPLETE LU DECOMPOSITION OF A WITH THE
-C       RECIPROCALS OF THE DIAGONAL ELEMENTS OF U.
-C     ILU AND JLU CORRESPONDS TO IA AND JA BUT FOR ALU.
-C     ID IS AN INTEGER ARRAY DIMENSIONED N.  IT CONTAINS
-C       INDICES TO THE DIAGONAL ELEMENTS OF U.
-C     V IS A REAL SCRATCH VECTOR OF LENGTH N.
-C     IE GIVES THE ROW NUMBER IN ERROR.
-C
-C     RETURN VALUE = .FALSE. IF ALL IS WELL, .TRUE. IF ERROR.
-C
-C NOTE: SLU0 SETS ARGUMENTS IC THROUGH V.
-C
-      ICBAD=0
-C ZERO COUNT OF ZERO DIAGONAL ELEMENTS IN U.
-C
-C FIRST CHECK STRUCTURE OF A AND BUILD COLUMN STRUCTURE
-      DO 10 I=1,N
-        IC(I)=0
-10    CONTINUE
-      DO 30 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        NODIAG=.TRUE.
-        DO 20 K=KS,KE
-          J=JA(K)
-          IF (J.LT.1.OR.J.GT.N) GO TO 210
-          IC(J)=IC(J)+1
-          IF (J.EQ.I) NODIAG=.FALSE.
-20      CONTINUE
-        IF (NODIAG) GO TO 210
-30    CONTINUE
-C MAKE IC INTO INDICES
-      KOLD=IC(1)
-      IC(1)=1
-      DO 40 I=1,N
-        KNEW=IC(I+1)
-        IF (KOLD.EQ.0) GO TO 210
-        IC(I+1)=IC(I)+KOLD
-        KOLD=KNEW
-40    CONTINUE
-C SET JC AND JCI FOR COLUMN STRUCTURE
-      DO 60 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        DO 50 K=KS,KE
-          J=JA(K)
-          L=IC(J)
-          IC(J)=L+1
-          JC(L)=I
-          JCI(L)=K
-50      CONTINUE
-60    CONTINUE
-C FIX UP IC
-      KOLD=IC(1)
-      IC(1)=1
-      DO 70 I=1,N
-        KNEW=IC(I+1)
-        IC(I+1)=KOLD
-        KOLD=KNEW
-70    CONTINUE
-C FIND SORTED ROW STRUCTURE FROM SORTED COLUMN STRUCTURE
-      NP=N+1
-      DO 80 I=1,NP
-        ILU(I)=IA(I)
-80    CONTINUE
-C MOVE ELEMENTS, SET JLU AND ID
-      DO 100 J=1,N
-        KS=IC(J)
-        KE=IC(J+1)-1
-        DO 90 K=KS,KE
-          I=JC(K)
-          L=ILU(I)
-          ILU(I)=L+1
-          JLU(L)=J
-          KK=JCI(K)
-          ALU(L)=A(KK)
-          IF (I.EQ.J) ID(J)=L
-90      CONTINUE
-100   CONTINUE
-C RESET ILU (COULD JUST USE IA)
-      DO 110 I=1,NP
-        ILU(I)=IA(I)
-110   CONTINUE
-C FINISHED WITH SORTED COLUMN AND ROW STRUCTURE
-C
-C DO LU DECOMPOSITION USING GAUSSIAN ELIMINATION
-      DO 120 I=1,N
-        V(I)=0.0
-120   CONTINUE
-      DO 200 IROW=1,N
-        I=ID(IROW)
-        PIVOT=ALU(I)
-        IF (PIVOT.NE.0.0) GO TO 140
-C THIS CASE MAKES THE ILU LESS ACCURATE
-        ICBAD=ICBAD+1
-        KS=ILU(IROW)
-        KE=ILU(IROW+1)-1
-        DO 130 K=KS,KE
-          PIVOT=PIVOT+ABS(ALU(K))
-130     CONTINUE
-        IF (PIVOT.EQ.0.0) GO TO 220
-140     PIVOT=1.0/PIVOT
-        ALU(I)=PIVOT
-        KKS=I+1
-        KKE=ILU(IROW+1)-1
-        IF (KKS.GT.KKE) GO TO 160
-        DO 150 K=KKS,KKE
-          J=JLU(K)
-          V(J)=ALU(K)
-150     CONTINUE
-C FIX L IN COLUMN IROW AND DO PARTIAL LU IN SUBMATRIX
-160     KS=IC(IROW)
-        KE=IC(IROW+1)-1
-        DO 190 K=KS,KE
-          I=JC(K)
-          IF (I.LE.IROW) GO TO 190
-          LS=ILU(I)
-          LE=ILU(I+1)-1
-          DO 180 L=LS,LE
-            J=JLU(L)
-            IF (J.LT.IROW) GO TO 180
-            IF (J.GT.IROW) GO TO 170
-            AMULT=ALU(L)*PIVOT
-            ALU(L)=AMULT
-            IF (AMULT.EQ.0.0) GO TO 190
-            GO TO 180
-170         IF (V(J).EQ.0.0) GO TO 180
-            ALU(L)=ALU(L)-AMULT*V(J)
-180       CONTINUE
-190     CONTINUE
-C RESET V
-        IF (KKS.GT.KKE) GO TO 200
-        DO 195 K=KKS,KKE
-          J=JLU(K)
-          V(J)=0.0
-195     CONTINUE
-200   CONTINUE
-C NORMAL RETURN
-      SLU0 = .FALSE.
-      RETURN
-C ERROR RETURNS
-210   IE=I
-      SLU0 = .TRUE.
-      RETURN
-C NORMAL RETURN
-220   IE=IROW
-      SLU0 = .TRUE.
-      RETURN
-      END
-      SUBROUTINE SNCG0(N,IA,JA,A,B,X,ILU,JLU,ID,ALU,R,P,S1,S2,
-     *  EPS,ITER,IE)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),B(N),X(N),ILU(*),JLU(*),ALU(*),ID(N),
-     *  R(N),P(N),S1(N),S2(N)
-C     NONSYMMETRIC CONJUGATE GRADIENT
-C WHERE:
-C     N,IA,JA,A,B, AND X ARE DESCRIBED IN SUBROUTINE SILUCG.
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW OF LU.
-C     JLU GIVES COLUMN NUMBER.
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U.
-C     ALU HAS NONZERO ELEMENTS OF LU MATRIX STORED BY ROW
-C       WITH RECIPROCALS OF DIAGONAL ELEMENTS OF U.
-C     R,P,S1, AND S2 ARE VECTORS OF LENGTH N USED IN THE
-C       ITERATIONS.
-C FOLLOWING PARAMETER ADDED BY J. THORNBURG, 17/MAY/85.
-C     EPS IS CONVERGENCE CRITERIA.  (DESCRIBED IN SUBROUTINE
-C       SILUCG).
-C     ITER IS MAX NUMBER OF ITERATIONS, OR 0 FOR "NO LIMIT".
-C     IE GIVES ACTUAL NUMBER OF ITERATIONS, NEGATIVE IF
-C       NO CONVERGENCE.
-C
-C R0=B-A*X0
-      CALL SMUL10(N,IA,JA,A,X,R)
-      DO 10 I=1,N
-        R(I)=B(I)-R(I)
-10    CONTINUE
-C P0=(UT*U)(-1)*AT*(L*LT)(-1)*R0
-C FIRST SOLVE L*LT*S1=R0
-      CALL SSUBL0(N,ILU,JLU,ID,ALU,R,S1)
-C TIMES TRANSPOSE OF A
-      CALL SMUL20(N,IA,JA,A,S1,S2)
-C THEN SOLVE UT*U*P=S2
-      CALL SSUBU0(N,ILU,JLU,ID,ALU,S2,P)
-      IE=0
-      RDOT = SGVV(R,S1,N)
-C LOOP BEGINS HERE
-20    CALL SMUL30(N,ILU,JLU,ID,ALU,P,S2)
-      PDOT = SGVV(P,S2,N)
-C
-      IF (PDOT.EQ.0.0) RETURN
-C
-      ALPHA=RDOT/PDOT
-C EQUATION 9PA                      ALPHA=(R,LINV*R)/(P,UT*U*P)
-      XMAX=1.0
-      XDIF=0.0
-      DO 30 I=1,N
-        AP=ALPHA*P(I)
-        X(I)=X(I)+AP
-C EQUATION 9PB                      X=X+ALPHA*P
-        AP=ABS(AP)
-        XX=ABS(X(I))
-        IF (AP.GT.XDIF) XDIF=AP
-        IF (XX.GT.XMAX) XMAX=XX
-30    CONTINUE
-      IE=IE+1
-C
-C CONVERGENCE TEST (CHANGED BY J. THORNBURG, 17/MAY/85, 4/AUG/89)
-C
-      IF ((EPS .GT. 0.0) .AND. (XDIF .LE. EPS * XMAX)) RETURN
-      IF ((EPS .LT. 0.0) .AND. (XMAX + XDIF/ABS(EPS) .EQ. XMAX))
-     *   RETURN
-C
-C EXCEEDED ITERATION LIMIT?
-C
-      IF ((ITER .NE. 0) .AND. (IE .GE. ITER)) GO TO 60
-      CALL SMUL10(N,IA,JA,A,P,S2)
-      DO 40 I=1,N
-        R(I)=R(I)-ALPHA*S2(I)
-C EQUATION 9PC                      R=R-ALPHA*A*P
-40    CONTINUE
-      CALL SSUBL0(N,ILU,JLU,ID,ALU,R,S1)
-C     CALL INPROD(R,S1,EDOT,RRDOT,N)
-      RRDOT = SGVV(R,S1,N)
-      BETA=RRDOT/RDOT
-C EQUATION 9PD                      BETA=(R+,LINV*R+)/(R,LINV*R)
-      RDOT=RRDOT
-      CALL SMUL20(N,IA,JA,A,S1,S2)
-      CALL SSUBU0(N,ILU,JLU,ID,ALU,S2,S1)
-      DO 50 I=1,N
-        P(I)=S1(I)+BETA*P(I)
-C EQUATION 9PE                      P=(UT*U)(-1)*AT*(L*LT)(-1)*R+BETA*P
-50    CONTINUE
-      GO TO 20
-60    IE=-IE
-      RETURN
-      END
-      SUBROUTINE SMUL10(N,IA,JA,A,B,X)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),B(N),X(N)
-C MULTIPLY A TIMES B TO GET X
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     IA GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW
-C     JA GIVES COLUMN NUMBER
-C     A CONTAINS THE NONZERO ELEMENTS OF THE NONSYMMETRIC
-C       MATRIX STORED BY ROW
-C     B IS THE VECTOR
-C     X IS THE PRODUCT (MUST BE DIFFERENT FROM B)
-C
-      DO 20 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        SUM=0.0
-        DO 10 K=KS,KE
-          J=JA(K)
-          SUM=SUM+A(K)*B(J)
-10      CONTINUE
-        X(I)=SUM
-20    CONTINUE
-      RETURN
-      END
-      SUBROUTINE SMUL20(N,IA,JA,A,B,X)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),B(N),X(N)
-C MULTIPLY TRANSPOSE OF A TIMES B TO GET X
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     IA GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW
-C     JA GIVES COLUMN NUMBER
-C     A CONTAINS THE NONZERO ELEMENTS OF THE NONSYMMETRIC
-C       MATRIX STORED BY ROW
-C     B IS THE VECTOR
-C     X IS THE PRODUCT (MUST BE DIFFERENT FROM B)
-C
-      DO 10 I=1,N
-        X(I)=0.0
-10    CONTINUE
-      DO 30 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        BB=B(I)
-        DO 20 K=KS,KE
-          J=JA(K)
-          X(J)=X(J)+A(K)*BB
-20      CONTINUE
-30    CONTINUE
-      RETURN
-      END
-      SUBROUTINE SMUL30(N,ILU,JLU,ID,ALU,B,X)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION ILU(*),JLU(*),ID(*),ALU(*),B(N),X(N)
-C MULTIPLY TRANSPOSE OF U TIMES U TIMES B TO GET X
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW OF LU
-C     JLU GIVES COLUMN NUMBER
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U
-C     ALU HAS NONZERO ELEMENTS OF LU MATRIX STORED BY ROW
-C       WITH RECIPROCALS OF DIAGONAL ELEMENTS
-C     B IS THE VECTOR
-C     X IS THE PRODUCT UT*U*B (X MUST BE DIFFERENT FROM B)
-C
-      DO 10 I=1,N
-        X(I)=0.0
-10    CONTINUE
-      DO 50 I=1,N
-        KS=ID(I)+1
-        KE=ILU(I+1)-1
-        DIAG=1.0/ALU(KS-1)
-        XX=DIAG*B(I)
-        IF (KS.GT.KE) GO TO 30
-        DO 20 K=KS,KE
-          J=JLU(K)
-          XX=XX+ALU(K)*B(J)
-20      CONTINUE
-30      X(I)=X(I)+DIAG*XX
-        IF (KS.GT.KE) GO TO 50
-        DO 40 K=KS,KE
-          J=JLU(K)
-          X(J)=X(J)+ALU(K)*XX
-40      CONTINUE
-50    CONTINUE
-      RETURN
-      END
-      SUBROUTINE SSUBU0(N,ILU,JLU,ID,ALU,B,X)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION ILU(*),JLU(*),ID(*),ALU(*),B(N),X(N)
-C DO FORWARD AND BACK SUBSTITUTION TO SOLVE UT*U*X=B
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW OF LU
-C     JLU GIVES COLUMN NUMBER
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U
-C     ALU HAS NONZERO ELMENTS OF LU MATRIX STORED BY ROW
-C       WITH RECIPROCALS OF DIAGONAL ELEMENTS OF U
-C     B IS THE RHS VECTOR
-C     X IS THE SOLUTION VECTOR
-C
-      NP=N+1
-      DO 10 I=1,N
-        X(I)=B(I)
-10    CONTINUE
-C FORWARD SUBSTITUTION
-      DO 30 I=1,N
-        KS=ID(I)+1
-        KE=ILU(I+1)-1
-        XX=X(I)*ALU(KS-1)
-        X(I)=XX
-        IF (KS.GT.KE) GO TO 30
-        DO 20 K=KS,KE
-          J=JLU(K)
-          X(J)=X(J)-ALU(K)*XX
-20      CONTINUE
-30    CONTINUE
-C BACK SUBSTITUTION
-      DO 60 II=1,N
-        I=NP-II
-        KS=ID(I)+1
-        KE=ILU(I+1)-1
-        SUM=0.0
-        IF (KS.GT.KE) GO TO 50
-        DO 40 K=KS,KE
-          J=JLU(K)
-          SUM=SUM+ALU(K)*X(J)
-40      CONTINUE
-50      X(I)=(X(I)-SUM)*ALU(KS-1)
-60    CONTINUE
-      RETURN
-      END
-      SUBROUTINE SSUBL0(N,ILU,JLU,ID,ALU,B,X)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION ILU(*),JLU(*),ID(*),ALU(*),B(N),X(N)
-C DO FORWARD AND BACK SUBSTITUTION TO SOLVE L*LT*X=B
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW LU
-C     JLU GIVES THE COLUMN NUMBER
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U
-C     ALU HAS NONZERO ELEMENTS OF LU MATRIX STORED BY ROW
-C       DIAGONAL ELEMENTS OF L ARE 1.0 AND NOT STORED
-C     B IS THE RHS VECTOR
-C     X IS THE SOLUTION VECTOR
-C
-      NP=N+1
-      DO 10 I=1,N
-        X(I)=B(I)
-10    CONTINUE
-C FORWARD SUBSTITUTION
-      DO 30 I=1,N
-        KS=ILU(I)
-        KE=ID(I)-1
-        IF (KS.GT.KE) GO TO 30
-        SUM=0.0
-        DO 20 K=KS,KE
-          J=JLU(K)
-          SUM=SUM+ALU(K)*X(J)
-20      CONTINUE
-        X(I)=X(I)-SUM
-30    CONTINUE
-C BACK SUBSTITUTION
-      DO 50 II=1,N
-        I=NP-II
-        KS=ILU(I)
-        KE=ID(I)-1
-        IF (KS.GT.KE) GO TO 50
-        XX=X(I)
-        IF (XX.EQ.0.0) GO TO 50
-        DO 40 K=KS,KE
-          J=JLU(K)
-          X(J)=X(J)-ALU(K)*XX
-40      CONTINUE
-50    CONTINUE
-      RETURN
-      END
-      REAL FUNCTION SGVV(V,W,N)
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION V(N),W(N)
-C     SUBROUTINE TO COMPUTE REAL VECTOR DOT PRODUCT.
-C
-      SUM = 0.0
-            DO 10 I = 1,N
-            SUM = SUM + V(I)*W(I)
-10          CONTINUE
-      SGVV = SUM
-      RETURN
-      END
-      LOGICAL FUNCTION DILUCG(N,IA,JA,A,B,X,ITEMP,RTEMP,EPS,ITER,IE)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-C
-C     INCOMPLETE LU DECOMPOSITION-CONJUGATE GRADIENT
-C     -          --               -         -
-C WHERE:
-C     |N| IS THE NUMBER OF EQUATIONS.  IF N < 0, ITEMP AND
-C       RTEMP CONTAIN THE ILU FROM A PREVIOUS CALL AND
-C       B AND X ARE THE NEW RHS AND INITIAL GUESS.
-C     IA IS AN INTEGER ARRAY DIMENSIONED |N|+1.  IA(I) IS THE
-C       INDEX INTO ARRAYS JA AND A OF THE FIRST NON-ZERO
-C       ELEMENT IN ROW I.  LET MAX=IA(|N|+1)-IA(1).
-C     JA IS AN INTEGER ARRAY DIMENSIONED MAX.  JA(K) GIVES
-C       THE COLUMN NUMBER OF A(K).
-C     A IS A DOUBLE PRECISION ARRAY DIMENSIONED MAX.  IT CONTAINS THE
-C       NONZERO ELEMENTS OF THE MATRIX STORED BY ROW.
-C     B CONTAINS THE RHS VECTOR.
-C     X IS A DOUBLE PRECISION ARRAY DIMENSIONED |N|.  ON ENTRY, IT CONTAINS
-C       AN INITIAL ESTIMATE; ON EXIT, THE SOLUTION.
-C     ITEMP IS AN INTEGER SCRATCH ARRAY DIMENSIONED 3*(|N|+MAX)+2.
-C     RTEMP IS A DOUBLE PRECISION SCRATCH ARRAY DIMENSIONED 4*|N|+MAX.
-C     EPS IS THE CONVERGENCE CRITERIA.  IT SPECIFIES THE RELATIVE
-C       ERROR ALLOWED IN THE SOLUTION.  TO BE PRECISE, CONVERGENCE
-C       IS DEEMED TO HAVE OCCURED WHEN THE INFINITY-NORM OF THE
-C       CHANGE IN THE SOLUTION IN ONE ITERATION IS .LE. EPS * THE
-C       INFINITY-NORM OF THE CURRENT SOLUTION.  HOWEVER, IF EPS
-C       .LT. 0.0D0, IT IS INTERNALLY SCALED BY THE MACHINE PRECISION,
-C       SO THAT, FOR EXAMPLE, EPS = -256.0 WILL ALLOW THE LAST TWO
-C       HEXADECIMAL DIGITS OF THE SOLUTION TO BE IN ERROR.
-C     ITER GIVES THE REQUESTED NUMBER OF ITERATIONS, OR IS 0 FOR "NO LIMIT".
-C     IE IS AN INTEGER VARIABLE.  FOR NORMAL RETURN, IT GIVES
-C       THE NUMBER OF ITERATIONS DONE (NEGATIVE IS NO CONVERGENCE).
-C       ON ERROR RETURN, THE ROW IN ERROR.
-C
-C THE FUNCTION RETURNS .FALSE. IF ALL'S WELL, .TRUE. ON AN ERROR RETURN.
-C
-C (MODIFIED TO RETURN LOGICAL VALUE, J. THORNBURG, 13/MAY/85.)
-C (MODIFIED TO ADD CONVERGENCE CRITERIA, J. THORNBURG, 17/MAY/85.)
-C (MODIFIED CONVERGENCE CRITERIA, J. THORNBURG, 4/AUG/89.)
-C (MODIFIED TO AVOID SKIP RETURNS, J. THORNBURG, 10/SEP/89.)
-C
-C REFERENCE:
-C     D.S.KERSHAW,"THE INCOMPLETE CHOLESKY-CONJUGATE GRADIENT
-C       METHOD FOR INTERATIVE SOLUTION OF LINEAR EQUATIONS",
-C       J.COMPUT.PHYS. JAN 1978 PP 43-65
-C
-      DIMENSION IA(*),JA(*),A(*),B(*),X(*),ITEMP(*),RTEMP(*)
-      LOGICAL DLU0
-      NP=IABS(N)
-      IE=0
-      IF (NP.EQ.0) GO TO 20
-C CALCULATE INDICES FOR BREAKING UP TEMPORARY ARRAYS.
-      N1=NP+1
-      MAX=IA(N1)-IA(1)
-      ILU=1
-      JLU=ILU+N1
-      ID=JLU+MAX
-      IC=ID+NP
-      JC=IC+N1
-      JCI=JC+MAX
-      IR=1
-      IP=IR+NP
-      IS1=IP+NP
-      IS2=IS1+NP
-      IALU=IS2+NP
-      IF (N.LT.0) GO TO 10
-C DO INCOMPLETE LU DECOMPOSITION
-      IF (DLU0(NP,IA,JA,A,ITEMP(IC),ITEMP(JC),ITEMP(JCI),RTEMP(IALU),
-     *    ITEMP(ILU),ITEMP(JLU),ITEMP(ID),RTEMP(IR),IE)) GOTO 20
-C AND DO CONJUGATE GRADIENT ITERATIONS
-C CALL MODIFIED TO ADD ADJUSTABLE CONVERGENCE CRITERIA EPS
-C - J. THORNBURG, 17/MAY/85.
-10    CALL DNCG0(NP,IA,JA,A,B,X,ITEMP(ILU),ITEMP(JLU),ITEMP(ID),
-     *  RTEMP(IALU),RTEMP(IR),RTEMP(IP),RTEMP(IS1),RTEMP(IS2),
-     *  EPS,ITER,IE)
-      DILUCG = .FALSE.
-      RETURN
-20    DILUCG = .TRUE.
-      RETURN
-      END
-      LOGICAL FUNCTION DLU0(N,IA,JA,A,IC,JC,JCI,ALU,ILU,JLU,ID,V,IE)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),IC(*),JC(*),JCI(*),
-     *  ALU(*),ILU(*),JLU(*),ID(N),V(N)
-      LOGICAL NODIAG
-      COMMON /ICBD00/ ICBAD
-C     INCOMPLETE LU DECOMPOSITION
-C WHERE:
-C     N,IA,JA, AND A ARE DESCRIBED IN SUBROUTINE ILUCG
-C     IC IS AN INTEGER ARRAY DIMENSIONED N+1, IC(J) GIVES THE
-C       INDEX OF THE FIRST NONZERO ELEMENT IN COLMN J IN
-C       ARRAY JC.
-C     JC IS AN INTEGER ARRAY WITH THE SAME DIMENSION AS A.
-C       JC(K) GIVES THE ROW NUMBER OF THE K'TH ELEMENT IN
-C       THE COLUMN STRUCTURE.
-C     JCI IS AN INTEGER ARRAY WITH THE SAME DIMENSION AS A.
-C       JCI(K) GIVES THE INDEX INTO ARRAY A OF THE K'TH ELEMENT
-C       OF THE COLUMN STRUCTURE.
-C     ALU HAS THE SAME DIMENSION AS A.  ON EXIT, IT WILL
-C       CONTAIN THE INCOMPLETE LU DECOMPOSITION OF A WITH THE
-C       RECIPROCALS OF THE DIAGONAL ELEMENTS OF U.
-C     ILU AND JLU CORRESPONDS TO IA AND JA BUT FOR ALU.
-C     ID IS AN INTEGER ARRAY DIMENSIONED N.  IT CONTAINS
-C       INDICES TO THE DIAGONAL ELEMENTS OF U.
-C     V IS A REAL SCRATCH VECTOR OF LENGTH N.
-C     IE GIVES THE ROW NUMBER IN ERROR.
-C
-C     RETURN VALUE = .FALSE. IF ALL IS WELL, .TRUE. IF ERROR.
-C
-C NOTE: DLU0 SETS ARGUMENTS IC THROUGH V.
-C
-      ICBAD=0
-C ZERO COUNT OF ZERO DIAGONAL ELEMENTS IN U.
-C
-C FIRST CHECK STRUCTURE OF A AND BUILD COLUMN STRUCTURE
-      DO 10 I=1,N
-        IC(I)=0
-10    CONTINUE
-      DO 30 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        NODIAG=.TRUE.
-        DO 20 K=KS,KE
-          J=JA(K)
-          IF (J.LT.1.OR.J.GT.N) GO TO 210
-          IC(J)=IC(J)+1
-          IF (J.EQ.I) NODIAG=.FALSE.
-20      CONTINUE
-        IF (NODIAG) GO TO 210
-30    CONTINUE
-C MAKE IC INTO INDICES
-      KOLD=IC(1)
-      IC(1)=1
-      DO 40 I=1,N
-        KNEW=IC(I+1)
-        IF (KOLD.EQ.0) GO TO 210
-        IC(I+1)=IC(I)+KOLD
-        KOLD=KNEW
-40    CONTINUE
-C SET JC AND JCI FOR COLUMN STRUCTURE
-      DO 60 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        DO 50 K=KS,KE
-          J=JA(K)
-          L=IC(J)
-          IC(J)=L+1
-          JC(L)=I
-          JCI(L)=K
-50      CONTINUE
-60    CONTINUE
-C FIX UP IC
-      KOLD=IC(1)
-      IC(1)=1
-      DO 70 I=1,N
-        KNEW=IC(I+1)
-        IC(I+1)=KOLD
-        KOLD=KNEW
-70    CONTINUE
-C FIND SORTED ROW STRUCTURE FROM SORTED COLUMN STRUCTURE
-      NP=N+1
-      DO 80 I=1,NP
-        ILU(I)=IA(I)
-80    CONTINUE
-C MOVE ELEMENTS, SET JLU AND ID
-      DO 100 J=1,N
-        KS=IC(J)
-        KE=IC(J+1)-1
-        DO 90 K=KS,KE
-          I=JC(K)
-          L=ILU(I)
-          ILU(I)=L+1
-          JLU(L)=J
-          KK=JCI(K)
-          ALU(L)=A(KK)
-          IF (I.EQ.J) ID(J)=L
-90      CONTINUE
-100   CONTINUE
-C RESET ILU (COULD JUST USE IA)
-      DO 110 I=1,NP
-        ILU(I)=IA(I)
-110   CONTINUE
-C FINISHED WITH SORTED COLUMN AND ROW STRUCTURE
-C
-C DO LU DECOMPOSITION USING GAUSSIAN ELIMINATION
-      DO 120 I=1,N
-        V(I)=0.0D0
-120   CONTINUE
-      DO 200 IROW=1,N
-        I=ID(IROW)
-        PIVOT=ALU(I)
-        IF (PIVOT.NE.0.0D0) GO TO 140
-C THIS CASE MAKES THE ILU LESS ACCURATE
-        ICBAD=ICBAD+1
-        KS=ILU(IROW)
-        KE=ILU(IROW+1)-1
-        DO 130 K=KS,KE
-          PIVOT=PIVOT+DABS(ALU(K))
-130     CONTINUE
-        IF (PIVOT.EQ.0.0D0) GO TO 220
-140     PIVOT=1.0D0/PIVOT
-        ALU(I)=PIVOT
-        KKS=I+1
-        KKE=ILU(IROW+1)-1
-        IF (KKS.GT.KKE) GO TO 160
-        DO 150 K=KKS,KKE
-          J=JLU(K)
-          V(J)=ALU(K)
-150     CONTINUE
-C FIX L IN COLUMN IROW AND DO PARTIAL LU IN SUBMATRIX
-160     KS=IC(IROW)
-        KE=IC(IROW+1)-1
-        DO 190 K=KS,KE
-          I=JC(K)
-          IF (I.LE.IROW) GO TO 190
-          LS=ILU(I)
-          LE=ILU(I+1)-1
-          DO 180 L=LS,LE
-            J=JLU(L)
-            IF (J.LT.IROW) GO TO 180
-            IF (J.GT.IROW) GO TO 170
-            AMULT=ALU(L)*PIVOT
-            ALU(L)=AMULT
-            IF (AMULT.EQ.0.0) GO TO 190
-            GO TO 180
-170         IF (V(J).EQ.0.0D0) GO TO 180
-            ALU(L)=ALU(L)-AMULT*V(J)
-180       CONTINUE
-190     CONTINUE
-C RESET V
-        IF (KKS.GT.KKE) GO TO 200
-        DO 195 K=KKS,KKE
-          J=JLU(K)
-          V(J)=0.0D0
-195     CONTINUE
-200   CONTINUE
-C NORMAL RETURN
-      DLU0 = .FALSE.
-      RETURN
-C ERROR RETURNS
-210   IE=I
-      DLU0 = .TRUE.
-      RETURN
-220   IE=IROW
-      DLU0 = .TRUE.
-      RETURN
-      END
-      SUBROUTINE DNCG0(N,IA,JA,A,B,X,ILU,JLU,ID,ALU,R,P,S1,S2,
-     *  EPS,ITER,IE)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),B(N),X(N),ILU(*),JLU(*),ALU(*),ID(N),
-     *  R(N),P(N),S1(N),S2(N)
-C     NONSYMMETRIC CONJUGATE GRADIENT
-C WHERE:
-C     N,IA,JA,A,B, AND X ARE DESCRIBED IN SUBROUTINE DILUCG.
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW OF LU.
-C     JLU GIVES COLUMN NUMBER.
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U.
-C     ALU HAS NONZERO ELEMENTS OF LU MATRIX STORED BY ROW
-C       WITH RECIPROCALS OF DIAGONAL ELEMENTS OF U.
-C     R,P,S1, AND S2 ARE VECTORS OF LENGTH N USED IN THE
-C       ITERATIONS.
-C FOLLOWING PARAMETER ADDED BY J. THORNBURG, 17/MAY/85.
-C     EPS IS CONVERGENCE CRITERIA.  (DESCRIBED IN SUBROUTINE
-C       DILUCG).
-C     ITER IS MAX NUMBER OF ITERATIONS, OR 0 FOR "NO LIMIT".
-C     IE GIVES ACTUAL NUMBER OF ITERATIONS, NEGATIVE IF
-C       NO CONVERGENCE.
-C
-C R0=B-A*X0
-      CALL DMUL10(N,IA,JA,A,X,R)
-      DO 10 I=1,N
-        R(I)=B(I)-R(I)
-10    CONTINUE
-C P0=(UT*U)(-1)*AT*(L*LT)(-1)*R0
-C FIRST SOLVE L*LT*S1=R0
-      CALL DSUBL0(N,ILU,JLU,ID,ALU,R,S1)
-C TIMES TRANSPOSE OF A
-      CALL DMUL20(N,IA,JA,A,S1,S2)
-C THEN SOLVE UT*U*P=S2
-      CALL DSUBU0(N,ILU,JLU,ID,ALU,S2,P)
-      IE=0
-C INPROD IS DOT PRODUCT ROUTINE IN *LIBRARY (UBC MATRIX P 28)
-C ALL CALLS ON IT COMMENTED OUT AND REPLACED WITH CALLS TO NEW
-C ROUTINE DGVV(...); SOURCE CODE FOR LATTER ADDED TO END OF
-C THIS FILE.  - J. THORNBURG, 10/MAY/85.
-C     CALL INPROD(R,S1,EDOT,RDOT,N)
-      RDOT = DGVV(R,S1,N)
-C LOOP BEGINS HERE
-20    CALL DMUL30(N,ILU,JLU,ID,ALU,P,S2)
-C     CALL INPROD(P,S2,EDOT,PDOT,N)
-      PDOT = DGVV(P,S2,N)
-C
-      IF (PDOT.EQ.0.0D0) RETURN
-C
-      ALPHA=RDOT/PDOT
-C EQUATION 9PA                      ALPHA=(R,LINV*R)/(P,UT*U*P)
-      XMAX=1.0D0
-      XDIF=0.0D0
-      DO 30 I=1,N
-        AP=ALPHA*P(I)
-        X(I)=X(I)+AP
-C EQUATION 9PB                      X=X+ALPHA*P
-        AP=DABS(AP)
-        XX=DABS(X(I))
-        IF (AP.GT.XDIF) XDIF=AP
-        IF (XX.GT.XMAX) XMAX=XX
-30    CONTINUE
-      IE=IE+1
-C
-C CONVERGENCE TEST (CHANGED BY J. THORNBURG, 17/MAY/85, 4/AUG/89)
-C
-      IF ((EPS .GT. 0.0D0) .AND. (XDIF .LE. EPS * XMAX)) RETURN
-      IF ((EPS .LT. 0.0D0) .AND. (XMAX + XDIF/DABS(EPS) .EQ. XMAX))
-     *   RETURN
-C
-C EXCEEDED ITERATION LIMIT?
-C
-      IF ((ITER .NE. 0) .AND. (IE .GE. ITER)) GO TO 60
-      CALL DMUL10(N,IA,JA,A,P,S2)
-      DO 40 I=1,N
-        R(I)=R(I)-ALPHA*S2(I)
-C EQUATION 9PC                      R=R-ALPHA*A*P
-40    CONTINUE
-      CALL DSUBL0(N,ILU,JLU,ID,ALU,R,S1)
-C     CALL INPROD(R,S1,EDOT,RRDOT,N)
-      RRDOT = DGVV(R,S1,N)
-      BETA=RRDOT/RDOT
-C EQUATION 9PD                      BETA=(R+,LINV*R+)/(R,LINV*R)
-      RDOT=RRDOT
-      CALL DMUL20(N,IA,JA,A,S1,S2)
-      CALL DSUBU0(N,ILU,JLU,ID,ALU,S2,S1)
-      DO 50 I=1,N
-        P(I)=S1(I)+BETA*P(I)
-C EQUATION 9PE                      P=(UT*U)(-1)*AT*(L*LT)(-1)*R+BETA*P
-50    CONTINUE
-      GO TO 20
-60    IE=-IE
-      RETURN
-      END
-      SUBROUTINE DMUL10(N,IA,JA,A,B,X)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),B(N),X(N)
-C MULTIPLY A TIMES B TO GET X
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     IA GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW
-C     JA GIVES COLUMN NUMBER
-C     A CONTAINS THE NONZERO ELEMENTS OF THE NONSYMMETRIC
-C       MATRIX STORED BY ROW
-C     B IS THE VECTOR
-C     X IS THE PRODUCT (MUST BE DIFFERENT FROM B)
-C
-      DO 20 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        SUM=0.0D0
-        DO 10 K=KS,KE
-          J=JA(K)
-          SUM=SUM+A(K)*B(J)
-10      CONTINUE
-        X(I)=SUM
-20    CONTINUE
-      RETURN
-      END
-      SUBROUTINE DMUL20(N,IA,JA,A,B,X)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION IA(*),JA(*),A(*),B(N),X(N)
-C MULTIPLY TRANSPOSE OF A TIMES B TO GET X
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     IA GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW
-C     JA GIVES COLUMN NUMBER
-C     A CONTAINS THE NONZERO ELEMENTS OF THE NONSYMMETRIC
-C       MATRIX STORED BY ROW
-C     B IS THE VECTOR
-C     X IS THE PRODUCT (MUST BE DIFFERENT FROM B)
-C
-      DO 10 I=1,N
-        X(I)=0.0D0
-10    CONTINUE
-      DO 30 I=1,N
-        KS=IA(I)
-        KE=IA(I+1)-1
-        BB=B(I)
-        DO 20 K=KS,KE
-          J=JA(K)
-          X(J)=X(J)+A(K)*BB
-20      CONTINUE
-30    CONTINUE
-      RETURN
-      END
-      SUBROUTINE DMUL30(N,ILU,JLU,ID,ALU,B,X)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION ILU(*),JLU(*),ID(*),ALU(*),B(N),X(N)
-C MULTIPLY TRANSPOSE OF U TIMES U TIMES B TO GET X
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW OF LU
-C     JLU GIVES COLUMN NUMBER
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U
-C     ALU HAS NONZERO ELEMENTS OF LU MATRIX STORED BY ROW
-C       WITH RECIPROCALS OF DIAGONAL ELEMENTS
-C     B IS THE VECTOR
-C     X IS THE PRODUCT UT*U*B (X MUST BE DIFFERENT FROM B)
-C
-      DO 10 I=1,N
-        X(I)=0.0D0
-10    CONTINUE
-      DO 50 I=1,N
-        KS=ID(I)+1
-        KE=ILU(I+1)-1
-        DIAG=1.0D0/ALU(KS-1)
-        XX=DIAG*B(I)
-        IF (KS.GT.KE) GO TO 30
-        DO 20 K=KS,KE
-          J=JLU(K)
-          XX=XX+ALU(K)*B(J)
-20      CONTINUE
-30      X(I)=X(I)+DIAG*XX
-        IF (KS.GT.KE) GO TO 50
-        DO 40 K=KS,KE
-          J=JLU(K)
-          X(J)=X(J)+ALU(K)*XX
-40      CONTINUE
-50    CONTINUE
-      RETURN
-      END
-      SUBROUTINE DSUBU0(N,ILU,JLU,ID,ALU,B,X)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION ILU(*),JLU(*),ID(*),ALU(*),B(N),X(N)
-C DO FORWARD AND BACK SUBSTITUTION TO SOLVE UT*U*X=B
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW OF LU
-C     JLU GIVES COLUMN NUMBER
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U
-C     ALU HAS NONZERO ELMENTS OF LU MATRIX STORED BY ROW
-C       WITH RECIPROCALS OF DIAGONAL ELEMENTS OF U
-C     B IS THE RHS VECTOR
-C     X IS THE SOLUTION VECTOR
-C
-      NP=N+1
-      DO 10 I=1,N
-        X(I)=B(I)
-10    CONTINUE
-C FORWARD SUBSTITUTION
-      DO 30 I=1,N
-        KS=ID(I)+1
-        KE=ILU(I+1)-1
-        XX=X(I)*ALU(KS-1)
-        X(I)=XX
-        IF (KS.GT.KE) GO TO 30
-        DO 20 K=KS,KE
-          J=JLU(K)
-          X(J)=X(J)-ALU(K)*XX
-20      CONTINUE
-30    CONTINUE
-C BACK SUBSTITUTION
-      DO 60 II=1,N
-        I=NP-II
-        KS=ID(I)+1
-        KE=ILU(I+1)-1
-        SUM=0.0D0
-        IF (KS.GT.KE) GO TO 50
-        DO 40 K=KS,KE
-          J=JLU(K)
-          SUM=SUM+ALU(K)*X(J)
-40      CONTINUE
-50      X(I)=(X(I)-SUM)*ALU(KS-1)
-60    CONTINUE
-      RETURN
-      END
-      SUBROUTINE DSUBL0(N,ILU,JLU,ID,ALU,B,X)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION ILU(*),JLU(*),ID(*),ALU(*),B(N),X(N)
-C DO FORWARD AND BACK SUBSTITUTION TO SOLVE L*LT*X=B
-C WHERE:
-C     N IS THE ORDER OF THE MATRIX
-C     ILU GIVES INDEX OF FIRST NONZERO ELEMENT IN ROW LU
-C     JLU GIVES THE COLUMN NUMBER
-C     ID GIVES INDEX OF DIAGONAL ELEMENT OF U
-C     ALU HAS NONZERO ELEMENTS OF LU MATRIX STORED BY ROW
-C       DIAGONAL ELEMENTS OF L ARE 1.0 AND NOT STORED
-C     B IS THE RHS VECTOR
-C     X IS THE SOLUTION VECTOR
-C
-      NP=N+1
-      DO 10 I=1,N
-        X(I)=B(I)
-10    CONTINUE
-C FORWARD SUBSTITUTION
-      DO 30 I=1,N
-        KS=ILU(I)
-        KE=ID(I)-1
-        IF (KS.GT.KE) GO TO 30
-        SUM=0.0D0
-        DO 20 K=KS,KE
-          J=JLU(K)
-          SUM=SUM+ALU(K)*X(J)
-20      CONTINUE
-        X(I)=X(I)-SUM
-30    CONTINUE
-C BACK SUBSTITUTION
-      DO 50 II=1,N
-        I=NP-II
-        KS=ILU(I)
-        KE=ID(I)-1
-        IF (KS.GT.KE) GO TO 50
-        XX=X(I)
-        IF (XX.EQ.0.0) GO TO 50
-        DO 40 K=KS,KE
-          J=JLU(K)
-          X(J)=X(J)-ALU(K)*XX
-40      CONTINUE
-50    CONTINUE
-      RETURN
-      END
-      DOUBLE PRECISION FUNCTION DGVV(V,W,N)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DIMENSION V(N),W(N)
-C     SUBROUTINE TO COMPUTE DOUBLE PRECISION VECTOR DOT PRODUCT.
-C
-      SUM = 0.0D0
-            DO 10 I = 1,N
-            SUM = SUM + V(I)*W(I)
-10          CONTINUE
-      DGVV = SUM
-      RETURN
-      END
diff --git a/src/elliptic/make.code.defn b/src/elliptic/make.code.defn
index 553d372..b00d006 100644
--- a/src/elliptic/make.code.defn
+++ b/src/elliptic/make.code.defn
@@ -4,7 +4,7 @@
 # Source files in this directory
 SRCS = Jacobian.cc \
        dense_Jacobian.cc lapack_wrapper.F77 \
-       row_sparse_Jacobian.cc ilucg_wrapper.F77 ilucg.f77
+       row_sparse_Jacobian.cc ilucg_wrapper.F77
 
 # Subdirectories containing source files
 SUBDIRS =