! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! Linear Algebra Data and Routines File ! ! Generated by KPP-2.2 symbolic chemistry Kinetics PreProcessor ! (http://www.cs.vt.edu/~asandu/Software/KPP) ! KPP is distributed under GPL, the general public licence ! (http://www.gnu.org/copyleft/gpl.html) ! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa ! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech ! With important contributions from: ! M. Damian, Villanova University, USA ! R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany ! ! File : gckpp_LinearAlgebra.f90 ! Time : Tue Mar 16 10:29:04 2010 ! Working directory : /mnt/lstr04/srv/home/c/ccarouge/KPP/geoschem_kppfiles/v8-03-01/Paulot_scheme ! Equation file : gckpp.kpp ! Output root filename : gckpp ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ MODULE gckpp_LinearAlgebra USE gckpp_Parameters USE gckpp_JacobianSP IMPLICIT NONE CONTAINS ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! SPARSE_UTIL - SPARSE utility functions ! Arguments : ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppDecomp( JVS, IER ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Sparse LU factorization ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: IER REAL(kind=dp) :: JVS(LU_NONZERO), W(NVAR), a INTEGER :: k, kk, j, jj a = 0. ! mz_rs_20050606 IER = 0 DO k=1,NVAR ! mz_rs_20050606: don't check if real value == 0 ! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN IER = k RETURN END IF DO kk = LU_CROW(k), LU_CROW(k+1)-1 W( LU_ICOL(kk) ) = JVS(kk) END DO DO kk = LU_CROW(k), LU_DIAG(k)-1 j = LU_ICOL(kk) a = -W(j) / JVS( LU_DIAG(j) ) W(j) = -a DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1 W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj) END DO END DO DO kk = LU_CROW(k), LU_CROW(k+1)-1 JVS(kk) = W( LU_ICOL(kk) ) END DO END DO END SUBROUTINE KppDecomp ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppDecompCmplx( JVS, IER ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Sparse LU factorization, complex ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: IER DOUBLE COMPLEX :: JVS(LU_NONZERO), W(NVAR), a REAL(kind=dp) :: b = 0.0 INTEGER :: k, kk, j, jj IER = 0 DO k=1,NVAR IF ( ABS(JVS(LU_DIAG(k))) < TINY(b) ) THEN IER = k RETURN END IF DO kk = LU_CROW(k), LU_CROW(k+1)-1 W( LU_ICOL(kk) ) = JVS(kk) END DO DO kk = LU_CROW(k), LU_DIAG(k)-1 j = LU_ICOL(kk) a = -W(j) / JVS( LU_DIAG(j) ) W(j) = -a DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1 W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj) END DO END DO DO kk = LU_CROW(k), LU_CROW(k+1)-1 JVS(kk) = W( LU_ICOL(kk) ) END DO END DO END SUBROUTINE KppDecompCmplx ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppDecompCmplxR( JVSR, JVSI, IER ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Sparse LU factorization, complex ! (Real and Imaginary parts are used instead of complex data type) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: IER REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO) REAL(kind=dp) :: WR(NVAR), WI(NVAR), ar, ai, den INTEGER :: k, kk, j, jj IER = 0 ar = 0.0 DO k=1,NVAR IF ( ( ABS(JVSR(LU_DIAG(k))) < TINY(ar) ) .AND. & ( ABS(JVSI(LU_DIAG(k))) < TINY(ar) ) ) THEN IER = k RETURN END IF DO kk = LU_CROW(k), LU_CROW(k+1)-1 WR( LU_ICOL(kk) ) = JVSR(kk) WI( LU_ICOL(kk) ) = JVSI(kk) END DO DO kk = LU_CROW(k), LU_DIAG(k)-1 j = LU_ICOL(kk) den = JVSR(LU_DIAG(j))**2 + JVSI(LU_DIAG(j))**2 ar = -(WR(j)*JVSR(LU_DIAG(j)) + WI(j)*JVSI(LU_DIAG(j)))/den ai = -(WI(j)*JVSR(LU_DIAG(j)) - WR(j)*JVSI(LU_DIAG(j)))/den WR(j) = -ar WI(j) = -ai DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1 WR( LU_ICOL(jj) ) = WR( LU_ICOL(jj) ) + ar*JVSR(jj) - ai*JVSI(jj) WI( LU_ICOL(jj) ) = WI( LU_ICOL(jj) ) + ar*JVSI(jj) + ai*JVSR(jj) END DO END DO DO kk = LU_CROW(k), LU_CROW(k+1)-1 JVSR(kk) = WR( LU_ICOL(kk) ) JVSI(kk) = WI( LU_ICOL(kk) ) END DO END DO END SUBROUTINE KppDecompCmplxR ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveIndirect( JVS, X ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Sparse solve subroutine using indirect addressing ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: i, j REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR), sum DO i=1,NVAR DO j = LU_CROW(i), LU_DIAG(i)-1 X(i) = X(i) - JVS(j)*X(LU_ICOL(j)); END DO END DO DO i=NVAR,1,-1 sum = X(i); DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1 sum = sum - JVS(j)*X(LU_ICOL(j)); END DO X(i) = sum/JVS(LU_DIAG(i)); END DO END SUBROUTINE KppSolveIndirect ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveTRIndirect( JVS, X ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Complex sparse solve transpose subroutine using indirect addressing ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: i, j REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR) DO i=1,NVAR X(i) = X(i)/JVS(LU_DIAG(i)) ! subtract all nonzero elements in row i of JVS from X DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1 X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i) END DO END DO DO i=NVAR, 1, -1 ! subtract all nonzero elements in row i of JVS from X DO j=LU_CROW(i),LU_DIAG(i)-1 X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i) END DO END DO END SUBROUTINE KppSolveTRIndirect ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveCmplx( JVS, X ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Complex sparse solve subroutine using indirect addressing ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: i, j DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR), sum DO i=1,NVAR DO j = LU_CROW(i), LU_DIAG(i)-1 X(i) = X(i) - JVS(j)*X(LU_ICOL(j)); END DO END DO DO i=NVAR,1,-1 sum = X(i); DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1 sum = sum - JVS(j)*X(LU_ICOL(j)); END DO X(i) = sum/JVS(LU_DIAG(i)); END DO END SUBROUTINE KppSolveCmplx ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveCmplxR( JVSR, JVSI, XR, XI ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Complex sparse solve subroutine using indirect addressing ! (Real and Imaginary parts are used instead of complex data type) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: i, j REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), sumr, sumi, den DO i=1,NVAR DO j = LU_CROW(i), LU_DIAG(i)-1 XR(i) = XR(i) - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j))) XI(i) = XI(i) - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j))) END DO END DO DO i=NVAR,1,-1 sumr = XR(i); sumi = XI(i) DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1 sumr = sumr - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j))) sumi = sumi - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j))) END DO den = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2 XR(i) = (sumr*JVSR(LU_DIAG(i)) + sumi*JVSI(LU_DIAG(i)))/den XI(i) = (sumi*JVSR(LU_DIAG(i)) - sumr*JVSI(LU_DIAG(i)))/den END DO END SUBROUTINE KppSolveCmplxR ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveTRCmplx( JVS, X ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Complex sparse solve transpose subroutine using indirect addressing ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: i, j DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR) DO i=1,NVAR X(i) = X(i)/JVS(LU_DIAG(i)) ! subtract all nonzero elements in row i of JVS from X DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1 X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i) END DO END DO DO i=NVAR, 1, -1 ! subtract all nonzero elements in row i of JVS from X DO j=LU_CROW(i),LU_DIAG(i)-1 X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i) END DO END DO END SUBROUTINE KppSolveTRCmplx ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveTRCmplxR( JVSR, JVSI, XR, XI ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Complex sparse solve transpose subroutine using indirect addressing ! (Real and Imaginary parts are used instead of complex data type) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE gckpp_Parameters USE gckpp_JacobianSP INTEGER :: i, j REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), den DO i=1,NVAR den = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2 XR(i) = (XR(i)*JVSR(LU_DIAG(i)) + XI(i)*JVSI(LU_DIAG(i)))/den XI(i) = (XI(i)*JVSR(LU_DIAG(i)) - XR(i)*JVSI(LU_DIAG(i)))/den ! subtract all nonzero elements in row i of JVS from X DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1 XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i)) XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i)) END DO END DO DO i=NVAR, 1, -1 ! subtract all nonzero elements in row i of JVS from X DO j=LU_CROW(i),LU_DIAG(i)-1 XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i)) XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i)) END DO END DO END SUBROUTINE KppSolveTRCmplxR ! ! Next few commented subroutines perform sparse big linear algebra ! !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !SUBROUTINE KppDecompBig( JVS, IP, IER ) !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !! Sparse LU factorization !! for the Runge Kutta (3n)x(3n) linear system !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! USE gckpp_Parameters ! USE gckpp_JacobianSP ! ! INTEGER :: IP3(3), IER, IP(3,NVAR) ! REAL(kind=dp) :: JVS(3,3,LU_NONZERO), W(3,3,NVAR), a(3,3), E(3,3) ! INTEGER :: k, kk, j, jj ! ! a = 0.0d0 ! IER = 0 ! DO k=1,NVAR ! DO kk = LU_CROW(k), LU_CROW(k+1)-1 ! W( 1:3,1:3,LU_ICOL(kk) ) = JVS(1:3,1:3,kk) ! END DO ! DO kk = LU_CROW(k), LU_DIAG(k)-1 ! j = LU_ICOL(kk) ! E(1:3,1:3) = JVS( 1:3,1:3,LU_DIAG(j) ) ! ! CALL DGETRF(3,3,E,3,IP3,IER) ! CALL FAC3(E,IP3,IER) ! IF ( IER /= 0 ) RETURN ! ! a = W(j) / JVS( LU_DIAG(j) ) ! a(1:3,1:3) = W( 1:3,1:3,j ) ! ! CALL DGETRS ('N',3,3,E,3,IP3,a,3,IER) ! CALL SOL3('N',E,IP3,a(1,1)) ! CALL SOL3('N',E,IP3,a(1,2)) ! CALL SOL3('N',E,IP3,a(1,3)) ! W(1:3,1:3,j) = a(1:3,1:3) ! DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1 ! W( 1:3,1:3,LU_ICOL(jj) ) = W( 1:3,1:3,LU_ICOL(jj) ) & ! - MATMUL( a(1:3,1:3) , JVS(1:3,1:3,jj) ) ! END DO ! END DO ! DO kk = LU_CROW(k), LU_CROW(k+1)-1 ! JVS(1:3,1:3,kk) = W( 1:3,1:3,LU_ICOL(kk) ) ! END DO ! END DO ! ! DO k=1,NVAR ! ! CALL WGEFA(JVS(1,1,LU_DIAG(k)),3,3,IP(1,k),IER) ! ! CALL DGETRF(3,3,JVS(1,1,LU_DIAG(k)),3,IP(1,k),IER) ! CALL FAC3(JVS(1,1,LU_DIAG(k)),IP(1,k),IER) ! IF ( IER /= 0 ) RETURN ! END DO ! !END SUBROUTINE KppDecompBig ! ! !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !SUBROUTINE KppSolveBig( JVS, IP, X ) !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !! Sparse solve subroutine using indirect addressing !! for the Runge Kutta (3n)x(3n) linear system !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! USE gckpp_Parameters ! USE gckpp_JacobianSP ! ! INTEGER :: i, j, k, m, IP3(3), IP(3,NVAR), IER ! REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR), sum(3) ! ! DO i=1,NVAR ! DO j = LU_CROW(i), LU_DIAG(i)-1 ! !X(1:3,i) = X(1:3,i) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j))); ! DO k=1,3 ! DO m=1,3 ! X(k,i) = X(k,i) - JVS(k,m,j)*X(m,LU_ICOL(j)) ! END DO ! END DO ! END DO ! END DO ! ! DO i=NVAR,1,-1 ! sum(1:3) = X(1:3,i); ! DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1 ! !sum(1:3) = sum(1:3) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j))); ! DO k=1,3 ! DO m=1,3 ! sum(k) = sum(k) - JVS(k,m,j)*X(m,LU_ICOL(j)) ! END DO ! END DO ! END DO ! ! X(i) = sum/JVS(LU_DIAG(i)); ! ! CALL DGETRS ('N',3,1,JVS(1:3,1:3,LU_DIAG(i)),3,IP(1,i),sum,3,0) ! ! CALL WGESL('N',JVS(1,1,LU_DIAG(i)),3,3,IP(1,i),sum) ! CALL SOL3('N',JVS(1,1,LU_DIAG(i)),IP(1,i),sum) ! X(1:3,i) = sum(1:3) ! END DO ! !END SUBROUTINE KppSolveBig ! ! !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !SUBROUTINE KppSolveBigTR( JVS, IP, X ) !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !! Big sparse transpose solve using indirect addressing !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! USE gckpp_Parameters ! USE gckpp_JacobianSP ! ! INTEGER :: i, j, k, m, IP(3,NVAR) ! REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR) ! ! DO i=1,NVAR ! ! X(i) = X(i)/JVS(LU_DIAG(i)) ! CALL SOL3('T',JVS(1,1,LU_DIAG(i)),IP(1,i),X(1,i)) ! DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1 ! !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) & ! ! - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) ) ! DO k=1,3 ! DO m=1,3 ! X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i) ! END DO ! END DO ! END DO ! END DO ! ! DO i=NVAR, 1, -1 ! DO j=LU_CROW(i),LU_DIAG(i)-1 ! !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) & ! ! - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) ) ! DO k=1,3 ! DO m=1,3 ! X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i) ! END DO ! END DO ! END DO ! END DO ! !END SUBROUTINE KppSolveBigTR ! ! ! !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !SUBROUTINE FAC3(A,IPVT,INFO) !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !! FAC3 FACTORS THE MATRIX A (3,3) BY !! GAUSS ELIMINATION WITH PARTIAL PIVOTING !! LINPACK - LIKE !! !! Remove comments to perform pivoting !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !! ! REAL(kind=dp) :: A(3,3) ! INTEGER :: IPVT(3),INFO !! INTEGER :: L !! REAL(kind=dp) :: t, dmax, da, TMP(3) ! REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0 ! ! info = 0 !! t = TINY(da) !! !! da = ABS(A(1,1)); L = 1 !! IF ( ABS(A(2,1))>da ) THEN !! da = ABS(A(2,1)); L = 2 !! IF ( ABS(A(3,1))>da ) THEN !! L = 3 !! END IF !! END IF !! IPVT(1) = L !! IF (L /=1 ) THEN !! TMP(1:3) = A(L,1:3) !! A(L,1:3) = A(1,1:3) !! A(1,1:3) = TMP(1:3) !! END IF !! IF (ABS(A(1,1)) < t) THEN !! info = 1 !! return !! END IF !! ! A(2,1) = A(2,1)/A(1,1) ! A(2,2) = A(2,2) - A(2,1)*A(1,2) ! A(2,3) = A(2,3) - A(2,1)*A(1,3) ! A(3,1) = A(3,1)/A(1,1) ! A(3,2) = A(3,2) - A(3,1)*A(1,2) ! A(3,3) = A(3,3) - A(3,1)*A(1,3) ! !! IPVT(2) = 2 !! IF (ABS(A(3,2))>ABS(A(2,2))) THEN !! IPVT(2) = 3 !! TMP(2:3) = A(3,2:3) !! A(3,2:3) = A(2,2:3) !! A(2,2:3) = TMP(2:3) !! END IF !! IF (ABS(A(2,2)) < t) THEN !! info = 1 !! return !! END IF !! ! A(3,2) = A(3,2)/A(2,2) ! A(3,3) = A(3,3) - A(3,2)*A(2,3) ! IPVT(3) = 3 ! !END SUBROUTINE FAC3 ! ! !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !SUBROUTINE SOL3(Trans,A,IPVT,b) !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !! SOL3 solves the system 3x3 !! A * x = b or trans(a) * x = b !! using the factors computed by WGEFA. !! !! Trans = 'N' to solve A*x = b , !! = 'T' to solve transpose(A)*x = b !! LINPACK - LIKE !! !! Remove comments to use pivoting !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! CHARACTER :: Trans ! REAL(kind=dp) :: a(3,3),b(3) ! INTEGER :: IPVT(3) !! INTEGER :: L !! REAL(kind=dp) :: TMP ! ! SELECT CASE (Trans) ! ! CASE ('n','N') ! Solve A * x = b ! !! Solve L*y = b !! L = IPVT(1) !! IF (L /= 1) THEN !! TMP = B(1); B(1) = B(L); B(L) = TMP !! END IF ! b(2) = b(2)-A(2,1)*b(1) ! b(3) = b(3)-A(3,1)*b(1) ! !! L = IPVT(2) !! IF (L /= 2) THEN !! TMP = B(2); B(2) = B(L); B(L) = TMP !! END IF ! b(3) = b(3)-A(3,2)*b(2) ! !! Solve U*x = y ! b(3) = b(3)/A(3,3) ! b(2) = (b(2)-A(2,3)*b(3))/A(2,2) ! b(1) = (b(1)-A(1,3)*b(3)-A(1,2)*b(2))/A(1,1) ! ! ! CASE ('t','T') ! Solve transpose(A) * x = b ! !! Solve transpose(U)*y = b ! b(1) = b(1)/A(1,1) ! b(2) = (b(2)-A(1,2)*b(1))/A(2,2) ! b(3) = (b(3)-A(1,3)*b(1)-A(2,3)*b(2))/A(3,3) ! !! Solve transpose(L)*x = y ! b(2) = b(2)-A(3,2)*b(3) !! L = ipvt(2) !! IF (L /= 2) THEN !! TMP = B(2); B(2) = B(L); B(L) = TMP !! END IF ! b(1) = b(1)-A(3,1)*b(3)-A(2,1)*b(2) !! L = ipvt(1) !! IF (L /= 1) THEN !! TMP = B(1); B(1) = B(L); B(L) = TMP !! END IF ! ! END SELECT ! !END SUBROUTINE SOL3 ! End of SPARSE_UTIL function ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! KppSolve - sparse back substitution ! Arguments : ! JVS - sparse Jacobian of variables ! X - Vector for variables ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolve ( JVS, X ) ! JVS - sparse Jacobian of variables REAL(kind=dp) :: JVS(LU_NONZERO) ! X - Vector for variables REAL(kind=dp) :: X(NVAR) X(75) = X(75)-JVS(284)*X(65) X(77) = X(77)-JVS(294)*X(76) X(84) = X(84)-JVS(331)*X(47)-JVS(332)*X(74) X(85) = X(85)-JVS(338)*X(57) X(86) = X(86)-JVS(348)*X(81) X(87) = X(87)-JVS(354)*X(78) X(90) = X(90)-JVS(402)*X(83)-JVS(403)*X(84) X(91) = X(91)-JVS(413)*X(70)-JVS(414)*X(76)-JVS(415)*X(77)-JVS(416)*X(81)-JVS(417)*X(84)-JVS(418)*X(85)-JVS(419)*X(86)& &-JVS(420)*X(90) X(92) = X(92)-JVS(448)*X(50) X(93) = X(93)-JVS(456)*X(82) X(94) = X(94)-JVS(465)*X(44)-JVS(466)*X(56) X(95) = X(95)-JVS(475)*X(79) X(96) = X(96)-JVS(483)*X(67) X(98) = X(98)-JVS(502)*X(44)-JVS(503)*X(57) X(99) = X(99)-JVS(512)*X(51)-JVS(513)*X(52)-JVS(514)*X(55)-JVS(515)*X(60)-JVS(516)*X(83)-JVS(517)*X(87)-JVS(518)*X(90) X(100) = X(100)-JVS(541)*X(73)-JVS(542)*X(76)-JVS(543)*X(83)-JVS(544)*X(84)-JVS(545)*X(93) X(101) = X(101)-JVS(560)*X(73) X(102) = X(102)-JVS(568)*X(64) X(103) = X(103)-JVS(576)*X(80) X(104) = X(104)-JVS(583)*X(82)-JVS(584)*X(103) X(105) = X(105)-JVS(594)*X(70) X(107) = X(107)-JVS(610)*X(80)-JVS(611)*X(82)-JVS(612)*X(103) X(108) = X(108)-JVS(620)*X(63)-JVS(621)*X(85)-JVS(622)*X(98) X(109) = X(109)-JVS(632)*X(69)-JVS(633)*X(86) X(110) = X(110)-JVS(643)*X(67)-JVS(644)*X(86) X(111) = X(111)-JVS(654)*X(76) X(112) = X(112)-JVS(662)*X(52)-JVS(663)*X(61) X(113) = X(113)-JVS(676)*X(45)-JVS(677)*X(58) X(114) = X(114)-JVS(688)*X(60)-JVS(689)*X(62)-JVS(690)*X(63)-JVS(691)*X(64)-JVS(692)*X(68)-JVS(693)*X(69)-JVS(694)& &*X(70)-JVS(695)*X(73)-JVS(696)*X(74)-JVS(697)*X(79)-JVS(698)*X(81)-JVS(699)*X(82)-JVS(700)*X(84)-JVS(701)& &*X(86)-JVS(702)*X(90)-JVS(703)*X(92)-JVS(704)*X(93)-JVS(705)*X(94)-JVS(706)*X(95)-JVS(707)*X(96)-JVS(708)& &*X(97)-JVS(709)*X(98)-JVS(710)*X(100)-JVS(711)*X(101)-JVS(712)*X(102)-JVS(713)*X(103)-JVS(714)*X(104)-JVS(715)& &*X(105)-JVS(716)*X(106)-JVS(717)*X(107)-JVS(718)*X(108)-JVS(719)*X(109)-JVS(720)*X(110)-JVS(721)*X(111)& &-JVS(722)*X(112)-JVS(723)*X(113) X(115) = X(115)-JVS(747)*X(58)-JVS(748)*X(61)-JVS(749)*X(69)-JVS(750)*X(86)-JVS(751)*X(106)-JVS(752)*X(109)-JVS(753)& &*X(110)-JVS(754)*X(112)-JVS(755)*X(113) X(116) = X(116)-JVS(772)*X(80)-JVS(773)*X(103)-JVS(774)*X(107) X(117) = X(117)-JVS(785)*X(71)-JVS(786)*X(81) X(118) = X(118)-JVS(795)*X(63)-JVS(796)*X(73)-JVS(797)*X(76)-JVS(798)*X(83)-JVS(799)*X(84)-JVS(800)*X(92)-JVS(801)& &*X(93)-JVS(802)*X(95)-JVS(803)*X(96)-JVS(804)*X(97)-JVS(805)*X(101)-JVS(806)*X(103)-JVS(807)*X(104)-JVS(808)& &*X(105)-JVS(809)*X(107)-JVS(810)*X(108)-JVS(811)*X(110)-JVS(812)*X(111) X(119) = X(119)-JVS(831)*X(69)-JVS(832)*X(70)-JVS(833)*X(76)-JVS(834)*X(77)-JVS(835)*X(83)-JVS(836)*X(84)-JVS(837)& &*X(92)-JVS(838)*X(95)-JVS(839)*X(102)-JVS(840)*X(103)-JVS(841)*X(104)-JVS(842)*X(105)-JVS(843)*X(107)-JVS(844)& &*X(108)-JVS(845)*X(109)-JVS(846)*X(111) X(120) = X(120)-JVS(864)*X(74)-JVS(865)*X(81)-JVS(866)*X(117) X(121) = X(121)-JVS(876)*X(74)-JVS(877)*X(81)-JVS(878)*X(117) X(122) = X(122)-JVS(888)*X(64)-JVS(889)*X(79)-JVS(890)*X(120)-JVS(891)*X(121) X(123) = X(123)-JVS(902)*X(74)-JVS(903)*X(81) X(124) = X(124)-JVS(912)*X(74)-JVS(913)*X(82)-JVS(914)*X(117)-JVS(915)*X(123) X(125) = X(125)-JVS(925)*X(117)-JVS(926)*X(123) X(126) = X(126)-JVS(935)*X(51)-JVS(936)*X(66)-JVS(937)*X(87)-JVS(938)*X(111)-JVS(939)*X(116)-JVS(940)*X(123)-JVS(941)& &*X(124)-JVS(942)*X(125) X(127) = X(127)-JVS(952)*X(56)-JVS(953)*X(66)-JVS(954)*X(67)-JVS(955)*X(71)-JVS(956)*X(80)-JVS(957)*X(87)-JVS(958)& &*X(88)-JVS(959)*X(94)-JVS(960)*X(97)-JVS(961)*X(98)-JVS(962)*X(101)-JVS(963)*X(102)-JVS(964)*X(103)-JVS(965)& &*X(105)-JVS(966)*X(106)-JVS(967)*X(107)-JVS(968)*X(109)-JVS(969)*X(110)-JVS(970)*X(111)-JVS(971)*X(116)& &-JVS(972)*X(117)-JVS(973)*X(120)-JVS(974)*X(121)-JVS(975)*X(122)-JVS(976)*X(123)-JVS(977)*X(124)-JVS(978)& &*X(125)-JVS(979)*X(126) X(128) = X(128)-JVS(991)*X(45)-JVS(992)*X(49)-JVS(993)*X(50)-JVS(994)*X(53)-JVS(995)*X(55)-JVS(996)*X(59)-JVS(997)& &*X(65)-JVS(998)*X(67)-JVS(999)*X(68)-JVS(1000)*X(71)-JVS(1001)*X(75)-JVS(1002)*X(79)-JVS(1003)*X(80)-JVS(1004)& &*X(82)-JVS(1005)*X(83)-JVS(1006)*X(84)-JVS(1007)*X(87)-JVS(1008)*X(92)-JVS(1009)*X(93)-JVS(1010)*X(94)& &-JVS(1011)*X(95)-JVS(1012)*X(96)-JVS(1013)*X(97)-JVS(1014)*X(98)-JVS(1015)*X(99)-JVS(1016)*X(100)-JVS(1017)& &*X(101)-JVS(1018)*X(102)-JVS(1019)*X(103)-JVS(1020)*X(104)-JVS(1021)*X(105)-JVS(1022)*X(106)-JVS(1023)*X(107)& &-JVS(1024)*X(108)-JVS(1025)*X(109)-JVS(1026)*X(110)-JVS(1027)*X(111)-JVS(1028)*X(112)-JVS(1029)*X(113)& &-JVS(1030)*X(114)-JVS(1031)*X(115)-JVS(1032)*X(116)-JVS(1033)*X(117)-JVS(1034)*X(118)-JVS(1035)*X(119)& &-JVS(1036)*X(120)-JVS(1037)*X(121)-JVS(1038)*X(122)-JVS(1039)*X(123)-JVS(1040)*X(124)-JVS(1041)*X(125)& &-JVS(1042)*X(126)-JVS(1043)*X(127) X(129) = X(129)-JVS(1054)*X(49)-JVS(1055)*X(51)-JVS(1056)*X(52)-JVS(1057)*X(55)-JVS(1058)*X(59)-JVS(1059)*X(60)& &-JVS(1060)*X(65)-JVS(1061)*X(75)-JVS(1062)*X(81)-JVS(1063)*X(86)-JVS(1064)*X(90)-JVS(1065)*X(92)-JVS(1066)& &*X(93)-JVS(1067)*X(95)-JVS(1068)*X(99)-JVS(1069)*X(100)-JVS(1070)*X(101)-JVS(1071)*X(102)-JVS(1072)*X(103)& &-JVS(1073)*X(104)-JVS(1074)*X(105)-JVS(1075)*X(107)-JVS(1076)*X(111)-JVS(1077)*X(114)-JVS(1078)*X(115)& &-JVS(1079)*X(116)-JVS(1080)*X(117)-JVS(1081)*X(118)-JVS(1082)*X(119)-JVS(1083)*X(120)-JVS(1084)*X(121)& &-JVS(1085)*X(122)-JVS(1086)*X(123)-JVS(1087)*X(124)-JVS(1088)*X(125)-JVS(1089)*X(126)-JVS(1090)*X(127)& &-JVS(1091)*X(128) X(130) = X(130)-JVS(1101)*X(49)-JVS(1102)*X(54)-JVS(1103)*X(63)-JVS(1104)*X(64)-JVS(1105)*X(65)-JVS(1106)*X(73)& &-JVS(1107)*X(75)-JVS(1108)*X(85)-JVS(1109)*X(93)-JVS(1110)*X(94)-JVS(1111)*X(95)-JVS(1112)*X(97)-JVS(1113)& &*X(98)-JVS(1114)*X(101)-JVS(1115)*X(102)-JVS(1116)*X(105)-JVS(1117)*X(106)-JVS(1118)*X(107)-JVS(1119)*X(108)& &-JVS(1120)*X(109)-JVS(1121)*X(110)-JVS(1122)*X(111)-JVS(1123)*X(112)-JVS(1124)*X(113)-JVS(1125)*X(115)& &-JVS(1126)*X(116)-JVS(1127)*X(117)-JVS(1128)*X(118)-JVS(1129)*X(119)-JVS(1130)*X(120)-JVS(1131)*X(121)& &-JVS(1132)*X(122)-JVS(1133)*X(123)-JVS(1134)*X(124)-JVS(1135)*X(125)-JVS(1136)*X(126)-JVS(1137)*X(127)& &-JVS(1138)*X(128)-JVS(1139)*X(129) X(131) = X(131)-JVS(1148)*X(71)-JVS(1149)*X(73)-JVS(1150)*X(77)-JVS(1151)*X(78)-JVS(1152)*X(101)-JVS(1153)*X(102)& &-JVS(1154)*X(103)-JVS(1155)*X(105)-JVS(1156)*X(106)-JVS(1157)*X(107)-JVS(1158)*X(111)-JVS(1159)*X(117)& &-JVS(1160)*X(120)-JVS(1161)*X(121)-JVS(1162)*X(122)-JVS(1163)*X(123)-JVS(1164)*X(124)-JVS(1165)*X(125)& &-JVS(1166)*X(126)-JVS(1167)*X(128)-JVS(1168)*X(129)-JVS(1169)*X(130) X(132) = X(132)-JVS(1177)*X(44)-JVS(1178)*X(46)-JVS(1179)*X(47)-JVS(1180)*X(48)-JVS(1181)*X(51)-JVS(1182)*X(52)& &-JVS(1183)*X(53)-JVS(1184)*X(54)-JVS(1185)*X(56)-JVS(1186)*X(57)-JVS(1187)*X(58)-JVS(1188)*X(59)-JVS(1189)& &*X(60)-JVS(1190)*X(61)-JVS(1191)*X(62)-JVS(1192)*X(63)-JVS(1193)*X(64)-JVS(1194)*X(66)-JVS(1195)*X(67)& &-JVS(1196)*X(68)-JVS(1197)*X(69)-JVS(1198)*X(70)-JVS(1199)*X(71)-JVS(1200)*X(72)-JVS(1201)*X(73)-JVS(1202)& &*X(74)-JVS(1203)*X(75)-JVS(1204)*X(76)-JVS(1205)*X(77)-JVS(1206)*X(78)-JVS(1207)*X(79)-JVS(1208)*X(80)& &-JVS(1209)*X(81)-JVS(1210)*X(82)-JVS(1211)*X(83)-JVS(1212)*X(84)-JVS(1213)*X(85)-JVS(1214)*X(86)-JVS(1215)& &*X(87)-JVS(1216)*X(88)-JVS(1217)*X(89)-JVS(1218)*X(90)-JVS(1219)*X(91)-JVS(1220)*X(92)-JVS(1221)*X(93)& &-JVS(1222)*X(94)-JVS(1223)*X(95)-JVS(1224)*X(96)-JVS(1225)*X(97)-JVS(1226)*X(98)-JVS(1227)*X(99)-JVS(1228)& &*X(100)-JVS(1229)*X(101)-JVS(1230)*X(102)-JVS(1231)*X(103)-JVS(1232)*X(104)-JVS(1233)*X(105)-JVS(1234)*X(106)& &-JVS(1235)*X(107)-JVS(1236)*X(108)-JVS(1237)*X(109)-JVS(1238)*X(110)-JVS(1239)*X(111)-JVS(1240)*X(112)& &-JVS(1241)*X(113)-JVS(1242)*X(114)-JVS(1243)*X(115)-JVS(1244)*X(116)-JVS(1245)*X(117)-JVS(1246)*X(118)& &-JVS(1247)*X(119)-JVS(1248)*X(120)-JVS(1249)*X(121)-JVS(1250)*X(122)-JVS(1251)*X(123)-JVS(1252)*X(124)& &-JVS(1253)*X(125)-JVS(1254)*X(126)-JVS(1255)*X(127)-JVS(1256)*X(128)-JVS(1257)*X(129)-JVS(1258)*X(130)& &-JVS(1259)*X(131) X(133) = X(133)-JVS(1266)*X(49)-JVS(1267)*X(54)-JVS(1268)*X(60)-JVS(1269)*X(62)-JVS(1270)*X(85)-JVS(1271)*X(86)& &-JVS(1272)*X(89)-JVS(1273)*X(92)-JVS(1274)*X(94)-JVS(1275)*X(95)-JVS(1276)*X(96)-JVS(1277)*X(97)-JVS(1278)& &*X(98)-JVS(1279)*X(101)-JVS(1280)*X(102)-JVS(1281)*X(103)-JVS(1282)*X(104)-JVS(1283)*X(105)-JVS(1284)*X(106)& &-JVS(1285)*X(107)-JVS(1286)*X(108)-JVS(1287)*X(109)-JVS(1288)*X(110)-JVS(1289)*X(111)-JVS(1290)*X(112)& &-JVS(1291)*X(113)-JVS(1292)*X(115)-JVS(1293)*X(116)-JVS(1294)*X(117)-JVS(1295)*X(118)-JVS(1296)*X(119)& &-JVS(1297)*X(120)-JVS(1298)*X(121)-JVS(1299)*X(122)-JVS(1300)*X(123)-JVS(1301)*X(124)-JVS(1302)*X(125)& &-JVS(1303)*X(126)-JVS(1304)*X(127)-JVS(1305)*X(128)-JVS(1306)*X(129)-JVS(1307)*X(130)-JVS(1308)*X(131)& &-JVS(1309)*X(132) X(134) = X(134)-JVS(1315)*X(53)-JVS(1316)*X(55)-JVS(1317)*X(75)-JVS(1318)*X(83)-JVS(1319)*X(84)-JVS(1320)*X(87)& &-JVS(1321)*X(92)-JVS(1322)*X(94)-JVS(1323)*X(95)-JVS(1324)*X(97)-JVS(1325)*X(98)-JVS(1326)*X(101)-JVS(1327)& &*X(102)-JVS(1328)*X(103)-JVS(1329)*X(104)-JVS(1330)*X(105)-JVS(1331)*X(106)-JVS(1332)*X(107)-JVS(1333)*X(108)& &-JVS(1334)*X(109)-JVS(1335)*X(110)-JVS(1336)*X(111)-JVS(1337)*X(112)-JVS(1338)*X(113)-JVS(1339)*X(116)& &-JVS(1340)*X(117)-JVS(1341)*X(120)-JVS(1342)*X(121)-JVS(1343)*X(122)-JVS(1344)*X(123)-JVS(1345)*X(124)& &-JVS(1346)*X(125)-JVS(1347)*X(126)-JVS(1348)*X(127)-JVS(1349)*X(128)-JVS(1350)*X(129)-JVS(1351)*X(130)& &-JVS(1352)*X(131)-JVS(1353)*X(132)-JVS(1354)*X(133) X(135) = X(135)-JVS(1359)*X(55)-JVS(1360)*X(78)-JVS(1361)*X(79)-JVS(1362)*X(81)-JVS(1363)*X(82)-JVS(1364)*X(86)& &-JVS(1365)*X(111)-JVS(1366)*X(113)-JVS(1367)*X(120)-JVS(1368)*X(121)-JVS(1369)*X(122)-JVS(1370)*X(123)& &-JVS(1371)*X(124)-JVS(1372)*X(125)-JVS(1373)*X(127)-JVS(1374)*X(128)-JVS(1375)*X(129)-JVS(1376)*X(130)& &-JVS(1377)*X(131)-JVS(1378)*X(132)-JVS(1379)*X(133)-JVS(1380)*X(134) X(136) = X(136)-JVS(1384)*X(46)-JVS(1385)*X(48)-JVS(1386)*X(56)-JVS(1387)*X(57)-JVS(1388)*X(58)-JVS(1389)*X(59)& &-JVS(1390)*X(60)-JVS(1391)*X(61)-JVS(1392)*X(62)-JVS(1393)*X(66)-JVS(1394)*X(67)-JVS(1395)*X(69)-JVS(1396)& &*X(70)-JVS(1397)*X(71)-JVS(1398)*X(72)-JVS(1399)*X(73)-JVS(1400)*X(74)-JVS(1401)*X(75)-JVS(1402)*X(76)& &-JVS(1403)*X(77)-JVS(1404)*X(78)-JVS(1405)*X(79)-JVS(1406)*X(80)-JVS(1407)*X(81)-JVS(1408)*X(82)-JVS(1409)& &*X(83)-JVS(1410)*X(84)-JVS(1411)*X(86)-JVS(1412)*X(87)-JVS(1413)*X(88)-JVS(1414)*X(90)-JVS(1415)*X(91)& &-JVS(1416)*X(92)-JVS(1417)*X(93)-JVS(1418)*X(94)-JVS(1419)*X(95)-JVS(1420)*X(96)-JVS(1421)*X(97)-JVS(1422)& &*X(98)-JVS(1423)*X(100)-JVS(1424)*X(101)-JVS(1425)*X(102)-JVS(1426)*X(103)-JVS(1427)*X(104)-JVS(1428)*X(105)& &-JVS(1429)*X(106)-JVS(1430)*X(107)-JVS(1431)*X(108)-JVS(1432)*X(109)-JVS(1433)*X(110)-JVS(1434)*X(111)& &-JVS(1435)*X(112)-JVS(1436)*X(113)-JVS(1437)*X(114)-JVS(1438)*X(115)-JVS(1439)*X(116)-JVS(1440)*X(117)& &-JVS(1441)*X(118)-JVS(1442)*X(119)-JVS(1443)*X(120)-JVS(1444)*X(121)-JVS(1445)*X(122)-JVS(1446)*X(123)& &-JVS(1447)*X(124)-JVS(1448)*X(125)-JVS(1449)*X(126)-JVS(1450)*X(127)-JVS(1451)*X(128)-JVS(1452)*X(129)& &-JVS(1453)*X(130)-JVS(1454)*X(131)-JVS(1455)*X(132)-JVS(1456)*X(133)-JVS(1457)*X(134)-JVS(1458)*X(135) X(137) = X(137)-JVS(1461)*X(71)-JVS(1462)*X(97)-JVS(1463)*X(106)-JVS(1464)*X(108)-JVS(1465)*X(116)-JVS(1466)*X(117)& &-JVS(1467)*X(121)-JVS(1468)*X(122)-JVS(1469)*X(123)-JVS(1470)*X(125)-JVS(1471)*X(126)-JVS(1472)*X(128)& &-JVS(1473)*X(129)-JVS(1474)*X(130)-JVS(1475)*X(131)-JVS(1476)*X(132)-JVS(1477)*X(133)-JVS(1478)*X(134)& &-JVS(1479)*X(135)-JVS(1480)*X(136) X(137) = X(137)/JVS(1481) X(136) = (X(136)-JVS(1460)*X(137))/(JVS(1459)) X(135) = (X(135)-JVS(1382)*X(136)-JVS(1383)*X(137))/(JVS(1381)) X(134) = (X(134)-JVS(1356)*X(135)-JVS(1357)*X(136)-JVS(1358)*X(137))/(JVS(1355)) X(133) = (X(133)-JVS(1311)*X(134)-JVS(1312)*X(135)-JVS(1313)*X(136)-JVS(1314)*X(137))/(JVS(1310)) X(132) = (X(132)-JVS(1261)*X(133)-JVS(1262)*X(134)-JVS(1263)*X(135)-JVS(1264)*X(136)-JVS(1265)*X(137))/(JVS(1260)) X(131) = (X(131)-JVS(1171)*X(132)-JVS(1172)*X(133)-JVS(1173)*X(134)-JVS(1174)*X(135)-JVS(1175)*X(136)-JVS(1176)& &*X(137))/(JVS(1170)) X(130) = (X(130)-JVS(1141)*X(131)-JVS(1142)*X(132)-JVS(1143)*X(133)-JVS(1144)*X(134)-JVS(1145)*X(135)-JVS(1146)*X(136)& &-JVS(1147)*X(137))/(JVS(1140)) X(129) = (X(129)-JVS(1093)*X(130)-JVS(1094)*X(131)-JVS(1095)*X(132)-JVS(1096)*X(133)-JVS(1097)*X(134)-JVS(1098)*X(135)& &-JVS(1099)*X(136)-JVS(1100)*X(137))/(JVS(1092)) X(128) = (X(128)-JVS(1045)*X(129)-JVS(1046)*X(130)-JVS(1047)*X(131)-JVS(1048)*X(132)-JVS(1049)*X(133)-JVS(1050)*X(134)& &-JVS(1051)*X(135)-JVS(1052)*X(136)-JVS(1053)*X(137))/(JVS(1044)) X(127) = (X(127)-JVS(981)*X(128)-JVS(982)*X(129)-JVS(983)*X(130)-JVS(984)*X(131)-JVS(985)*X(132)-JVS(986)*X(133)& &-JVS(987)*X(134)-JVS(988)*X(135)-JVS(989)*X(136)-JVS(990)*X(137))/(JVS(980)) X(126) = (X(126)-JVS(944)*X(129)-JVS(945)*X(130)-JVS(946)*X(132)-JVS(947)*X(133)-JVS(948)*X(134)-JVS(949)*X(135)& &-JVS(950)*X(136)-JVS(951)*X(137))/(JVS(943)) X(125) = (X(125)-JVS(928)*X(129)-JVS(929)*X(130)-JVS(930)*X(132)-JVS(931)*X(133)-JVS(932)*X(134)-JVS(933)*X(135)& &-JVS(934)*X(136))/(JVS(927)) X(124) = (X(124)-JVS(917)*X(125)-JVS(918)*X(129)-JVS(919)*X(130)-JVS(920)*X(132)-JVS(921)*X(133)-JVS(922)*X(134)& &-JVS(923)*X(135)-JVS(924)*X(136))/(JVS(916)) X(123) = (X(123)-JVS(905)*X(129)-JVS(906)*X(130)-JVS(907)*X(132)-JVS(908)*X(133)-JVS(909)*X(134)-JVS(910)*X(135)& &-JVS(911)*X(136))/(JVS(904)) X(122) = (X(122)-JVS(893)*X(123)-JVS(894)*X(128)-JVS(895)*X(129)-JVS(896)*X(130)-JVS(897)*X(132)-JVS(898)*X(133)& &-JVS(899)*X(134)-JVS(900)*X(135)-JVS(901)*X(136))/(JVS(892)) X(121) = (X(121)-JVS(880)*X(123)-JVS(881)*X(129)-JVS(882)*X(130)-JVS(883)*X(132)-JVS(884)*X(133)-JVS(885)*X(134)& &-JVS(886)*X(135)-JVS(887)*X(136))/(JVS(879)) X(120) = (X(120)-JVS(868)*X(123)-JVS(869)*X(129)-JVS(870)*X(130)-JVS(871)*X(132)-JVS(872)*X(133)-JVS(873)*X(134)& &-JVS(874)*X(135)-JVS(875)*X(136))/(JVS(867)) X(119) = (X(119)-JVS(848)*X(120)-JVS(849)*X(121)-JVS(850)*X(122)-JVS(851)*X(123)-JVS(852)*X(124)-JVS(853)*X(125)& &-JVS(854)*X(126)-JVS(855)*X(128)-JVS(856)*X(129)-JVS(857)*X(130)-JVS(858)*X(132)-JVS(859)*X(133)-JVS(860)& &*X(134)-JVS(861)*X(135)-JVS(862)*X(136)-JVS(863)*X(137))/(JVS(847)) X(118) = (X(118)-JVS(814)*X(119)-JVS(815)*X(120)-JVS(816)*X(121)-JVS(817)*X(122)-JVS(818)*X(123)-JVS(819)*X(124)& &-JVS(820)*X(125)-JVS(821)*X(126)-JVS(822)*X(128)-JVS(823)*X(129)-JVS(824)*X(130)-JVS(825)*X(132)-JVS(826)& &*X(133)-JVS(827)*X(134)-JVS(828)*X(135)-JVS(829)*X(136)-JVS(830)*X(137))/(JVS(813)) X(117) = (X(117)-JVS(788)*X(129)-JVS(789)*X(130)-JVS(790)*X(132)-JVS(791)*X(133)-JVS(792)*X(134)-JVS(793)*X(135)& &-JVS(794)*X(136))/(JVS(787)) X(116) = (X(116)-JVS(776)*X(123)-JVS(777)*X(125)-JVS(778)*X(130)-JVS(779)*X(132)-JVS(780)*X(133)-JVS(781)*X(134)& &-JVS(782)*X(135)-JVS(783)*X(136)-JVS(784)*X(137))/(JVS(775)) X(115) = (X(115)-JVS(757)*X(116)-JVS(758)*X(118)-JVS(759)*X(120)-JVS(760)*X(126)-JVS(761)*X(127)-JVS(762)*X(128)& &-JVS(763)*X(129)-JVS(764)*X(130)-JVS(765)*X(131)-JVS(766)*X(132)-JVS(767)*X(133)-JVS(768)*X(134)-JVS(769)& &*X(135)-JVS(770)*X(136)-JVS(771)*X(137))/(JVS(756)) X(114) = (X(114)-JVS(725)*X(115)-JVS(726)*X(116)-JVS(727)*X(117)-JVS(728)*X(119)-JVS(729)*X(120)-JVS(730)*X(121)& &-JVS(731)*X(122)-JVS(732)*X(123)-JVS(733)*X(124)-JVS(734)*X(125)-JVS(735)*X(126)-JVS(736)*X(127)-JVS(737)& &*X(128)-JVS(738)*X(129)-JVS(739)*X(130)-JVS(740)*X(131)-JVS(741)*X(132)-JVS(742)*X(133)-JVS(743)*X(134)& &-JVS(744)*X(135)-JVS(745)*X(136)-JVS(746)*X(137))/(JVS(724)) X(113) = (X(113)-JVS(679)*X(127)-JVS(680)*X(128)-JVS(681)*X(129)-JVS(682)*X(130)-JVS(683)*X(131)-JVS(684)*X(132)& &-JVS(685)*X(133)-JVS(686)*X(134)-JVS(687)*X(136))/(JVS(678)) X(112) = (X(112)-JVS(665)*X(113)-JVS(666)*X(126)-JVS(667)*X(127)-JVS(668)*X(129)-JVS(669)*X(130)-JVS(670)*X(131)& &-JVS(671)*X(132)-JVS(672)*X(133)-JVS(673)*X(134)-JVS(674)*X(136)-JVS(675)*X(137))/(JVS(664)) X(111) = (X(111)-JVS(656)*X(124)-JVS(657)*X(130)-JVS(658)*X(132)-JVS(659)*X(133)-JVS(660)*X(134)-JVS(661)*X(136))& &/(JVS(655)) X(110) = (X(110)-JVS(646)*X(120)-JVS(647)*X(129)-JVS(648)*X(130)-JVS(649)*X(132)-JVS(650)*X(133)-JVS(651)*X(134)& &-JVS(652)*X(135)-JVS(653)*X(136))/(JVS(645)) X(109) = (X(109)-JVS(635)*X(120)-JVS(636)*X(129)-JVS(637)*X(130)-JVS(638)*X(132)-JVS(639)*X(133)-JVS(640)*X(134)& &-JVS(641)*X(135)-JVS(642)*X(136))/(JVS(634)) X(108) = (X(108)-JVS(624)*X(122)-JVS(625)*X(126)-JVS(626)*X(130)-JVS(627)*X(132)-JVS(628)*X(133)-JVS(629)*X(134)& &-JVS(630)*X(136)-JVS(631)*X(137))/(JVS(623)) X(107) = (X(107)-JVS(614)*X(123)-JVS(615)*X(125)-JVS(616)*X(132)-JVS(617)*X(134)-JVS(618)*X(135)-JVS(619)*X(136))& &/(JVS(613)) X(106) = (X(106)-JVS(603)*X(129)-JVS(604)*X(130)-JVS(605)*X(131)-JVS(606)*X(132)-JVS(607)*X(133)-JVS(608)*X(134)& &-JVS(609)*X(136))/(JVS(602)) X(105) = (X(105)-JVS(596)*X(121)-JVS(597)*X(130)-JVS(598)*X(132)-JVS(599)*X(133)-JVS(600)*X(134)-JVS(601)*X(136))& &/(JVS(595)) X(104) = (X(104)-JVS(586)*X(105)-JVS(587)*X(107)-JVS(588)*X(123)-JVS(589)*X(125)-JVS(590)*X(132)-JVS(591)*X(134)& &-JVS(592)*X(135)-JVS(593)*X(136))/(JVS(585)) X(103) = (X(103)-JVS(578)*X(107)-JVS(579)*X(125)-JVS(580)*X(132)-JVS(581)*X(134)-JVS(582)*X(136))/(JVS(577)) X(102) = (X(102)-JVS(570)*X(122)-JVS(571)*X(130)-JVS(572)*X(132)-JVS(573)*X(133)-JVS(574)*X(134)-JVS(575)*X(136))& &/(JVS(569)) X(101) = (X(101)-JVS(562)*X(120)-JVS(563)*X(130)-JVS(564)*X(132)-JVS(565)*X(133)-JVS(566)*X(134)-JVS(567)*X(136))& &/(JVS(561)) X(100) = (X(100)-JVS(547)*X(101)-JVS(548)*X(103)-JVS(549)*X(107)-JVS(550)*X(111)-JVS(551)*X(123)-JVS(552)*X(124)& &-JVS(553)*X(125)-JVS(554)*X(130)-JVS(555)*X(132)-JVS(556)*X(133)-JVS(557)*X(134)-JVS(558)*X(135)-JVS(559)& &*X(136))/(JVS(546)) X(99) = (X(99)-JVS(520)*X(100)-JVS(521)*X(102)-JVS(522)*X(103)-JVS(523)*X(107)-JVS(524)*X(111)-JVS(525)*X(114)& &-JVS(526)*X(115)-JVS(527)*X(117)-JVS(528)*X(118)-JVS(529)*X(121)-JVS(530)*X(123)-JVS(531)*X(127)-JVS(532)& &*X(128)-JVS(533)*X(129)-JVS(534)*X(130)-JVS(535)*X(131)-JVS(536)*X(132)-JVS(537)*X(133)-JVS(538)*X(134)& &-JVS(539)*X(135)-JVS(540)*X(136))/(JVS(519)) X(98) = (X(98)-JVS(505)*X(126)-JVS(506)*X(130)-JVS(507)*X(132)-JVS(508)*X(133)-JVS(509)*X(134)-JVS(510)*X(136)& &-JVS(511)*X(137))/(JVS(504)) X(97) = (X(97)-JVS(496)*X(121)-JVS(497)*X(129)-JVS(498)*X(130)-JVS(499)*X(133)-JVS(500)*X(134)-JVS(501)*X(136))& &/(JVS(495)) X(96) = (X(96)-JVS(485)*X(97)-JVS(486)*X(103)-JVS(487)*X(110)-JVS(488)*X(125)-JVS(489)*X(130)-JVS(490)*X(132)-JVS(491)& &*X(133)-JVS(492)*X(134)-JVS(493)*X(135)-JVS(494)*X(136))/(JVS(484)) X(95) = (X(95)-JVS(477)*X(122)-JVS(478)*X(128)-JVS(479)*X(132)-JVS(480)*X(134)-JVS(481)*X(135)-JVS(482)*X(136))& &/(JVS(476)) X(94) = (X(94)-JVS(468)*X(126)-JVS(469)*X(130)-JVS(470)*X(132)-JVS(471)*X(133)-JVS(472)*X(134)-JVS(473)*X(136)& &-JVS(474)*X(137))/(JVS(467)) X(93) = (X(93)-JVS(458)*X(101)-JVS(459)*X(107)-JVS(460)*X(123)-JVS(461)*X(132)-JVS(462)*X(134)-JVS(463)*X(135)& &-JVS(464)*X(136))/(JVS(457)) X(92) = (X(92)-JVS(450)*X(95)-JVS(451)*X(104)-JVS(452)*X(128)-JVS(453)*X(132)-JVS(454)*X(134)-JVS(455)*X(136))& &/(JVS(449)) X(91) = (X(91)-JVS(422)*X(96)-JVS(423)*X(97)-JVS(424)*X(98)-JVS(425)*X(100)-JVS(426)*X(104)-JVS(427)*X(105)-JVS(428)& &*X(111)-JVS(429)*X(114)-JVS(430)*X(115)-JVS(431)*X(118)-JVS(432)*X(119)-JVS(433)*X(120)-JVS(434)*X(121)& &-JVS(435)*X(122)-JVS(436)*X(123)-JVS(437)*X(124)-JVS(438)*X(126)-JVS(439)*X(127)-JVS(440)*X(129)-JVS(441)& &*X(130)-JVS(442)*X(132)-JVS(443)*X(133)-JVS(444)*X(134)-JVS(445)*X(135)-JVS(446)*X(136)-JVS(447)*X(137))& &/(JVS(421)) X(90) = (X(90)-JVS(405)*X(100)-JVS(406)*X(111)-JVS(407)*X(123)-JVS(408)*X(129)-JVS(409)*X(130)-JVS(410)*X(132)& &-JVS(411)*X(134)-JVS(412)*X(136))/(JVS(404)) X(89) = (X(89)-JVS(379)*X(92)-JVS(380)*X(94)-JVS(381)*X(97)-JVS(382)*X(98)-JVS(383)*X(101)-JVS(384)*X(102)-JVS(385)& &*X(105)-JVS(386)*X(106)-JVS(387)*X(108)-JVS(388)*X(109)-JVS(389)*X(110)-JVS(390)*X(111)-JVS(391)*X(112)& &-JVS(392)*X(116)-JVS(393)*X(117)-JVS(394)*X(119)-JVS(395)*X(123)-JVS(396)*X(126)-JVS(397)*X(130)-JVS(398)& &*X(132)-JVS(399)*X(133)-JVS(400)*X(134)-JVS(401)*X(136))/(JVS(378)) X(88) = (X(88)-JVS(363)*X(94)-JVS(364)*X(98)-JVS(365)*X(101)-JVS(366)*X(102)-JVS(367)*X(105)-JVS(368)*X(106)-JVS(369)& &*X(109)-JVS(370)*X(110)-JVS(371)*X(111)-JVS(372)*X(116)-JVS(373)*X(117)-JVS(374)*X(123)-JVS(375)*X(126)& &-JVS(376)*X(132)-JVS(377)*X(133))/(JVS(362)) X(87) = (X(87)-JVS(356)*X(111)-JVS(357)*X(130)-JVS(358)*X(132)-JVS(359)*X(134)-JVS(360)*X(135)-JVS(361)*X(136))& &/(JVS(355)) X(86) = (X(86)-JVS(350)*X(120)-JVS(351)*X(129)-JVS(352)*X(132)-JVS(353)*X(135))/(JVS(349)) X(85) = (X(85)-JVS(340)*X(98)-JVS(341)*X(126)-JVS(342)*X(130)-JVS(343)*X(132)-JVS(344)*X(133)-JVS(345)*X(134)-JVS(346)& &*X(136)-JVS(347)*X(137))/(JVS(339)) X(84) = (X(84)-JVS(334)*X(123)-JVS(335)*X(132)-JVS(336)*X(134)-JVS(337)*X(136))/(JVS(333)) X(83) = (X(83)-JVS(327)*X(123)-JVS(328)*X(130)-JVS(329)*X(134)-JVS(330)*X(136))/(JVS(326)) X(82) = (X(82)-JVS(322)*X(123)-JVS(323)*X(132)-JVS(324)*X(134)-JVS(325)*X(135))/(JVS(321)) X(81) = (X(81)-JVS(318)*X(129)-JVS(319)*X(132)-JVS(320)*X(135))/(JVS(317)) X(80) = (X(80)-JVS(313)*X(103)-JVS(314)*X(107)-JVS(315)*X(132)-JVS(316)*X(136))/(JVS(312)) X(79) = (X(79)-JVS(308)*X(122)-JVS(309)*X(128)-JVS(310)*X(132)-JVS(311)*X(135))/(JVS(307)) X(78) = (X(78)-JVS(302)*X(111)-JVS(303)*X(130)-JVS(304)*X(132)-JVS(305)*X(134)-JVS(306)*X(135))/(JVS(301)) X(77) = (X(77)-JVS(296)*X(111)-JVS(297)*X(130)-JVS(298)*X(132)-JVS(299)*X(134)-JVS(300)*X(136))/(JVS(295)) X(76) = (X(76)-JVS(291)*X(111)-JVS(292)*X(132)-JVS(293)*X(136))/(JVS(290)) X(75) = (X(75)-JVS(286)*X(95)-JVS(287)*X(128)-JVS(288)*X(134)-JVS(289)*X(136))/(JVS(285)) X(74) = (X(74)-JVS(281)*X(123)-JVS(282)*X(132)-JVS(283)*X(136))/(JVS(280)) X(73) = (X(73)-JVS(277)*X(101)-JVS(278)*X(132)-JVS(279)*X(136))/(JVS(276)) X(72) = (X(72)-JVS(265)*X(84)-JVS(266)*X(92)-JVS(267)*X(93)-JVS(268)*X(100)-JVS(269)*X(103)-JVS(270)*X(119)-JVS(271)& &*X(124)-JVS(272)*X(132)-JVS(273)*X(134)-JVS(274)*X(135)-JVS(275)*X(136))/(JVS(264)) X(71) = (X(71)-JVS(261)*X(117)-JVS(262)*X(132)-JVS(263)*X(136))/(JVS(260)) X(70) = (X(70)-JVS(257)*X(105)-JVS(258)*X(132)-JVS(259)*X(136))/(JVS(256)) X(69) = (X(69)-JVS(253)*X(109)-JVS(254)*X(132)-JVS(255)*X(136))/(JVS(252)) X(68) = (X(68)-JVS(246)*X(103)-JVS(247)*X(125)-JVS(248)*X(132)-JVS(249)*X(134)-JVS(250)*X(135)-JVS(251)*X(136))& &/(JVS(245)) X(67) = (X(67)-JVS(242)*X(110)-JVS(243)*X(132)-JVS(244)*X(136))/(JVS(241)) X(66) = (X(66)-JVS(237)*X(87)-JVS(238)*X(126)-JVS(239)*X(132)-JVS(240)*X(136))/(JVS(236)) X(65) = (X(65)-JVS(232)*X(75)-JVS(233)*X(95)-JVS(234)*X(128)-JVS(235)*X(134))/(JVS(231)) X(64) = (X(64)-JVS(228)*X(122)-JVS(229)*X(132)-JVS(230)*X(136))/(JVS(227)) X(63) = (X(63)-JVS(224)*X(108)-JVS(225)*X(132)-JVS(226)*X(136))/(JVS(223)) X(62) = (X(62)-JVS(220)*X(132)-JVS(221)*X(133)-JVS(222)*X(136))/(JVS(219)) X(61) = (X(61)-JVS(216)*X(112)-JVS(217)*X(132)-JVS(218)*X(136))/(JVS(215)) X(60) = (X(60)-JVS(213)*X(129)-JVS(214)*X(132))/(JVS(212)) X(59) = (X(59)-JVS(209)*X(128)-JVS(210)*X(132)-JVS(211)*X(136))/(JVS(208)) X(58) = (X(58)-JVS(205)*X(113)-JVS(206)*X(132)-JVS(207)*X(136))/(JVS(204)) X(57) = (X(57)-JVS(201)*X(98)-JVS(202)*X(132)-JVS(203)*X(136))/(JVS(200)) X(56) = (X(56)-JVS(197)*X(94)-JVS(198)*X(132)-JVS(199)*X(136))/(JVS(196)) X(55) = (X(55)-JVS(194)*X(128)-JVS(195)*X(129))/(JVS(193)) X(54) = (X(54)-JVS(190)*X(130)-JVS(191)*X(132)-JVS(192)*X(136))/(JVS(189)) X(53) = (X(53)-JVS(186)*X(128)-JVS(187)*X(132)-JVS(188)*X(134))/(JVS(185)) X(52) = (X(52)-JVS(183)*X(129)-JVS(184)*X(132))/(JVS(182)) X(51) = (X(51)-JVS(180)*X(129)-JVS(181)*X(132))/(JVS(179)) X(50) = (X(50)-JVS(175)*X(92)-JVS(176)*X(95)-JVS(177)*X(128)-JVS(178)*X(134))/(JVS(174)) X(49) = (X(49)-JVS(172)*X(128)-JVS(173)*X(130))/(JVS(171)) X(48) = (X(48)-JVS(168)*X(60)-JVS(169)*X(129)-JVS(170)*X(132))/(JVS(167)) X(47) = (X(47)-JVS(165)*X(74)-JVS(166)*X(132))/(JVS(164)) X(46) = (X(46)-JVS(162)*X(132)-JVS(163)*X(136))/(JVS(161)) X(45) = (X(45)-JVS(159)*X(113)-JVS(160)*X(128))/(JVS(158)) X(44) = (X(44)-JVS(157)*X(132))/(JVS(156)) X(43) = (X(43)-JVS(135)*X(45)-JVS(136)*X(46)-JVS(137)*X(47)-JVS(138)*X(49)-JVS(139)*X(50)-JVS(140)*X(55)-JVS(141)& &*X(65)-JVS(142)*X(72)-JVS(143)*X(74)-JVS(144)*X(79)-JVS(145)*X(82)-JVS(146)*X(89)-JVS(147)*X(93)-JVS(148)*X(96)& &-JVS(149)*X(99)-JVS(150)*X(104)-JVS(151)*X(114)-JVS(152)*X(125)-JVS(153)*X(128)-JVS(154)*X(135)-JVS(155)& &*X(137))/(JVS(134)) X(42) = (X(42)-JVS(108)*X(65)-JVS(109)*X(68)-JVS(110)*X(70)-JVS(111)*X(72)-JVS(112)*X(75)-JVS(113)*X(78)-JVS(114)& &*X(81)-JVS(115)*X(84)-JVS(116)*X(87)-JVS(117)*X(89)-JVS(118)*X(91)-JVS(119)*X(92)-JVS(120)*X(95)-JVS(121)*X(97)& &-JVS(122)*X(100)-JVS(123)*X(102)-JVS(124)*X(119)-JVS(125)*X(121)-JVS(126)*X(122)-JVS(127)*X(123)-JVS(128)& &*X(130)-JVS(129)*X(132)-JVS(130)*X(133)-JVS(131)*X(134)-JVS(132)*X(135)-JVS(133)*X(136))/(JVS(107)) X(41) = (X(41)-JVS(106)*X(65))/(JVS(105)) X(40) = (X(40)-JVS(104)*X(89))/(JVS(103)) X(39) = (X(39)-JVS(102)*X(72))/(JVS(101)) X(38) = (X(38)-JVS(100)*X(96))/(JVS(99)) X(37) = (X(37)-JVS(98)*X(93))/(JVS(97)) X(36) = (X(36)-JVS(96)*X(104))/(JVS(95)) X(35) = (X(35)-JVS(94)*X(47))/(JVS(93)) X(34) = (X(34)-JVS(92)*X(74))/(JVS(91)) X(33) = (X(33)-JVS(90)*X(82))/(JVS(89)) X(32) = (X(32)-JVS(88)*X(125))/(JVS(87)) X(31) = (X(31)-JVS(85)*X(93)-JVS(86)*X(132))/(JVS(84)) X(30) = (X(30)-JVS(83)*X(132))/(JVS(82)) X(29) = (X(29)-JVS(78)*X(30)-JVS(79)*X(132)-JVS(80)*X(134)-JVS(81)*X(136))/(JVS(77)) X(28) = (X(28)-JVS(75)*X(29)-JVS(76)*X(134))/(JVS(74)) X(27) = (X(27)-JVS(72)*X(29)-JVS(73)*X(136))/(JVS(71)) X(26) = (X(26)-JVS(70)*X(132))/(JVS(69)) X(25) = (X(25)-JVS(65)*X(26)-JVS(66)*X(132)-JVS(67)*X(134)-JVS(68)*X(136))/(JVS(64)) X(24) = (X(24)-JVS(62)*X(25)-JVS(63)*X(134))/(JVS(61)) X(23) = (X(23)-JVS(59)*X(25)-JVS(60)*X(136))/(JVS(58)) X(22) = (X(22)-JVS(57)*X(132))/(JVS(56)) X(21) = (X(21)-JVS(52)*X(22)-JVS(53)*X(132)-JVS(54)*X(134)-JVS(55)*X(136))/(JVS(51)) X(20) = (X(20)-JVS(49)*X(21)-JVS(50)*X(134))/(JVS(48)) X(19) = (X(19)-JVS(46)*X(21)-JVS(47)*X(136))/(JVS(45)) X(18) = (X(18)-JVS(44)*X(132))/(JVS(43)) X(17) = (X(17)-JVS(39)*X(18)-JVS(40)*X(132)-JVS(41)*X(134)-JVS(42)*X(136))/(JVS(38)) X(16) = (X(16)-JVS(36)*X(17)-JVS(37)*X(134))/(JVS(35)) X(15) = (X(15)-JVS(33)*X(17)-JVS(34)*X(136))/(JVS(32)) X(14) = (X(14)-JVS(30)*X(81)-JVS(31)*X(132))/(JVS(29)) X(13) = (X(13)-JVS(27)*X(60)-JVS(28)*X(132))/(JVS(26)) X(12) = (X(12)-JVS(24)*X(48)-JVS(25)*X(132))/(JVS(23)) X(11) = (X(11)-JVS(22)*X(137))/(JVS(21)) X(10) = (X(10)-JVS(20)*X(45))/(JVS(19)) X(9) = (X(9)-JVS(18)*X(50))/(JVS(17)) X(8) = (X(8)-JVS(16)*X(79))/(JVS(15)) X(7) = (X(7)-JVS(14)*X(49))/(JVS(13)) X(6) = (X(6)-JVS(12)*X(135))/(JVS(11)) X(5) = (X(5)-JVS(10)*X(128))/(JVS(9)) X(4) = (X(4)-JVS(8)*X(55))/(JVS(7)) X(3) = (X(3)-JVS(6)*X(99))/(JVS(5)) X(2) = (X(2)-JVS(4)*X(46))/(JVS(3)) X(1) = (X(1)-JVS(2)*X(114))/(JVS(1)) END SUBROUTINE KppSolve ! End of KppSolve function ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! KppSolveTR - sparse, transposed back substitution ! Arguments : ! JVS - sparse Jacobian of variables ! X - Vector for variables ! XX - Vector for output variables ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppSolveTR ( JVS, X, XX ) ! JVS - sparse Jacobian of variables REAL(kind=dp) :: JVS(LU_NONZERO) ! X - Vector for variables REAL(kind=dp) :: X(NVAR) ! XX - Vector for output variables REAL(kind=dp) :: XX(NVAR) XX(1) = X(1)/JVS(1) XX(2) = X(2)/JVS(3) XX(3) = X(3)/JVS(5) XX(4) = X(4)/JVS(7) XX(5) = X(5)/JVS(9) XX(6) = X(6)/JVS(11) XX(7) = X(7)/JVS(13) XX(8) = X(8)/JVS(15) XX(9) = X(9)/JVS(17) XX(10) = X(10)/JVS(19) XX(11) = X(11)/JVS(21) XX(12) = X(12)/JVS(23) XX(13) = X(13)/JVS(26) XX(14) = X(14)/JVS(29) XX(15) = X(15)/JVS(32) XX(16) = X(16)/JVS(35) XX(17) = (X(17)-JVS(33)*XX(15)-JVS(36)*XX(16))/(JVS(38)) XX(18) = (X(18)-JVS(39)*XX(17))/(JVS(43)) XX(19) = X(19)/JVS(45) XX(20) = X(20)/JVS(48) XX(21) = (X(21)-JVS(46)*XX(19)-JVS(49)*XX(20))/(JVS(51)) XX(22) = (X(22)-JVS(52)*XX(21))/(JVS(56)) XX(23) = X(23)/JVS(58) XX(24) = X(24)/JVS(61) XX(25) = (X(25)-JVS(59)*XX(23)-JVS(62)*XX(24))/(JVS(64)) XX(26) = (X(26)-JVS(65)*XX(25))/(JVS(69)) XX(27) = X(27)/JVS(71) XX(28) = X(28)/JVS(74) XX(29) = (X(29)-JVS(72)*XX(27)-JVS(75)*XX(28))/(JVS(77)) XX(30) = (X(30)-JVS(78)*XX(29))/(JVS(82)) XX(31) = X(31)/JVS(84) XX(32) = X(32)/JVS(87) XX(33) = X(33)/JVS(89) XX(34) = X(34)/JVS(91) XX(35) = X(35)/JVS(93) XX(36) = X(36)/JVS(95) XX(37) = X(37)/JVS(97) XX(38) = X(38)/JVS(99) XX(39) = X(39)/JVS(101) XX(40) = X(40)/JVS(103) XX(41) = X(41)/JVS(105) XX(42) = X(42)/JVS(107) XX(43) = X(43)/JVS(134) XX(44) = X(44)/JVS(156) XX(45) = (X(45)-JVS(20)*XX(10)-JVS(135)*XX(43))/(JVS(158)) XX(46) = (X(46)-JVS(4)*XX(2)-JVS(136)*XX(43))/(JVS(161)) XX(47) = (X(47)-JVS(94)*XX(35)-JVS(137)*XX(43))/(JVS(164)) XX(48) = (X(48)-JVS(24)*XX(12))/(JVS(167)) XX(49) = (X(49)-JVS(14)*XX(7)-JVS(138)*XX(43))/(JVS(171)) XX(50) = (X(50)-JVS(18)*XX(9)-JVS(139)*XX(43))/(JVS(174)) XX(51) = X(51)/JVS(179) XX(52) = X(52)/JVS(182) XX(53) = X(53)/JVS(185) XX(54) = X(54)/JVS(189) XX(55) = (X(55)-JVS(8)*XX(4)-JVS(140)*XX(43))/(JVS(193)) XX(56) = X(56)/JVS(196) XX(57) = X(57)/JVS(200) XX(58) = X(58)/JVS(204) XX(59) = X(59)/JVS(208) XX(60) = (X(60)-JVS(27)*XX(13)-JVS(168)*XX(48))/(JVS(212)) XX(61) = X(61)/JVS(215) XX(62) = X(62)/JVS(219) XX(63) = X(63)/JVS(223) XX(64) = X(64)/JVS(227) XX(65) = (X(65)-JVS(106)*XX(41)-JVS(108)*XX(42)-JVS(141)*XX(43))/(JVS(231)) XX(66) = X(66)/JVS(236) XX(67) = X(67)/JVS(241) XX(68) = (X(68)-JVS(109)*XX(42))/(JVS(245)) XX(69) = X(69)/JVS(252) XX(70) = (X(70)-JVS(110)*XX(42))/(JVS(256)) XX(71) = X(71)/JVS(260) XX(72) = (X(72)-JVS(102)*XX(39)-JVS(111)*XX(42)-JVS(142)*XX(43))/(JVS(264)) XX(73) = X(73)/JVS(276) XX(74) = (X(74)-JVS(92)*XX(34)-JVS(143)*XX(43)-JVS(165)*XX(47))/(JVS(280)) XX(75) = (X(75)-JVS(112)*XX(42)-JVS(232)*XX(65))/(JVS(285)) XX(76) = X(76)/JVS(290) XX(77) = X(77)/JVS(295) XX(78) = (X(78)-JVS(113)*XX(42))/(JVS(301)) XX(79) = (X(79)-JVS(16)*XX(8)-JVS(144)*XX(43))/(JVS(307)) XX(80) = X(80)/JVS(312) XX(81) = (X(81)-JVS(30)*XX(14)-JVS(114)*XX(42))/(JVS(317)) XX(82) = (X(82)-JVS(90)*XX(33)-JVS(145)*XX(43))/(JVS(321)) XX(83) = X(83)/JVS(326) XX(84) = (X(84)-JVS(115)*XX(42)-JVS(265)*XX(72))/(JVS(333)) XX(85) = X(85)/JVS(339) XX(86) = X(86)/JVS(349) XX(87) = (X(87)-JVS(116)*XX(42)-JVS(237)*XX(66))/(JVS(355)) XX(88) = X(88)/JVS(362) XX(89) = (X(89)-JVS(104)*XX(40)-JVS(117)*XX(42)-JVS(146)*XX(43))/(JVS(378)) XX(90) = X(90)/JVS(404) XX(91) = (X(91)-JVS(118)*XX(42))/(JVS(421)) XX(92) = (X(92)-JVS(119)*XX(42)-JVS(175)*XX(50)-JVS(266)*XX(72)-JVS(379)*XX(89))/(JVS(449)) XX(93) = (X(93)-JVS(85)*XX(31)-JVS(98)*XX(37)-JVS(147)*XX(43)-JVS(267)*XX(72))/(JVS(457)) XX(94) = (X(94)-JVS(197)*XX(56)-JVS(363)*XX(88)-JVS(380)*XX(89))/(JVS(467)) XX(95) = (X(95)-JVS(120)*XX(42)-JVS(176)*XX(50)-JVS(233)*XX(65)-JVS(286)*XX(75)-JVS(450)*XX(92))/(JVS(476)) XX(96) = (X(96)-JVS(100)*XX(38)-JVS(148)*XX(43)-JVS(422)*XX(91))/(JVS(484)) XX(97) = (X(97)-JVS(121)*XX(42)-JVS(381)*XX(89)-JVS(423)*XX(91)-JVS(485)*XX(96))/(JVS(495)) XX(98) = (X(98)-JVS(201)*XX(57)-JVS(340)*XX(85)-JVS(364)*XX(88)-JVS(382)*XX(89)-JVS(424)*XX(91))/(JVS(504)) XX(99) = (X(99)-JVS(6)*XX(3)-JVS(149)*XX(43))/(JVS(519)) XX(100) = (X(100)-JVS(122)*XX(42)-JVS(268)*XX(72)-JVS(405)*XX(90)-JVS(425)*XX(91)-JVS(520)*XX(99))/(JVS(546)) XX(101) = (X(101)-JVS(277)*XX(73)-JVS(365)*XX(88)-JVS(383)*XX(89)-JVS(458)*XX(93)-JVS(547)*XX(100))/(JVS(561)) XX(102) = (X(102)-JVS(123)*XX(42)-JVS(366)*XX(88)-JVS(384)*XX(89)-JVS(521)*XX(99))/(JVS(569)) XX(103) = (X(103)-JVS(246)*XX(68)-JVS(269)*XX(72)-JVS(313)*XX(80)-JVS(486)*XX(96)-JVS(522)*XX(99)-JVS(548)*XX(100))& &/(JVS(577)) XX(104) = (X(104)-JVS(96)*XX(36)-JVS(150)*XX(43)-JVS(426)*XX(91)-JVS(451)*XX(92))/(JVS(585)) XX(105) = (X(105)-JVS(257)*XX(70)-JVS(367)*XX(88)-JVS(385)*XX(89)-JVS(427)*XX(91)-JVS(586)*XX(104))/(JVS(595)) XX(106) = (X(106)-JVS(368)*XX(88)-JVS(386)*XX(89))/(JVS(602)) XX(107) = (X(107)-JVS(314)*XX(80)-JVS(459)*XX(93)-JVS(523)*XX(99)-JVS(549)*XX(100)-JVS(578)*XX(103)-JVS(587)*XX(104))& &/(JVS(613)) XX(108) = (X(108)-JVS(224)*XX(63)-JVS(387)*XX(89))/(JVS(623)) XX(109) = (X(109)-JVS(253)*XX(69)-JVS(369)*XX(88)-JVS(388)*XX(89))/(JVS(634)) XX(110) = (X(110)-JVS(242)*XX(67)-JVS(370)*XX(88)-JVS(389)*XX(89)-JVS(487)*XX(96))/(JVS(645)) XX(111) = (X(111)-JVS(291)*XX(76)-JVS(296)*XX(77)-JVS(302)*XX(78)-JVS(356)*XX(87)-JVS(371)*XX(88)-JVS(390)*XX(89)& &-JVS(406)*XX(90)-JVS(428)*XX(91)-JVS(524)*XX(99)-JVS(550)*XX(100))/(JVS(655)) XX(112) = (X(112)-JVS(216)*XX(61)-JVS(391)*XX(89))/(JVS(664)) XX(113) = (X(113)-JVS(159)*XX(45)-JVS(205)*XX(58)-JVS(665)*XX(112))/(JVS(678)) XX(114) = (X(114)-JVS(2)*XX(1)-JVS(151)*XX(43)-JVS(429)*XX(91)-JVS(525)*XX(99))/(JVS(724)) XX(115) = (X(115)-JVS(430)*XX(91)-JVS(526)*XX(99)-JVS(725)*XX(114))/(JVS(756)) XX(116) = (X(116)-JVS(372)*XX(88)-JVS(392)*XX(89)-JVS(726)*XX(114)-JVS(757)*XX(115))/(JVS(775)) XX(117) = (X(117)-JVS(261)*XX(71)-JVS(373)*XX(88)-JVS(393)*XX(89)-JVS(527)*XX(99)-JVS(727)*XX(114))/(JVS(787)) XX(118) = (X(118)-JVS(431)*XX(91)-JVS(528)*XX(99)-JVS(758)*XX(115))/(JVS(813)) XX(119) = (X(119)-JVS(124)*XX(42)-JVS(270)*XX(72)-JVS(394)*XX(89)-JVS(432)*XX(91)-JVS(728)*XX(114)-JVS(814)*XX(118))& &/(JVS(847)) XX(120) = (X(120)-JVS(350)*XX(86)-JVS(433)*XX(91)-JVS(562)*XX(101)-JVS(635)*XX(109)-JVS(646)*XX(110)-JVS(729)*XX(114)& &-JVS(759)*XX(115)-JVS(815)*XX(118)-JVS(848)*XX(119))/(JVS(867)) XX(121) = (X(121)-JVS(125)*XX(42)-JVS(434)*XX(91)-JVS(496)*XX(97)-JVS(529)*XX(99)-JVS(596)*XX(105)-JVS(730)*XX(114)& &-JVS(816)*XX(118)-JVS(849)*XX(119))/(JVS(879)) XX(122) = (X(122)-JVS(126)*XX(42)-JVS(228)*XX(64)-JVS(308)*XX(79)-JVS(435)*XX(91)-JVS(477)*XX(95)-JVS(570)*XX(102)& &-JVS(624)*XX(108)-JVS(731)*XX(114)-JVS(817)*XX(118)-JVS(850)*XX(119))/(JVS(892)) XX(123) = (X(123)-JVS(127)*XX(42)-JVS(281)*XX(74)-JVS(322)*XX(82)-JVS(327)*XX(83)-JVS(334)*XX(84)-JVS(374)*XX(88)& &-JVS(395)*XX(89)-JVS(407)*XX(90)-JVS(436)*XX(91)-JVS(460)*XX(93)-JVS(530)*XX(99)-JVS(551)*XX(100)-JVS(588)& &*XX(104)-JVS(614)*XX(107)-JVS(732)*XX(114)-JVS(776)*XX(116)-JVS(818)*XX(118)-JVS(851)*XX(119)-JVS(868)& &*XX(120)-JVS(880)*XX(121)-JVS(893)*XX(122))/(JVS(904)) XX(124) = (X(124)-JVS(271)*XX(72)-JVS(437)*XX(91)-JVS(552)*XX(100)-JVS(656)*XX(111)-JVS(733)*XX(114)-JVS(819)*XX(118)& &-JVS(852)*XX(119))/(JVS(916)) XX(125) = (X(125)-JVS(88)*XX(32)-JVS(152)*XX(43)-JVS(247)*XX(68)-JVS(488)*XX(96)-JVS(553)*XX(100)-JVS(579)*XX(103)& &-JVS(589)*XX(104)-JVS(615)*XX(107)-JVS(734)*XX(114)-JVS(777)*XX(116)-JVS(820)*XX(118)-JVS(853)*XX(119)& &-JVS(917)*XX(124))/(JVS(927)) XX(126) = (X(126)-JVS(238)*XX(66)-JVS(341)*XX(85)-JVS(375)*XX(88)-JVS(396)*XX(89)-JVS(438)*XX(91)-JVS(468)*XX(94)& &-JVS(505)*XX(98)-JVS(625)*XX(108)-JVS(666)*XX(112)-JVS(735)*XX(114)-JVS(760)*XX(115)-JVS(821)*XX(118)& &-JVS(854)*XX(119))/(JVS(943)) XX(127) = (X(127)-JVS(439)*XX(91)-JVS(531)*XX(99)-JVS(667)*XX(112)-JVS(679)*XX(113)-JVS(736)*XX(114)-JVS(761)*XX(115))& &/(JVS(980)) XX(128) = (X(128)-JVS(10)*XX(5)-JVS(153)*XX(43)-JVS(160)*XX(45)-JVS(172)*XX(49)-JVS(177)*XX(50)-JVS(186)*XX(53)& &-JVS(194)*XX(55)-JVS(209)*XX(59)-JVS(234)*XX(65)-JVS(287)*XX(75)-JVS(309)*XX(79)-JVS(452)*XX(92)-JVS(478)& &*XX(95)-JVS(532)*XX(99)-JVS(680)*XX(113)-JVS(737)*XX(114)-JVS(762)*XX(115)-JVS(822)*XX(118)-JVS(855)*XX(119)& &-JVS(894)*XX(122)-JVS(981)*XX(127))/(JVS(1044)) XX(129) = (X(129)-JVS(169)*XX(48)-JVS(180)*XX(51)-JVS(183)*XX(52)-JVS(195)*XX(55)-JVS(213)*XX(60)-JVS(318)*XX(81)& &-JVS(351)*XX(86)-JVS(408)*XX(90)-JVS(440)*XX(91)-JVS(497)*XX(97)-JVS(533)*XX(99)-JVS(603)*XX(106)-JVS(636)& &*XX(109)-JVS(647)*XX(110)-JVS(668)*XX(112)-JVS(681)*XX(113)-JVS(738)*XX(114)-JVS(763)*XX(115)-JVS(788)& &*XX(117)-JVS(823)*XX(118)-JVS(856)*XX(119)-JVS(869)*XX(120)-JVS(881)*XX(121)-JVS(895)*XX(122)-JVS(905)& &*XX(123)-JVS(918)*XX(124)-JVS(928)*XX(125)-JVS(944)*XX(126)-JVS(982)*XX(127)-JVS(1045)*XX(128))/(JVS(1092)) XX(130) = (X(130)-JVS(128)*XX(42)-JVS(173)*XX(49)-JVS(190)*XX(54)-JVS(297)*XX(77)-JVS(303)*XX(78)-JVS(328)*XX(83)& &-JVS(342)*XX(85)-JVS(357)*XX(87)-JVS(397)*XX(89)-JVS(409)*XX(90)-JVS(441)*XX(91)-JVS(469)*XX(94)-JVS(489)& &*XX(96)-JVS(498)*XX(97)-JVS(506)*XX(98)-JVS(534)*XX(99)-JVS(554)*XX(100)-JVS(563)*XX(101)-JVS(571)*XX(102)& &-JVS(597)*XX(105)-JVS(604)*XX(106)-JVS(626)*XX(108)-JVS(637)*XX(109)-JVS(648)*XX(110)-JVS(657)*XX(111)& &-JVS(669)*XX(112)-JVS(682)*XX(113)-JVS(739)*XX(114)-JVS(764)*XX(115)-JVS(778)*XX(116)-JVS(789)*XX(117)& &-JVS(824)*XX(118)-JVS(857)*XX(119)-JVS(870)*XX(120)-JVS(882)*XX(121)-JVS(896)*XX(122)-JVS(906)*XX(123)& &-JVS(919)*XX(124)-JVS(929)*XX(125)-JVS(945)*XX(126)-JVS(983)*XX(127)-JVS(1046)*XX(128)-JVS(1093)*XX(129))& &/(JVS(1140)) XX(131) = (X(131)-JVS(535)*XX(99)-JVS(605)*XX(106)-JVS(670)*XX(112)-JVS(683)*XX(113)-JVS(740)*XX(114)-JVS(765)*XX(115)& &-JVS(984)*XX(127)-JVS(1047)*XX(128)-JVS(1094)*XX(129)-JVS(1141)*XX(130))/(JVS(1170)) XX(132) = (X(132)-JVS(25)*XX(12)-JVS(28)*XX(13)-JVS(31)*XX(14)-JVS(40)*XX(17)-JVS(44)*XX(18)-JVS(53)*XX(21)-JVS(57)& &*XX(22)-JVS(66)*XX(25)-JVS(70)*XX(26)-JVS(79)*XX(29)-JVS(83)*XX(30)-JVS(86)*XX(31)-JVS(129)*XX(42)-JVS(157)& &*XX(44)-JVS(162)*XX(46)-JVS(166)*XX(47)-JVS(170)*XX(48)-JVS(181)*XX(51)-JVS(184)*XX(52)-JVS(187)*XX(53)& &-JVS(191)*XX(54)-JVS(198)*XX(56)-JVS(202)*XX(57)-JVS(206)*XX(58)-JVS(210)*XX(59)-JVS(214)*XX(60)-JVS(217)& &*XX(61)-JVS(220)*XX(62)-JVS(225)*XX(63)-JVS(229)*XX(64)-JVS(239)*XX(66)-JVS(243)*XX(67)-JVS(248)*XX(68)& &-JVS(254)*XX(69)-JVS(258)*XX(70)-JVS(262)*XX(71)-JVS(272)*XX(72)-JVS(278)*XX(73)-JVS(282)*XX(74)-JVS(292)& &*XX(76)-JVS(298)*XX(77)-JVS(304)*XX(78)-JVS(310)*XX(79)-JVS(315)*XX(80)-JVS(319)*XX(81)-JVS(323)*XX(82)& &-JVS(335)*XX(84)-JVS(343)*XX(85)-JVS(352)*XX(86)-JVS(358)*XX(87)-JVS(376)*XX(88)-JVS(398)*XX(89)-JVS(410)& &*XX(90)-JVS(442)*XX(91)-JVS(453)*XX(92)-JVS(461)*XX(93)-JVS(470)*XX(94)-JVS(479)*XX(95)-JVS(490)*XX(96)& &-JVS(507)*XX(98)-JVS(536)*XX(99)-JVS(555)*XX(100)-JVS(564)*XX(101)-JVS(572)*XX(102)-JVS(580)*XX(103)-JVS(590)& &*XX(104)-JVS(598)*XX(105)-JVS(606)*XX(106)-JVS(616)*XX(107)-JVS(627)*XX(108)-JVS(638)*XX(109)-JVS(649)& &*XX(110)-JVS(658)*XX(111)-JVS(671)*XX(112)-JVS(684)*XX(113)-JVS(741)*XX(114)-JVS(766)*XX(115)-JVS(779)& &*XX(116)-JVS(790)*XX(117)-JVS(825)*XX(118)-JVS(858)*XX(119)-JVS(871)*XX(120)-JVS(883)*XX(121)-JVS(897)& &*XX(122)-JVS(907)*XX(123)-JVS(920)*XX(124)-JVS(930)*XX(125)-JVS(946)*XX(126)-JVS(985)*XX(127)-JVS(1048)& &*XX(128)-JVS(1095)*XX(129)-JVS(1142)*XX(130)-JVS(1171)*XX(131))/(JVS(1260)) XX(133) = (X(133)-JVS(130)*XX(42)-JVS(221)*XX(62)-JVS(344)*XX(85)-JVS(377)*XX(88)-JVS(399)*XX(89)-JVS(443)*XX(91)& &-JVS(471)*XX(94)-JVS(491)*XX(96)-JVS(499)*XX(97)-JVS(508)*XX(98)-JVS(537)*XX(99)-JVS(556)*XX(100)-JVS(565)& &*XX(101)-JVS(573)*XX(102)-JVS(599)*XX(105)-JVS(607)*XX(106)-JVS(628)*XX(108)-JVS(639)*XX(109)-JVS(650)& &*XX(110)-JVS(659)*XX(111)-JVS(672)*XX(112)-JVS(685)*XX(113)-JVS(742)*XX(114)-JVS(767)*XX(115)-JVS(780)& &*XX(116)-JVS(791)*XX(117)-JVS(826)*XX(118)-JVS(859)*XX(119)-JVS(872)*XX(120)-JVS(884)*XX(121)-JVS(898)& &*XX(122)-JVS(908)*XX(123)-JVS(921)*XX(124)-JVS(931)*XX(125)-JVS(947)*XX(126)-JVS(986)*XX(127)-JVS(1049)& &*XX(128)-JVS(1096)*XX(129)-JVS(1143)*XX(130)-JVS(1172)*XX(131)-JVS(1261)*XX(132))/(JVS(1310)) XX(134) = (X(134)-JVS(37)*XX(16)-JVS(41)*XX(17)-JVS(50)*XX(20)-JVS(54)*XX(21)-JVS(63)*XX(24)-JVS(67)*XX(25)-JVS(76)& &*XX(28)-JVS(80)*XX(29)-JVS(131)*XX(42)-JVS(178)*XX(50)-JVS(188)*XX(53)-JVS(235)*XX(65)-JVS(249)*XX(68)& &-JVS(273)*XX(72)-JVS(288)*XX(75)-JVS(299)*XX(77)-JVS(305)*XX(78)-JVS(324)*XX(82)-JVS(329)*XX(83)-JVS(336)& &*XX(84)-JVS(345)*XX(85)-JVS(359)*XX(87)-JVS(400)*XX(89)-JVS(411)*XX(90)-JVS(444)*XX(91)-JVS(454)*XX(92)& &-JVS(462)*XX(93)-JVS(472)*XX(94)-JVS(480)*XX(95)-JVS(492)*XX(96)-JVS(500)*XX(97)-JVS(509)*XX(98)-JVS(538)& &*XX(99)-JVS(557)*XX(100)-JVS(566)*XX(101)-JVS(574)*XX(102)-JVS(581)*XX(103)-JVS(591)*XX(104)-JVS(600)*XX(105)& &-JVS(608)*XX(106)-JVS(617)*XX(107)-JVS(629)*XX(108)-JVS(640)*XX(109)-JVS(651)*XX(110)-JVS(660)*XX(111)& &-JVS(673)*XX(112)-JVS(686)*XX(113)-JVS(743)*XX(114)-JVS(768)*XX(115)-JVS(781)*XX(116)-JVS(792)*XX(117)& &-JVS(827)*XX(118)-JVS(860)*XX(119)-JVS(873)*XX(120)-JVS(885)*XX(121)-JVS(899)*XX(122)-JVS(909)*XX(123)& &-JVS(922)*XX(124)-JVS(932)*XX(125)-JVS(948)*XX(126)-JVS(987)*XX(127)-JVS(1050)*XX(128)-JVS(1097)*XX(129)& &-JVS(1144)*XX(130)-JVS(1173)*XX(131)-JVS(1262)*XX(132)-JVS(1311)*XX(133))/(JVS(1355)) XX(135) = (X(135)-JVS(12)*XX(6)-JVS(132)*XX(42)-JVS(154)*XX(43)-JVS(250)*XX(68)-JVS(274)*XX(72)-JVS(306)*XX(78)& &-JVS(311)*XX(79)-JVS(320)*XX(81)-JVS(325)*XX(82)-JVS(353)*XX(86)-JVS(360)*XX(87)-JVS(445)*XX(91)-JVS(463)& &*XX(93)-JVS(481)*XX(95)-JVS(493)*XX(96)-JVS(539)*XX(99)-JVS(558)*XX(100)-JVS(592)*XX(104)-JVS(618)*XX(107)& &-JVS(641)*XX(109)-JVS(652)*XX(110)-JVS(744)*XX(114)-JVS(769)*XX(115)-JVS(782)*XX(116)-JVS(793)*XX(117)& &-JVS(828)*XX(118)-JVS(861)*XX(119)-JVS(874)*XX(120)-JVS(886)*XX(121)-JVS(900)*XX(122)-JVS(910)*XX(123)& &-JVS(923)*XX(124)-JVS(933)*XX(125)-JVS(949)*XX(126)-JVS(988)*XX(127)-JVS(1051)*XX(128)-JVS(1098)*XX(129)& &-JVS(1145)*XX(130)-JVS(1174)*XX(131)-JVS(1263)*XX(132)-JVS(1312)*XX(133)-JVS(1356)*XX(134))/(JVS(1381)) XX(136) = (X(136)-JVS(34)*XX(15)-JVS(42)*XX(17)-JVS(47)*XX(19)-JVS(55)*XX(21)-JVS(60)*XX(23)-JVS(68)*XX(25)-JVS(73)& &*XX(27)-JVS(81)*XX(29)-JVS(133)*XX(42)-JVS(163)*XX(46)-JVS(192)*XX(54)-JVS(199)*XX(56)-JVS(203)*XX(57)& &-JVS(207)*XX(58)-JVS(211)*XX(59)-JVS(218)*XX(61)-JVS(222)*XX(62)-JVS(226)*XX(63)-JVS(230)*XX(64)-JVS(240)& &*XX(66)-JVS(244)*XX(67)-JVS(251)*XX(68)-JVS(255)*XX(69)-JVS(259)*XX(70)-JVS(263)*XX(71)-JVS(275)*XX(72)& &-JVS(279)*XX(73)-JVS(283)*XX(74)-JVS(289)*XX(75)-JVS(293)*XX(76)-JVS(300)*XX(77)-JVS(316)*XX(80)-JVS(330)& &*XX(83)-JVS(337)*XX(84)-JVS(346)*XX(85)-JVS(361)*XX(87)-JVS(401)*XX(89)-JVS(412)*XX(90)-JVS(446)*XX(91)& &-JVS(455)*XX(92)-JVS(464)*XX(93)-JVS(473)*XX(94)-JVS(482)*XX(95)-JVS(494)*XX(96)-JVS(501)*XX(97)-JVS(510)& &*XX(98)-JVS(540)*XX(99)-JVS(559)*XX(100)-JVS(567)*XX(101)-JVS(575)*XX(102)-JVS(582)*XX(103)-JVS(593)*XX(104)& &-JVS(601)*XX(105)-JVS(609)*XX(106)-JVS(619)*XX(107)-JVS(630)*XX(108)-JVS(642)*XX(109)-JVS(653)*XX(110)& &-JVS(661)*XX(111)-JVS(674)*XX(112)-JVS(687)*XX(113)-JVS(745)*XX(114)-JVS(770)*XX(115)-JVS(783)*XX(116)& &-JVS(794)*XX(117)-JVS(829)*XX(118)-JVS(862)*XX(119)-JVS(875)*XX(120)-JVS(887)*XX(121)-JVS(901)*XX(122)& &-JVS(911)*XX(123)-JVS(924)*XX(124)-JVS(934)*XX(125)-JVS(950)*XX(126)-JVS(989)*XX(127)-JVS(1052)*XX(128)& &-JVS(1099)*XX(129)-JVS(1146)*XX(130)-JVS(1175)*XX(131)-JVS(1264)*XX(132)-JVS(1313)*XX(133)-JVS(1357)*XX(134)& &-JVS(1382)*XX(135))/(JVS(1459)) XX(137) = (X(137)-JVS(22)*XX(11)-JVS(155)*XX(43)-JVS(347)*XX(85)-JVS(447)*XX(91)-JVS(474)*XX(94)-JVS(511)*XX(98)& &-JVS(631)*XX(108)-JVS(675)*XX(112)-JVS(746)*XX(114)-JVS(771)*XX(115)-JVS(784)*XX(116)-JVS(830)*XX(118)& &-JVS(863)*XX(119)-JVS(951)*XX(126)-JVS(990)*XX(127)-JVS(1053)*XX(128)-JVS(1100)*XX(129)-JVS(1147)*XX(130)& &-JVS(1176)*XX(131)-JVS(1265)*XX(132)-JVS(1314)*XX(133)-JVS(1358)*XX(134)-JVS(1383)*XX(135)-JVS(1460)*XX(136))& &/(JVS(1481)) XX(137) = XX(137) XX(136) = XX(136)-JVS(1480)*XX(137) XX(135) = XX(135)-JVS(1458)*XX(136)-JVS(1479)*XX(137) XX(134) = XX(134)-JVS(1380)*XX(135)-JVS(1457)*XX(136)-JVS(1478)*XX(137) XX(133) = XX(133)-JVS(1354)*XX(134)-JVS(1379)*XX(135)-JVS(1456)*XX(136)-JVS(1477)*XX(137) XX(132) = XX(132)-JVS(1309)*XX(133)-JVS(1353)*XX(134)-JVS(1378)*XX(135)-JVS(1455)*XX(136)-JVS(1476)*XX(137) XX(131) = XX(131)-JVS(1259)*XX(132)-JVS(1308)*XX(133)-JVS(1352)*XX(134)-JVS(1377)*XX(135)-JVS(1454)*XX(136)-JVS(1475)& &*XX(137) XX(130) = XX(130)-JVS(1169)*XX(131)-JVS(1258)*XX(132)-JVS(1307)*XX(133)-JVS(1351)*XX(134)-JVS(1376)*XX(135)-JVS(1453)& &*XX(136)-JVS(1474)*XX(137) XX(129) = XX(129)-JVS(1139)*XX(130)-JVS(1168)*XX(131)-JVS(1257)*XX(132)-JVS(1306)*XX(133)-JVS(1350)*XX(134)-JVS(1375)& &*XX(135)-JVS(1452)*XX(136)-JVS(1473)*XX(137) XX(128) = XX(128)-JVS(1091)*XX(129)-JVS(1138)*XX(130)-JVS(1167)*XX(131)-JVS(1256)*XX(132)-JVS(1305)*XX(133)-JVS(1349)& &*XX(134)-JVS(1374)*XX(135)-JVS(1451)*XX(136)-JVS(1472)*XX(137) XX(127) = XX(127)-JVS(1043)*XX(128)-JVS(1090)*XX(129)-JVS(1137)*XX(130)-JVS(1255)*XX(132)-JVS(1304)*XX(133)-JVS(1348)& &*XX(134)-JVS(1373)*XX(135)-JVS(1450)*XX(136) XX(126) = XX(126)-JVS(979)*XX(127)-JVS(1042)*XX(128)-JVS(1089)*XX(129)-JVS(1136)*XX(130)-JVS(1166)*XX(131)-JVS(1254)& &*XX(132)-JVS(1303)*XX(133)-JVS(1347)*XX(134)-JVS(1449)*XX(136)-JVS(1471)*XX(137) XX(125) = XX(125)-JVS(942)*XX(126)-JVS(978)*XX(127)-JVS(1041)*XX(128)-JVS(1088)*XX(129)-JVS(1135)*XX(130)-JVS(1165)& &*XX(131)-JVS(1253)*XX(132)-JVS(1302)*XX(133)-JVS(1346)*XX(134)-JVS(1372)*XX(135)-JVS(1448)*XX(136)-JVS(1470)& &*XX(137) XX(124) = XX(124)-JVS(941)*XX(126)-JVS(977)*XX(127)-JVS(1040)*XX(128)-JVS(1087)*XX(129)-JVS(1134)*XX(130)-JVS(1164)& &*XX(131)-JVS(1252)*XX(132)-JVS(1301)*XX(133)-JVS(1345)*XX(134)-JVS(1371)*XX(135)-JVS(1447)*XX(136) XX(123) = XX(123)-JVS(915)*XX(124)-JVS(926)*XX(125)-JVS(940)*XX(126)-JVS(976)*XX(127)-JVS(1039)*XX(128)-JVS(1086)& &*XX(129)-JVS(1133)*XX(130)-JVS(1163)*XX(131)-JVS(1251)*XX(132)-JVS(1300)*XX(133)-JVS(1344)*XX(134)-JVS(1370)& &*XX(135)-JVS(1446)*XX(136)-JVS(1469)*XX(137) XX(122) = XX(122)-JVS(975)*XX(127)-JVS(1038)*XX(128)-JVS(1085)*XX(129)-JVS(1132)*XX(130)-JVS(1162)*XX(131)-JVS(1250)& &*XX(132)-JVS(1299)*XX(133)-JVS(1343)*XX(134)-JVS(1369)*XX(135)-JVS(1445)*XX(136)-JVS(1468)*XX(137) XX(121) = XX(121)-JVS(891)*XX(122)-JVS(974)*XX(127)-JVS(1037)*XX(128)-JVS(1084)*XX(129)-JVS(1131)*XX(130)-JVS(1161)& &*XX(131)-JVS(1249)*XX(132)-JVS(1298)*XX(133)-JVS(1342)*XX(134)-JVS(1368)*XX(135)-JVS(1444)*XX(136)-JVS(1467)& &*XX(137) XX(120) = XX(120)-JVS(890)*XX(122)-JVS(973)*XX(127)-JVS(1036)*XX(128)-JVS(1083)*XX(129)-JVS(1130)*XX(130)-JVS(1160)& &*XX(131)-JVS(1248)*XX(132)-JVS(1297)*XX(133)-JVS(1341)*XX(134)-JVS(1367)*XX(135)-JVS(1443)*XX(136) XX(119) = XX(119)-JVS(1035)*XX(128)-JVS(1082)*XX(129)-JVS(1129)*XX(130)-JVS(1247)*XX(132)-JVS(1296)*XX(133)-JVS(1442)& &*XX(136) XX(118) = XX(118)-JVS(1034)*XX(128)-JVS(1081)*XX(129)-JVS(1128)*XX(130)-JVS(1246)*XX(132)-JVS(1295)*XX(133)-JVS(1441)& &*XX(136) XX(117) = XX(117)-JVS(866)*XX(120)-JVS(878)*XX(121)-JVS(914)*XX(124)-JVS(925)*XX(125)-JVS(972)*XX(127)-JVS(1033)& &*XX(128)-JVS(1080)*XX(129)-JVS(1127)*XX(130)-JVS(1159)*XX(131)-JVS(1245)*XX(132)-JVS(1294)*XX(133)-JVS(1340)& &*XX(134)-JVS(1440)*XX(136)-JVS(1466)*XX(137) XX(116) = XX(116)-JVS(939)*XX(126)-JVS(971)*XX(127)-JVS(1032)*XX(128)-JVS(1079)*XX(129)-JVS(1126)*XX(130)-JVS(1244)& &*XX(132)-JVS(1293)*XX(133)-JVS(1339)*XX(134)-JVS(1439)*XX(136)-JVS(1465)*XX(137) XX(115) = XX(115)-JVS(1031)*XX(128)-JVS(1078)*XX(129)-JVS(1125)*XX(130)-JVS(1243)*XX(132)-JVS(1292)*XX(133)-JVS(1438)& &*XX(136) XX(114) = XX(114)-JVS(1030)*XX(128)-JVS(1077)*XX(129)-JVS(1242)*XX(132)-JVS(1437)*XX(136) XX(113) = XX(113)-JVS(723)*XX(114)-JVS(755)*XX(115)-JVS(1029)*XX(128)-JVS(1124)*XX(130)-JVS(1241)*XX(132)-JVS(1291)& &*XX(133)-JVS(1338)*XX(134)-JVS(1366)*XX(135)-JVS(1436)*XX(136) XX(112) = XX(112)-JVS(722)*XX(114)-JVS(754)*XX(115)-JVS(1028)*XX(128)-JVS(1123)*XX(130)-JVS(1240)*XX(132)-JVS(1290)& &*XX(133)-JVS(1337)*XX(134)-JVS(1435)*XX(136) XX(111) = XX(111)-JVS(721)*XX(114)-JVS(812)*XX(118)-JVS(846)*XX(119)-JVS(938)*XX(126)-JVS(970)*XX(127)-JVS(1027)& &*XX(128)-JVS(1076)*XX(129)-JVS(1122)*XX(130)-JVS(1158)*XX(131)-JVS(1239)*XX(132)-JVS(1289)*XX(133)-JVS(1336)& &*XX(134)-JVS(1365)*XX(135)-JVS(1434)*XX(136) XX(110) = XX(110)-JVS(720)*XX(114)-JVS(753)*XX(115)-JVS(811)*XX(118)-JVS(969)*XX(127)-JVS(1026)*XX(128)-JVS(1121)& &*XX(130)-JVS(1238)*XX(132)-JVS(1288)*XX(133)-JVS(1335)*XX(134)-JVS(1433)*XX(136) XX(109) = XX(109)-JVS(719)*XX(114)-JVS(752)*XX(115)-JVS(845)*XX(119)-JVS(968)*XX(127)-JVS(1025)*XX(128)-JVS(1120)& &*XX(130)-JVS(1237)*XX(132)-JVS(1287)*XX(133)-JVS(1334)*XX(134)-JVS(1432)*XX(136) XX(108) = XX(108)-JVS(718)*XX(114)-JVS(810)*XX(118)-JVS(844)*XX(119)-JVS(1024)*XX(128)-JVS(1119)*XX(130)-JVS(1236)& &*XX(132)-JVS(1286)*XX(133)-JVS(1333)*XX(134)-JVS(1431)*XX(136)-JVS(1464)*XX(137) XX(107) = XX(107)-JVS(717)*XX(114)-JVS(774)*XX(116)-JVS(809)*XX(118)-JVS(843)*XX(119)-JVS(967)*XX(127)-JVS(1023)& &*XX(128)-JVS(1075)*XX(129)-JVS(1118)*XX(130)-JVS(1157)*XX(131)-JVS(1235)*XX(132)-JVS(1285)*XX(133)-JVS(1332)& &*XX(134)-JVS(1430)*XX(136) XX(106) = XX(106)-JVS(716)*XX(114)-JVS(751)*XX(115)-JVS(966)*XX(127)-JVS(1022)*XX(128)-JVS(1117)*XX(130)-JVS(1156)& &*XX(131)-JVS(1234)*XX(132)-JVS(1284)*XX(133)-JVS(1331)*XX(134)-JVS(1429)*XX(136)-JVS(1463)*XX(137) XX(105) = XX(105)-JVS(715)*XX(114)-JVS(808)*XX(118)-JVS(842)*XX(119)-JVS(965)*XX(127)-JVS(1021)*XX(128)-JVS(1074)& &*XX(129)-JVS(1116)*XX(130)-JVS(1155)*XX(131)-JVS(1233)*XX(132)-JVS(1283)*XX(133)-JVS(1330)*XX(134)-JVS(1428)& &*XX(136) XX(104) = XX(104)-JVS(714)*XX(114)-JVS(807)*XX(118)-JVS(841)*XX(119)-JVS(1020)*XX(128)-JVS(1073)*XX(129)-JVS(1232)& &*XX(132)-JVS(1282)*XX(133)-JVS(1329)*XX(134)-JVS(1427)*XX(136) XX(103) = XX(103)-JVS(584)*XX(104)-JVS(612)*XX(107)-JVS(713)*XX(114)-JVS(773)*XX(116)-JVS(806)*XX(118)-JVS(840)& &*XX(119)-JVS(964)*XX(127)-JVS(1019)*XX(128)-JVS(1072)*XX(129)-JVS(1154)*XX(131)-JVS(1231)*XX(132)-JVS(1281)& &*XX(133)-JVS(1328)*XX(134)-JVS(1426)*XX(136) XX(102) = XX(102)-JVS(712)*XX(114)-JVS(839)*XX(119)-JVS(963)*XX(127)-JVS(1018)*XX(128)-JVS(1071)*XX(129)-JVS(1115)& &*XX(130)-JVS(1153)*XX(131)-JVS(1230)*XX(132)-JVS(1280)*XX(133)-JVS(1327)*XX(134)-JVS(1425)*XX(136) XX(101) = XX(101)-JVS(711)*XX(114)-JVS(805)*XX(118)-JVS(962)*XX(127)-JVS(1017)*XX(128)-JVS(1070)*XX(129)-JVS(1114)& &*XX(130)-JVS(1152)*XX(131)-JVS(1229)*XX(132)-JVS(1279)*XX(133)-JVS(1326)*XX(134)-JVS(1424)*XX(136) XX(100) = XX(100)-JVS(710)*XX(114)-JVS(1016)*XX(128)-JVS(1069)*XX(129)-JVS(1228)*XX(132)-JVS(1423)*XX(136) XX(99) = XX(99)-JVS(1015)*XX(128)-JVS(1068)*XX(129)-JVS(1227)*XX(132) XX(98) = XX(98)-JVS(622)*XX(108)-JVS(709)*XX(114)-JVS(961)*XX(127)-JVS(1014)*XX(128)-JVS(1113)*XX(130)-JVS(1226)& &*XX(132)-JVS(1278)*XX(133)-JVS(1325)*XX(134)-JVS(1422)*XX(136) XX(97) = XX(97)-JVS(708)*XX(114)-JVS(804)*XX(118)-JVS(960)*XX(127)-JVS(1013)*XX(128)-JVS(1112)*XX(130)-JVS(1225)& &*XX(132)-JVS(1277)*XX(133)-JVS(1324)*XX(134)-JVS(1421)*XX(136)-JVS(1462)*XX(137) XX(96) = XX(96)-JVS(707)*XX(114)-JVS(803)*XX(118)-JVS(1012)*XX(128)-JVS(1224)*XX(132)-JVS(1276)*XX(133)-JVS(1420)& &*XX(136) XX(95) = XX(95)-JVS(706)*XX(114)-JVS(802)*XX(118)-JVS(838)*XX(119)-JVS(1011)*XX(128)-JVS(1067)*XX(129)-JVS(1111)& &*XX(130)-JVS(1223)*XX(132)-JVS(1275)*XX(133)-JVS(1323)*XX(134)-JVS(1419)*XX(136) XX(94) = XX(94)-JVS(705)*XX(114)-JVS(959)*XX(127)-JVS(1010)*XX(128)-JVS(1110)*XX(130)-JVS(1222)*XX(132)-JVS(1274)& &*XX(133)-JVS(1322)*XX(134)-JVS(1418)*XX(136) XX(93) = XX(93)-JVS(545)*XX(100)-JVS(704)*XX(114)-JVS(801)*XX(118)-JVS(1009)*XX(128)-JVS(1066)*XX(129)-JVS(1109)& &*XX(130)-JVS(1221)*XX(132)-JVS(1417)*XX(136) XX(92) = XX(92)-JVS(703)*XX(114)-JVS(800)*XX(118)-JVS(837)*XX(119)-JVS(1008)*XX(128)-JVS(1065)*XX(129)-JVS(1220)& &*XX(132)-JVS(1273)*XX(133)-JVS(1321)*XX(134)-JVS(1416)*XX(136) XX(91) = XX(91)-JVS(1219)*XX(132)-JVS(1415)*XX(136) XX(90) = XX(90)-JVS(420)*XX(91)-JVS(518)*XX(99)-JVS(702)*XX(114)-JVS(1064)*XX(129)-JVS(1218)*XX(132)-JVS(1414)*XX(136) XX(89) = XX(89)-JVS(1217)*XX(132)-JVS(1272)*XX(133) XX(88) = XX(88)-JVS(958)*XX(127)-JVS(1216)*XX(132)-JVS(1413)*XX(136) XX(87) = XX(87)-JVS(517)*XX(99)-JVS(937)*XX(126)-JVS(957)*XX(127)-JVS(1007)*XX(128)-JVS(1215)*XX(132)-JVS(1320)& &*XX(134)-JVS(1412)*XX(136) XX(86) = XX(86)-JVS(419)*XX(91)-JVS(633)*XX(109)-JVS(644)*XX(110)-JVS(701)*XX(114)-JVS(750)*XX(115)-JVS(1063)*XX(129)& &-JVS(1214)*XX(132)-JVS(1271)*XX(133)-JVS(1364)*XX(135)-JVS(1411)*XX(136) XX(85) = XX(85)-JVS(418)*XX(91)-JVS(621)*XX(108)-JVS(1108)*XX(130)-JVS(1213)*XX(132)-JVS(1270)*XX(133) XX(84) = XX(84)-JVS(403)*XX(90)-JVS(417)*XX(91)-JVS(544)*XX(100)-JVS(700)*XX(114)-JVS(799)*XX(118)-JVS(836)*XX(119)& &-JVS(1006)*XX(128)-JVS(1212)*XX(132)-JVS(1319)*XX(134)-JVS(1410)*XX(136) XX(83) = XX(83)-JVS(402)*XX(90)-JVS(516)*XX(99)-JVS(543)*XX(100)-JVS(798)*XX(118)-JVS(835)*XX(119)-JVS(1005)*XX(128)& &-JVS(1211)*XX(132)-JVS(1318)*XX(134)-JVS(1409)*XX(136) XX(82) = XX(82)-JVS(456)*XX(93)-JVS(583)*XX(104)-JVS(611)*XX(107)-JVS(699)*XX(114)-JVS(913)*XX(124)-JVS(1004)*XX(128)& &-JVS(1210)*XX(132)-JVS(1363)*XX(135)-JVS(1408)*XX(136) XX(81) = XX(81)-JVS(348)*XX(86)-JVS(416)*XX(91)-JVS(698)*XX(114)-JVS(786)*XX(117)-JVS(865)*XX(120)-JVS(877)*XX(121)& &-JVS(903)*XX(123)-JVS(1062)*XX(129)-JVS(1209)*XX(132)-JVS(1362)*XX(135)-JVS(1407)*XX(136) XX(80) = XX(80)-JVS(576)*XX(103)-JVS(610)*XX(107)-JVS(772)*XX(116)-JVS(956)*XX(127)-JVS(1003)*XX(128)-JVS(1208)& &*XX(132)-JVS(1406)*XX(136) XX(79) = XX(79)-JVS(475)*XX(95)-JVS(697)*XX(114)-JVS(889)*XX(122)-JVS(1002)*XX(128)-JVS(1207)*XX(132)-JVS(1361)& &*XX(135)-JVS(1405)*XX(136) XX(78) = XX(78)-JVS(354)*XX(87)-JVS(1151)*XX(131)-JVS(1206)*XX(132)-JVS(1360)*XX(135)-JVS(1404)*XX(136) XX(77) = XX(77)-JVS(415)*XX(91)-JVS(834)*XX(119)-JVS(1150)*XX(131)-JVS(1205)*XX(132)-JVS(1403)*XX(136) XX(76) = XX(76)-JVS(294)*XX(77)-JVS(414)*XX(91)-JVS(542)*XX(100)-JVS(654)*XX(111)-JVS(797)*XX(118)-JVS(833)*XX(119)& &-JVS(1204)*XX(132)-JVS(1402)*XX(136) XX(75) = XX(75)-JVS(1001)*XX(128)-JVS(1061)*XX(129)-JVS(1107)*XX(130)-JVS(1203)*XX(132)-JVS(1317)*XX(134)-JVS(1401)& &*XX(136) XX(74) = XX(74)-JVS(332)*XX(84)-JVS(696)*XX(114)-JVS(864)*XX(120)-JVS(876)*XX(121)-JVS(902)*XX(123)-JVS(912)*XX(124)& &-JVS(1202)*XX(132)-JVS(1400)*XX(136) XX(73) = XX(73)-JVS(541)*XX(100)-JVS(560)*XX(101)-JVS(695)*XX(114)-JVS(796)*XX(118)-JVS(1106)*XX(130)-JVS(1149)& &*XX(131)-JVS(1201)*XX(132)-JVS(1399)*XX(136) XX(72) = XX(72)-JVS(1200)*XX(132)-JVS(1398)*XX(136) XX(71) = XX(71)-JVS(785)*XX(117)-JVS(955)*XX(127)-JVS(1000)*XX(128)-JVS(1148)*XX(131)-JVS(1199)*XX(132)-JVS(1397)& &*XX(136)-JVS(1461)*XX(137) XX(70) = XX(70)-JVS(413)*XX(91)-JVS(594)*XX(105)-JVS(694)*XX(114)-JVS(832)*XX(119)-JVS(1198)*XX(132)-JVS(1396)*XX(136) XX(69) = XX(69)-JVS(632)*XX(109)-JVS(693)*XX(114)-JVS(749)*XX(115)-JVS(831)*XX(119)-JVS(1197)*XX(132)-JVS(1395)& &*XX(136) XX(68) = XX(68)-JVS(692)*XX(114)-JVS(999)*XX(128)-JVS(1196)*XX(132) XX(67) = XX(67)-JVS(483)*XX(96)-JVS(643)*XX(110)-JVS(954)*XX(127)-JVS(998)*XX(128)-JVS(1195)*XX(132)-JVS(1394)*XX(136) XX(66) = XX(66)-JVS(936)*XX(126)-JVS(953)*XX(127)-JVS(1194)*XX(132)-JVS(1393)*XX(136) XX(65) = XX(65)-JVS(284)*XX(75)-JVS(997)*XX(128)-JVS(1060)*XX(129)-JVS(1105)*XX(130) XX(64) = XX(64)-JVS(568)*XX(102)-JVS(691)*XX(114)-JVS(888)*XX(122)-JVS(1104)*XX(130)-JVS(1193)*XX(132) XX(63) = XX(63)-JVS(620)*XX(108)-JVS(690)*XX(114)-JVS(795)*XX(118)-JVS(1103)*XX(130)-JVS(1192)*XX(132) XX(62) = XX(62)-JVS(689)*XX(114)-JVS(1191)*XX(132)-JVS(1269)*XX(133)-JVS(1392)*XX(136) XX(61) = XX(61)-JVS(663)*XX(112)-JVS(748)*XX(115)-JVS(1190)*XX(132)-JVS(1391)*XX(136) XX(60) = XX(60)-JVS(515)*XX(99)-JVS(688)*XX(114)-JVS(1059)*XX(129)-JVS(1189)*XX(132)-JVS(1268)*XX(133)-JVS(1390)& &*XX(136) XX(59) = XX(59)-JVS(996)*XX(128)-JVS(1058)*XX(129)-JVS(1188)*XX(132)-JVS(1389)*XX(136) XX(58) = XX(58)-JVS(677)*XX(113)-JVS(747)*XX(115)-JVS(1187)*XX(132)-JVS(1388)*XX(136) XX(57) = XX(57)-JVS(338)*XX(85)-JVS(503)*XX(98)-JVS(1186)*XX(132)-JVS(1387)*XX(136) XX(56) = XX(56)-JVS(466)*XX(94)-JVS(952)*XX(127)-JVS(1185)*XX(132)-JVS(1386)*XX(136) XX(55) = XX(55)-JVS(514)*XX(99)-JVS(995)*XX(128)-JVS(1057)*XX(129)-JVS(1316)*XX(134)-JVS(1359)*XX(135) XX(54) = XX(54)-JVS(1102)*XX(130)-JVS(1184)*XX(132)-JVS(1267)*XX(133) XX(53) = XX(53)-JVS(994)*XX(128)-JVS(1183)*XX(132)-JVS(1315)*XX(134) XX(52) = XX(52)-JVS(513)*XX(99)-JVS(662)*XX(112)-JVS(1056)*XX(129)-JVS(1182)*XX(132) XX(51) = XX(51)-JVS(512)*XX(99)-JVS(935)*XX(126)-JVS(1055)*XX(129)-JVS(1181)*XX(132) XX(50) = XX(50)-JVS(448)*XX(92)-JVS(993)*XX(128) XX(49) = XX(49)-JVS(992)*XX(128)-JVS(1054)*XX(129)-JVS(1101)*XX(130)-JVS(1266)*XX(133) XX(48) = XX(48)-JVS(1180)*XX(132)-JVS(1385)*XX(136) XX(47) = XX(47)-JVS(331)*XX(84)-JVS(1179)*XX(132) XX(46) = XX(46)-JVS(1178)*XX(132)-JVS(1384)*XX(136) XX(45) = XX(45)-JVS(676)*XX(113)-JVS(991)*XX(128) XX(44) = XX(44)-JVS(465)*XX(94)-JVS(502)*XX(98)-JVS(1177)*XX(132) XX(43) = XX(43) XX(42) = XX(42) XX(41) = XX(41) XX(40) = XX(40) XX(39) = XX(39) XX(38) = XX(38) XX(37) = XX(37) XX(36) = XX(36) XX(35) = XX(35) XX(34) = XX(34) XX(33) = XX(33) XX(32) = XX(32) XX(31) = XX(31) XX(30) = XX(30) XX(29) = XX(29) XX(28) = XX(28) XX(27) = XX(27) XX(26) = XX(26) XX(25) = XX(25) XX(24) = XX(24) XX(23) = XX(23) XX(22) = XX(22) XX(21) = XX(21) XX(20) = XX(20) XX(19) = XX(19) XX(18) = XX(18) XX(17) = XX(17) XX(16) = XX(16) XX(15) = XX(15) XX(14) = XX(14) XX(13) = XX(13) XX(12) = XX(12) XX(11) = XX(11) XX(10) = XX(10) XX(9) = XX(9) XX(8) = XX(8) XX(7) = XX(7) XX(6) = XX(6) XX(5) = XX(5) XX(4) = XX(4) XX(3) = XX(3) XX(2) = XX(2) XX(1) = XX(1) END SUBROUTINE KppSolveTR ! End of KppSolveTR function ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! BLAS_UTIL - BLAS-LIKE utility functions ! Arguments : ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !-------------------------------------------------------------- ! ! BLAS/LAPACK-like subroutines used by the integration algorithms ! It is recommended to replace them by calls to the optimized ! BLAS/LAPACK library for your machine ! ! (C) Adrian Sandu, Aug. 2004 ! Virginia Polytechnic Institute and State University !-------------------------------------------------------------- !-------------------------------------------------------------- SUBROUTINE WCOPY(N,X,incX,Y,incY) !-------------------------------------------------------------- ! copies a vector, x, to a vector, y: y <- x ! only for incX=incY=1 ! after BLAS ! replace this by the function from the optimized BLAS implementation: ! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1) !-------------------------------------------------------------- ! USE gckpp_Precision INTEGER :: i,incX,incY,M,MP1,N REAL(kind=dp) :: X(N),Y(N) IF (N.LE.0) RETURN M = MOD(N,8) IF( M .NE. 0 ) THEN DO i = 1,M Y(i) = X(i) END DO IF( N .LT. 8 ) RETURN END IF MP1 = M+1 DO i = MP1,N,8 Y(i) = X(i) Y(i + 1) = X(i + 1) Y(i + 2) = X(i + 2) Y(i + 3) = X(i + 3) Y(i + 4) = X(i + 4) Y(i + 5) = X(i + 5) Y(i + 6) = X(i + 6) Y(i + 7) = X(i + 7) END DO END SUBROUTINE WCOPY !-------------------------------------------------------------- SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY) !-------------------------------------------------------------- ! constant times a vector plus a vector: y <- y + Alpha*x ! only for incX=incY=1 ! after BLAS ! replace this by the function from the optimized BLAS implementation: ! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1) !-------------------------------------------------------------- INTEGER :: i,incX,incY,M,MP1,N REAL(kind=dp) :: X(N),Y(N),Alpha REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp IF (Alpha .EQ. ZERO) RETURN IF (N .LE. 0) RETURN M = MOD(N,4) IF( M .NE. 0 ) THEN DO i = 1,M Y(i) = Y(i) + Alpha*X(i) END DO IF( N .LT. 4 ) RETURN END IF MP1 = M + 1 DO i = MP1,N,4 Y(i) = Y(i) + Alpha*X(i) Y(i + 1) = Y(i + 1) + Alpha*X(i + 1) Y(i + 2) = Y(i + 2) + Alpha*X(i + 2) Y(i + 3) = Y(i + 3) + Alpha*X(i + 3) END DO END SUBROUTINE WAXPY !-------------------------------------------------------------- SUBROUTINE WSCAL(N,Alpha,X,incX) !-------------------------------------------------------------- ! constant times a vector: x(1:N) <- Alpha*x(1:N) ! only for incX=incY=1 ! after BLAS ! replace this by the function from the optimized BLAS implementation: ! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1) !-------------------------------------------------------------- INTEGER :: i,incX,M,MP1,N REAL(kind=dp) :: X(N),Alpha REAL(kind=dp), PARAMETER :: ZERO=0.0_dp, ONE=1.0_dp IF (Alpha .EQ. ONE) RETURN IF (N .LE. 0) RETURN M = MOD(N,5) IF( M .NE. 0 ) THEN IF (Alpha .EQ. (-ONE)) THEN DO i = 1,M X(i) = -X(i) END DO ELSEIF (Alpha .EQ. ZERO) THEN DO i = 1,M X(i) = ZERO END DO ELSE DO i = 1,M X(i) = Alpha*X(i) END DO END IF IF( N .LT. 5 ) RETURN END IF MP1 = M + 1 IF (Alpha .EQ. (-ONE)) THEN DO i = MP1,N,5 X(i) = -X(i) X(i + 1) = -X(i + 1) X(i + 2) = -X(i + 2) X(i + 3) = -X(i + 3) X(i + 4) = -X(i + 4) END DO ELSEIF (Alpha .EQ. ZERO) THEN DO i = MP1,N,5 X(i) = ZERO X(i + 1) = ZERO X(i + 2) = ZERO X(i + 3) = ZERO X(i + 4) = ZERO END DO ELSE DO i = MP1,N,5 X(i) = Alpha*X(i) X(i + 1) = Alpha*X(i + 1) X(i + 2) = Alpha*X(i + 2) X(i + 3) = Alpha*X(i + 3) X(i + 4) = Alpha*X(i + 4) END DO END IF END SUBROUTINE WSCAL !-------------------------------------------------------------- REAL(kind=dp) FUNCTION WLAMCH( C ) !-------------------------------------------------------------- ! returns epsilon machine ! after LAPACK ! replace this by the function from the optimized LAPACK implementation: ! CALL SLAMCH('E') or CALL DLAMCH('E') !-------------------------------------------------------------- ! USE gckpp_Precision CHARACTER :: C INTEGER :: i REAL(kind=dp), SAVE :: Eps REAL(kind=dp) :: Suma REAL(kind=dp), PARAMETER :: ONE=1.0_dp, HALF=0.5_dp LOGICAL, SAVE :: First=.TRUE. !$OMP THREADPRIVATE( Eps, First ) IF (First) THEN First = .FALSE. Eps = HALF**(16) DO i = 17, 80 Eps = Eps*HALF CALL WLAMCH_ADD(ONE,Eps,Suma) IF (Suma.LE.ONE) GOTO 10 END DO PRINT*,'ERROR IN WLAMCH. EPS < ',Eps RETURN 10 Eps = Eps*2 i = i-1 END IF WLAMCH = Eps END FUNCTION WLAMCH SUBROUTINE WLAMCH_ADD( A, B, Suma ) ! USE gckpp_Precision REAL(kind=dp) A, B, Suma Suma = A + B END SUBROUTINE WLAMCH_ADD !-------------------------------------------------------------- !-------------------------------------------------------------- SUBROUTINE SET2ZERO(N,Y) !-------------------------------------------------------------- ! copies zeros into the vector y: y <- 0 ! after BLAS !-------------------------------------------------------------- INTEGER :: i,M,MP1,N REAL(kind=dp) :: Y(N) REAL(kind=dp), PARAMETER :: ZERO = 0.0d0 IF (N.LE.0) RETURN M = MOD(N,8) IF( M .NE. 0 ) THEN DO i = 1,M Y(i) = ZERO END DO IF( N .LT. 8 ) RETURN END IF MP1 = M+1 DO i = MP1,N,8 Y(i) = ZERO Y(i + 1) = ZERO Y(i + 2) = ZERO Y(i + 3) = ZERO Y(i + 4) = ZERO Y(i + 5) = ZERO Y(i + 6) = ZERO Y(i + 7) = ZERO END DO END SUBROUTINE SET2ZERO !-------------------------------------------------------------- REAL(kind=dp) FUNCTION WDOT (N, DX, incX, DY, incY) !-------------------------------------------------------------- ! dot produce: wdot = x(1:N)*y(1:N) ! only for incX=incY=1 ! after BLAS ! replace this by the function from the optimized BLAS implementation: ! CALL SDOT(N,X,1,Y,1) or CALL DDOT(N,X,1,Y,1) !-------------------------------------------------------------- ! USE messy_mecca_kpp_Precision !-------------------------------------------------------------- IMPLICIT NONE INTEGER :: N, incX, incY REAL(kind=dp) :: DX(N), DY(N) INTEGER :: i, IX, IY, M, MP1, NS WDOT = 0.0D0 IF (N .LE. 0) RETURN IF (incX .EQ. incY) IF (incX-1) 5,20,60 ! ! Code for unequal or nonpositive increments. ! 5 IX = 1 IY = 1 IF (incX .LT. 0) IX = (-N+1)*incX + 1 IF (incY .LT. 0) IY = (-N+1)*incY + 1 DO i = 1,N WDOT = WDOT + DX(IX)*DY(IY) IX = IX + incX IY = IY + incY END DO RETURN ! ! Code for both increments equal to 1. ! ! Clean-up loop so remaining vector length is a multiple of 5. ! 20 M = MOD(N,5) IF (M .EQ. 0) GO TO 40 DO i = 1,M WDOT = WDOT + DX(i)*DY(i) END DO IF (N .LT. 5) RETURN 40 MP1 = M + 1 DO i = MP1,N,5 WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) + & DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4) END DO RETURN ! ! Code for equal, positive, non-unit increments. ! 60 NS = N*incX DO i = 1,NS,incX WDOT = WDOT + DX(i)*DY(i) END DO END FUNCTION WDOT !-------------------------------------------------------------- SUBROUTINE WADD(N,X,Y,Z) !-------------------------------------------------------------- ! adds two vectors: z <- x + y ! BLAS - like !-------------------------------------------------------------- ! USE gckpp_Precision INTEGER :: i, M, MP1, N REAL(kind=dp) :: X(N),Y(N),Z(N) IF (N.LE.0) RETURN M = MOD(N,5) IF( M /= 0 ) THEN DO i = 1,M Z(i) = X(i) + Y(i) END DO IF( N < 5 ) RETURN END IF MP1 = M+1 DO i = MP1,N,5 Z(i) = X(i) + Y(i) Z(i + 1) = X(i + 1) + Y(i + 1) Z(i + 2) = X(i + 2) + Y(i + 2) Z(i + 3) = X(i + 3) + Y(i + 3) Z(i + 4) = X(i + 4) + Y(i + 4) END DO END SUBROUTINE WADD !-------------------------------------------------------------- SUBROUTINE WGEFA(N,A,Ipvt,info) !-------------------------------------------------------------- ! WGEFA FACTORS THE MATRIX A (N,N) BY ! GAUSS ELIMINATION WITH PARTIAL PIVOTING ! LINPACK - LIKE !-------------------------------------------------------------- ! INTEGER :: N,Ipvt(N),info REAL(kind=dp) :: A(N,N) REAL(kind=dp) :: t, dmax, da INTEGER :: j,k,l REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0 info = 0 size: IF (n > 1) THEN col: DO k = 1, n-1 ! find l = pivot index ! l = idamax(n-k+1,A(k,k),1) + k - 1 l = k; dmax = abs(A(k,k)) DO j = k+1,n da = ABS(A(j,k)) IF (da > dmax) THEN l = j; dmax = da END IF END DO Ipvt(k) = l ! zero pivot implies this column already triangularized IF (ABS(A(l,k)) < TINY(ZERO)) THEN info = k return ELSE IF (l /= k) THEN t = A(l,k); A(l,k) = A(k,k); A(k,k) = t END IF t = -ONE/A(k,k) CALL WSCAL(n-k,t,A(k+1,k),1) DO j = k+1, n t = A(l,j) IF (l /= k) THEN A(l,j) = A(k,j); A(k,j) = t END IF CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) END DO END IF END DO col END IF size Ipvt(N) = N IF (ABS(A(N,N)) == ZERO) info = N END SUBROUTINE WGEFA !-------------------------------------------------------------- SUBROUTINE WGESL(Trans,N,A,Ipvt,b) !-------------------------------------------------------------- ! WGESL solves the system ! a * x = b or trans(a) * x = b ! using the factors computed by WGEFA. ! ! Trans = 'N' to solve A*x = b , ! = 'T' to solve transpose(A)*x = b ! LINPACK - LIKE !-------------------------------------------------------------- INTEGER :: N,Ipvt(N) CHARACTER :: trans REAL(kind=dp) :: A(N,N),b(N) REAL(kind=dp) :: t INTEGER :: k,kb,l SELECT CASE (Trans) CASE ('n','N') ! Solve A * x = b ! first solve L*y = b IF (n >= 2) THEN DO k = 1, n-1 l = Ipvt(k) t = b(l) IF (l /= k) THEN b(l) = b(k) b(k) = t END IF CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) END DO END IF ! now solve U*x = y DO kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) CALL WAXPY(k-1,t,a(1,k),1,b(1),1) END DO CASE ('t','T') ! Solve transpose(A) * x = b ! first solve trans(U)*y = b DO k = 1, n t = WDOT(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) END DO ! now solve trans(L)*x = y IF (n >= 2) THEN DO kb = 1, n-1 k = n - kb b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1) l = Ipvt(k) IF (l /= k) THEN t = b(l); b(l) = b(k); b(k) = t END IF END DO END IF END SELECT END SUBROUTINE WGESL ! End of BLAS_UTIL function ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ END MODULE gckpp_LinearAlgebra