MODULE cbm_LinearAlgebra USE cbm_Parameters USE cbm_JacobianSP IMPLICIT NONE CONTAINS ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE KppDecomp( JVS, IER ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Sparse LU factorization ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters USE cbm_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 cbm_Parameters ! USE cbm_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 cbm_Parameters ! USE cbm_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 cbm_Parameters ! USE cbm_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(11) = X(11)-JVS(38)*X(5)-JVS(39)*X(7) X(12) = X(12)-JVS(43)*X(6) X(14) = X(14)-JVS(55)*X(5)-JVS(56)*X(7)-JVS(57)*X(11) X(15) = X(15)-JVS(62)*X(7) X(16) = X(16)-JVS(68)*X(15) X(18) = X(18)-JVS(85)*X(5)-JVS(86)*X(7)-JVS(87)*X(13)-JVS(88)*X(14)-JVS(89)*X(15)-JVS(90)*X(17) X(19) = X(19)-JVS(105)*X(11)-JVS(106)*X(14) X(20) = X(20)-JVS(112)*X(7)-JVS(113)*X(13) X(21) = X(21)-JVS(122)*X(17)-JVS(123)*X(19) X(23) = X(23)-JVS(140)*X(22) X(24) = X(24)-JVS(146)*X(13)-JVS(147)*X(17)-JVS(148)*X(19)-JVS(149)*X(20)-JVS(150)*X(22)-JVS(151)*X(23) X(25) = X(25)-JVS(159)*X(17)-JVS(160)*X(19)-JVS(161)*X(22)-JVS(162)*X(23) X(26) = X(26)-JVS(170)*X(3)-JVS(171)*X(4)-JVS(172)*X(6)-JVS(173)*X(9)-JVS(174)*X(10)-JVS(175)*X(11)-JVS(176)*X(13)& &-JVS(177)*X(14)-JVS(178)*X(18)-JVS(179)*X(19)-JVS(180)*X(20)-JVS(181)*X(22)-JVS(182)*X(23)-JVS(183)*X(24)& &-JVS(184)*X(25) X(27) = X(27)-JVS(192)*X(1)-JVS(193)*X(2)-JVS(194)*X(5)-JVS(195)*X(7)-JVS(196)*X(9)-JVS(197)*X(10)-JVS(198)*X(12)& &-JVS(199)*X(14)-JVS(200)*X(15)-JVS(201)*X(16)-JVS(202)*X(17)-JVS(203)*X(19)-JVS(204)*X(20)-JVS(205)*X(21)& &-JVS(206)*X(22)-JVS(207)*X(23)-JVS(208)*X(24)-JVS(209)*X(25)-JVS(210)*X(26) X(28) = X(28)-JVS(217)*X(2)-JVS(218)*X(5)-JVS(219)*X(7)-JVS(220)*X(10)-JVS(221)*X(11)-JVS(222)*X(13)-JVS(223)*X(14)& &-JVS(224)*X(15)-JVS(225)*X(16)-JVS(226)*X(17)-JVS(227)*X(19)-JVS(228)*X(20)-JVS(229)*X(21)-JVS(230)*X(22)& &-JVS(231)*X(23)-JVS(232)*X(24)-JVS(233)*X(25)-JVS(234)*X(26)-JVS(235)*X(27) X(29) = X(29)-JVS(241)*X(1)-JVS(242)*X(17)-JVS(243)*X(21)-JVS(244)*X(22)-JVS(245)*X(23)-JVS(246)*X(24)-JVS(247)*X(25)& &-JVS(248)*X(26)-JVS(249)*X(27)-JVS(250)*X(28) X(30) = X(30)-JVS(255)*X(6)-JVS(256)*X(12)-JVS(257)*X(14)-JVS(258)*X(21)-JVS(259)*X(22)-JVS(260)*X(23)-JVS(261)*X(24)& &-JVS(262)*X(25)-JVS(263)*X(26)-JVS(264)*X(27)-JVS(265)*X(28)-JVS(266)*X(29) X(31) = X(31)-JVS(270)*X(8)-JVS(271)*X(9)-JVS(272)*X(11)-JVS(273)*X(13)-JVS(274)*X(18)-JVS(275)*X(19)-JVS(276)*X(20)& &-JVS(277)*X(22)-JVS(278)*X(23)-JVS(279)*X(24)-JVS(280)*X(25)-JVS(281)*X(26)-JVS(282)*X(27)-JVS(283)*X(28)& &-JVS(284)*X(29)-JVS(285)*X(30) X(32) = X(32)-JVS(288)*X(3)-JVS(289)*X(15)-JVS(290)*X(19)-JVS(291)*X(22)-JVS(292)*X(24)-JVS(293)*X(25)-JVS(294)*X(26)& &-JVS(295)*X(27)-JVS(296)*X(28)-JVS(297)*X(29)-JVS(298)*X(30)-JVS(299)*X(31) X(32) = X(32)/JVS(300) X(31) = (X(31)-JVS(287)*X(32))/(JVS(286)) X(30) = (X(30)-JVS(268)*X(31)-JVS(269)*X(32))/(JVS(267)) X(29) = (X(29)-JVS(252)*X(30)-JVS(253)*X(31)-JVS(254)*X(32))/(JVS(251)) X(28) = (X(28)-JVS(237)*X(29)-JVS(238)*X(30)-JVS(239)*X(31)-JVS(240)*X(32))/(JVS(236)) X(27) = (X(27)-JVS(212)*X(28)-JVS(213)*X(29)-JVS(214)*X(30)-JVS(215)*X(31)-JVS(216)*X(32))/(JVS(211)) X(26) = (X(26)-JVS(186)*X(27)-JVS(187)*X(28)-JVS(188)*X(29)-JVS(189)*X(30)-JVS(190)*X(31)-JVS(191)*X(32))/(JVS(185)) X(25) = (X(25)-JVS(164)*X(26)-JVS(165)*X(27)-JVS(166)*X(28)-JVS(167)*X(29)-JVS(168)*X(30)-JVS(169)*X(31))/(JVS(163)) X(24) = (X(24)-JVS(153)*X(25)-JVS(154)*X(26)-JVS(155)*X(27)-JVS(156)*X(29)-JVS(157)*X(30)-JVS(158)*X(31))/(JVS(152)) X(23) = (X(23)-JVS(142)*X(25)-JVS(143)*X(27)-JVS(144)*X(29)-JVS(145)*X(30))/(JVS(141)) X(22) = (X(22)-JVS(136)*X(25)-JVS(137)*X(27)-JVS(138)*X(29)-JVS(139)*X(30))/(JVS(135)) X(21) = (X(21)-JVS(125)*X(22)-JVS(126)*X(23)-JVS(127)*X(24)-JVS(128)*X(25)-JVS(129)*X(27)-JVS(130)*X(28)-JVS(131)& &*X(29)-JVS(132)*X(30)-JVS(133)*X(31)-JVS(134)*X(32))/(JVS(124)) X(20) = (X(20)-JVS(115)*X(22)-JVS(116)*X(23)-JVS(117)*X(25)-JVS(118)*X(26)-JVS(119)*X(27)-JVS(120)*X(29)-JVS(121)& &*X(30))/(JVS(114)) X(19) = (X(19)-JVS(108)*X(25)-JVS(109)*X(27)-JVS(110)*X(30)-JVS(111)*X(31))/(JVS(107)) X(18) = (X(18)-JVS(92)*X(19)-JVS(93)*X(20)-JVS(94)*X(22)-JVS(95)*X(23)-JVS(96)*X(24)-JVS(97)*X(25)-JVS(98)*X(26)& &-JVS(99)*X(27)-JVS(100)*X(28)-JVS(101)*X(29)-JVS(102)*X(30)-JVS(103)*X(31)-JVS(104)*X(32))/(JVS(91)) X(17) = (X(17)-JVS(81)*X(22)-JVS(82)*X(25)-JVS(83)*X(27)-JVS(84)*X(29))/(JVS(80)) X(16) = (X(16)-JVS(70)*X(17)-JVS(71)*X(19)-JVS(72)*X(21)-JVS(73)*X(22)-JVS(74)*X(23)-JVS(75)*X(24)-JVS(76)*X(25)& &-JVS(77)*X(27)-JVS(78)*X(29)-JVS(79)*X(30))/(JVS(69)) X(15) = (X(15)-JVS(64)*X(19)-JVS(65)*X(22)-JVS(66)*X(25)-JVS(67)*X(27))/(JVS(63)) X(14) = (X(14)-JVS(59)*X(27)-JVS(60)*X(30)-JVS(61)*X(31))/(JVS(58)) X(13) = (X(13)-JVS(52)*X(20)-JVS(53)*X(26)-JVS(54)*X(27))/(JVS(51)) X(12) = (X(12)-JVS(45)*X(14)-JVS(46)*X(21)-JVS(47)*X(24)-JVS(48)*X(26)-JVS(49)*X(27)-JVS(50)*X(30))/(JVS(44)) X(11) = (X(11)-JVS(41)*X(27)-JVS(42)*X(31))/(JVS(40)) X(10) = (X(10)-JVS(35)*X(26)-JVS(36)*X(27)-JVS(37)*X(28))/(JVS(34)) X(9) = (X(9)-JVS(31)*X(26)-JVS(32)*X(27)-JVS(33)*X(31))/(JVS(30)) X(8) = (X(8)-JVS(22)*X(13)-JVS(23)*X(20)-JVS(24)*X(22)-JVS(25)*X(23)-JVS(26)*X(27)-JVS(27)*X(29)-JVS(28)*X(30)-JVS(29)& &*X(31))/(JVS(21)) X(7) = (X(7)-JVS(20)*X(27))/(JVS(19)) X(6) = (X(6)-JVS(17)*X(26)-JVS(18)*X(30))/(JVS(16)) X(5) = (X(5)-JVS(15)*X(27))/(JVS(14)) X(4) = (X(4)-JVS(10)*X(14)-JVS(11)*X(26)-JVS(12)*X(27)-JVS(13)*X(30))/(JVS(9)) X(3) = (X(3)-JVS(7)*X(26)-JVS(8)*X(32))/(JVS(6)) X(2) = (X(2)-JVS(4)*X(27)-JVS(5)*X(28))/(JVS(3)) X(1) = (X(1)-JVS(2)*X(25))/(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(6) XX(4) = X(4)/JVS(9) XX(5) = X(5)/JVS(14) XX(6) = X(6)/JVS(16) XX(7) = X(7)/JVS(19) XX(8) = X(8)/JVS(21) XX(9) = X(9)/JVS(30) XX(10) = X(10)/JVS(34) XX(11) = X(11)/JVS(40) XX(12) = X(12)/JVS(44) XX(13) = (X(13)-JVS(22)*XX(8))/(JVS(51)) XX(14) = (X(14)-JVS(10)*XX(4)-JVS(45)*XX(12))/(JVS(58)) XX(15) = X(15)/JVS(63) XX(16) = X(16)/JVS(69) XX(17) = (X(17)-JVS(70)*XX(16))/(JVS(80)) XX(18) = X(18)/JVS(91) XX(19) = (X(19)-JVS(64)*XX(15)-JVS(71)*XX(16)-JVS(92)*XX(18))/(JVS(107)) XX(20) = (X(20)-JVS(23)*XX(8)-JVS(52)*XX(13)-JVS(93)*XX(18))/(JVS(114)) XX(21) = (X(21)-JVS(46)*XX(12)-JVS(72)*XX(16))/(JVS(124)) XX(22) = (X(22)-JVS(24)*XX(8)-JVS(65)*XX(15)-JVS(73)*XX(16)-JVS(81)*XX(17)-JVS(94)*XX(18)-JVS(115)*XX(20)-JVS(125)& &*XX(21))/(JVS(135)) XX(23) = (X(23)-JVS(25)*XX(8)-JVS(74)*XX(16)-JVS(95)*XX(18)-JVS(116)*XX(20)-JVS(126)*XX(21))/(JVS(141)) XX(24) = (X(24)-JVS(47)*XX(12)-JVS(75)*XX(16)-JVS(96)*XX(18)-JVS(127)*XX(21))/(JVS(152)) XX(25) = (X(25)-JVS(2)*XX(1)-JVS(66)*XX(15)-JVS(76)*XX(16)-JVS(82)*XX(17)-JVS(97)*XX(18)-JVS(108)*XX(19)-JVS(117)& &*XX(20)-JVS(128)*XX(21)-JVS(136)*XX(22)-JVS(142)*XX(23)-JVS(153)*XX(24))/(JVS(163)) XX(26) = (X(26)-JVS(7)*XX(3)-JVS(11)*XX(4)-JVS(17)*XX(6)-JVS(31)*XX(9)-JVS(35)*XX(10)-JVS(48)*XX(12)-JVS(53)*XX(13)& &-JVS(98)*XX(18)-JVS(118)*XX(20)-JVS(154)*XX(24)-JVS(164)*XX(25))/(JVS(185)) XX(27) = (X(27)-JVS(4)*XX(2)-JVS(12)*XX(4)-JVS(15)*XX(5)-JVS(20)*XX(7)-JVS(26)*XX(8)-JVS(32)*XX(9)-JVS(36)*XX(10)& &-JVS(41)*XX(11)-JVS(49)*XX(12)-JVS(54)*XX(13)-JVS(59)*XX(14)-JVS(67)*XX(15)-JVS(77)*XX(16)-JVS(83)*XX(17)& &-JVS(99)*XX(18)-JVS(109)*XX(19)-JVS(119)*XX(20)-JVS(129)*XX(21)-JVS(137)*XX(22)-JVS(143)*XX(23)-JVS(155)& &*XX(24)-JVS(165)*XX(25)-JVS(186)*XX(26))/(JVS(211)) XX(28) = (X(28)-JVS(5)*XX(2)-JVS(37)*XX(10)-JVS(100)*XX(18)-JVS(130)*XX(21)-JVS(166)*XX(25)-JVS(187)*XX(26)-JVS(212)& &*XX(27))/(JVS(236)) XX(29) = (X(29)-JVS(27)*XX(8)-JVS(78)*XX(16)-JVS(84)*XX(17)-JVS(101)*XX(18)-JVS(120)*XX(20)-JVS(131)*XX(21)-JVS(138)& &*XX(22)-JVS(144)*XX(23)-JVS(156)*XX(24)-JVS(167)*XX(25)-JVS(188)*XX(26)-JVS(213)*XX(27)-JVS(237)*XX(28))& &/(JVS(251)) XX(30) = (X(30)-JVS(13)*XX(4)-JVS(18)*XX(6)-JVS(28)*XX(8)-JVS(50)*XX(12)-JVS(60)*XX(14)-JVS(79)*XX(16)-JVS(102)*XX(18)& &-JVS(110)*XX(19)-JVS(121)*XX(20)-JVS(132)*XX(21)-JVS(139)*XX(22)-JVS(145)*XX(23)-JVS(157)*XX(24)-JVS(168)& &*XX(25)-JVS(189)*XX(26)-JVS(214)*XX(27)-JVS(238)*XX(28)-JVS(252)*XX(29))/(JVS(267)) XX(31) = (X(31)-JVS(29)*XX(8)-JVS(33)*XX(9)-JVS(42)*XX(11)-JVS(61)*XX(14)-JVS(103)*XX(18)-JVS(111)*XX(19)-JVS(133)& &*XX(21)-JVS(158)*XX(24)-JVS(169)*XX(25)-JVS(190)*XX(26)-JVS(215)*XX(27)-JVS(239)*XX(28)-JVS(253)*XX(29)& &-JVS(268)*XX(30))/(JVS(286)) XX(32) = (X(32)-JVS(8)*XX(3)-JVS(104)*XX(18)-JVS(134)*XX(21)-JVS(191)*XX(26)-JVS(216)*XX(27)-JVS(240)*XX(28)-JVS(254)& &*XX(29)-JVS(269)*XX(30)-JVS(287)*XX(31))/(JVS(300)) XX(32) = XX(32) XX(31) = XX(31)-JVS(299)*XX(32) XX(30) = XX(30)-JVS(285)*XX(31)-JVS(298)*XX(32) XX(29) = XX(29)-JVS(266)*XX(30)-JVS(284)*XX(31)-JVS(297)*XX(32) XX(28) = XX(28)-JVS(250)*XX(29)-JVS(265)*XX(30)-JVS(283)*XX(31)-JVS(296)*XX(32) XX(27) = XX(27)-JVS(235)*XX(28)-JVS(249)*XX(29)-JVS(264)*XX(30)-JVS(282)*XX(31)-JVS(295)*XX(32) XX(26) = XX(26)-JVS(210)*XX(27)-JVS(234)*XX(28)-JVS(248)*XX(29)-JVS(263)*XX(30)-JVS(281)*XX(31)-JVS(294)*XX(32) XX(25) = XX(25)-JVS(184)*XX(26)-JVS(209)*XX(27)-JVS(233)*XX(28)-JVS(247)*XX(29)-JVS(262)*XX(30)-JVS(280)*XX(31)& &-JVS(293)*XX(32) XX(24) = XX(24)-JVS(183)*XX(26)-JVS(208)*XX(27)-JVS(232)*XX(28)-JVS(246)*XX(29)-JVS(261)*XX(30)-JVS(279)*XX(31)& &-JVS(292)*XX(32) XX(23) = XX(23)-JVS(151)*XX(24)-JVS(162)*XX(25)-JVS(182)*XX(26)-JVS(207)*XX(27)-JVS(231)*XX(28)-JVS(245)*XX(29)& &-JVS(260)*XX(30)-JVS(278)*XX(31) XX(22) = XX(22)-JVS(140)*XX(23)-JVS(150)*XX(24)-JVS(161)*XX(25)-JVS(181)*XX(26)-JVS(206)*XX(27)-JVS(230)*XX(28)& &-JVS(244)*XX(29)-JVS(259)*XX(30)-JVS(277)*XX(31)-JVS(291)*XX(32) XX(21) = XX(21)-JVS(205)*XX(27)-JVS(229)*XX(28)-JVS(243)*XX(29)-JVS(258)*XX(30) XX(20) = XX(20)-JVS(149)*XX(24)-JVS(180)*XX(26)-JVS(204)*XX(27)-JVS(228)*XX(28)-JVS(276)*XX(31) XX(19) = XX(19)-JVS(123)*XX(21)-JVS(148)*XX(24)-JVS(160)*XX(25)-JVS(179)*XX(26)-JVS(203)*XX(27)-JVS(227)*XX(28)& &-JVS(275)*XX(31)-JVS(290)*XX(32) XX(18) = XX(18)-JVS(178)*XX(26)-JVS(274)*XX(31) XX(17) = XX(17)-JVS(90)*XX(18)-JVS(122)*XX(21)-JVS(147)*XX(24)-JVS(159)*XX(25)-JVS(202)*XX(27)-JVS(226)*XX(28)& &-JVS(242)*XX(29) XX(16) = XX(16)-JVS(201)*XX(27)-JVS(225)*XX(28) XX(15) = XX(15)-JVS(68)*XX(16)-JVS(89)*XX(18)-JVS(200)*XX(27)-JVS(224)*XX(28)-JVS(289)*XX(32) XX(14) = XX(14)-JVS(88)*XX(18)-JVS(106)*XX(19)-JVS(177)*XX(26)-JVS(199)*XX(27)-JVS(223)*XX(28)-JVS(257)*XX(30) XX(13) = XX(13)-JVS(87)*XX(18)-JVS(113)*XX(20)-JVS(146)*XX(24)-JVS(176)*XX(26)-JVS(222)*XX(28)-JVS(273)*XX(31) XX(12) = XX(12)-JVS(198)*XX(27)-JVS(256)*XX(30) XX(11) = XX(11)-JVS(57)*XX(14)-JVS(105)*XX(19)-JVS(175)*XX(26)-JVS(221)*XX(28)-JVS(272)*XX(31) XX(10) = XX(10)-JVS(174)*XX(26)-JVS(197)*XX(27)-JVS(220)*XX(28) XX(9) = XX(9)-JVS(173)*XX(26)-JVS(196)*XX(27)-JVS(271)*XX(31) XX(8) = XX(8)-JVS(270)*XX(31) XX(7) = XX(7)-JVS(39)*XX(11)-JVS(56)*XX(14)-JVS(62)*XX(15)-JVS(86)*XX(18)-JVS(112)*XX(20)-JVS(195)*XX(27)-JVS(219)& &*XX(28) XX(6) = XX(6)-JVS(43)*XX(12)-JVS(172)*XX(26)-JVS(255)*XX(30) XX(5) = XX(5)-JVS(38)*XX(11)-JVS(55)*XX(14)-JVS(85)*XX(18)-JVS(194)*XX(27)-JVS(218)*XX(28) XX(4) = XX(4)-JVS(171)*XX(26) XX(3) = XX(3)-JVS(170)*XX(26)-JVS(288)*XX(32) XX(2) = XX(2)-JVS(193)*XX(27)-JVS(217)*XX(28) XX(1) = XX(1)-JVS(192)*XX(27)-JVS(241)*XX(29) 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 cbm_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 cbm_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 cbm_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 cbm_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 cbm_LinearAlgebra