ILUSOLVDS Subroutine

      SUBROUTINE ILUSOLVDS(NN, Q, LENQ, ROWPOSP, COLPOSP, B)
C
C     Computes Q^{-1}*B -- returns result as B.
C---------------------------------------------------------------------
C
C     Input variables:
C
C       Q        triangular factors of preconditioning matrix, stored
C                in the general sparse row scheme.
C       ROWPOSP  indices of row-starts within B.
C       COLPOSP  column indices for matrix B stored in the general
C                sparse row scheme.
C       NN       logical row dimension of Q.
C       LENQ     number of data entries in Q.
C       B        right hand side -- should have dimension NN.
C
C     Output variables:
C
C       B        solution of Qx = B.
C
C---------------------------------------------------------------------
C
      USE REAL_PRECISION
      INTEGER, INTENT(IN) ::  NN, LENQ, ROWPOSP(NN+1), COLPOSP(LENQ)
      REAL (KIND=R8), INTENT(IN) :: Q(LENQ) 
      REAL (KIND=R8), INTENT(IN OUT) :: B(NN)
C 
C LOCAL VARIABLES
      INTEGER:: DIAG(NN), I, K, J
C------------------------------------------------ COMPUTE B = INV(L)*B
      B(1) = B(1)/Q(1)
      DIAG(1) = 1
      DO  I = 2, NN
        K = ROWPOSP(I)
 42     J = COLPOSP(K)
        IF (J .LT. I) THEN
          B(I) = B(I) - Q(K)*B(J)
          K = K + 1
          GOTO 42
        ELSE 
          DIAG(I) = K
          B(I) = B(I)/Q(K)
        ENDIF
      END DO
C------------------------------------------------ COMPUTE B = INV(U)*B
      DO  I = NN-1, 1, -1
        DO  K = DIAG(I)+1, ROWPOSP(I+1)-1
          B(I) = B(I) - Q(K)*B(COLPOSP(K))
        END DO
      END DO
      RETURN
      END SUBROUTINE ILUSOLVDS