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

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