1631 lines
71 KiB
Fortran
1631 lines
71 KiB
Fortran
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
!
|
|
! 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_adj_LinearAlgebra.f90
|
|
! Time : Tue May 14 19:43:54 2013
|
|
! Working directory : /home/daven/kpp-2.2.1/GC_KPP
|
|
! Equation file : gckpp_adj.kpp
|
|
! Output root filename : gckpp_adj
|
|
!
|
|
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
|
|
|
|
MODULE gckpp_adj_LinearAlgebra
|
|
|
|
USE gckpp_adj_Parameters
|
|
USE gckpp_adj_JacobianSP
|
|
|
|
IMPLICIT NONE
|
|
|
|
CONTAINS
|
|
|
|
|
|
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
!
|
|
! SPARSE_UTIL - SPARSE utility functions
|
|
! Arguments :
|
|
!
|
|
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
|
|
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
SUBROUTINE KppDecomp( JVS, IER )
|
|
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
! Sparse LU factorization
|
|
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
USE gckpp_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
USE gckpp_adj_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_adj_Parameters
|
|
! USE gckpp_adj_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_adj_Parameters
|
|
! USE gckpp_adj_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_adj_Parameters
|
|
! USE gckpp_adj_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(47) = X(47)-JVS(167)*X(40)-JVS(168)*X(41)-JVS(169)*X(46)
|
|
X(48) = X(48)-JVS(190)*X(46)
|
|
X(49) = X(49)-JVS(196)*X(32)
|
|
X(50) = X(50)-JVS(206)*X(41)-JVS(207)*X(42)
|
|
X(52) = X(52)-JVS(225)*X(16)-JVS(226)*X(31)
|
|
X(53) = X(53)-JVS(235)*X(16)-JVS(236)*X(32)
|
|
X(56) = X(56)-JVS(259)*X(22)-JVS(260)*X(23)-JVS(261)*X(25)-JVS(262)*X(34)-JVS(263)*X(51)
|
|
X(58) = X(58)-JVS(291)*X(45)-JVS(292)*X(57)
|
|
X(59) = X(59)-JVS(301)*X(40)
|
|
X(61) = X(61)-JVS(317)*X(40)-JVS(318)*X(41)-JVS(319)*X(44)-JVS(320)*X(58)-JVS(321)*X(59)
|
|
X(62) = X(62)-JVS(337)*X(49)-JVS(338)*X(53)
|
|
X(63) = X(63)-JVS(347)*X(38)-JVS(348)*X(48)
|
|
X(64) = X(64)-JVS(358)*X(42)
|
|
X(65) = X(65)-JVS(366)*X(43)-JVS(367)*X(51)-JVS(368)*X(55)
|
|
X(66) = X(66)-JVS(377)*X(41)-JVS(378)*X(58)
|
|
X(67) = X(67)-JVS(387)*X(39)-JVS(388)*X(46)
|
|
X(68) = X(68)-JVS(397)*X(30)-JVS(398)*X(31)-JVS(399)*X(38)-JVS(400)*X(39)-JVS(401)*X(42)-JVS(402)*X(43)-JVS(403)*X(51)&
|
|
&-JVS(404)*X(52)-JVS(405)*X(54)-JVS(406)*X(55)-JVS(407)*X(63)-JVS(408)*X(64)-JVS(409)*X(65)-JVS(410)*X(67)
|
|
X(69) = X(69)-JVS(424)*X(26)-JVS(425)*X(27)-JVS(426)*X(28)-JVS(427)*X(34)-JVS(428)*X(36)-JVS(429)*X(37)-JVS(430)*X(40)&
|
|
&-JVS(431)*X(42)-JVS(432)*X(44)-JVS(433)*X(45)-JVS(434)*X(46)-JVS(435)*X(48)-JVS(436)*X(50)-JVS(437)*X(51)&
|
|
&-JVS(438)*X(52)-JVS(439)*X(53)-JVS(440)*X(54)-JVS(441)*X(55)-JVS(442)*X(57)-JVS(443)*X(58)-JVS(444)*X(59)&
|
|
&-JVS(445)*X(60)-JVS(446)*X(61)-JVS(447)*X(62)-JVS(448)*X(63)-JVS(449)*X(64)-JVS(450)*X(65)-JVS(451)*X(66)&
|
|
&-JVS(452)*X(67)
|
|
X(70) = X(70)-JVS(474)*X(37)-JVS(475)*X(48)
|
|
X(71) = X(71)-JVS(485)*X(33)-JVS(486)*X(35)-JVS(487)*X(37)-JVS(488)*X(48)-JVS(489)*X(54)-JVS(490)*X(60)-JVS(491)*X(63)&
|
|
&-JVS(492)*X(70)
|
|
X(72) = X(72)-JVS(508)*X(22)-JVS(509)*X(30)-JVS(510)*X(54)
|
|
X(73) = X(73)-JVS(519)*X(54)-JVS(520)*X(60)-JVS(521)*X(62)-JVS(522)*X(72)
|
|
X(74) = X(74)-JVS(531)*X(23)-JVS(532)*X(35)-JVS(533)*X(68)-JVS(534)*X(70)-JVS(535)*X(72)-JVS(536)*X(73)
|
|
X(75) = X(75)-JVS(549)*X(41)-JVS(550)*X(42)-JVS(551)*X(51)-JVS(552)*X(55)-JVS(553)*X(58)-JVS(554)*X(59)-JVS(555)*X(60)&
|
|
&-JVS(556)*X(61)-JVS(557)*X(62)-JVS(558)*X(64)-JVS(559)*X(65)-JVS(560)*X(66)-JVS(561)*X(70)-JVS(562)*X(72)&
|
|
&-JVS(563)*X(73)
|
|
X(76) = X(76)-JVS(578)*X(57)-JVS(579)*X(59)-JVS(580)*X(60)-JVS(581)*X(62)-JVS(582)*X(64)-JVS(583)*X(66)-JVS(584)*X(72)&
|
|
&-JVS(585)*X(73)
|
|
X(77) = X(77)-JVS(597)*X(45)-JVS(598)*X(46)-JVS(599)*X(57)-JVS(600)*X(67)
|
|
X(78) = X(78)-JVS(610)*X(26)-JVS(611)*X(44)-JVS(612)*X(58)-JVS(613)*X(77)
|
|
X(79) = X(79)-JVS(625)*X(45)-JVS(626)*X(46)-JVS(627)*X(57)
|
|
X(80) = X(80)-JVS(636)*X(45)-JVS(637)*X(46)-JVS(638)*X(57)-JVS(639)*X(67)-JVS(640)*X(79)
|
|
X(81) = X(81)-JVS(649)*X(18)-JVS(650)*X(33)-JVS(651)*X(68)-JVS(652)*X(70)-JVS(653)*X(72)-JVS(654)*X(73)-JVS(655)*X(76)&
|
|
&-JVS(656)*X(77)-JVS(657)*X(79)-JVS(658)*X(80)
|
|
X(82) = X(82)-JVS(668)*X(18)-JVS(669)*X(19)-JVS(670)*X(21)-JVS(671)*X(24)-JVS(672)*X(25)-JVS(673)*X(29)-JVS(674)*X(38)&
|
|
&-JVS(675)*X(39)-JVS(676)*X(43)-JVS(677)*X(44)-JVS(678)*X(51)-JVS(679)*X(52)-JVS(680)*X(53)-JVS(681)*X(54)&
|
|
&-JVS(682)*X(55)-JVS(683)*X(56)-JVS(684)*X(57)-JVS(685)*X(59)-JVS(686)*X(60)-JVS(687)*X(62)-JVS(688)*X(63)&
|
|
&-JVS(689)*X(64)-JVS(690)*X(65)-JVS(691)*X(66)-JVS(692)*X(67)-JVS(693)*X(68)-JVS(694)*X(69)-JVS(695)*X(70)&
|
|
&-JVS(696)*X(71)-JVS(697)*X(72)-JVS(698)*X(73)-JVS(699)*X(74)-JVS(700)*X(75)-JVS(701)*X(76)-JVS(702)*X(77)&
|
|
&-JVS(703)*X(78)-JVS(704)*X(79)-JVS(705)*X(80)-JVS(706)*X(81)
|
|
X(83) = X(83)-JVS(716)*X(16)-JVS(717)*X(17)-JVS(718)*X(20)-JVS(719)*X(22)-JVS(720)*X(23)-JVS(721)*X(24)-JVS(722)*X(26)&
|
|
&-JVS(723)*X(27)-JVS(724)*X(28)-JVS(725)*X(29)-JVS(726)*X(30)-JVS(727)*X(31)-JVS(728)*X(32)-JVS(729)*X(33)&
|
|
&-JVS(730)*X(34)-JVS(731)*X(35)-JVS(732)*X(36)-JVS(733)*X(37)-JVS(734)*X(38)-JVS(735)*X(39)-JVS(736)*X(40)&
|
|
&-JVS(737)*X(41)-JVS(738)*X(42)-JVS(739)*X(43)-JVS(740)*X(44)-JVS(741)*X(45)-JVS(742)*X(46)-JVS(743)*X(47)&
|
|
&-JVS(744)*X(48)-JVS(745)*X(49)-JVS(746)*X(50)-JVS(747)*X(51)-JVS(748)*X(52)-JVS(749)*X(53)-JVS(750)*X(55)&
|
|
&-JVS(751)*X(56)-JVS(752)*X(57)-JVS(753)*X(58)-JVS(754)*X(59)-JVS(755)*X(61)-JVS(756)*X(62)-JVS(757)*X(63)&
|
|
&-JVS(758)*X(64)-JVS(759)*X(65)-JVS(760)*X(66)-JVS(761)*X(67)-JVS(762)*X(68)-JVS(763)*X(69)-JVS(764)*X(70)&
|
|
&-JVS(765)*X(71)-JVS(766)*X(72)-JVS(767)*X(73)-JVS(768)*X(74)-JVS(769)*X(75)-JVS(770)*X(76)-JVS(771)*X(77)&
|
|
&-JVS(772)*X(78)-JVS(773)*X(79)-JVS(774)*X(80)-JVS(775)*X(81)-JVS(776)*X(82)
|
|
X(84) = X(84)-JVS(785)*X(17)-JVS(786)*X(20)-JVS(787)*X(28)-JVS(788)*X(29)-JVS(789)*X(30)-JVS(790)*X(31)-JVS(791)*X(32)&
|
|
&-JVS(792)*X(33)-JVS(793)*X(34)-JVS(794)*X(35)-JVS(795)*X(36)-JVS(796)*X(37)-JVS(797)*X(38)-JVS(798)*X(39)&
|
|
&-JVS(799)*X(40)-JVS(800)*X(41)-JVS(801)*X(42)-JVS(802)*X(43)-JVS(803)*X(44)-JVS(804)*X(45)-JVS(805)*X(46)&
|
|
&-JVS(806)*X(47)-JVS(807)*X(48)-JVS(808)*X(49)-JVS(809)*X(50)-JVS(810)*X(51)-JVS(811)*X(52)-JVS(812)*X(53)&
|
|
&-JVS(813)*X(54)-JVS(814)*X(55)-JVS(815)*X(57)-JVS(816)*X(58)-JVS(817)*X(59)-JVS(818)*X(60)-JVS(819)*X(61)&
|
|
&-JVS(820)*X(62)-JVS(821)*X(63)-JVS(822)*X(64)-JVS(823)*X(65)-JVS(824)*X(66)-JVS(825)*X(67)-JVS(826)*X(68)&
|
|
&-JVS(827)*X(69)-JVS(828)*X(70)-JVS(829)*X(71)-JVS(830)*X(72)-JVS(831)*X(73)-JVS(832)*X(74)-JVS(833)*X(75)&
|
|
&-JVS(834)*X(76)-JVS(835)*X(77)-JVS(836)*X(78)-JVS(837)*X(79)-JVS(838)*X(80)-JVS(839)*X(81)-JVS(840)*X(82)&
|
|
&-JVS(841)*X(83)
|
|
X(85) = X(85)-JVS(849)*X(24)-JVS(850)*X(25)-JVS(851)*X(51)-JVS(852)*X(52)-JVS(853)*X(53)-JVS(854)*X(54)-JVS(855)*X(55)&
|
|
&-JVS(856)*X(57)-JVS(857)*X(59)-JVS(858)*X(60)-JVS(859)*X(62)-JVS(860)*X(63)-JVS(861)*X(64)-JVS(862)*X(65)&
|
|
&-JVS(863)*X(66)-JVS(864)*X(67)-JVS(865)*X(70)-JVS(866)*X(72)-JVS(867)*X(73)-JVS(868)*X(74)-JVS(869)*X(76)&
|
|
&-JVS(870)*X(77)-JVS(871)*X(78)-JVS(872)*X(79)-JVS(873)*X(80)-JVS(874)*X(81)-JVS(875)*X(82)-JVS(876)*X(83)&
|
|
&-JVS(877)*X(84)
|
|
X(86) = X(86)-JVS(884)*X(21)-JVS(885)*X(26)-JVS(886)*X(27)-JVS(887)*X(42)-JVS(888)*X(49)-JVS(889)*X(51)-JVS(890)*X(52)&
|
|
&-JVS(891)*X(53)-JVS(892)*X(54)-JVS(893)*X(55)-JVS(894)*X(57)-JVS(895)*X(59)-JVS(896)*X(60)-JVS(897)*X(61)&
|
|
&-JVS(898)*X(62)-JVS(899)*X(63)-JVS(900)*X(64)-JVS(901)*X(65)-JVS(902)*X(66)-JVS(903)*X(67)-JVS(904)*X(70)&
|
|
&-JVS(905)*X(71)-JVS(906)*X(72)-JVS(907)*X(73)-JVS(908)*X(74)-JVS(909)*X(75)-JVS(910)*X(76)-JVS(911)*X(77)&
|
|
&-JVS(912)*X(78)-JVS(913)*X(79)-JVS(914)*X(80)-JVS(915)*X(81)-JVS(916)*X(82)-JVS(917)*X(83)-JVS(918)*X(84)&
|
|
&-JVS(919)*X(85)
|
|
X(87) = X(87)-JVS(925)*X(21)-JVS(926)*X(22)-JVS(927)*X(23)-JVS(928)*X(25)-JVS(929)*X(29)-JVS(930)*X(34)-JVS(931)*X(46)&
|
|
&-JVS(932)*X(48)-JVS(933)*X(56)-JVS(934)*X(57)-JVS(935)*X(59)-JVS(936)*X(64)-JVS(937)*X(65)-JVS(938)*X(66)&
|
|
&-JVS(939)*X(67)-JVS(940)*X(68)-JVS(941)*X(69)-JVS(942)*X(70)-JVS(943)*X(71)-JVS(944)*X(72)-JVS(945)*X(73)&
|
|
&-JVS(946)*X(74)-JVS(947)*X(75)-JVS(948)*X(76)-JVS(949)*X(77)-JVS(950)*X(78)-JVS(951)*X(79)-JVS(952)*X(80)&
|
|
&-JVS(953)*X(81)-JVS(954)*X(82)-JVS(955)*X(83)-JVS(956)*X(84)-JVS(957)*X(85)-JVS(958)*X(86)
|
|
X(88) = X(88)-JVS(963)*X(19)-JVS(964)*X(36)-JVS(965)*X(50)-JVS(966)*X(51)-JVS(967)*X(58)-JVS(968)*X(64)-JVS(969)*X(65)&
|
|
&-JVS(970)*X(66)-JVS(971)*X(77)-JVS(972)*X(79)-JVS(973)*X(80)-JVS(974)*X(82)-JVS(975)*X(83)-JVS(976)*X(84)&
|
|
&-JVS(977)*X(85)-JVS(978)*X(86)-JVS(979)*X(87)
|
|
X(89) = X(89)-JVS(983)*X(25)-JVS(984)*X(44)-JVS(985)*X(46)-JVS(986)*X(48)-JVS(987)*X(58)-JVS(988)*X(77)-JVS(989)*X(78)&
|
|
&-JVS(990)*X(79)-JVS(991)*X(80)-JVS(992)*X(81)-JVS(993)*X(82)-JVS(994)*X(83)-JVS(995)*X(84)-JVS(996)*X(85)&
|
|
&-JVS(997)*X(86)-JVS(998)*X(87)-JVS(999)*X(88)
|
|
X(90) = X(90)-JVS(1002)*X(21)-JVS(1003)*X(27)-JVS(1004)*X(28)-JVS(1005)*X(34)-JVS(1006)*X(48)-JVS(1007)*X(49)&
|
|
&-JVS(1008)*X(51)-JVS(1009)*X(52)-JVS(1010)*X(53)-JVS(1011)*X(54)-JVS(1012)*X(55)-JVS(1013)*X(57)-JVS(1014)&
|
|
&*X(59)-JVS(1015)*X(60)-JVS(1016)*X(62)-JVS(1017)*X(63)-JVS(1018)*X(64)-JVS(1019)*X(65)-JVS(1020)*X(66)&
|
|
&-JVS(1021)*X(67)-JVS(1022)*X(70)-JVS(1023)*X(71)-JVS(1024)*X(72)-JVS(1025)*X(73)-JVS(1026)*X(74)-JVS(1027)&
|
|
&*X(75)-JVS(1028)*X(76)-JVS(1029)*X(77)-JVS(1030)*X(78)-JVS(1031)*X(79)-JVS(1032)*X(80)-JVS(1033)*X(81)&
|
|
&-JVS(1034)*X(82)-JVS(1035)*X(83)-JVS(1036)*X(84)-JVS(1037)*X(85)-JVS(1038)*X(86)-JVS(1039)*X(87)-JVS(1040)&
|
|
&*X(88)-JVS(1041)*X(89)
|
|
X(90) = X(90)/JVS(1042)
|
|
X(89) = (X(89)-JVS(1001)*X(90))/(JVS(1000))
|
|
X(88) = (X(88)-JVS(981)*X(89)-JVS(982)*X(90))/(JVS(980))
|
|
X(87) = (X(87)-JVS(960)*X(88)-JVS(961)*X(89)-JVS(962)*X(90))/(JVS(959))
|
|
X(86) = (X(86)-JVS(921)*X(87)-JVS(922)*X(88)-JVS(923)*X(89)-JVS(924)*X(90))/(JVS(920))
|
|
X(85) = (X(85)-JVS(879)*X(86)-JVS(880)*X(87)-JVS(881)*X(88)-JVS(882)*X(89)-JVS(883)*X(90))/(JVS(878))
|
|
X(84) = (X(84)-JVS(843)*X(85)-JVS(844)*X(86)-JVS(845)*X(87)-JVS(846)*X(88)-JVS(847)*X(89)-JVS(848)*X(90))/(JVS(842))
|
|
X(83) = (X(83)-JVS(778)*X(84)-JVS(779)*X(85)-JVS(780)*X(86)-JVS(781)*X(87)-JVS(782)*X(88)-JVS(783)*X(89)-JVS(784)&
|
|
&*X(90))/(JVS(777))
|
|
X(82) = (X(82)-JVS(708)*X(83)-JVS(709)*X(84)-JVS(710)*X(85)-JVS(711)*X(86)-JVS(712)*X(87)-JVS(713)*X(88)-JVS(714)&
|
|
&*X(89)-JVS(715)*X(90))/(JVS(707))
|
|
X(81) = (X(81)-JVS(660)*X(82)-JVS(661)*X(83)-JVS(662)*X(84)-JVS(663)*X(85)-JVS(664)*X(86)-JVS(665)*X(87)-JVS(666)&
|
|
&*X(89)-JVS(667)*X(90))/(JVS(659))
|
|
X(80) = (X(80)-JVS(642)*X(83)-JVS(643)*X(84)-JVS(644)*X(85)-JVS(645)*X(86)-JVS(646)*X(87)-JVS(647)*X(89)-JVS(648)&
|
|
&*X(90))/(JVS(641))
|
|
X(79) = (X(79)-JVS(629)*X(83)-JVS(630)*X(84)-JVS(631)*X(85)-JVS(632)*X(86)-JVS(633)*X(87)-JVS(634)*X(89)-JVS(635)&
|
|
&*X(90))/(JVS(628))
|
|
X(78) = (X(78)-JVS(615)*X(79)-JVS(616)*X(80)-JVS(617)*X(82)-JVS(618)*X(83)-JVS(619)*X(84)-JVS(620)*X(85)-JVS(621)&
|
|
&*X(86)-JVS(622)*X(87)-JVS(623)*X(89)-JVS(624)*X(90))/(JVS(614))
|
|
X(77) = (X(77)-JVS(602)*X(79)-JVS(603)*X(83)-JVS(604)*X(84)-JVS(605)*X(85)-JVS(606)*X(86)-JVS(607)*X(87)-JVS(608)&
|
|
&*X(89)-JVS(609)*X(90))/(JVS(601))
|
|
X(76) = (X(76)-JVS(587)*X(77)-JVS(588)*X(79)-JVS(589)*X(80)-JVS(590)*X(83)-JVS(591)*X(84)-JVS(592)*X(85)-JVS(593)&
|
|
&*X(86)-JVS(594)*X(87)-JVS(595)*X(89)-JVS(596)*X(90))/(JVS(586))
|
|
X(75) = (X(75)-JVS(565)*X(76)-JVS(566)*X(77)-JVS(567)*X(78)-JVS(568)*X(79)-JVS(569)*X(80)-JVS(570)*X(82)-JVS(571)&
|
|
&*X(83)-JVS(572)*X(84)-JVS(573)*X(85)-JVS(574)*X(86)-JVS(575)*X(87)-JVS(576)*X(89)-JVS(577)*X(90))/(JVS(564))
|
|
X(74) = (X(74)-JVS(538)*X(76)-JVS(539)*X(77)-JVS(540)*X(80)-JVS(541)*X(81)-JVS(542)*X(83)-JVS(543)*X(84)-JVS(544)&
|
|
&*X(85)-JVS(545)*X(86)-JVS(546)*X(87)-JVS(547)*X(89)-JVS(548)*X(90))/(JVS(537))
|
|
X(73) = (X(73)-JVS(524)*X(76)-JVS(525)*X(83)-JVS(526)*X(84)-JVS(527)*X(85)-JVS(528)*X(86)-JVS(529)*X(87)-JVS(530)&
|
|
&*X(90))/(JVS(523))
|
|
X(72) = (X(72)-JVS(512)*X(73)-JVS(513)*X(83)-JVS(514)*X(84)-JVS(515)*X(85)-JVS(516)*X(86)-JVS(517)*X(87)-JVS(518)&
|
|
&*X(90))/(JVS(511))
|
|
X(71) = (X(71)-JVS(494)*X(72)-JVS(495)*X(73)-JVS(496)*X(74)-JVS(497)*X(75)-JVS(498)*X(76)-JVS(499)*X(77)-JVS(500)&
|
|
&*X(81)-JVS(501)*X(83)-JVS(502)*X(84)-JVS(503)*X(85)-JVS(504)*X(86)-JVS(505)*X(87)-JVS(506)*X(89)-JVS(507)&
|
|
&*X(90))/(JVS(493))
|
|
X(70) = (X(70)-JVS(477)*X(77)-JVS(478)*X(83)-JVS(479)*X(84)-JVS(480)*X(85)-JVS(481)*X(86)-JVS(482)*X(87)-JVS(483)&
|
|
&*X(89)-JVS(484)*X(90))/(JVS(476))
|
|
X(69) = (X(69)-JVS(454)*X(70)-JVS(455)*X(71)-JVS(456)*X(72)-JVS(457)*X(73)-JVS(458)*X(74)-JVS(459)*X(76)-JVS(460)&
|
|
&*X(77)-JVS(461)*X(78)-JVS(462)*X(79)-JVS(463)*X(80)-JVS(464)*X(81)-JVS(465)*X(82)-JVS(466)*X(83)-JVS(467)*X(84)&
|
|
&-JVS(468)*X(85)-JVS(469)*X(86)-JVS(470)*X(87)-JVS(471)*X(88)-JVS(472)*X(89)-JVS(473)*X(90))/(JVS(453))
|
|
X(68) = (X(68)-JVS(412)*X(70)-JVS(413)*X(72)-JVS(414)*X(73)-JVS(415)*X(77)-JVS(416)*X(80)-JVS(417)*X(83)-JVS(418)&
|
|
&*X(84)-JVS(419)*X(85)-JVS(420)*X(86)-JVS(421)*X(87)-JVS(422)*X(89)-JVS(423)*X(90))/(JVS(411))
|
|
X(67) = (X(67)-JVS(390)*X(83)-JVS(391)*X(84)-JVS(392)*X(85)-JVS(393)*X(86)-JVS(394)*X(87)-JVS(395)*X(89)-JVS(396)&
|
|
&*X(90))/(JVS(389))
|
|
X(66) = (X(66)-JVS(380)*X(79)-JVS(381)*X(83)-JVS(382)*X(84)-JVS(383)*X(85)-JVS(384)*X(86)-JVS(385)*X(89)-JVS(386)&
|
|
&*X(90))/(JVS(379))
|
|
X(65) = (X(65)-JVS(370)*X(80)-JVS(371)*X(83)-JVS(372)*X(84)-JVS(373)*X(85)-JVS(374)*X(86)-JVS(375)*X(87)-JVS(376)&
|
|
&*X(90))/(JVS(369))
|
|
X(64) = (X(64)-JVS(360)*X(77)-JVS(361)*X(83)-JVS(362)*X(84)-JVS(363)*X(85)-JVS(364)*X(86)-JVS(365)*X(90))/(JVS(359))
|
|
X(63) = (X(63)-JVS(350)*X(77)-JVS(351)*X(83)-JVS(352)*X(84)-JVS(353)*X(85)-JVS(354)*X(86)-JVS(355)*X(87)-JVS(356)&
|
|
&*X(89)-JVS(357)*X(90))/(JVS(349))
|
|
X(62) = (X(62)-JVS(340)*X(72)-JVS(341)*X(73)-JVS(342)*X(83)-JVS(343)*X(84)-JVS(344)*X(85)-JVS(345)*X(86)-JVS(346)&
|
|
&*X(90))/(JVS(339))
|
|
X(61) = (X(61)-JVS(323)*X(62)-JVS(324)*X(65)-JVS(325)*X(66)-JVS(326)*X(70)-JVS(327)*X(78)-JVS(328)*X(79)-JVS(329)&
|
|
&*X(80)-JVS(330)*X(82)-JVS(331)*X(83)-JVS(332)*X(84)-JVS(333)*X(85)-JVS(334)*X(86)-JVS(335)*X(89)-JVS(336)&
|
|
&*X(90))/(JVS(322))
|
|
X(60) = (X(60)-JVS(310)*X(76)-JVS(311)*X(83)-JVS(312)*X(84)-JVS(313)*X(85)-JVS(314)*X(86)-JVS(315)*X(87)-JVS(316)&
|
|
&*X(90))/(JVS(309))
|
|
X(59) = (X(59)-JVS(303)*X(80)-JVS(304)*X(83)-JVS(305)*X(84)-JVS(306)*X(85)-JVS(307)*X(86)-JVS(308)*X(90))/(JVS(302))
|
|
X(58) = (X(58)-JVS(294)*X(79)-JVS(295)*X(83)-JVS(296)*X(84)-JVS(297)*X(85)-JVS(298)*X(86)-JVS(299)*X(89)-JVS(300)&
|
|
&*X(90))/(JVS(293))
|
|
X(57) = (X(57)-JVS(286)*X(79)-JVS(287)*X(84)-JVS(288)*X(85)-JVS(289)*X(86)-JVS(290)*X(90))/(JVS(285))
|
|
X(56) = (X(56)-JVS(265)*X(57)-JVS(266)*X(59)-JVS(267)*X(64)-JVS(268)*X(65)-JVS(269)*X(66)-JVS(270)*X(67)-JVS(271)&
|
|
&*X(68)-JVS(272)*X(69)-JVS(273)*X(71)-JVS(274)*X(75)-JVS(275)*X(76)-JVS(276)*X(79)-JVS(277)*X(80)-JVS(278)*X(82)&
|
|
&-JVS(279)*X(83)-JVS(280)*X(84)-JVS(281)*X(85)-JVS(282)*X(86)-JVS(283)*X(87)-JVS(284)*X(90))/(JVS(264))
|
|
X(55) = (X(55)-JVS(253)*X(80)-JVS(254)*X(84)-JVS(255)*X(85)-JVS(256)*X(86)-JVS(257)*X(87)-JVS(258)*X(90))/(JVS(252))
|
|
X(54) = (X(54)-JVS(246)*X(73)-JVS(247)*X(83)-JVS(248)*X(84)-JVS(249)*X(85)-JVS(250)*X(86)-JVS(251)*X(90))/(JVS(245))
|
|
X(53) = (X(53)-JVS(238)*X(72)-JVS(239)*X(73)-JVS(240)*X(83)-JVS(241)*X(84)-JVS(242)*X(85)-JVS(243)*X(86)-JVS(244)&
|
|
&*X(90))/(JVS(237))
|
|
X(52) = (X(52)-JVS(228)*X(72)-JVS(229)*X(73)-JVS(230)*X(83)-JVS(231)*X(84)-JVS(232)*X(85)-JVS(233)*X(86)-JVS(234)&
|
|
&*X(90))/(JVS(227))
|
|
X(51) = (X(51)-JVS(221)*X(84)-JVS(222)*X(85)-JVS(223)*X(86)-JVS(224)*X(90))/(JVS(220))
|
|
X(50) = (X(50)-JVS(209)*X(51)-JVS(210)*X(58)-JVS(211)*X(64)-JVS(212)*X(65)-JVS(213)*X(66)-JVS(214)*X(83)-JVS(215)&
|
|
&*X(84)-JVS(216)*X(85)-JVS(217)*X(86)-JVS(218)*X(89)-JVS(219)*X(90))/(JVS(208))
|
|
X(49) = (X(49)-JVS(198)*X(53)-JVS(199)*X(72)-JVS(200)*X(73)-JVS(201)*X(83)-JVS(202)*X(84)-JVS(203)*X(85)-JVS(204)&
|
|
&*X(86)-JVS(205)*X(90))/(JVS(197))
|
|
X(48) = (X(48)-JVS(192)*X(77)-JVS(193)*X(83)-JVS(194)*X(87)-JVS(195)*X(89))/(JVS(191))
|
|
X(47) = (X(47)-JVS(171)*X(48)-JVS(172)*X(49)-JVS(173)*X(50)-JVS(174)*X(58)-JVS(175)*X(59)-JVS(176)*X(66)-JVS(177)&
|
|
&*X(68)-JVS(178)*X(69)-JVS(179)*X(71)-JVS(180)*X(75)-JVS(181)*X(77)-JVS(182)*X(80)-JVS(183)*X(83)-JVS(184)*X(84)&
|
|
&-JVS(185)*X(85)-JVS(186)*X(86)-JVS(187)*X(87)-JVS(188)*X(89)-JVS(189)*X(90))/(JVS(170))
|
|
X(46) = (X(46)-JVS(164)*X(83)-JVS(165)*X(87)-JVS(166)*X(89))/(JVS(163))
|
|
X(45) = (X(45)-JVS(159)*X(57)-JVS(160)*X(79)-JVS(161)*X(83)-JVS(162)*X(84))/(JVS(158))
|
|
X(44) = (X(44)-JVS(154)*X(78)-JVS(155)*X(82)-JVS(156)*X(83)-JVS(157)*X(89))/(JVS(153))
|
|
X(43) = (X(43)-JVS(148)*X(51)-JVS(149)*X(55)-JVS(150)*X(65)-JVS(151)*X(83)-JVS(152)*X(84))/(JVS(147))
|
|
X(42) = (X(42)-JVS(144)*X(64)-JVS(145)*X(83)-JVS(146)*X(84))/(JVS(143))
|
|
X(41) = (X(41)-JVS(140)*X(66)-JVS(141)*X(83)-JVS(142)*X(84))/(JVS(139))
|
|
X(40) = (X(40)-JVS(136)*X(59)-JVS(137)*X(83)-JVS(138)*X(84))/(JVS(135))
|
|
X(39) = (X(39)-JVS(132)*X(67)-JVS(133)*X(83)-JVS(134)*X(84))/(JVS(131))
|
|
X(38) = (X(38)-JVS(128)*X(63)-JVS(129)*X(83)-JVS(130)*X(84))/(JVS(127))
|
|
X(37) = (X(37)-JVS(124)*X(70)-JVS(125)*X(83)-JVS(126)*X(84))/(JVS(123))
|
|
X(36) = (X(36)-JVS(120)*X(83)-JVS(121)*X(84)-JVS(122)*X(88))/(JVS(119))
|
|
X(35) = (X(35)-JVS(116)*X(74)-JVS(117)*X(83)-JVS(118)*X(84))/(JVS(115))
|
|
X(34) = (X(34)-JVS(113)*X(83)-JVS(114)*X(87))/(JVS(112))
|
|
X(33) = (X(33)-JVS(109)*X(81)-JVS(110)*X(83)-JVS(111)*X(84))/(JVS(108))
|
|
X(32) = (X(32)-JVS(105)*X(53)-JVS(106)*X(83)-JVS(107)*X(84))/(JVS(104))
|
|
X(31) = (X(31)-JVS(101)*X(52)-JVS(102)*X(83)-JVS(103)*X(84))/(JVS(100))
|
|
X(30) = (X(30)-JVS(97)*X(72)-JVS(98)*X(83)-JVS(99)*X(84))/(JVS(96))
|
|
X(29) = (X(29)-JVS(93)*X(82)-JVS(94)*X(83)-JVS(95)*X(84))/(JVS(92))
|
|
X(28) = (X(28)-JVS(89)*X(83)-JVS(90)*X(84)-JVS(91)*X(90))/(JVS(88))
|
|
X(27) = (X(27)-JVS(85)*X(83)-JVS(86)*X(84)-JVS(87)*X(86))/(JVS(84))
|
|
X(26) = (X(26)-JVS(81)*X(78)-JVS(82)*X(83)-JVS(83)*X(84))/(JVS(80))
|
|
X(25) = (X(25)-JVS(78)*X(82)-JVS(79)*X(87))/(JVS(77))
|
|
X(24) = (X(24)-JVS(74)*X(82)-JVS(75)*X(83)-JVS(76)*X(85))/(JVS(73))
|
|
X(23) = (X(23)-JVS(71)*X(83)-JVS(72)*X(87))/(JVS(70))
|
|
X(22) = (X(22)-JVS(68)*X(83)-JVS(69)*X(87))/(JVS(67))
|
|
X(21) = (X(21)-JVS(65)*X(82)-JVS(66)*X(86))/(JVS(64))
|
|
X(20) = (X(20)-JVS(61)*X(34)-JVS(62)*X(83)-JVS(63)*X(87))/(JVS(60))
|
|
X(19) = (X(19)-JVS(58)*X(82)-JVS(59)*X(88))/(JVS(57))
|
|
X(18) = (X(18)-JVS(55)*X(81)-JVS(56)*X(82))/(JVS(54))
|
|
X(17) = (X(17)-JVS(52)*X(83)-JVS(53)*X(84))/(JVS(51))
|
|
X(16) = (X(16)-JVS(50)*X(83))/(JVS(49))
|
|
X(15) = (X(15)-JVS(47)*X(46)-JVS(48)*X(83))/(JVS(46))
|
|
X(14) = (X(14)-JVS(36)*X(17)-JVS(37)*X(18)-JVS(38)*X(21)-JVS(39)*X(25)-JVS(40)*X(44)-JVS(41)*X(56)-JVS(42)*X(69)&
|
|
&-JVS(43)*X(73)-JVS(44)*X(82)-JVS(45)*X(89))/(JVS(35))
|
|
X(13) = (X(13)-JVS(28)*X(46)-JVS(29)*X(47)-JVS(30)*X(80)-JVS(31)*X(83)-JVS(32)*X(85)-JVS(33)*X(86)-JVS(34)*X(89))&
|
|
&/(JVS(27))
|
|
X(12) = (X(12)-JVS(25)*X(34)-JVS(26)*X(83))/(JVS(24))
|
|
X(11) = (X(11)-JVS(22)*X(20)-JVS(23)*X(83))/(JVS(21))
|
|
X(10) = (X(10)-JVS(20)*X(73))/(JVS(19))
|
|
X(9) = (X(9)-JVS(18)*X(18))/(JVS(17))
|
|
X(8) = (X(8)-JVS(16)*X(44))/(JVS(15))
|
|
X(7) = (X(7)-JVS(14)*X(21))/(JVS(13))
|
|
X(6) = (X(6)-JVS(12)*X(89))/(JVS(11))
|
|
X(5) = (X(5)-JVS(10)*X(82))/(JVS(9))
|
|
X(4) = (X(4)-JVS(8)*X(25))/(JVS(7))
|
|
X(3) = (X(3)-JVS(6)*X(56))/(JVS(5))
|
|
X(2) = (X(2)-JVS(4)*X(17))/(JVS(3))
|
|
X(1) = (X(1)-JVS(2)*X(69))/(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(24)
|
|
XX(13) = X(13)/JVS(27)
|
|
XX(14) = X(14)/JVS(35)
|
|
XX(15) = X(15)/JVS(46)
|
|
XX(16) = X(16)/JVS(49)
|
|
XX(17) = (X(17)-JVS(4)*XX(2)-JVS(36)*XX(14))/(JVS(51))
|
|
XX(18) = (X(18)-JVS(18)*XX(9)-JVS(37)*XX(14))/(JVS(54))
|
|
XX(19) = X(19)/JVS(57)
|
|
XX(20) = (X(20)-JVS(22)*XX(11))/(JVS(60))
|
|
XX(21) = (X(21)-JVS(14)*XX(7)-JVS(38)*XX(14))/(JVS(64))
|
|
XX(22) = X(22)/JVS(67)
|
|
XX(23) = X(23)/JVS(70)
|
|
XX(24) = X(24)/JVS(73)
|
|
XX(25) = (X(25)-JVS(8)*XX(4)-JVS(39)*XX(14))/(JVS(77))
|
|
XX(26) = X(26)/JVS(80)
|
|
XX(27) = X(27)/JVS(84)
|
|
XX(28) = X(28)/JVS(88)
|
|
XX(29) = X(29)/JVS(92)
|
|
XX(30) = X(30)/JVS(96)
|
|
XX(31) = X(31)/JVS(100)
|
|
XX(32) = X(32)/JVS(104)
|
|
XX(33) = X(33)/JVS(108)
|
|
XX(34) = (X(34)-JVS(25)*XX(12)-JVS(61)*XX(20))/(JVS(112))
|
|
XX(35) = X(35)/JVS(115)
|
|
XX(36) = X(36)/JVS(119)
|
|
XX(37) = X(37)/JVS(123)
|
|
XX(38) = X(38)/JVS(127)
|
|
XX(39) = X(39)/JVS(131)
|
|
XX(40) = X(40)/JVS(135)
|
|
XX(41) = X(41)/JVS(139)
|
|
XX(42) = X(42)/JVS(143)
|
|
XX(43) = X(43)/JVS(147)
|
|
XX(44) = (X(44)-JVS(16)*XX(8)-JVS(40)*XX(14))/(JVS(153))
|
|
XX(45) = X(45)/JVS(158)
|
|
XX(46) = (X(46)-JVS(28)*XX(13)-JVS(47)*XX(15))/(JVS(163))
|
|
XX(47) = (X(47)-JVS(29)*XX(13))/(JVS(170))
|
|
XX(48) = (X(48)-JVS(171)*XX(47))/(JVS(191))
|
|
XX(49) = (X(49)-JVS(172)*XX(47))/(JVS(197))
|
|
XX(50) = (X(50)-JVS(173)*XX(47))/(JVS(208))
|
|
XX(51) = (X(51)-JVS(148)*XX(43)-JVS(209)*XX(50))/(JVS(220))
|
|
XX(52) = (X(52)-JVS(101)*XX(31))/(JVS(227))
|
|
XX(53) = (X(53)-JVS(105)*XX(32)-JVS(198)*XX(49))/(JVS(237))
|
|
XX(54) = X(54)/JVS(245)
|
|
XX(55) = (X(55)-JVS(149)*XX(43))/(JVS(252))
|
|
XX(56) = (X(56)-JVS(6)*XX(3)-JVS(41)*XX(14))/(JVS(264))
|
|
XX(57) = (X(57)-JVS(159)*XX(45)-JVS(265)*XX(56))/(JVS(285))
|
|
XX(58) = (X(58)-JVS(174)*XX(47)-JVS(210)*XX(50))/(JVS(293))
|
|
XX(59) = (X(59)-JVS(136)*XX(40)-JVS(175)*XX(47)-JVS(266)*XX(56))/(JVS(302))
|
|
XX(60) = X(60)/JVS(309)
|
|
XX(61) = X(61)/JVS(322)
|
|
XX(62) = (X(62)-JVS(323)*XX(61))/(JVS(339))
|
|
XX(63) = (X(63)-JVS(128)*XX(38))/(JVS(349))
|
|
XX(64) = (X(64)-JVS(144)*XX(42)-JVS(211)*XX(50)-JVS(267)*XX(56))/(JVS(359))
|
|
XX(65) = (X(65)-JVS(150)*XX(43)-JVS(212)*XX(50)-JVS(268)*XX(56)-JVS(324)*XX(61))/(JVS(369))
|
|
XX(66) = (X(66)-JVS(140)*XX(41)-JVS(176)*XX(47)-JVS(213)*XX(50)-JVS(269)*XX(56)-JVS(325)*XX(61))/(JVS(379))
|
|
XX(67) = (X(67)-JVS(132)*XX(39)-JVS(270)*XX(56))/(JVS(389))
|
|
XX(68) = (X(68)-JVS(177)*XX(47)-JVS(271)*XX(56))/(JVS(411))
|
|
XX(69) = (X(69)-JVS(2)*XX(1)-JVS(42)*XX(14)-JVS(178)*XX(47)-JVS(272)*XX(56))/(JVS(453))
|
|
XX(70) = (X(70)-JVS(124)*XX(37)-JVS(326)*XX(61)-JVS(412)*XX(68)-JVS(454)*XX(69))/(JVS(476))
|
|
XX(71) = (X(71)-JVS(179)*XX(47)-JVS(273)*XX(56)-JVS(455)*XX(69))/(JVS(493))
|
|
XX(72) = (X(72)-JVS(97)*XX(30)-JVS(199)*XX(49)-JVS(228)*XX(52)-JVS(238)*XX(53)-JVS(340)*XX(62)-JVS(413)*XX(68)&
|
|
&-JVS(456)*XX(69)-JVS(494)*XX(71))/(JVS(511))
|
|
XX(73) = (X(73)-JVS(20)*XX(10)-JVS(43)*XX(14)-JVS(200)*XX(49)-JVS(229)*XX(52)-JVS(239)*XX(53)-JVS(246)*XX(54)-JVS(341)&
|
|
&*XX(62)-JVS(414)*XX(68)-JVS(457)*XX(69)-JVS(495)*XX(71)-JVS(512)*XX(72))/(JVS(523))
|
|
XX(74) = (X(74)-JVS(116)*XX(35)-JVS(458)*XX(69)-JVS(496)*XX(71))/(JVS(537))
|
|
XX(75) = (X(75)-JVS(180)*XX(47)-JVS(274)*XX(56)-JVS(497)*XX(71))/(JVS(564))
|
|
XX(76) = (X(76)-JVS(275)*XX(56)-JVS(310)*XX(60)-JVS(459)*XX(69)-JVS(498)*XX(71)-JVS(524)*XX(73)-JVS(538)*XX(74)&
|
|
&-JVS(565)*XX(75))/(JVS(586))
|
|
XX(77) = (X(77)-JVS(181)*XX(47)-JVS(192)*XX(48)-JVS(350)*XX(63)-JVS(360)*XX(64)-JVS(415)*XX(68)-JVS(460)*XX(69)&
|
|
&-JVS(477)*XX(70)-JVS(499)*XX(71)-JVS(539)*XX(74)-JVS(566)*XX(75)-JVS(587)*XX(76))/(JVS(601))
|
|
XX(78) = (X(78)-JVS(81)*XX(26)-JVS(154)*XX(44)-JVS(327)*XX(61)-JVS(461)*XX(69)-JVS(567)*XX(75))/(JVS(614))
|
|
XX(79) = (X(79)-JVS(160)*XX(45)-JVS(276)*XX(56)-JVS(286)*XX(57)-JVS(294)*XX(58)-JVS(328)*XX(61)-JVS(380)*XX(66)&
|
|
&-JVS(462)*XX(69)-JVS(568)*XX(75)-JVS(588)*XX(76)-JVS(602)*XX(77)-JVS(615)*XX(78))/(JVS(628))
|
|
XX(80) = (X(80)-JVS(30)*XX(13)-JVS(182)*XX(47)-JVS(253)*XX(55)-JVS(277)*XX(56)-JVS(303)*XX(59)-JVS(329)*XX(61)&
|
|
&-JVS(370)*XX(65)-JVS(416)*XX(68)-JVS(463)*XX(69)-JVS(540)*XX(74)-JVS(569)*XX(75)-JVS(589)*XX(76)-JVS(616)&
|
|
&*XX(78))/(JVS(641))
|
|
XX(81) = (X(81)-JVS(55)*XX(18)-JVS(109)*XX(33)-JVS(464)*XX(69)-JVS(500)*XX(71)-JVS(541)*XX(74))/(JVS(659))
|
|
XX(82) = (X(82)-JVS(10)*XX(5)-JVS(44)*XX(14)-JVS(56)*XX(18)-JVS(58)*XX(19)-JVS(65)*XX(21)-JVS(74)*XX(24)-JVS(78)&
|
|
&*XX(25)-JVS(93)*XX(29)-JVS(155)*XX(44)-JVS(278)*XX(56)-JVS(330)*XX(61)-JVS(465)*XX(69)-JVS(570)*XX(75)&
|
|
&-JVS(617)*XX(78)-JVS(660)*XX(81))/(JVS(707))
|
|
XX(83) = (X(83)-JVS(23)*XX(11)-JVS(26)*XX(12)-JVS(31)*XX(13)-JVS(48)*XX(15)-JVS(50)*XX(16)-JVS(52)*XX(17)-JVS(62)&
|
|
&*XX(20)-JVS(68)*XX(22)-JVS(71)*XX(23)-JVS(75)*XX(24)-JVS(82)*XX(26)-JVS(85)*XX(27)-JVS(89)*XX(28)-JVS(94)&
|
|
&*XX(29)-JVS(98)*XX(30)-JVS(102)*XX(31)-JVS(106)*XX(32)-JVS(110)*XX(33)-JVS(113)*XX(34)-JVS(117)*XX(35)&
|
|
&-JVS(120)*XX(36)-JVS(125)*XX(37)-JVS(129)*XX(38)-JVS(133)*XX(39)-JVS(137)*XX(40)-JVS(141)*XX(41)-JVS(145)&
|
|
&*XX(42)-JVS(151)*XX(43)-JVS(156)*XX(44)-JVS(161)*XX(45)-JVS(164)*XX(46)-JVS(183)*XX(47)-JVS(193)*XX(48)&
|
|
&-JVS(201)*XX(49)-JVS(214)*XX(50)-JVS(230)*XX(52)-JVS(240)*XX(53)-JVS(247)*XX(54)-JVS(279)*XX(56)-JVS(295)&
|
|
&*XX(58)-JVS(304)*XX(59)-JVS(311)*XX(60)-JVS(331)*XX(61)-JVS(342)*XX(62)-JVS(351)*XX(63)-JVS(361)*XX(64)&
|
|
&-JVS(371)*XX(65)-JVS(381)*XX(66)-JVS(390)*XX(67)-JVS(417)*XX(68)-JVS(466)*XX(69)-JVS(478)*XX(70)-JVS(501)&
|
|
&*XX(71)-JVS(513)*XX(72)-JVS(525)*XX(73)-JVS(542)*XX(74)-JVS(571)*XX(75)-JVS(590)*XX(76)-JVS(603)*XX(77)&
|
|
&-JVS(618)*XX(78)-JVS(629)*XX(79)-JVS(642)*XX(80)-JVS(661)*XX(81)-JVS(708)*XX(82))/(JVS(777))
|
|
XX(84) = (X(84)-JVS(53)*XX(17)-JVS(83)*XX(26)-JVS(86)*XX(27)-JVS(90)*XX(28)-JVS(95)*XX(29)-JVS(99)*XX(30)-JVS(103)&
|
|
&*XX(31)-JVS(107)*XX(32)-JVS(111)*XX(33)-JVS(118)*XX(35)-JVS(121)*XX(36)-JVS(126)*XX(37)-JVS(130)*XX(38)&
|
|
&-JVS(134)*XX(39)-JVS(138)*XX(40)-JVS(142)*XX(41)-JVS(146)*XX(42)-JVS(152)*XX(43)-JVS(162)*XX(45)-JVS(184)&
|
|
&*XX(47)-JVS(202)*XX(49)-JVS(215)*XX(50)-JVS(221)*XX(51)-JVS(231)*XX(52)-JVS(241)*XX(53)-JVS(248)*XX(54)&
|
|
&-JVS(254)*XX(55)-JVS(280)*XX(56)-JVS(287)*XX(57)-JVS(296)*XX(58)-JVS(305)*XX(59)-JVS(312)*XX(60)-JVS(332)&
|
|
&*XX(61)-JVS(343)*XX(62)-JVS(352)*XX(63)-JVS(362)*XX(64)-JVS(372)*XX(65)-JVS(382)*XX(66)-JVS(391)*XX(67)&
|
|
&-JVS(418)*XX(68)-JVS(467)*XX(69)-JVS(479)*XX(70)-JVS(502)*XX(71)-JVS(514)*XX(72)-JVS(526)*XX(73)-JVS(543)&
|
|
&*XX(74)-JVS(572)*XX(75)-JVS(591)*XX(76)-JVS(604)*XX(77)-JVS(619)*XX(78)-JVS(630)*XX(79)-JVS(643)*XX(80)&
|
|
&-JVS(662)*XX(81)-JVS(709)*XX(82)-JVS(778)*XX(83))/(JVS(842))
|
|
XX(85) = (X(85)-JVS(32)*XX(13)-JVS(76)*XX(24)-JVS(185)*XX(47)-JVS(203)*XX(49)-JVS(216)*XX(50)-JVS(222)*XX(51)-JVS(232)&
|
|
&*XX(52)-JVS(242)*XX(53)-JVS(249)*XX(54)-JVS(255)*XX(55)-JVS(281)*XX(56)-JVS(288)*XX(57)-JVS(297)*XX(58)&
|
|
&-JVS(306)*XX(59)-JVS(313)*XX(60)-JVS(333)*XX(61)-JVS(344)*XX(62)-JVS(353)*XX(63)-JVS(363)*XX(64)-JVS(373)&
|
|
&*XX(65)-JVS(383)*XX(66)-JVS(392)*XX(67)-JVS(419)*XX(68)-JVS(468)*XX(69)-JVS(480)*XX(70)-JVS(503)*XX(71)&
|
|
&-JVS(515)*XX(72)-JVS(527)*XX(73)-JVS(544)*XX(74)-JVS(573)*XX(75)-JVS(592)*XX(76)-JVS(605)*XX(77)-JVS(620)&
|
|
&*XX(78)-JVS(631)*XX(79)-JVS(644)*XX(80)-JVS(663)*XX(81)-JVS(710)*XX(82)-JVS(779)*XX(83)-JVS(843)*XX(84))&
|
|
&/(JVS(878))
|
|
XX(86) = (X(86)-JVS(33)*XX(13)-JVS(66)*XX(21)-JVS(87)*XX(27)-JVS(186)*XX(47)-JVS(204)*XX(49)-JVS(217)*XX(50)-JVS(223)&
|
|
&*XX(51)-JVS(233)*XX(52)-JVS(243)*XX(53)-JVS(250)*XX(54)-JVS(256)*XX(55)-JVS(282)*XX(56)-JVS(289)*XX(57)&
|
|
&-JVS(298)*XX(58)-JVS(307)*XX(59)-JVS(314)*XX(60)-JVS(334)*XX(61)-JVS(345)*XX(62)-JVS(354)*XX(63)-JVS(364)&
|
|
&*XX(64)-JVS(374)*XX(65)-JVS(384)*XX(66)-JVS(393)*XX(67)-JVS(420)*XX(68)-JVS(469)*XX(69)-JVS(481)*XX(70)&
|
|
&-JVS(504)*XX(71)-JVS(516)*XX(72)-JVS(528)*XX(73)-JVS(545)*XX(74)-JVS(574)*XX(75)-JVS(593)*XX(76)-JVS(606)&
|
|
&*XX(77)-JVS(621)*XX(78)-JVS(632)*XX(79)-JVS(645)*XX(80)-JVS(664)*XX(81)-JVS(711)*XX(82)-JVS(780)*XX(83)&
|
|
&-JVS(844)*XX(84)-JVS(879)*XX(85))/(JVS(920))
|
|
XX(87) = (X(87)-JVS(63)*XX(20)-JVS(69)*XX(22)-JVS(72)*XX(23)-JVS(79)*XX(25)-JVS(114)*XX(34)-JVS(165)*XX(46)-JVS(187)&
|
|
&*XX(47)-JVS(194)*XX(48)-JVS(257)*XX(55)-JVS(283)*XX(56)-JVS(315)*XX(60)-JVS(355)*XX(63)-JVS(375)*XX(65)&
|
|
&-JVS(394)*XX(67)-JVS(421)*XX(68)-JVS(470)*XX(69)-JVS(482)*XX(70)-JVS(505)*XX(71)-JVS(517)*XX(72)-JVS(529)&
|
|
&*XX(73)-JVS(546)*XX(74)-JVS(575)*XX(75)-JVS(594)*XX(76)-JVS(607)*XX(77)-JVS(622)*XX(78)-JVS(633)*XX(79)&
|
|
&-JVS(646)*XX(80)-JVS(665)*XX(81)-JVS(712)*XX(82)-JVS(781)*XX(83)-JVS(845)*XX(84)-JVS(880)*XX(85)-JVS(921)&
|
|
&*XX(86))/(JVS(959))
|
|
XX(88) = (X(88)-JVS(59)*XX(19)-JVS(122)*XX(36)-JVS(471)*XX(69)-JVS(713)*XX(82)-JVS(782)*XX(83)-JVS(846)*XX(84)&
|
|
&-JVS(881)*XX(85)-JVS(922)*XX(86)-JVS(960)*XX(87))/(JVS(980))
|
|
XX(89) = (X(89)-JVS(12)*XX(6)-JVS(34)*XX(13)-JVS(45)*XX(14)-JVS(157)*XX(44)-JVS(166)*XX(46)-JVS(188)*XX(47)-JVS(195)&
|
|
&*XX(48)-JVS(218)*XX(50)-JVS(299)*XX(58)-JVS(335)*XX(61)-JVS(356)*XX(63)-JVS(385)*XX(66)-JVS(395)*XX(67)&
|
|
&-JVS(422)*XX(68)-JVS(472)*XX(69)-JVS(483)*XX(70)-JVS(506)*XX(71)-JVS(547)*XX(74)-JVS(576)*XX(75)-JVS(595)&
|
|
&*XX(76)-JVS(608)*XX(77)-JVS(623)*XX(78)-JVS(634)*XX(79)-JVS(647)*XX(80)-JVS(666)*XX(81)-JVS(714)*XX(82)&
|
|
&-JVS(783)*XX(83)-JVS(847)*XX(84)-JVS(882)*XX(85)-JVS(923)*XX(86)-JVS(961)*XX(87)-JVS(981)*XX(88))/(JVS(1000))
|
|
XX(90) = (X(90)-JVS(91)*XX(28)-JVS(189)*XX(47)-JVS(205)*XX(49)-JVS(219)*XX(50)-JVS(224)*XX(51)-JVS(234)*XX(52)&
|
|
&-JVS(244)*XX(53)-JVS(251)*XX(54)-JVS(258)*XX(55)-JVS(284)*XX(56)-JVS(290)*XX(57)-JVS(300)*XX(58)-JVS(308)&
|
|
&*XX(59)-JVS(316)*XX(60)-JVS(336)*XX(61)-JVS(346)*XX(62)-JVS(357)*XX(63)-JVS(365)*XX(64)-JVS(376)*XX(65)&
|
|
&-JVS(386)*XX(66)-JVS(396)*XX(67)-JVS(423)*XX(68)-JVS(473)*XX(69)-JVS(484)*XX(70)-JVS(507)*XX(71)-JVS(518)&
|
|
&*XX(72)-JVS(530)*XX(73)-JVS(548)*XX(74)-JVS(577)*XX(75)-JVS(596)*XX(76)-JVS(609)*XX(77)-JVS(624)*XX(78)&
|
|
&-JVS(635)*XX(79)-JVS(648)*XX(80)-JVS(667)*XX(81)-JVS(715)*XX(82)-JVS(784)*XX(83)-JVS(848)*XX(84)-JVS(883)&
|
|
&*XX(85)-JVS(924)*XX(86)-JVS(962)*XX(87)-JVS(982)*XX(88)-JVS(1001)*XX(89))/(JVS(1042))
|
|
XX(90) = XX(90)
|
|
XX(89) = XX(89)-JVS(1041)*XX(90)
|
|
XX(88) = XX(88)-JVS(999)*XX(89)-JVS(1040)*XX(90)
|
|
XX(87) = XX(87)-JVS(979)*XX(88)-JVS(998)*XX(89)-JVS(1039)*XX(90)
|
|
XX(86) = XX(86)-JVS(958)*XX(87)-JVS(978)*XX(88)-JVS(997)*XX(89)-JVS(1038)*XX(90)
|
|
XX(85) = XX(85)-JVS(919)*XX(86)-JVS(957)*XX(87)-JVS(977)*XX(88)-JVS(996)*XX(89)-JVS(1037)*XX(90)
|
|
XX(84) = XX(84)-JVS(877)*XX(85)-JVS(918)*XX(86)-JVS(956)*XX(87)-JVS(976)*XX(88)-JVS(995)*XX(89)-JVS(1036)*XX(90)
|
|
XX(83) = XX(83)-JVS(841)*XX(84)-JVS(876)*XX(85)-JVS(917)*XX(86)-JVS(955)*XX(87)-JVS(975)*XX(88)-JVS(994)*XX(89)&
|
|
&-JVS(1035)*XX(90)
|
|
XX(82) = XX(82)-JVS(776)*XX(83)-JVS(840)*XX(84)-JVS(875)*XX(85)-JVS(916)*XX(86)-JVS(954)*XX(87)-JVS(974)*XX(88)&
|
|
&-JVS(993)*XX(89)-JVS(1034)*XX(90)
|
|
XX(81) = XX(81)-JVS(706)*XX(82)-JVS(775)*XX(83)-JVS(839)*XX(84)-JVS(874)*XX(85)-JVS(915)*XX(86)-JVS(953)*XX(87)&
|
|
&-JVS(992)*XX(89)-JVS(1033)*XX(90)
|
|
XX(80) = XX(80)-JVS(658)*XX(81)-JVS(705)*XX(82)-JVS(774)*XX(83)-JVS(838)*XX(84)-JVS(873)*XX(85)-JVS(914)*XX(86)&
|
|
&-JVS(952)*XX(87)-JVS(973)*XX(88)-JVS(991)*XX(89)-JVS(1032)*XX(90)
|
|
XX(79) = XX(79)-JVS(640)*XX(80)-JVS(657)*XX(81)-JVS(704)*XX(82)-JVS(773)*XX(83)-JVS(837)*XX(84)-JVS(872)*XX(85)&
|
|
&-JVS(913)*XX(86)-JVS(951)*XX(87)-JVS(972)*XX(88)-JVS(990)*XX(89)-JVS(1031)*XX(90)
|
|
XX(78) = XX(78)-JVS(703)*XX(82)-JVS(772)*XX(83)-JVS(836)*XX(84)-JVS(871)*XX(85)-JVS(912)*XX(86)-JVS(950)*XX(87)&
|
|
&-JVS(989)*XX(89)-JVS(1030)*XX(90)
|
|
XX(77) = XX(77)-JVS(613)*XX(78)-JVS(656)*XX(81)-JVS(702)*XX(82)-JVS(771)*XX(83)-JVS(835)*XX(84)-JVS(870)*XX(85)&
|
|
&-JVS(911)*XX(86)-JVS(949)*XX(87)-JVS(971)*XX(88)-JVS(988)*XX(89)-JVS(1029)*XX(90)
|
|
XX(76) = XX(76)-JVS(655)*XX(81)-JVS(701)*XX(82)-JVS(770)*XX(83)-JVS(834)*XX(84)-JVS(869)*XX(85)-JVS(910)*XX(86)&
|
|
&-JVS(948)*XX(87)-JVS(1028)*XX(90)
|
|
XX(75) = XX(75)-JVS(700)*XX(82)-JVS(769)*XX(83)-JVS(833)*XX(84)-JVS(909)*XX(86)-JVS(947)*XX(87)-JVS(1027)*XX(90)
|
|
XX(74) = XX(74)-JVS(699)*XX(82)-JVS(768)*XX(83)-JVS(832)*XX(84)-JVS(868)*XX(85)-JVS(908)*XX(86)-JVS(946)*XX(87)&
|
|
&-JVS(1026)*XX(90)
|
|
XX(73) = XX(73)-JVS(536)*XX(74)-JVS(563)*XX(75)-JVS(585)*XX(76)-JVS(654)*XX(81)-JVS(698)*XX(82)-JVS(767)*XX(83)&
|
|
&-JVS(831)*XX(84)-JVS(867)*XX(85)-JVS(907)*XX(86)-JVS(945)*XX(87)-JVS(1025)*XX(90)
|
|
XX(72) = XX(72)-JVS(522)*XX(73)-JVS(535)*XX(74)-JVS(562)*XX(75)-JVS(584)*XX(76)-JVS(653)*XX(81)-JVS(697)*XX(82)&
|
|
&-JVS(766)*XX(83)-JVS(830)*XX(84)-JVS(866)*XX(85)-JVS(906)*XX(86)-JVS(944)*XX(87)-JVS(1024)*XX(90)
|
|
XX(71) = XX(71)-JVS(696)*XX(82)-JVS(765)*XX(83)-JVS(829)*XX(84)-JVS(905)*XX(86)-JVS(943)*XX(87)-JVS(1023)*XX(90)
|
|
XX(70) = XX(70)-JVS(492)*XX(71)-JVS(534)*XX(74)-JVS(561)*XX(75)-JVS(652)*XX(81)-JVS(695)*XX(82)-JVS(764)*XX(83)&
|
|
&-JVS(828)*XX(84)-JVS(865)*XX(85)-JVS(904)*XX(86)-JVS(942)*XX(87)-JVS(1022)*XX(90)
|
|
XX(69) = XX(69)-JVS(694)*XX(82)-JVS(763)*XX(83)-JVS(827)*XX(84)-JVS(941)*XX(87)
|
|
XX(68) = XX(68)-JVS(533)*XX(74)-JVS(651)*XX(81)-JVS(693)*XX(82)-JVS(762)*XX(83)-JVS(826)*XX(84)-JVS(940)*XX(87)
|
|
XX(67) = XX(67)-JVS(410)*XX(68)-JVS(452)*XX(69)-JVS(600)*XX(77)-JVS(639)*XX(80)-JVS(692)*XX(82)-JVS(761)*XX(83)&
|
|
&-JVS(825)*XX(84)-JVS(864)*XX(85)-JVS(903)*XX(86)-JVS(939)*XX(87)-JVS(1021)*XX(90)
|
|
XX(66) = XX(66)-JVS(451)*XX(69)-JVS(560)*XX(75)-JVS(583)*XX(76)-JVS(691)*XX(82)-JVS(760)*XX(83)-JVS(824)*XX(84)&
|
|
&-JVS(863)*XX(85)-JVS(902)*XX(86)-JVS(938)*XX(87)-JVS(970)*XX(88)-JVS(1020)*XX(90)
|
|
XX(65) = XX(65)-JVS(409)*XX(68)-JVS(450)*XX(69)-JVS(559)*XX(75)-JVS(690)*XX(82)-JVS(759)*XX(83)-JVS(823)*XX(84)&
|
|
&-JVS(862)*XX(85)-JVS(901)*XX(86)-JVS(937)*XX(87)-JVS(969)*XX(88)-JVS(1019)*XX(90)
|
|
XX(64) = XX(64)-JVS(408)*XX(68)-JVS(449)*XX(69)-JVS(558)*XX(75)-JVS(582)*XX(76)-JVS(689)*XX(82)-JVS(758)*XX(83)&
|
|
&-JVS(822)*XX(84)-JVS(861)*XX(85)-JVS(900)*XX(86)-JVS(936)*XX(87)-JVS(968)*XX(88)-JVS(1018)*XX(90)
|
|
XX(63) = XX(63)-JVS(407)*XX(68)-JVS(448)*XX(69)-JVS(491)*XX(71)-JVS(688)*XX(82)-JVS(757)*XX(83)-JVS(821)*XX(84)&
|
|
&-JVS(860)*XX(85)-JVS(899)*XX(86)-JVS(1017)*XX(90)
|
|
XX(62) = XX(62)-JVS(447)*XX(69)-JVS(521)*XX(73)-JVS(557)*XX(75)-JVS(581)*XX(76)-JVS(687)*XX(82)-JVS(756)*XX(83)&
|
|
&-JVS(820)*XX(84)-JVS(859)*XX(85)-JVS(898)*XX(86)-JVS(1016)*XX(90)
|
|
XX(61) = XX(61)-JVS(446)*XX(69)-JVS(556)*XX(75)-JVS(755)*XX(83)-JVS(819)*XX(84)-JVS(897)*XX(86)
|
|
XX(60) = XX(60)-JVS(445)*XX(69)-JVS(490)*XX(71)-JVS(520)*XX(73)-JVS(555)*XX(75)-JVS(580)*XX(76)-JVS(686)*XX(82)&
|
|
&-JVS(818)*XX(84)-JVS(858)*XX(85)-JVS(896)*XX(86)-JVS(1015)*XX(90)
|
|
XX(59) = XX(59)-JVS(321)*XX(61)-JVS(444)*XX(69)-JVS(554)*XX(75)-JVS(579)*XX(76)-JVS(685)*XX(82)-JVS(754)*XX(83)&
|
|
&-JVS(817)*XX(84)-JVS(857)*XX(85)-JVS(895)*XX(86)-JVS(935)*XX(87)-JVS(1014)*XX(90)
|
|
XX(58) = XX(58)-JVS(320)*XX(61)-JVS(378)*XX(66)-JVS(443)*XX(69)-JVS(553)*XX(75)-JVS(612)*XX(78)-JVS(753)*XX(83)&
|
|
&-JVS(816)*XX(84)-JVS(967)*XX(88)-JVS(987)*XX(89)
|
|
XX(57) = XX(57)-JVS(292)*XX(58)-JVS(442)*XX(69)-JVS(578)*XX(76)-JVS(599)*XX(77)-JVS(627)*XX(79)-JVS(638)*XX(80)&
|
|
&-JVS(684)*XX(82)-JVS(752)*XX(83)-JVS(815)*XX(84)-JVS(856)*XX(85)-JVS(894)*XX(86)-JVS(934)*XX(87)-JVS(1013)&
|
|
&*XX(90)
|
|
XX(56) = XX(56)-JVS(683)*XX(82)-JVS(751)*XX(83)-JVS(933)*XX(87)
|
|
XX(55) = XX(55)-JVS(368)*XX(65)-JVS(406)*XX(68)-JVS(441)*XX(69)-JVS(552)*XX(75)-JVS(682)*XX(82)-JVS(750)*XX(83)&
|
|
&-JVS(814)*XX(84)-JVS(855)*XX(85)-JVS(893)*XX(86)-JVS(1012)*XX(90)
|
|
XX(54) = XX(54)-JVS(405)*XX(68)-JVS(440)*XX(69)-JVS(489)*XX(71)-JVS(510)*XX(72)-JVS(519)*XX(73)-JVS(681)*XX(82)&
|
|
&-JVS(813)*XX(84)-JVS(854)*XX(85)-JVS(892)*XX(86)-JVS(1011)*XX(90)
|
|
XX(53) = XX(53)-JVS(338)*XX(62)-JVS(439)*XX(69)-JVS(680)*XX(82)-JVS(749)*XX(83)-JVS(812)*XX(84)-JVS(853)*XX(85)&
|
|
&-JVS(891)*XX(86)-JVS(1010)*XX(90)
|
|
XX(52) = XX(52)-JVS(404)*XX(68)-JVS(438)*XX(69)-JVS(679)*XX(82)-JVS(748)*XX(83)-JVS(811)*XX(84)-JVS(852)*XX(85)&
|
|
&-JVS(890)*XX(86)-JVS(1009)*XX(90)
|
|
XX(51) = XX(51)-JVS(263)*XX(56)-JVS(367)*XX(65)-JVS(403)*XX(68)-JVS(437)*XX(69)-JVS(551)*XX(75)-JVS(678)*XX(82)&
|
|
&-JVS(747)*XX(83)-JVS(810)*XX(84)-JVS(851)*XX(85)-JVS(889)*XX(86)-JVS(966)*XX(88)-JVS(1008)*XX(90)
|
|
XX(50) = XX(50)-JVS(436)*XX(69)-JVS(746)*XX(83)-JVS(809)*XX(84)-JVS(965)*XX(88)
|
|
XX(49) = XX(49)-JVS(337)*XX(62)-JVS(745)*XX(83)-JVS(808)*XX(84)-JVS(888)*XX(86)-JVS(1007)*XX(90)
|
|
XX(48) = XX(48)-JVS(348)*XX(63)-JVS(435)*XX(69)-JVS(475)*XX(70)-JVS(488)*XX(71)-JVS(744)*XX(83)-JVS(807)*XX(84)&
|
|
&-JVS(932)*XX(87)-JVS(986)*XX(89)-JVS(1006)*XX(90)
|
|
XX(47) = XX(47)-JVS(743)*XX(83)-JVS(806)*XX(84)
|
|
XX(46) = XX(46)-JVS(169)*XX(47)-JVS(190)*XX(48)-JVS(388)*XX(67)-JVS(434)*XX(69)-JVS(598)*XX(77)-JVS(626)*XX(79)&
|
|
&-JVS(637)*XX(80)-JVS(742)*XX(83)-JVS(805)*XX(84)-JVS(931)*XX(87)-JVS(985)*XX(89)
|
|
XX(45) = XX(45)-JVS(291)*XX(58)-JVS(433)*XX(69)-JVS(597)*XX(77)-JVS(625)*XX(79)-JVS(636)*XX(80)-JVS(741)*XX(83)&
|
|
&-JVS(804)*XX(84)
|
|
XX(44) = XX(44)-JVS(319)*XX(61)-JVS(432)*XX(69)-JVS(611)*XX(78)-JVS(677)*XX(82)-JVS(740)*XX(83)-JVS(803)*XX(84)&
|
|
&-JVS(984)*XX(89)
|
|
XX(43) = XX(43)-JVS(366)*XX(65)-JVS(402)*XX(68)-JVS(676)*XX(82)-JVS(739)*XX(83)-JVS(802)*XX(84)
|
|
XX(42) = XX(42)-JVS(207)*XX(50)-JVS(358)*XX(64)-JVS(401)*XX(68)-JVS(431)*XX(69)-JVS(550)*XX(75)-JVS(738)*XX(83)&
|
|
&-JVS(801)*XX(84)-JVS(887)*XX(86)
|
|
XX(41) = XX(41)-JVS(168)*XX(47)-JVS(206)*XX(50)-JVS(318)*XX(61)-JVS(377)*XX(66)-JVS(549)*XX(75)-JVS(737)*XX(83)&
|
|
&-JVS(800)*XX(84)
|
|
XX(40) = XX(40)-JVS(167)*XX(47)-JVS(301)*XX(59)-JVS(317)*XX(61)-JVS(430)*XX(69)-JVS(736)*XX(83)-JVS(799)*XX(84)
|
|
XX(39) = XX(39)-JVS(387)*XX(67)-JVS(400)*XX(68)-JVS(675)*XX(82)-JVS(735)*XX(83)-JVS(798)*XX(84)
|
|
XX(38) = XX(38)-JVS(347)*XX(63)-JVS(399)*XX(68)-JVS(674)*XX(82)-JVS(734)*XX(83)-JVS(797)*XX(84)
|
|
XX(37) = XX(37)-JVS(429)*XX(69)-JVS(474)*XX(70)-JVS(487)*XX(71)-JVS(733)*XX(83)-JVS(796)*XX(84)
|
|
XX(36) = XX(36)-JVS(428)*XX(69)-JVS(732)*XX(83)-JVS(795)*XX(84)-JVS(964)*XX(88)
|
|
XX(35) = XX(35)-JVS(486)*XX(71)-JVS(532)*XX(74)-JVS(731)*XX(83)-JVS(794)*XX(84)
|
|
XX(34) = XX(34)-JVS(262)*XX(56)-JVS(427)*XX(69)-JVS(730)*XX(83)-JVS(793)*XX(84)-JVS(930)*XX(87)-JVS(1005)*XX(90)
|
|
XX(33) = XX(33)-JVS(485)*XX(71)-JVS(650)*XX(81)-JVS(729)*XX(83)-JVS(792)*XX(84)
|
|
XX(32) = XX(32)-JVS(196)*XX(49)-JVS(236)*XX(53)-JVS(728)*XX(83)-JVS(791)*XX(84)
|
|
XX(31) = XX(31)-JVS(226)*XX(52)-JVS(398)*XX(68)-JVS(727)*XX(83)-JVS(790)*XX(84)
|
|
XX(30) = XX(30)-JVS(397)*XX(68)-JVS(509)*XX(72)-JVS(726)*XX(83)-JVS(789)*XX(84)
|
|
XX(29) = XX(29)-JVS(673)*XX(82)-JVS(725)*XX(83)-JVS(788)*XX(84)-JVS(929)*XX(87)
|
|
XX(28) = XX(28)-JVS(426)*XX(69)-JVS(724)*XX(83)-JVS(787)*XX(84)-JVS(1004)*XX(90)
|
|
XX(27) = XX(27)-JVS(425)*XX(69)-JVS(723)*XX(83)-JVS(886)*XX(86)-JVS(1003)*XX(90)
|
|
XX(26) = XX(26)-JVS(424)*XX(69)-JVS(610)*XX(78)-JVS(722)*XX(83)-JVS(885)*XX(86)
|
|
XX(25) = XX(25)-JVS(261)*XX(56)-JVS(672)*XX(82)-JVS(850)*XX(85)-JVS(928)*XX(87)-JVS(983)*XX(89)
|
|
XX(24) = XX(24)-JVS(671)*XX(82)-JVS(721)*XX(83)-JVS(849)*XX(85)
|
|
XX(23) = XX(23)-JVS(260)*XX(56)-JVS(531)*XX(74)-JVS(720)*XX(83)-JVS(927)*XX(87)
|
|
XX(22) = XX(22)-JVS(259)*XX(56)-JVS(508)*XX(72)-JVS(719)*XX(83)-JVS(926)*XX(87)
|
|
XX(21) = XX(21)-JVS(670)*XX(82)-JVS(884)*XX(86)-JVS(925)*XX(87)-JVS(1002)*XX(90)
|
|
XX(20) = XX(20)-JVS(718)*XX(83)-JVS(786)*XX(84)
|
|
XX(19) = XX(19)-JVS(669)*XX(82)-JVS(963)*XX(88)
|
|
XX(18) = XX(18)-JVS(649)*XX(81)-JVS(668)*XX(82)
|
|
XX(17) = XX(17)-JVS(717)*XX(83)-JVS(785)*XX(84)
|
|
XX(16) = XX(16)-JVS(225)*XX(52)-JVS(235)*XX(53)-JVS(716)*XX(83)
|
|
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_adj_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_adj_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.
|
|
|
|
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_adj_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_adj_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_adj_LinearAlgebra
|
|
|