Mercurial > octave-nkf
view libcruft/odessa/odessa_sprim.f @ 5083:8386cf9811ee ss-2-1-63
[project @ 2004-11-17 15:49:26 by jwe]
author | jwe |
---|---|
date | Wed, 17 Nov 2004 15:51:27 +0000 |
parents | 7a37caf6ed43 |
children |
line wrap: on
line source
SUBROUTINE ODESSA_SPRIME (NEQ, Y, YH, NYH, NROW, NCOL, WM, IWM, 1 EWT, SAVF, FTEM, DFDP, PAR, F, JAC, DF, PJAC, PDF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION NEQ(*), Y(*), YH(NROW,NCOL,*), WM(*), IWM(*), 1 EWT(*), SAVF(*), FTEM(*), DFDP(NROW,*), PAR(*) EXTERNAL F, JAC, DF, PJAC, PDF PARAMETER (ONE=1.0D0,ZERO=0.0D0) COMMON /ODE001/ ROWND, ROWNS(173), 1 RDUM1(37),EL0, H, RDUM2(6), 2 IOWND1(14), IOWNS(4), 3 IDUM1(3), IERPJ, IDUM2(6), 4 MITER, IDUM3(4), N, IDUM4(5) COMMON /ODE002/ RDUM3(3), 1 IOWND2(3), IDUM5, NSV, IDUM6, NSPE, IDUM7, IERSP, JOPT, IDUM8 C----------------------------------------------------------------------- C ODESSA_SPRIME IS CALLED BY ODESSA TO INITIALIZE THE YH ARRAY. IT IS ALSO C CALLED BY ODESSA_STODE TO REEVALUATE FIRST ORDER DERIVATIVES WHEN KFLAG C .LE. -3. ODESSA_SPRIME COMPUTES THE FIRST DERIVATIVES OF THE SENSITIVITY C COEFFICIENTS WITH RESPECT TO THE INDEPENDENT VARIABLE T... C C ODESSA_SPRIME = D(DY/DP)/DT = JAC*DY/DP + DF/DP C WHERE JAC = JACOBIAN MATRIX C DY/DP = SENSITIVITY MATRIX C DF/DP = INHOMOGENEITY MATRIX C THIS ROUTINE USES THE COMMON VARIABLES EL0, H, IERPJ, MITER, N, C NSV, NSPE, IERSP, JOPT C----------------------------------------------------------------------- C CALL ODESSA_PREPJ WITH JOPT = 1. C IF MITER = 2 OR 5, EL0 IS TEMPORARILY SET TO -1.0 AND H IS C TEMPORARILY SET TO 1.0D0. C----------------------------------------------------------------------- NSPE = NSPE + 1 JOPT = 1 IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 10 HTEMP = H ETEMP = EL0 H = ONE EL0 = -ONE 10 CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, FTEM, 1 PAR, F, JAC, JOPT) IF (IERPJ .NE. 0) GO TO 300 JOPT = 0 IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 20 H = HTEMP EL0 = ETEMP C----------------------------------------------------------------------- C CALL ODESSA_PREPDF AND LOAD DFDP(*,JPAR). C----------------------------------------------------------------------- 20 DO 30 J = 2,NSV JPAR = J - 1 CALL PDF (NEQ, Y, WM, SAVF, FTEM, DFDP(1,JPAR), PAR, 1 F, DF, JPAR) 30 CONTINUE C----------------------------------------------------------------------- C COMPUTE JAC*DY/DP AND STORE RESULTS IN YH(*,*,2). C----------------------------------------------------------------------- GO TO (40,40,310,100,100) MITER C THE JACOBIAN IS FULL.------------------------------------------------ C FOR EACH ROW OF THE JACOBIAN.. 40 DO 70 IROW = 1,N C AND EACH COLUMN OF THE SENSITIVITY MATRIX.. DO 60 J = 2,NSV SUM = ZERO C TAKE THE VECTOR DOT PRODUCT.. DO 50 I = 1,N IPD = IROW + N*(I-1) + 2 SUM = SUM + WM(IPD)*YH(I,J,1) 50 CONTINUE YH(IROW,J,2) = SUM 60 CONTINUE 70 CONTINUE GO TO 200 C THE JACOBIAN IS BANDED.----------------------------------------------- 100 ML = IWM(1) MU = IWM(2) ICOUNT = 1 MBAND = ML + MU + 1 MEBAND = MBAND + ML NMU = N - MU ML1 = ML + 1 C FOR EACH ROW OF THE JACOBIAN.. DO 160 IROW = 1,N IF (IROW .GT. ML1) GO TO 110 IPD = MBAND + IROW + 1 IYH = 1 LBAND = MU + IROW GO TO 120 110 ICOUNT = ICOUNT + 1 IPD = ICOUNT*MEBAND + 2 IYH = IYH + 1 LBAND = LBAND - 1 IF (IROW .LE. NMU) LBAND = MBAND C AND EACH COLUMN OF THE SENSITIVITY MATRIX.. 120 DO 150 J = 2,NSV SUM = ZERO I1 = IPD I2 = IYH C TAKE THE VECTOR DOT PRODUCT. DO 140 I = 1,LBAND SUM = SUM + WM(I1)*YH(I2,J,1) I1 = I1 + MEBAND - 1 I2 = I2 + 1 140 CONTINUE YH(IROW,J,2) = SUM 150 CONTINUE 160 CONTINUE C----------------------------------------------------------------------- C ADD THE INHOMOGENEITY TERM, I.E., ADD DFDP(*,JPAR) TO YH(*,JPAR+1,2). C----------------------------------------------------------------------- 200 DO 220 J = 2,NSV JPAR = J - 1 DO 210 I = 1,N YH(I,J,2) = YH(I,J,2) + DFDP(I,JPAR) 210 CONTINUE 220 CONTINUE RETURN C----------------------------------------------------------------------- C ERROR RETURNS. C----------------------------------------------------------------------- 300 IERSP = -1 RETURN 310 IERSP = -2 RETURN C------------------------END OF SUBROUTINE ODESSA_SPRIME----------------------- END