chem/KPP/kpp/kpp-2.1/util/sutil.f

   1
   2 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   3       SUBROUTINE KppDecomp( JVS, IER )
   4 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   5 C        Sparse LU factorization
   6 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   7
   8       INCLUDE 'KPP_ROOT_Parameters.h'
   9       INCLUDE 'KPP_ROOT_Sparse.h'
  10
  11       INTEGER  IER
  12       KPP_REAL JVS(KPP_LU_NONZERO), W(KPP_NVAR)
  13       INTEGER  k, kk, j, jj
  14       KPP_REAL a
  15
  16       IER = 0
  17       DO k=1,NVAR
  18         IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
  19             IER = k
  20             RETURN
  21         END IF
  22         DO kk = LU_CROW(k), LU_CROW(k+1)-1
  23               W( LU_ICOL(kk) ) = JVS(kk)
  24         END DO
  25         DO kk = LU_CROW(k), LU_DIAG(k)-1
  26             j = LU_ICOL(kk)
  27             a = -W(j) / JVS( LU_DIAG(j) )
  28             W(j) = -a
  29             DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
  30                W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
  31             END DO
  32          END DO
  33          DO kk = LU_CROW(k), LU_CROW(k+1)-1
  34             JVS(kk) = W( LU_ICOL(kk) )
  35          END DO
  36       END DO
  37       RETURN
  38       END
  39
  40 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  41       SUBROUTINE KppDecompCmplx( JVS, IER )
  42 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  43 C        Sparse LU factorization, complex
  44 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  45
  46       INCLUDE 'KPP_ROOT_Parameters.h'
  47       INCLUDE 'KPP_ROOT_Sparse.h'
  48
  49       INTEGER  IER
  50       DOUBLE COMPLEX JVS(KPP_LU_NONZERO), W(KPP_NVAR), a
  51       INTEGER  k, kk, j, jj
  52
  53       IER = 0
  54       DO k=1,NVAR
  55         IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
  56             IER = k
  57             RETURN
  58         END IF
  59         DO kk = LU_CROW(k), LU_CROW(k+1)-1
  60               W( LU_ICOL(kk) ) = JVS(kk)
  61         END DO
  62         DO kk = LU_CROW(k), LU_DIAG(k)-1
  63             j = LU_ICOL(kk)
  64             a = -W(j) / JVS( LU_DIAG(j) )
  65             W(j) = -a
  66             DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
  67                W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
  68             END DO
  69          END DO
  70          DO kk = LU_CROW(k), LU_CROW(k+1)-1
  71             JVS(kk) = W( LU_ICOL(kk) )
  72          END DO
  73       END DO
  74       RETURN
  75       END
  76
  77 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  78       SUBROUTINE KppSolveIndirect( JVS, X )
  79 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  80 C        Sparse solve subroutine using indirect addressing
  81 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  82       INCLUDE 'KPP_ROOT_Parameters.h'
  83       INCLUDE 'KPP_ROOT_Sparse.h'
  84
  85       INTEGER i, j
  86       KPP_REAL JVS(KPP_LU_NONZERO), X(KPP_NVAR), sum
  87
  88       DO i=1,NVAR
  89          DO j = LU_CROW(i), LU_DIAG(i)-1
  90              X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
  91          END DO
  92       END DO
  93
  94       DO i=NVAR,1,-1
  95         sum = X(i);
  96         DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
  97           sum = sum - JVS(j)*X(LU_ICOL(j));
  98         END DO
  99         X(i) = sum/JVS(LU_DIAG(i));
 100       END DO
 101
 102       RETURN
 103       END
 104
 105 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 106       SUBROUTINE KppSolveCmplx( JVS, X )
 107 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 108 C        Complex sparse solve subroutine using indirect addressing
 109 C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 110       INCLUDE 'KPP_ROOT_Parameters.h'
 111       INCLUDE 'KPP_ROOT_Sparse.h'
 112
 113       INTEGER i, j
 114       DOUBLE COMPLEX JVS(KPP_LU_NONZERO), X(KPP_NVAR), sum
 115
 116       DO i=1,NVAR
 117          DO j = LU_CROW(i), LU_DIAG(i)-1
 118              X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
 119          END DO
 120       END DO
 121
 122       DO i=NVAR,1,-1
 123         sum = X(i);
 124         DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
 125           sum = sum - JVS(j)*X(LU_ICOL(j));
 126         END DO
 127         X(i) = sum/JVS(LU_DIAG(i));
 128       END DO
 129
 130       RETURN
 131       END