MODULE Sparse_Gaussian

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)


CONTAINS

SUBROUTINE dspslv (n, a, ia, ja, b, r, c, MAX, x, itemp, rtemp, ierr)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-07-14  Time: 20:38:15
 
!-----------------------------------------------------------------------
!          SOLUTION OF DOUBLE PRECISION SPARSE EQUATIONS
!-----------------------------------------------------------------------
!  DSPSLV CALLS DNSPIV WHICH USES SPARSE GAUSSIAN ELIMINATION WITH COLUMN
!  INTERCHANGES TO SOLVE THE LINEAR SYSTEM A X = B.  THE ELIMINATION PHASE
!  PERFORMS ROW OPERATIONS ON A AND B TO OBTAIN A UNIT UPPER TRIANGULAR
!  MATRIX U AND A VECTOR Y.  THE SOLUTION PHASE SOLVES U X = Y.

!  INPUT ARGUMENTS---

!  N      INTEGER NUMBER OF EQUATIONS AND UNKNOWNS

!  A      DOUBLE PRECISION ARRAY WITH ONE ENTRY PER NONZERO IN A,
!         CONTAINING THE ACTUAL NONZEROS. (SEE THE MATRIX STORAGE
!         DESCRIPTION BELOW)

!  IA     INTEGER ARRAY OF N+1 ENTRIES CONTAINING ROW POINTERS TO A
!         (SEE MATRIX STORAGE DESCRIPTION BELOW)

!  JA     INTEGER ARRAY WITH ONE ENTRY PER NONZERO IN A, CONTAINING COLUMN
!         NUMBERS OF THE NONZEROS OF A.  (SEE MATRIX STORAGE
!         DESCRIPTION BELOW)

!  B      DOUBLE PRECISION ARRAY OF N ENTRIES CONTAINING THE RIGHT
!         HAND SIDE DATA

!  R      INTEGER ARRAY OF N ENTRIES SPECIFYING THE ORDER OF THE
!         ROWS OF A (I.E., THE ELIMINATION ORDER FOR THE EQUATIONS)

!  C      INTEGER ARRAY OF N ENTRIES SPECIFYING THE ORDER OF THE
!         COLUMNS OF A.  C IS ALSO AN OUTPUT ARGUMENT

!  MAX    INTEGER NUMBER SPECIFYING MAXIMUM NUMBER OF OFF-DIAGONAL
!         NONZERO ENTRIES OF U WHICH MAY BE STORED

!  ITEMP  INTEGER ARRAY OF 3*N + MAX + 2 ENTRIES, FOR INTERNAL USE

!  RTEMP  DOUBLE PRECISION ARRAY OF N + MAX ENTRIES FOR INTERNAL USE


!  OUTPUT ARGUMENTS---

!  C      INTEGER ARRAY OF N ENTRIES SPECIFYING THE ORDER OF THE
!         COLUMNS OF U.  C IS ALSO AN INPUT ARGUMENT

!  X      DOUBLE PRECISION ARRAY OF N ENTRIES CONTAINING THE SOLUTION VECTOR

!  IERR   INTEGER NUMBER WHICH INDICATES ERROR CONDITIONS OR THE ACTUAL
!         NUMBER OF OFF-DIAGONAL ENTRIES IN U (FOR SUCCESSFUL COMPLETION)

!         IERR VALUES ARE---

!         0 < IERR             SUCCESSFUL COMPLETION.  IERR=MAX(1,M)
!                              WHERE M IS THE NUMBER OF OFF-DIAGONAL
!                              NONZERO ENTRIES OF U.

!         IERR = 0             ERROR.  N IS LESS THAN OR EQUAL TO 0

!         -N <= IERR < 0       ERROR.  ROW NUMBER IABS(IERR) OF A IS IS NULL

!         -2*N <= IERR < -N    ERROR.  ROW NUMBER IABS(IERR+N) HAS A
!                              DUPLICATE ENTRY

!         -3*N <= IERR < -2*N  ERROR.  ROW NUMBER IABS(IERR+2*N)
!                              HAS A ZERO PIVOT

!         -4*N <= IERR < -3*N  ERROR.  ROW NUMBER IABS(IERR+3*N)
!                              EXCEEDS STORAGE


!  STORAGE OF SPARSE MATRICES---

!  THE SPARSE MATRIX A IS STORED USING THREE ARRAYS IA, JA, AND A.
!  THE ARRAY A CONTAINS THE NONZEROS OF THE MATRIX ROW-BY-ROW, NOT
!  NECESSARILY IN ORDER OF INCREASING COLUMN NUMBER.  THE ARRAY JA CONTAINS
!  THE COLUMN NUMBERS CORRESPONDING TO THE NONZEROS STORED IN THE ARRAY A
!  (I.E., IF THE NONZERO STORED IN A(K) IS IN COLUMN J, THEN JA(K) = J).
!  THE ARRAY IA CONTAINS POINTERS TO THE ROWS OF NONZEROS/COLUMN INDICES IN
!  THE ARRAY A/JA (I.E., A(IA(I))/JA(IA(I)) IS THE FIRST ENTRY FOR ROW I IN
!  THE ARRAY A/JA).  IA(N+1) IS SET SO THAT IA(N+1) - IA(1) = THE NUMBER OF
!  NONZERO ELEMENTS IN A.
!------------------------

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: a(*)
INTEGER, INTENT(IN OUT)    :: ia(*)
INTEGER, INTENT(IN OUT)    :: ja(*)
REAL (dp), INTENT(IN OUT)  :: b(n)
INTEGER, INTENT(IN OUT)    :: r(n)
INTEGER, INTENT(IN OUT)    :: c(n)
INTEGER, INTENT(IN OUT)    :: MAX
REAL (dp), INTENT(IN OUT)  :: x(n)
INTEGER, INTENT(OUT)       :: itemp(*)
REAL (dp), INTENT(IN OUT)  :: rtemp(*)
INTEGER, INTENT(OUT)       :: ierr

INTEGER :: iu, ju, k, l, u, y, p

ierr = 0
IF (n <= 0) RETURN

!  SET INDICES TO DIVIDE TEMPORARY STORAGE FOR DNSPIV

y = 1
u = y + n
p = n + 1
iu = p + n + 1
ju = iu + n + 1

!  COMPUTE THE INVERSE PERMUTATION OF C

DO  k = 1,n
  l = c(k)
  itemp(l) = k
END DO

!  CALL DNSPIV TO PERFORM COMPUTATIONS

CALL dnspiv (n, ia, ja, a, b, MAX, r, c, itemp(1), x, rtemp(y), itemp(p),  &
             itemp(iu), itemp(ju), rtemp(u), ierr)
IF (ierr == 0) ierr = 1
RETURN
END SUBROUTINE dspslv



SUBROUTINE dnspiv(n, ia, ja, a, b, MAX, r, c, ic, x, y, p, iu, ju, u, ierr)
 
!  DNSPIV USES SPARSE GAUSSIAN ELIMINATION WITH COLUMN INTERCHANGES TO SOLVE
!  THE LINEAR SYSTEM A X = B.  THE ELIMINATION PHASE PERFORMS ROW OPERATIONS
!  ON A AND B TO OBTAIN A UNIT UPPER TRIANGULAR MATRIX U AND A VECTOR Y.
!  THE SOLUTION PHASE SOLVES U X = Y.


!  SEE DSPSLV FOR DESCRIPTIONS OF ALL INPUT AND OUTPUT ARGUMENTS
!  OTHER THAN THOSE DESCRIBED BELOW

!  IC  INTEGER ARRAY OF N ENTRIES WHICH IS THE INVERSE OF C
!      (I.E., IC(C(I)) = I). IC IS BOTH AN INPUT AND OUTPUT ARGUMENT.

!  INPUT ARGUMENTS (USED INTERNALLY ONLY)---

!  Y   DOUBLE PRECISION ARRAY OF N ENTRIES USED TO COMPUTE THE UPDATED
!      RIGHT HAND SIDE

!  P   INTEGER ARRAY OF N+1 ENTRIES USED FOR A LINKED LIST.
!      P(N+1) IS THE LIST HEADER, AND THE ENTRY FOLLOWING P(K) IS IN P(P(K)).
!      THUS, P(N+1) IS THE FIRST DATA ITEM, P(P(N+1)) IS THE SECOND, ETC.
!      A POINTER OF N+1 MARKS THE END OF THE LIST

!  IU  INTEGER ARRAY OF N+1 ENTRIES USED FOR ROW POINTERS TO U
!      (SEE MATRIX STORAGE DESCRIPTION BELOW)

!  JU  INTEGER ARRAY OF MAX ENTRIES USED FOR COLUMN NUMBERS OF THE NONZEROS IN
!      THE STRICT UPPER TRIANGLE OF U.  (SEE MATRIX STORAGE DESCRIPTION BELOW)

!  U   DOUBLE PRECISION ARRAY OF MAX ENTRIES USED FOR THE ACTUAL NONZEROS IN
!      THE STRICT UPPER TRIANGLE OF U. (SEE MATRIX STORAGE DESCRIPTION BELOW)


!  STORAGE OF SPARSE MATRICES---

!  THE SPARSE MATRIX A IS STORED USING THREE ARRAYS IA, JA, AND A.
!  THE ARRAY A CONTAINS THE NONZEROS OF THE MATRIX ROW-BY-ROW, NOT
!  NECESSARILY IN ORDER OF INCREASING COLUMN NUMBER.  THE ARRAY JA CONTAINS
!  THE COLUMN NUMBERS CORRESPONDING TO THE NONZEROS STORED IN THE ARRAY A
!  (I.E., IF THE NONZERO STORED IN A(K) IS IN COLUMN J, THEN JA(K) = J).
!  THE ARRAY IA CONTAINS POINTERS TO THE ROWS OF NONZEROS/COLUMN INDICES IN
!  THE ARRAY A/JA (I.E., A(IA(I))/JA(IA(I)) IS THE FIRST ENTRY FOR ROW I IN
!  THE ARRAY A/JA). IA(N+1) IS SET SO THAT IA(N+1) - IA(1) = THE NUMBER OF
!  NONZEROS IN A.  IU, JU, AND U ARE USED IN A SIMILAR WAY TO STORE THE STRICT
!  UPPER TRIANGLE OF U, EXCEPT THAT JU ACTUALLY CONTAINS C(J) INSTEAD OF J


INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(IN)        :: ia(*)
INTEGER, INTENT(IN)        :: ja(*)
REAL (dp), INTENT(IN)      :: a(*)
REAL (dp), INTENT(IN)      :: b(n)
INTEGER, INTENT(IN)        :: MAX
INTEGER, INTENT(IN)        :: r(n)
INTEGER, INTENT(IN OUT)    :: c(n)
INTEGER, INTENT(IN OUT)    :: ic(n)
REAL (dp), INTENT(OUT)     :: x(n)
REAL (dp), INTENT(IN OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: p(*)
INTEGER, INTENT(OUT)       :: iu(*)
INTEGER, INTENT(IN OUT)    :: ju(MAX)
REAL (dp), INTENT(IN OUT)  :: u(MAX)
INTEGER, INTENT(OUT)       :: ierr

REAL (dp)  :: dk, lki, xpv, xpvmax, yk

INTEGER    :: ck, i, j, jaj, jmin, jmax, juj, juptr, k, maxc, maxcl, nzcnt,  &
              pk, ppk, pv, v, vi, vj, vk

!  INITIALIZE WORK STORAGE AND POINTERS TO JU

x(1:n) = 0.d0
iu(1) = 1
juptr = 0

!  PERFORM SYMBOLIC AND NUMERIC FACTORIZATION ROW BY ROW
!  VK (VI,VJ) IS THE GRAPH VERTEX FOR ROW K (I,J) OF U

DO  k = 1, n
  
!  INITIALIZE LINKED LIST AND FREE STORAGE FOR THIS ROW
!  THE R(K)-TH ROW OF A BECOMES THE K-TH ROW OF U.
  
  p(n+1) = n + 1
  vk = r(k)
  
!  SET UP ADJACENCY LIST FOR VK, ORDERED IN CURRENT COLUMN ORDER OF U.
!  THE LOOP INDEX GOES DOWNWARD TO EXPLOIT ANY COLUMNS FROM A
!  IN CORRECT RELATIVE ORDER
  
  jmin = ia(vk)
  jmax = ia(vk+1) - 1
  IF (jmin > jmax) GO TO 170
  j = jmax
  20   jaj = ja(j)
  vj = ic(jaj)
  
!  STORE A(K,J) IN WORK VECTOR
  
  x(vj) = a(j)
!  THIS CODE INSERTS VJ INTO ADJACENCY LIST OF VK
  ppk = n + 1
  30   pk = ppk
  ppk = p(pk)
  IF (ppk-vj < 0.0) THEN
    GO TO    30
  ELSE IF (ppk-vj == 0.0) THEN
    GO TO   180
  END IF
  p(vj) = ppk
  p(pk) = vj
  j = j - 1
  IF (j >= jmin) GO TO 20
  
!  THE FOLLOWING CODE COMPUTES THE K-TH ROW OF U
  
  vi = n + 1
  yk = b(vk)
  50   vi = p(vi)
  IF (vi < k) THEN
    
!  VI LT VK -- PROCESS THE L(K,I) ELEMENT AND MERGE THE
!  ADJACENCY OF VI WITH THE ORDERED ADJACENCY OF VK
    
    lki = -x(vi)
    x(vi) = 0.d0
    
!  ADJUST RIGHT HAND SIDE TO REFLECT ELIMINATION
    
    yk = yk + lki * y(vi)
    ppk = vi
    jmin = iu(vi)
    jmax = iu(vi+1) - 1
    IF (jmin > jmax) GO TO 50
    DO  j = jmin, jmax
      juj = ju(j)
      vj = ic(juj)
      
!  IF VJ IS ALREADY IN THE ADJACENCY OF VK, SKIP THE INSERTION
      
      IF (x(vj) == 0.d0) THEN
        
!  INSERT VJ IN ADJACENCY LIST OF VK.
!  RESET PPK TO VI IF WE HAVE PASSED THE CORRECT INSERTION SPOT.
!  (THIS HAPPENS WHEN THE ADJACENCY OF
!  VI IS NOT IN CURRENT COLUMN ORDER DUE TO PIVOTING.)
        
        IF (vj-ppk < 0.0) THEN
          GO TO    60
        ELSE IF (vj-ppk == 0.0) THEN
          GO TO    90
        ELSE
          GO TO    70
        END IF
        60  ppk = vi
        70  pk = ppk
        ppk = p(pk)
        IF (ppk-vj < 0.0) THEN
          GO TO    70
        ELSE IF (ppk-vj == 0.0) THEN
          GO TO    90
        END IF
        p(vj) = ppk
        p(pk) = vj
        ppk = vj
      END IF
      
!  COMPUTE L(K,J) = L(K,J) - L(K,I)*U(I,J) FOR L(K,I) NONZERO
!  COMPUTE U*(K,J) = U*(K,J) - L(K,I)*U(I,J) FOR U(K,J) NONZERO
!  (U*(K,J) = U(K,J)*D(K,K))
      
      90  x(vj) = x(vj) + lki * u(j)
    END DO
    GO TO 50
  END IF
  
!  PIVOT--INTERCHANGE LARGEST ENTRY OF K-TH ROW OF U WITH THE DIAGONAL ENTRY.
  
!  FIND LARGEST ENTRY, COUNTING OFF-DIAGONAL NONZEROS
  
  IF (vi > n) GO TO 190
  xpvmax = ABS(x(vi))
  maxc = vi
  nzcnt = 0
  pv = vi
  110   v = pv
  pv = p(pv)
  IF (pv <= n) THEN
    nzcnt = nzcnt + 1
    xpv = ABS(x(pv))
    IF (xpv <= xpvmax) GO TO 110
    xpvmax = xpv
    maxc = pv
    maxcl = v
    GO TO 110
  END IF
  IF (xpvmax == 0.d0) GO TO 190
  
!  IF VI = K, THEN THERE IS AN ENTRY FOR DIAGONAL
!  WHICH MUST BE DELETED.  OTHERWISE, DELETE THE
!  ENTRY WHICH WILL BECOME THE DIAGONAL ENTRY
  
  IF (vi /= k) THEN
    IF (vi /= maxc) THEN
      p(maxcl) = p(maxc)
      GO TO 120
    END IF
  END IF
  vi = p(vi)
  
!  COMPUTE D(K) = 1/L(K,K) AND PERFORM INTERCHANGE.
  
  120 dk = 1.d0 / x(maxc)
  x(maxc) = x(k)
  i = c(k)
  c(k) = c(maxc)
  c(maxc) = i
  ck = c(k)
  ic(ck) = k
  ic(i) = maxc
  x(k) = 0.d0
  
!  UPDATE RIGHT HAND SIDE.
  
  y(k) = yk * dk
  
!  COMPUTE VALUE FOR IU(K+1) AND CHECK FOR STORAGE OVERFLOW
  
  iu(k+1) = iu(k) + nzcnt
  IF (iu(k+1) > MAX+1) GO TO 200
  
!  MOVE COLUMN INDICES FROM LINKED LIST TO JU.
!  COLUMNS ARE STORED IN CURRENT ORDER WITH ORIGINAL
!  COLUMN NUMBER (C(J)) STORED FOR CURRENT COLUMN J
  
  IF (vi <= n) THEN
    j = vi
    130  juptr = juptr + 1
    ju(juptr) = c(j)
    u(juptr) = x(j) * dk
    x(j) = 0.d0
    j = p(j)
    IF (j <= n) GO TO 130
  END IF
END DO

!  BACKSOLVE U X = Y, AND REORDER X TO CORRESPOND WITH A

k = n
DO  i = 1, n
  yk = y(k)
  jmin = iu(k)
  jmax = iu(k+1) - 1
  IF (jmin <= jmax) THEN
    DO  j = jmin, jmax
      juj = ju(j)
      juj = ic(juj)
      yk = yk - u(j) * y(juj)
    END DO
  END IF
  y(k) = yk
  ck = c(k)
  x(ck) = yk
  k = k - 1
END DO

!  RETURN WITH IERR = NUMBER OF OFF-DIAGONAL NONZEROS IN U

ierr = iu(n+1) - iu(1)
RETURN

!  ERROR RETURNS

!  ROW K OF A IS NULL

170 ierr = -k
RETURN

!  ROW K OF A HAS A DUPLICATE ENTRY

180 ierr = -(n+k)
RETURN

!  ZERO PIVOT IN ROW K

190 ierr = -(2*n+k)
RETURN

!  STORAGE FOR U EXCEEDED ON ROW K

200 ierr = -(3*n+k)
RETURN
END SUBROUTINE dnspiv

END MODULE Sparse_Gaussian