1C********************************************************** 2 SUBROUTINE HEM5(NQ,NV,NU,NL,FPROB, 3 & T,Q,V,U,A,RLAM,TEND,H, 4 & RTOL,ATOL,ITOL, 5 & SOLOUT,IOUT, 6 & WK,LWK,IWK,LIWK,IDID) 7C----------------------------------------------------------------------- 8C NUMERICAL SOLUTION OF A CONSTRAINED MECHANICAL SYSTEM 9C 10C q' = T(q,t)v 11C M(t,q)v' = f(t,q,v,u) - L(q,v,u,t)*lamda 12C 0 = H(t,q)v + k(t,q) 13C u' = d(q,v,u,lambda,t) 14C 15C 16C THE LOCAL ERROR ESTIMATION AND THE STEP SIZE CONTROL IS BASED ON 17C EMBEDDED FORMULAS OF ORDERS 5 AND 3. THIS METHOD IS PROVIDED WITH 18C TWO DENSE OUTPUT OF ORDER 4. 19C 20C !!! VERSION WITH MA28 !!! 21C 22C AUTHOR: V. BRASEY 23C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES 24C CH-1211 GENEVE 24, SWITZERLAND 25C E-MAIL: BRASEY@DIVSUN.UNIGE.CH, HAIRER@DIVSUN.UNIGE.CH 26C 27C DESCRIPTION OF THIS CODE IS GIVEN IN: 28C V. BRASEY , HALF-EXPLICIT METHOD FOR SEMI-EXPLICIT DIFFERENTIAL- 29C ALGEBRAIC EQUATIONS OF INDEX 2. THESIS, 1994 30C UNIV. DE GENEVE, SECT. DE MATHEMATIQUES. 31C (SEE ALSO ENCLOSED USER-GUIDE) 32C 33C VERSION OF JUNE 8, 1994 34C 35C MANY THANKS TO CH. ENGSTLER FOR CORRECTIONS AND SUGGESTIONS 36C 37C DIMENSIONS 38C ---------- 39C 40C NQ (SIZE OF POSITION VECTOR) 41C NV (SIZE OF VELOCITY VECTOR, NQ>=NV) 42C NL (SIZE OF LAGRANGE MULTIPLIER VECTOR) 43C NU (SIZE OF EXTERNAL DYNAMIC VECTOR) 44C IWK(1) = LDG (SPARSE ALGEBRA: NON ZERO ENTRIES) 45C IWK(2) = LDF 46C 47C 48C SOPHISTICATED SETTING OF PARAMETERS 49C ----------------------------------- 50C SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK 51C WELL. THEY MAY BE DEFINED BY SETTING WK(1),..,WK(8) 52C AS WELL AS IWK(11) .. IWK(16) DIFFERENT FROM ZERO. 53C FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: 54C 55C IWK(11) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. 56C THE DEFAULT VALUE (FOR IWK(11)=0) IS 100000. 57C 58C IWK(12) SWITCH FOR A PROJECTION TO ENSURE CONSISTENT INITIAL VALUE 59C FOR IWK(12)=1 AN INITIAL PROJECTION IS PERFORMED. 60C NO PROJECTION IS DONE IF IWK(12)=0. 61C THE DEFAULT VALUE FOR IWK(12) IS 0. 62C 63C IWK(13) FOR IWK(13).GT.0 IT IS THE NUMBER OF STEPS BETWEEN 64C TWO PROJECTIONS ON THE MANIFOLD DEFINED BY 0 = g(q,t). 65C FOR IWK(13).LE.0 NO PROECTION IS PERFORMED. 66C THE DEFAULT VALUE FOR IWK(13) IS 0. 67C 68C IWK(14) MODE (=0: FULL LINEAR ALGEBRA WITH DEC, =1: IDEM WITH FL, 69C =2: FULL LINEAR ALGEBRA WITH DGETRF, =3: FL 70C =4: SPARSE, =5: IDEM WITH FL) 71C 72C IWK(15) IACC (=1: COMPUTE THE ACCELERATION) 73C 74C IWK(16) IGIIN (=1: COMPUTE NUMERICALLY GII) 75C 76C 77C IWK(21->29) IPAR 78C IPAR(1) = IWK(21) = NMRC (SIZE OF A BLOCK OF AM) 79C IPAR(2) = IWK(22) = NBLK (NUMBER OF BLOCK OF AM) 80C IPAR(3) = IWK(23) = NPGP (0 IF GP AS THE SAME PATTERN AS PREVIOUS CALL) 81C IPAR(4) = IWK(24) = NPFL (0 IF FL AS THE SAME PATTERN AS PREVIOUS CALL) 82C IPAR(5) = IWK(25) = IS (SIZE OF INTEGER WORK SPACE FOR MA28 (MIN 13*NM)) 83C IPAR(6) = IWK(26) = IXS (SIZE OF REAL WORK SPACE FOR MA28 (MIN NM+4*NZA)) 84C IPAR(7) = IWK(27) = PREVL 85C IPAR(8) = IWK(28) = IO 86C IPAR(9) = FLAG TO INDICATE IF UMDFAC HAS BEEN CALLED AT LEAST ONCE 87C 88C IWK(31->38) ISTAT 89C ISTAT(1) = IWK(31) = NSTEP 90C ISTAT(2) = IWK(32) = NACCPT 91C ISTAT(3) = IWK(33) = NREJCT 92C ISTAT(4) = IWK(34) = NFCN 93C ISTAT(5) = IWK(35) = NGQCN 94C ISTAT(6) = IWK(36) = NAMAT 95C ISTAT(7) = IWK(37) = NDEC 96C ISTAT(8) = IWK(38) = NSOL 97C --- NFCN NUMBER OF f - EVALUATIONS 98C --- NGCN NUMBER OF g - EVALUATIONS 99C --- NSTEP NUMBER OF ALL COMPUTED STEPS 100C --- NACCPT NUMBER OF ACCEPTED STEPS 101C --- NREJCT NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP 102C HAS BEEN ACCEPTED) 103C --- NDEC NUMBER OF LU-DECOMPOSITIONS 104C --- NSOL NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS 105C IWK(41->41+NM-1) = IP 106C IWK(29->29+NM+2*LDG-1) = INDG(1..2*LDG) 107C IWK(29+NM+2*LDG..) = INDG1(1..2*LDG) 108C IWK(29+NM+4*LDG..) = INDG2(1..2*LDG) 109C IWK(29+NM+6*LDG..) = INDGD(1..2*LDG) 110C IWK(29+NM+8*LDG..) = INDFL(1..2*LDF) 111C IWK(29+NM+4*LDG+2*LDF..) = INDAM(1..2*NZA) 112C IUMF(1) = IWK(29+NM+4*LDG+2*LDF+2*NZA..) = IKEEP(1..5*NM) 113C IUMF(IIK) = IWK(29+NM+4*LDG+2*LDF+4*NZA+5*NM..) = IW(1..8*NM) 114C IWK(ICONT) = IDOWK 115C 116C WK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. 117C 118C WK(2) THE SAFETY FACTOR IN STEP SIZE PREDICTION, 119C DEFAULT 0.85D0. 120C 121C WK(3), WK(4) PARAMETERS FOR STEP SIZE SELECTION 122C THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION 123C WK(3) <= HNEW/HOLD <= WK(4). 124C DEFAULT VALUES: WK(3)=0.2D0, WK(4)=10.D0 125C 126C WK(6) MAXIMAL STEP SIZE, DEFAULT TEND-T. 127C 128C WK(7) = BETA, DEFAULT 0.D0 129C 130C WK(8) = ALPHA, DEFAULT 1/5 131C 132C---------------------------------------------------------------------- 133C XLA(1..NZA) = AVALUE(1..4*NZA) 134C XUMF(IXW) = XLA(1+2*NZA..) = W(NM) 135C LXLA = 4*NZA +NM 136C 137C----------------------------------------------------------------------- 138C FPROB 139C ----- 140C 141C 0 -> UDOT 142C 1 -> f,M,QDOT,(FL) 143C 2 -> QDOT 144C 3 -> gt 145C 4 -> g 146C 5 -> f,(GII) 147C 6 -> GQ,gt 148C 7 -> (GII),f,M 149C 8 -> f,M,(FL) 150C 9 -> M 151C 10 -> gt,GQ,M,QDOT,(FL) 152C 11 -> gt,GQ,M 153C-------------------------------------------------------------------------- 154C OUTPUT PARAMETERS 155C ----------------- 156C T T-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED 157C (AFTER SUCCESSFUL RETURN T=TEND). 158C 159C Q(N) NUMERICAL APPROXIMATION OF POSITION VECTOR AT T. 160C 161C V(N) NUMERICAL APPROXIMATION OF VELOCITY VECTOR AT T. 162C 163C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP. 164C 165C IDID REPORTS ON SUCCESSFULNESS UPON RETURN: 166C IDID= 1 COMPUTATION SUCCESSFUL, 167C IDID=-1 INPUT IS NOT CONSISTENT, 168C IDID=-2 LARGER NMAX IS NEEDED, 169C IDID=-3 STEP SIZE BECOMES TOO SMALL, 170C IDID=-4 MATRIX IS SINGULAR. 171C IDID=-5 INITIAL PROJECTION: NO CONVERGENCE C----------------------------------------------------------------------- 172C *** *** *** *** *** *** *** *** *** *** *** *** *** 173C DECLARATIONS 174C *** *** *** *** *** *** *** *** *** *** *** *** *** 175 IMPLICIT REAL*8 (A-H,O-Z) 176 DIMENSION Q(NQ),V(NV),U(NU),A(NV),RLAM(NL) 177 DIMENSION ATOL(1),RTOL(1),WK(LWK),IWK(LIWK) 178 LOGICAL ARRET 179 EXTERNAL FPROB,SOLOUT 180C *** *** *** *** *** *** *** 181C SETTING THE PARAMETERS 182C *** *** *** *** *** *** *** 183 ARRET=.FALSE. 184 MODE=IWK(14) 185C -------- MODE, THE CHOICE OF LINEAR ALGEBRA 186 IF ((IWK(14).LE.-1).OR.(IWK(14).GE.6)) THEN 187 WRITE(6,*)' WRONG CHOICE OF IWK(14) (MODE):' 188 & ,IWK(14) 189 ARRET=.TRUE. 190 ELSE 191 MODE=IWK(14) 192 END IF 193 LDG = IWK(1) 194 LDF = IWK(2) 195 NM=NV+NL 196 NMRC = IWK(21) 197 NBLK = IWK(22) 198 NZA = NBLK*NMRC**2+LDF+LDG 199 IF (MODE.EQ.4) LDF=LDG 200 IF (MODE.GE.4) THEN 201 IF (LDG.LE.0) THEN 202 WRITE(6,*)' IWK(1) (LDG) MUST BE POSITIVE' 203 ARRET=.TRUE. 204 END IF 205 IF (LDF.LE.0) THEN 206 WRITE(6,*)' IWK(2) (LDF) MUST BE POSITIVE' 207 ARRET=.TRUE. 208 END IF 209 IF ((NMRC.EQ.0).OR.(NBLK.EQ.0)) THEN 210 WRITE(6,*)' IWK(21) (NMRC) AND IWK(22) (NBLK) MUST 211 & BE POSITIVE' 212 ARRET=.TRUE. 213 END IF 214 IF (IWK(28).EQ.0) IWK(28)=6 215 IF (IWK(25).LT.13*NM) THEN 216 WRITE(6,*)' INTEGER WORK SPACE (IWK(25)) 217 & FOR MA28 TOO SMALL' 218 IWK(25)=13*NM 219 END IF 220 IF (IWK(26).LT.(4*NZA+NM)) THEN 221 WRITE(6,*)' REAL WORK SPACE (IWK(26)) 222 & FOR MA28 TOO SMALL' 223 IWK(26)=4*NZA+NM 224 END IF 225 END IF 226 IF (MODE.LE.3) THEN 227 LDG=NL 228 LDF=NV 229 LDA=NM 230 NDIM=NV 231 MDIM=NL 232 NMDIM=NM 233 ELSE 234 LDA = NBLK*NMRC 235 NDIM=1 236 MDIM=1 237 NMDIM=NMRC 238 END IF 239C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- 240 IF(IWK(11).EQ.0)THEN 241 NMAX=100000 242 ELSE 243 NMAX=IWK(11) 244 END IF 245C -------- IPCIV SWITCH FOR INITIAL PROECTION 246 IF (IWK(12).EQ.1) THEN 247 IPCIV=1 248 ELSE 249 IPCIV=0 250 ENDIF 251C -------- IGII AND IACC 252 IF (IWK(15).EQ.1) THEN 253 IACC=1 254 ELSE 255 IACC=0 256 ENDIF 257 IF (IWK(16).EQ.1) THEN 258 IGII=1 259 ELSE 260 IGII=0 261 ENDIF 262C -------- IPRO NUMBERS OF STEPS BETWEEN TWO PROECTIONS 263 IPROJ = IWK(13) 264C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 265 IF(WK(1).EQ.0.D0)THEN 266 UROUND=1.D-16 267 ELSE 268 UROUND=WK(1) 269 IF(UROUND.LE.1.D-35.OR.UROUND.GE.1.D0)THEN 270 WRITE(6,*)' WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:' 271 & ,WK(1) 272 ARRET=.TRUE. 273 END IF 274 END IF 275C ------- SAFETY FACTOR ------------- 276 IF(WK(2).EQ.0.D0)THEN 277 SAFE=0.9D0 278 ELSE 279 SAFE=WK(2) 280 IF(SAFE.GE.1.D0.OR.SAFE.LE.1.D-4)THEN 281 WRITE(6,*)' CURIOUS INPUT FOR SAFETY FACTOR WK(2)=' 282 & ,WK(2) 283 ARRET=.TRUE. 284 END IF 285 END IF 286C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION 287 IF(WK(3).EQ.0.D0)THEN 288 F1=0.2D0 289 ELSE 290 F1=WK(3) 291 END IF 292 IF(WK(4).EQ.0.D0)THEN 293 F2=10.D0 294 ELSE 295 F2=WK(4) 296 END IF 297C -------- MAXIMAL STEP SIZE 298 IF(WK(6).EQ.0.D0)THEN 299 HMAX=TEND-T 300 ELSE 301 HMAX=WK(6) 302 END IF 303C -------- GUSTAFFSON STRATEGIE 304 BETA = WK(7) 305 IF(WK(8).EQ.0.D0)THEN 306 ALPHA = 1/5.d0 307 ELSE 308 ALPHA=WK(8) 309 END IF 310C -- WK SPACE 311 LXUMF = NM 312 LIUMF = 13*NM+4*nza 313 LIPAR = 9 314 LISTA = 9 315 LRDO = 8+6*(NQ+NV+NU)+NM+2*LDG*NDIM+LDA*NMDIM+NL+ 316 & LDF*MDIM+NZA+LXUMF 317 LIDO = 10+NM+LIPAR+LIUMF+4*LDG+2*LDF+2*NZA 318C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WK ----- 319 IQ1 = 11 320 IQ2 = IQ1+NQ 321 IQ3 = IQ2+NQ 322 IQ4 = IQ3+NQ 323 IQ5 = IQ4+NQ 324 IQ6 = IQ5+NQ 325 IQ7 = IQ6+NQ 326 IQDOT1 = IQ7+NQ 327 IQDOT2 = IQDOT1+NQ 328 IQDOT3 = IQDOT2+NQ 329 IQDOT4 = IQDOT3+NQ 330 IQDOT5 = IQDOT4+NQ 331 IQDOT6 = IQDOT5+NQ 332 IQDOT7 = IQDOT6+NQ 333 IQDOT = IQDOT7+NQ 334 IV1 = IQDOT+NQ 335 IV2 = IV1+NV 336 IV3 = IV2+NV 337 IV4 = IV3+NV 338 IV5 = IV4+NV 339 IV6 = IV5+NV 340 IV7 = IV6+NV 341 IU1 = IV7+NV 342 IU2 = IU1+NU 343 IU3 = IU2+NU 344 IU4 = IU3+NU 345 IU5 = IU4+NU 346 IU6 = IU5+NU 347 IU7 = IU6+NU 348 IVP1 = IU7+NU 349 IVP2 = IVP1+NV 350 IVP3 = IVP2+NV 351 IVP4 = IVP3+NV 352 IVP5 = IVP4+NV 353 IVP6 = IVP5+NV 354 IVP7 = IVP6+NV 355 IXL = IVP7+NV 356 IUDOT1 = IXL+NL 357 IUDOT2 = IUDOT1+NU 358 IUDOT3 = IUDOT2+NU 359 IUDOT4 = IUDOT3+NU 360 IUDOT5 = IUDOT4+NU 361 IUDOT6 = IUDOT5+NU 362 IUDOT7 = IUDOT6+NU 363 IUDOT = IUDOT7+NU 364 IGQ = IUDOT+NU 365 IGQ0 = IGQ+LDG*NDIM 366 IGQ1 = IGQ0+LDG*NDIM 367 IB = IGQ1+LDG*NDIM 368 IX0 = IB+LDA*NMDIM 369 ITEMP = IX0+NM 370 IAM = ITEMP +NV 371 IGT = IAM+LDA*NMDIM 372 IG = IGT+NL 373 IFL = IG+NL 374 IAV = IFL+LDF*MDIM 375 IXUMF = IAV+2*NZA 376 IBD=IXUMF+LXUMF 377 ICQ = IBD+LDA*NMDIM 378 ICV = ICQ+5*NQ 379 ICU = ICV+5*NV 380 IGD = ICU+5*NU 381 IGTD = IGD+LDG*NDIM 382 IXD = IGTD+NL 383 IQD = IXD+NM 384 IVD = IQD+NQ 385 IUD = IVD+NV 386 IVP = IUD+NU 387 IDO = IVP+NV 388C ------ TOTAL STORAGE REQUIREMENT ----------- 389 IS=IDO+LRDO 390 IF(IS.GT.LWK)THEN 391 WRITE(6,*)' INSUFFICIENT STORAGE FOR WK, MIN. LWK=',IS 392 ARRET=.TRUE. 393 END IF 394C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN IWK ----- 395 IPA = 21 396 ISTAT = 31 397 IIP = 41 398 ING = 41+NM 399 ING0 = ING+2*LDG 400 ING1 = ING0+2*LDG 401 INGD = ING1+2*LDG 402 INFL = INGD+2*LDG 403 INAM = INFL+2*LDF 404 INUMF =INAM+2*NZA 405 ICONT = INUMF+LIUMF 406C ------ TOTAL STORAGE REQUIREMENT ----------- 407 IS=ICONT+LIDO 408 IF(IS.GT.LIWK)THEN 409 WRITE(6,*)' INSUFFICIENT STORAGE FOR IWK, MIN. LIWK=',IS 410 ARRET=.TRUE. 411 END IF 412C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 413 IF (ARRET) THEN 414 IDID=-1 415 RETURN 416 END IF 417C -------- CALL TO CORE INTEGRATOR ------------ 418 CALL HCOR(NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM,NZA,LRDO,LIDO,LIPAR, 419 & LISTA,LIUMF,LXUMF,LDG,LDF,LDA,MODE,NMAX,IPCIV,IPROJ,IOUT, 420 & IGII,IACC,ITOL,IWK(IPA),IWK(ISTAT),IWK(IIP),IWK(ING), 421 & IWK(ING0),IWK(ING1),IWK(INGD),IWK(INFL),IWK(INAM),IWK(INUMF), 422 & IWK(ICONT),FPROB,SOLOUT,UROUND,T,TEND,H,HMAX,RTOL,ATOL,SAFE, 423 & ALPHA,BETA,F1,F2,Q,V,U,A,RLAM,WK(IQ1),WK(IQ2),WK(IQ3), 424 & WK(IQ4),WK(IQ5),WK(IQ6),WK(IQ7),WK(IQDOT1),WK(IQDOT2), 425 & WK(IQDOT3),WK(IQDOT4),WK(IQDOT5),WK(IQDOT6), 426 & WK(IQDOT7),WK(IQDOT),WK(IV1),WK(IV2),WK(IV3),WK(IV4), 427 & WK(IV5),WK(IV6),WK(IV7),WK(IU1),WK(IU2),WK(IU3), 428 & WK(IU4),WK(IU5),WK(IU6),WK(IU7),WK(IVP1),WK(IVP2), 429 & WK(IVP3),WK(IVP4),WK(IVP5),WK(IVP6),WK(IVP7),WK(IXL), 430 & WK(IUDOT1),WK(IUDOT2),WK(IUDOT3),WK(IUDOT4),WK(IUDOT5), 431 & WK(IUDOT6),WK(IUDOT7),WK(IUDOT),WK(IGQ),WK(IGQ0), 432 & WK(IGQ1),WK(IB),WK(IX0),WK(ITEMP),WK(IAM), 433 & WK(IGT),WK(IG),WK(IFL),WK(IAV),WK(IXUMF),WK(IBD), 434 & WK(ICQ),WK(ICV),WK(ICU),WK(IGD),WK(IGTD),WK(IXD), 435 & WK(IQD),WK(IVD),WK(IUD),WK(IVP),WK(IDO)) 436C ----------- RETURN ----------- 437 RETURN 438 END 439C 440C 441C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- 442C 443 SUBROUTINE HCOR(NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM,NZA,LRDO, 444 & LIDO,LIPAR,LISTA,LIUMF,LXUMF,LDG,LDF,LDA,MODE,NMAX, 445 & IPCIV,IPROJ,IOUT,IGII,IACC,ITOL,IPAR,ISTAT,IP,INDG, 446 & INDG0,INDG1,INDGD,INDFL,INDA,IUMF,IDOWK,FPROB,SOLOUT, 447 & UROUND,T,TEND,H,HMAX,RTOL,ATOL,SAFE,ALPHA,BETA, 448 & FAC1,FAC2,Q,V,U,A,RLAM,Q1,Q2,Q3,Q4,Q5,Q6,Q7, 449 & QDOT1,QDOT2,QDOT3,QDOT4,QDOT5,QDOT6,QDOT7,QDOT,V1, 450 & V2,V3,V4,V5,V6,V7,U1,U2,U3,U4,U5,U6,U7,VP1,VP2, 451 & VP3,VP4,VP5,VP6,VP7,XL,UDOT1,UDOT2,UDOT3,UDOT4, 452 & UDOT5,UDOT6,UDOT7,UDOT,GQ,GQ0,GQ1,B,X0,TEMP, 453 & AM,GT,G,FL,AVALUE,XUMF,BD,CONTQ,CONTV,CONTU, 454 & GD,GTD,XD,QD,VD,UD,VP,DOWK) 455C ---------------------------------------------------------- 456C CORE INTEGRATOR FOR HEM5 457C PARAMETERS SAME AS IN HEM5 WITH WORKSPACE ADDED 458C ---------------------------------------------------------- 459C DECLARATIONS 460C ---------------------------------------------------------- 461 IMPLICIT REAL*8 (A-H,O-Z) 462 DIMENSION Q(NQ),Q1(NQ),Q2(NQ),Q3(NQ),Q4(NQ),Q5(NQ) 463 DIMENSION Q6(NQ),Q7(NQ),V(NV),V1(NV),V2(NV),V3(NV) 464 DIMENSION V4(NV),V5(NV),V6(NV),V7(NV),VP1(NV),VP2(NV) 465 DIMENSION VP3(NV),VP4(NV),VP5(NV),VP6(NV),VP7(NV) 466 DIMENSION GQ(LDG,NDIM),GQ0(LDG,NDIM),GQ1(LDG,NDIM) 467 DIMENSION B(LDA,NMDIM),AM(LDA,NMDIM),INDG(2*LDG) 468 DIMENSION X0(NM),IP(NM),TEMP(NV),AVALUE(2*NZA),IDOWK(LIDO) 469 DIMENSION GT(NL),G(NL),B10(4),CONTQ(5*NQ),CONTV(5*NV) 470 DIMENSION GTD(NL),XD(NM),QD(NQ),VD(NV),VP(NV),DOWK(LRDO) 471 DIMENSION ATOL(1),RTOL(1),XUMF(LXUMF),FL(LDF,MDIM) 472 DIMENSION GD(LDG,NDIM),BD(LDA,NMDIM),IPAR(LIPAR) 473 DIMENSION INDG0(2*LDG),INDG1(2*LDG),INDGD(2*LDG),A(NV) 474 DIMENSION INDFL(2*LDF),INDA(2*NZA),IUMF(LIUMF),ISTAT(LISTA) 475 DIMENSION UDOT1(NU),UDOT2(NU),UDOT3(NU),UDOT4(NU) 476 DIMENSION UDOT5(NU),UDOT6(NU),UDOT7(NU),UDOT(NU),RLAM(NL) 477 DIMENSION QDOT1(NQ),QDOT2(NQ),QDOT3(NQ),QDOT4(NQ) 478 DIMENSION QDOT5(NQ),QDOT6(NQ),QDOT7(NQ),QDOT(NQ) 479 DIMENSION U(NU),U1(NU),U2(NU),U3(NU),U4(NU),U5(NU) 480 DIMENSION U6(NU),U7(NU),CONTU(5*NU),UD(NU),XL(NL) 481 LOGICAL REJECT,ACCEPT,LAST 482 EXTERNAL FPROB,SOLOUT 483C *** *** *** *** *** *** *** 484C INITIALISATIONS 485C *** *** *** *** *** *** *** 486 CALL COEF(C2,C3,C4,C5,C6,C7,C10,A21,A31,A32, 487 & A41,A42,A43,A51,A52,A53,A54,A61,A62,A63,A64,A65, 488 & A71,A72,A73,A74,A75,A76,A81,A82,A83,A84,A85,A86, 489 & A87,A106,A107,A109,B1,B2,B3,B4,B5,B6,B7,B8, 490 & D14,D15,D16,D17,D18,D19,D24,D25,D26,D27,D28,D29, 491 & D34,D35,D36,D37,D38,D39,D44,D45,D46,D47,D48,D49) 492 B10(1)=D19 493 B10(2)=D29 494 B10(3)=D39 495 B10(4)=D49 496 DO I=1,9 497 ISTAT(I)=0 498 END DO 499 NSTEP=0 500 NQ2=2*NQ 501 NQ3=3*NQ 502 NQ4=4*NQ 503 NV2=2*NV 504 NV3=3*NV 505 NV4=4*NV 506 NU2=2*NU 507 NU3=3*NU 508 NU4=4*NU 509 NQVU=NQ+NV+NU 510 NQVU2=2*NQVU 511 NQVU3=3*NQVU 512 NQVU4=4*NQVU 513 NBLK = IPAR(2) 514 NMRC = IPAR(1) 515cc 516CC ipar(3)=1 517cc 518 IPAR(9)=1 519C -- NBS : COUNTS THE STEPS BETWEEN THE PROJECTIONS 520 NBS = 0 521 LAST=.FALSE. 522 ACCEPT = .TRUE. 523 REJECT=.FALSE. 524 FACOLD=1.D0 525 FACC1=1.D0/FAC1 526 FACC2=1.D0/FAC2 527 POSNEG=SIGN(1.D0,TEND-T) 528 IRTRN=1 529 HMAX=ABS(HMAX) 530 H=MIN(MAX(1.D-4,ABS(H)),HMAX) 531 H=SIGN(H,POSNEG) 532 TOLD=T 533 CALL FPROB(10,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 534 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 535 & INDFL(LDF+1),T,Q,V,U,XL, 536 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT1,UDOT1,B) 537 ISTAT(5)=1 538 ISTAT(6)=1 539C --- PROJECTION TO ENSURE CONSISTENT INITIAL VALUES 540 IF ((IPCIV.EQ.1).OR.(IOUT.EQ.2)) THEN 541 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA, 542 & NZA,LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG, 543 & INDG,INDFL,INDA,IP,IUMF,XUMF, 544 & B,GQ,GQ,FL,AVALUE,IER) 545 IF (IER.NE.0) GOTO 176 546 IF (IPCIV.EQ.1) THEN 547 CALL APROJ(MODE,NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM, 548 * NBLK,NMRC,LDG,LDF,LDA,NZA,LIPAR,LIUMF,LXUMF, 549 * LISTA,ISTAT,IPAR,IUMF,INDA,INDG,INDFL,FL,XUMF, 550 * AVALUE,T,FPROB,Q,Q,Q2,QD,V,V,V2,U,XL,UDOT,G,GT, 551 * GQ,AM,B,X0,IP,ATOL,RTOL,ITOL,ACCEPT) 552 ISTAT(6) = ISTAT(6)+2 553 ISTAT(5) = ISTAT(5)+1 554 IF (.NOT.(ACCEPT)) GOTO 179 555 END IF 556 IF (IOUT.EQ.2) THEN 557 CALL FPROB(5,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 558 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 559 & INDFL(LDF+1),T,Q,V,U,XL, 560 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 561 IF (IGII.EQ.1) CALL GIINUM(MODE,NQ,NV,NU,NL,NDIM, 562 & MDIM,NMDIM,NM, 563 & LDG,LDF,LDA,NBLK,NMRC,IPAR(3),IPAR(4),INDG,INDGD, 564 & INDFL,T,UROUND,FPROB,Q,QD,QDOT1,V,U,GQ,GD,GT,GTD, 565 & FL,AM,X0(NV+1),G) 566 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF, 567 & IP,NZA,AVALUE,XUMF,B,X0,VP,XL) 568 ISTAT(4)=1 569 ISTAT(9)=1 570 END IF 571 ENDIF 572 IF(IOUT.GE.1) THEN 573 DOWK(1)=TOLD 574 DOWK(2)=H 575 CALL SOLOUT(MODE,ISTAT(2)+1,NQ,NV,NU,NL,LDG,LDF,LDA, 576 & LRDO,LIDO,FPROB,Q,V,U,DOWK,IDOWK) 577 END IF 578C -- FIRST STAGE : V1 AND Q1 CONTAIN THE INITIAL VALUES 579 DO I=1,NV 580 Q1(I) = Q(I) 581 V1(I) = V(I) 582 END DO 583 DO I=NV+1,NQ 584 Q1(I) = Q(I) 585 END DO 586 DO I=1,NU 587 U1(I) = U(I) 588 END DO 589C --- BASIC INTEGRATION STEP 590 909 IF(NSTEP.GT.NMAX) GOTO 178 591 IF (0.1D0*ABS(H).LE.ABS(T)*UROUND) GOTO 177 592 IF((T+H-TEND)*POSNEG.GT.0.D0) THEN 593 H=TEND-T 594 LAST=.TRUE. 595 END IF 596 ISTAT(1)=ISTAT(1)+1 597 ACCEPT = .TRUE. 598C -- 2d STAGE 599 CALL FPROB(8,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 600 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 601 & INDFL(LDF+1),T,Q1,V1,U1,XL, 602 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT1,UDOT1,B) 603 DO I=1,NQ 604 Q2(I) = Q1(I) + H*A21*QDOT1(I) 605 END DO 606 TCH = T+C2*H 607 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 608 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 609 & INDFL(LDF+1),TCH,Q2,V2,U2,XL, 610 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT1,UDOT1,B) 611 FAC = -1.D0/(H*A21) 612 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG1, 613 & GQ1,FAC,V1,GT,X0(NV+1)) 614 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA, 615 & NZA,LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG,INDG1, 616 & INDFL,INDA,IP,IUMF,XUMF,B,GQ,GQ1,FL,AVALUE,IER) 617 IF (IER.NE.0) GOTO 176 618 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF, 619 & IP,NZA,AVALUE,XUMF,B,X0,VP1,XL) 620cc 621 ipar(3)=0 622cc 623 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 624 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 625 & INDFL(LDF+1),T,Q1,V1,U1,XL, 626 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT1,UDOT1,B) 627 DO I=1,NU 628 U2(I) = U1(I) + H*A21*UDOT1(I) 629 END DO 630 DO I=1,NV 631 V2(I) = V1(I) + H*A21*VP1(I) 632 TEMP(I) = V1(I) + H*A31*VP1(I) 633 END DO 634C - II.3D STAGE 635 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 636 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 637 & INDFL(LDF+1),TCH,Q2,V2,U2,XL, 638 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT2,UDOT2,B) 639 DO I=1,NQ 640 Q3(I) = Q1(I) + H*(A31*QDOT1(I)+A32*QDOT2(I)) 641 END DO 642 TCH = T+C3*H 643 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 644 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 645 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 646 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT3,UDOT3,B) 647 FAC = -1.D0/(H*A32) 648 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG0, 649 & GQ0,FAC,TEMP,GT,X0(NV+1)) 650 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA, 651 & NZA,LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG1,INDG0, 652 & INDFL,INDA,IP,IUMF,XUMF,B,GQ1,GQ0,FL,AVALUE,IER) 653 IF (IER.NE.0) GOTO 176 654 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF, 655 & IP,NZA,AVALUE,XUMF,B,X0,VP2,XL) 656 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 657 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 658 & INDFL(LDF+1),T+C2*H,Q2,V2,U2,XL, 659 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT2,UDOT2,B) 660 DO I=1,NU 661 U3(I) = U1(I) + H*(A31*UDOT1(I)+A32*UDOT2(I)) 662 END DO 663 DO I=1,NV 664 V3(I) = TEMP(I) + H*A32*VP2(I) 665 TEMP(I) = V1(I) + H*(A41*VP1(I)+A42*VP2(I)) 666 END DO 667C - II. 4TH STAGE 668 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 669 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 670 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 671 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT3,UDOT3,B) 672 DO I=1,NQ 673 Q4(I) = Q1(I) + H*(A41*QDOT1(I)+A42*QDOT2(I)+A43*QDOT3(I)) 674 END DO 675 TCH = T+C4*H 676 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 677 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 678 & INDFL(LDF+1),TCH,Q4,V4,U4,XL, 679 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT4,UDOT4,B) 680 FAC = -1.D0/(H*A43) 681 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG1, 682 & GQ1,FAC,TEMP,GT,X0(NV+1)) 683 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 684 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG0,INDG1,INDFL, 685 & INDA,IP,IUMF,XUMF,B,GQ0,GQ1,FL,AVALUE,IER) 686 IF (IER.NE.0) GOTO 176 687 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF, 688 & IP,NZA,AVALUE,XUMF,B,X0,VP3,XL) 689 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 690 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 691 & INDFL(LDF+1),T+C3*H,Q3,V3,U3,XL, 692 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT3,UDOT3,B) 693 DO I=1,NU 694 U4(I) = U1(I) + H*(A41*UDOT1(I)+A42*UDOT2(I)+A43*UDOT3(I)) 695 END DO 696 DO I=1,NV 697 V4(I) = TEMP(I) + H*A43*VP3(I) 698 TEMP(I) = V1(I) + H*(A51*VP1(I)+A52*VP2(I) 699 & +A53*VP3(I)) 700 701 END DO 702C - II. 5TH STAGE 703 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 704 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 705 & INDFL(LDF+1),TCH,Q4,V4,U4,XL, 706 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT4,UDOT4,B) 707 DO I=1,NQ 708 Q5(I) = Q1(I)+H*(A51*QDOT1(I)+A52*QDOT2(I) 709 & +A53*QDOT3(I)+A54*QDOT4(I)) 710 END DO 711 TCH = T+C5*H 712 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 713 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 714 & INDFL(LDF+1),TCH,Q5,V5,U5,XL, 715 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT5,UDOT5,B) 716 FAC = -1.D0/(H*A54) 717 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG0, 718 & GQ0,FAC,TEMP,GT,X0(NV+1)) 719 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 720 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG1,INDG0, 721 & INDFL,INDA,IP,IUMF,XUMF,B,GQ1,GQ0,FL,AVALUE,IER) 722 IF (IER.NE.0) GOTO 176 723 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF, 724 & IP,NZA,AVALUE,XUMF,B,X0,VP4,XL) 725 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 726 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 727 & INDFL(LDF+1),T+C4*H,Q4,V4,U4,XL, 728 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT4,UDOT4,B) 729 DO I=1,NU 730 U5(I) = U1(I)+H*(A51*UDOT1(I)+A52*UDOT2(I) 731 & +A53*UDOT3(I)+A54*UDOT4(I)) 732 END DO 733 DO I=1,NV 734 V5(I) = TEMP(I) + H*A54*VP4(I) 735 TEMP(I) = V1(I) + H*(A61*VP1(I)+A62*VP2(I) 736 & +A63*VP3(I)+A64*VP4(I)) 737 END DO 738C - II. 6TH STAGE 739 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 740 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 741 & INDFL(LDF+1),TCH,Q5,V5,U5,XL, 742 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT5,UDOT5,B) 743 DO I=1,NQ 744 Q6(I) = Q1(I) + H*(A61*QDOT1(I)+A62*QDOT2(I)+ 745 & A63*QDOT3(I)+A64*QDOT4(I)+A65*QDOT5(I)) 746 END DO 747 TCH = T+C6*H 748 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 749 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 750 & INDFL(LDF+1),TCH,Q6,V6,U6,XL, 751 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT6,UDOT6,B) 752 FAC = -1.D0/(H*A65) 753 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG1, 754 & GQ1,FAC,TEMP,GT,X0(NV+1)) 755 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 756 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG0,INDG1, 757 & INDFL,INDA,IP,IUMF,XUMF,B,GQ0,GQ1,FL,AVALUE,IER) 758 IF (IER.NE.0) GOTO 176 759 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 760 & NZA,AVALUE,XUMF,B,X0,VP5,XL) 761 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 762 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 763 & INDFL(LDF+1),T+C5*H,Q5,V5,U5,XL, 764 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT5,UDOT5,B) 765 DO I=1,NU 766 U6(I) = U1(I)+H*(A61*UDOT1(I)+A62*UDOT2(I)+ 767 & A63*UDOT3(I)+A64*UDOT4(I)+A65*UDOT5(I)) 768 END DO 769 DO I=1,NV 770 V6(I) = TEMP(I) + H*A65*VP5(I) 771 TEMP(I) = V1(I) + H*(A71*VP1(I)+A72*VP2(I) 772 & +A73*VP3(I)+A74*VP4(I)+A75*VP5(I)) 773 END DO 774C - II. 7TH STAGE 775 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 776 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 777 & INDFL(LDF+1),TCH,Q6,V6,U6,XL, 778 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT6,UDOT6,B) 779 DO I=1,NQ 780 Q7(I) = Q1(I) + H*(A71*QDOT1(I)+A72*QDOT2(I)+ 781 & A73*QDOT3(I)+A74*QDOT4(I)+A75*QDOT5(I)+ 782 & A76*QDOT6(I)) 783 END DO 784 TCH = T+C7*H 785 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 786 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 787 & INDFL(LDF+1),TCH,Q7,V7,U7,XL, 788 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT7,UDOT7,B) 789 FAC = -1.D0/(H*A76) 790 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG0, 791 & GQ0,FAC,TEMP,GT,X0(NV+1)) 792 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 793 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG1,INDG0, 794 & INDFL,INDA,IP,IUMF,XUMF,B,GQ1,GQ0,FL,AVALUE,IER) 795 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 796 & NZA,AVALUE,XUMF,B,X0,VP6,XL) 797 IF (IER.NE.0) GOTO 176 798 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 799 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 800 & INDFL(LDF+1),T+C6*H,Q6,V6,U6,XL, 801 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT6,UDOT6,B) 802 DO I=1,NU 803 U7(I) = U1(I)+ H*(A71*UDOT1(I)+A72*UDOT2(I)+ 804 & A73*UDOT3(I)+A74*UDOT4(I)+A75*UDOT5(I)+ 805 & A76*UDOT6(I)) 806 END DO 807 DO I=1,NV 808 V7(I) = TEMP(I) + H*A76*VP6(I) 809 TEMP(I) = V1(I) + H*(A81*VP1(I)+A82*VP2(I) 810 & +A83*VP3(I)+A84*VP4(I)+A85*VP5(I)+A86*VP6(I)) 811 END DO 812C - II. 8TH STAGE 813 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 814 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 815 & INDFL(LDF+1),TCH,Q7,V7,U7,XL, 816 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT7,UDOT7,B) 817 DO I=1,NQ 818 Q2(I) = Q1(I) + H*(A81*QDOT1(I)+A82*QDOT2(I)+ 819 & A83*QDOT3(I)+A84*QDOT4(I)+A85*QDOT5(I) 820 & +A86*QDOT6(I)+A87*QDOT7(I)) 821 END DO 822 TH = T+H 823 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 824 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 825 & INDFL(LDF+1),TH,Q2,V2,U2,XL, 826 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT2,UDOT2,B) 827 FAC = -1.D0/(H*A87) 828 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG1, 829 & GQ1,FAC,TEMP,GT,X0(NV+1)) 830 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 831 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG0,INDG1, 832 & INDFL,INDA,IP,IUMF,XUMF,B,GQ0,GQ1,FL,AVALUE,IER) 833 IF (IER.NE.0) GOTO 176 834 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 835 & NZA,AVALUE,XUMF,B,X0,VP7,XL) 836 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 837 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 838 & INDFL(LDF+1),T+C7*H,Q7,V7,U7,XL, 839 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT7,UDOT7,B) 840 DO I=1,NU 841 U2(I) = U1(I)+ H*(A81*UDOT1(I)+A82*UDOT2(I)+ 842 & A83*UDOT3(I)+A84*UDOT4(I)+A85*UDOT5(I) 843 & +A86*UDOT6(I)+A87*UDOT7(I)) 844 END DO 845 DO I=1,NV 846 V2(I) = TEMP(I) + H*A87*VP7(I) 847 TEMP(I) = V1(I) + H*(B4*VP4(I)+B5*VP5(I) 848 & +B6*VP6(I)+B7*VP7(I)) 849 END DO 850C - II. LAST STAGE 851 CALL FPROB(1,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 852 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 853 & INDFL(LDF+1),TH,Q2,V2,U2,XL, 854 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT2,UDOT2,B) 855 DO I=1,NQ 856 Q(I) = Q1(I) +H*(B4*QDOT4(I)+B5*QDOT5(I)+B6*QDOT6(I) 857 & +B7*QDOT7(I)+B8*QDOT2(I)) 858 END DO 859 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 860 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 861 & INDFL(LDF+1),TH,Q,V,U,XL, 862 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 863 FAC = -1.D0/(H*B8) 864 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG0, 865 & GQ0,FAC,TEMP,GT,X0(NV+1)) 866 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 867 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG1,INDG0, 868 & INDFL,INDA,IP,IUMF,XUMF,B,GQ1,GQ0,FL,AVALUE,IER) 869 IF (IER.NE.0) GOTO 176 870 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 871 & NZA,AVALUE,XUMF,B,X0,VP2,XL) 872 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 873 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 874 & INDFL(LDF+1),TH,Q2,V2,U2,XL, 875 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT2,UDOT2,B) 876 DO I=1,NU 877 U(I) = U1(I)+H*(B4*UDOT4(I)+B5*UDOT5(I)+B6*UDOT6(I) 878 & +B7*UDOT7(I)+B8*UDOT2(I)) 879 END DO 880 DO I=1,NV 881 V(I) = TEMP(I) + H*B8*VP2(I) 882 END DO 883C -- STATISTICS 884 ISTAT(4) = ISTAT(4)+8 885 ISTAT(5) = ISTAT(5)+8 886 ISTAT(6) = ISTAT(6)+8 887 ISTAT(8) = ISTAT(8)+8 888C --- ERROR ESTIMATION 889 ERR1 = 0.D0 890 ERR2 = 0.D0 891 DO I=1,NV 892 IF (ITOL.EQ.0) THEN 893 SK1 = ATOL(1)+RTOL(1)*MAX(DABS(Q(I)),ABS(Q1(I))) 894 SK2 = ATOL(1)+RTOL(1)*MAX(DABS(V(I)),ABS(V1(I))) 895 ELSE 896 SK1 = ATOL(I)+RTOL(I)*MAX(ABS(Q(I)),ABS(Q1(I))) 897 SK2 = ATOL(NQ+I)+RTOL(NQ+I)*MAX(ABS(V(I)),ABS(V1(I))) 898 END IF 899 ERR1 = ERR1+((Q(I)-Q2(I))/SK1)**2+ 900 * ((V(I)-V2(I))/SK2)**2 901 QII= Q1(I)+H*(2.5D0*QDOT7(I)-1.5D0*QDOT2(I)) 902 VII= V1(I)+H*(2.5D0*VP7(I)-1.5D0*VP2(I)) 903 ERR2 = ERR2+((Q(I)-QII)/SK1)**2+ 904 * ((V(I)-VII)/SK2)**2 905 END DO 906 DO I=NV+1,NQ 907 IF (ITOL.EQ.0) THEN 908 SK1 = ATOL(1)+RTOL(1)*MAX(DABS(Q(I)),ABS(Q1(I))) 909 ELSE 910 SK1 = ATOL(I)+RTOL(I)*MAX(ABS(Q(I)),ABS(Q1(I))) 911 END IF 912 ERR1 = ERR1+((Q(I)-Q2(I))/SK1)**2 913 QII= Q1(I)+H*(2.5D0*QDOT7(I)-1.5D0*QDOT2(I)) 914 ERR2 = ERR2+((Q(I)-QII)/SK1)**2 915 END DO 916 NQV=NQ+NV 917 DO I=1,NU 918 IF (ITOL.EQ.0) THEN 919 SK1 = ATOL(1)+RTOL(1)*MAX(DABS(U(I)),ABS(U1(I))) 920 ELSE 921 SK1 = ATOL(NQV+I)+RTOL(NQV+I)*MAX(ABS(U(I)),ABS(U1(I))) 922 END IF 923 ERR1 = ERR1+((U(I)-U2(I))/SK1)**2 924 UII= U1(I)+H*(2.5D0*UDOT7(I)-1.5D0*UDOT2(I)) 925 ERR2 = ERR2+((U(I)-UII)/SK1)**2 926 END DO 927 ERR1 = SQRT(ERR1/(NQV+NU)) 928 ERR2 = SQRT(ERR2/(NQV+NU)) 929 ERR =ERR1+0.01D0*ERR2 930c 931 EPS=1.D-13 932 IF (ERR.LT.EPS) THEN 933 HNEW=1.2D0*H 934 GOTO 222 935 END IF 936 ERR=(ERR1**2)/ERR 937C --- COMPUTATION OF HNEW WITH GUSTAFSSON STABILIZATION 938 FAC11=ERR**ALPHA 939 FAC=FAC11/FACOLD**BETA 940C --- WE REQUIRE FAC1 <= HNEW/H <= FAC2 941 FAC=MAX(FACC2,MIN(FACC1,FAC/SAFE)) 942 HNEW=H/FAC 943 IF (ERR.GT.1.D0) ACCEPT=.FALSE. 944c 945 222 CONTINUE 946c 947 IF (ACCEPT) THEN 948 NBS = NBS+1 949 IF (NBS.EQ.IPROJ) THEN 950 NBS = 0 951C -- PROJECTION ON MANIFOLD DEFINED BY g(q,t)=0 952 CALL APROJ(MODE,NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM,NBLK, 953 * NMRC,LDG,LDF,LDA,NZA,LIPAR,LIUMF,LXUMF,LISTA, 954 * ISTAT,IPAR,IUMF,INDA,INDG0,INDFL,FL,XUMF, 955 * AVALUE,T,FPROB,Q,Q1,Q2,QD,V,V1,V2,U,XL,UDOT, 956 * G,GT,GQ0,AM,B,X0,IP,ATOL,RTOL,ITOL,ACCEPT) 957 ISTAT(6) = ISTAT(6)+2 958 ISTAT(5) = ISTAT(5)+1 959 END IF 960 IF (ACCEPT) THEN 961C --- STEP IS ACCEPTED 962 FACOLD=MAX(ERR,1.D-4) 963 ISTAT(2)=ISTAT(2)+1 964 CALL FPROB(2,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 965 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 966 & INDFL(LDF+1),TH,Q,V,U,XL, 967 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 968 IF ((IOUT.EQ.1).OR.(IOUT.EQ.2)) THEN 969C --- COMPUTE INTERMEDIATE SUM FOR DENSE OUTPUT ONLY 970 DO I=1,NV 971 TEMP(I) = V1(I) + H*(A106*VP6(I)+A107*VP7(I)) 972 END DO 973 CALL FPROB(8,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 974 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 975 & INDFL(LDF+1),TH,Q,V,U,XL, 976 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 977 DO I=1,NQ 978 Q3(I) = Q1(I) +H*(A106*QDOT6(I)+A107*QDOT7(I) 979 & +A109*QDOT(I)) 980 END DO 981 TCH = T+C10*H 982 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 983 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 984 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 985 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 986 FAC = -1.D0/(H*A109) 987 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG1, 988 & GQ1,FAC,TEMP,GT,X0(NV+1)) 989 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 990 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG0,INDG1, 991 & INDFL,INDA,IP,IUMF,XUMF,B,GQ0,GQ1,FL,AVALUE,IER) 992 IF (IER.NE.0) GOTO 176 993 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF, 994 & IP,NZA,AVALUE,XUMF,B,X0,VP1,XL) 995 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 996 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 997 & INDFL(LDF+1),TH,Q,V,U,XL, 998 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 999 DO I=1,NU 1000 U3(I) = U1(I)+H*(A106*UDOT6(I)+A107*UDOT7(I) 1001 & +A109*UDOT(I)) 1002 END DO 1003 DO I=1,NV 1004 V3(I) = TEMP(I) + H*A109*VP1(I) 1005 END DO 1006 CALL FPROB(2,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1007 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 1008 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 1009 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT3,UDOT3,B) 1010 DO I=1,NQ 1011 CONTQ(I) = Q1(I) 1012 CONTQ(I+NQ)= D14*QDOT4(I)+D15*QDOT5(I)+D16*QDOT6(I)+ 1013 & D17*QDOT7(I)+D18*QDOT2(I)+D19*QDOT3(I) 1014 CONTQ(I+NQ2)= D24*QDOT4(I)+D25*QDOT5(I)+D26*QDOT6(I)+ 1015 & D27*QDOT7(I)+D28*QDOT2(I)+D29*QDOT3(I) 1016 CONTQ(I+NQ3)= D34*QDOT4(I)+D35*QDOT5(I)+D36*QDOT6(I)+ 1017 & D37*QDOT7(I)+D38*QDOT2(I)+D39*QDOT3(I) 1018 CONTQ(I+NQ4)= D44*QDOT4(I)+D45*QDOT5(I)+D46*QDOT6(I)+ 1019 & D47*QDOT7(I)+D48*QDOT2(I)+D49*QDOT3(I) 1020 END DO 1021 DO I=1,NV 1022 CONTV(I)= V1(I) 1023 CONTV(I+NV)= D14*VP4(I)+D15*VP5(I)+D16*VP6(I) 1024 & +D17*VP7(I)+D18*VP2(I) 1025 CONTV(I+NV2)= D24*VP4(I)+D25*VP5(I)+D26*VP6(I) 1026 & +D27*VP7(I)+D28*VP2(I) 1027 CONTV(I+NV3)= D34*VP4(I)+D35*VP5(I)+D36*VP6(I) 1028 & +D37*VP7(I)+D38*VP2(I) 1029 CONTV(I+NV4)= D44*VP4(I)+D45*VP5(I)+D46*VP6(I) 1030 & +D47*VP7(I)+D48*VP2(I) 1031 END DO 1032C-- STATISTICS 1033 ISTAT(4)=ISTAT(4)+1 1034 ISTAT(5)=ISTAT(5)+1 1035 ISTAT(6)=ISTAT(6)+1 1036 ISTAT(8)=ISTAT(8)+1 1037 END IF 1038 IF (IOUT.EQ.1) THEN 1039 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1040 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 1041 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 1042 & G,GQ0,X0(1),X0(NV+1),GT,FL,QDOT3,UDOT3,B) 1043 DO I=1,NU 1044 CONTU(I)= U1(I) 1045 CONTU(I+NU)= D14*UDOT4(I)+D15*UDOT5(I)+D16*UDOT6(I) 1046 & +D17*UDOT7(I)+D18*UDOT2(I)+D19*UDOT3(I) 1047 CONTU(I+NU2)= D24*UDOT4(I)+D25*UDOT5(I)+D26*UDOT6(I) 1048 & +D27*UDOT7(I)+D28*UDOT2(I)+D29*UDOT3(I) 1049 CONTU(I+NU3)= D34*UDOT4(I)+D35*UDOT5(I)+D36*UDOT6(I) 1050 & +D37*UDOT7(I)+D38*UDOT2(I)+D39*UDOT3(I) 1051 CONTU(I+NU4)= D44*UDOT4(I)+D45*UDOT5(I)+D46*UDOT6(I) 1052 & +D47*UDOT7(I)+D48*UDOT2(I)+D49*UDOT3(I) 1053 END DO 1054 CALL CPCT(MODE,NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM,LDG,LDF, 1055 & LDA,NZA,LIPAR,LXUMF,LIUMF,LRDO,LIDO,IUMF,IPAR, 1056 & ISTAT,INDG1,INDFL,IP,IDOWK,H,T,TCH,B10,CONTQ, 1057 & CONTV,CONTU,Q3,V3,U3,GQ1,FL,AVALUE,XUMF,DOWK) 1058 CALL SOLOUT(MODE,ISTAT(2)+1,NQ,NV,NU,NL,LDG,LDF,LDA,LRDO, 1059 & LIDO,FPROB,Q,V,U,DOWK,IDOWK) 1060 END IF 1061C --- COMPUTE THE ACCELERATION 1062 IF((IACC.EQ.1).OR.(IOUT.EQ.2)) THEN 1063 DO I=1,NU 1064 CONTU(I)= U1(I) 1065 CONTU(I+NU)= D14*UDOT4(I)+D15*UDOT5(I)+D16*UDOT6(I) 1066 & +D17*UDOT7(I)+D18*UDOT2(I) 1067 CONTU(I+NU2)= D24*UDOT4(I)+D25*UDOT5(I)+D26*UDOT6(I) 1068 & +D27*UDOT7(I)+D28*UDOT2(I) 1069 CONTU(I+NU3)= D34*UDOT4(I)+D35*UDOT5(I)+D36*UDOT6(I) 1070 & +D37*UDOT7(I)+D38*UDOT2(I) 1071 CONTU(I+NU4)= D44*UDOT4(I)+D45*UDOT5(I)+D46*UDOT6(I) 1072 & +D47*UDOT7(I)+D48*UDOT2(I) 1073 END DO 1074 CALL FPROB(7,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1075 & IPAR(3),IPAR(4),INDG1(1),INDG1(LDG+1),INDFL(1), 1076 & INDFL(LDF+1),TH,Q,V,U,XL, 1077 & G,GQ1,X0(1),X0(NV+1),GT,FL,QDOT1,UDOT1,B) 1078 IF (IGII.EQ.1) THEN 1079 CALL GIINUM(MODE,NQ,NV,NU,NL,NDIM,MDIM,NMDIM,NM, 1080 & LDG,LDF,LDA,NBLK,NMRC,IPAR(3),IPAR(4),INDG0, 1081 & INDGD,INDFL,TH,UROUND,FPROB,Q,QD,QDOT,V,U,GQ0, 1082 & GD,GT,GTD,FL,AM,X0(NV+1),G) 1083 END IF 1084 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 1085 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG0,INDG0, 1086 & INDFL,INDA,IP,IUMF,XUMF,B,GQ0,GQ0,FL,AVALUE,IER) 1087 IF (IER.NE.0) GOTO 176 1088 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 1089 & NZA,AVALUE,XUMF,B,X0,VP1,XL) 1090 ISTAT(4)=ISTAT(4)+1 1091 ISTAT(6)=ISTAT(6)+1 1092 ISTAT(8)=ISTAT(8)+1 1093 END IF 1094 IF (IACC.EQ.1) THEN 1095 DO I=1,NV 1096 A(I)=VP1(I) 1097 END DO 1098 DO I=1,NL 1099 RLAM(I)=XL(I) 1100 END DO 1101 END IF 1102 IF (IOUT.EQ.2) THEN 1103C -- COMPUTE Q ,V AND U AT S=THETA 1104 S=0.540179418d0 1105 TSH=T+S*H 1106 BS=S*(B10(1)+S*(B10(2)+S*(B10(3)+S*B10(4)))) 1107 DO J=1,NV 1108 QD(J)=CONTQ(J)+H*S*(CONTQ(J+NQ)+S*(CONTQ(J+NQ2)+ 1109 & S*(CONTQ(J+NQ3)+S*CONTQ(J+NQ4)))) 1110 VD(J)=CONTV(J)+H*S*(CONTV(J+NV)+S*(CONTV(J+NV2)+ 1111 & S*(CONTV(J+NV3)+S*CONTV(J+NV4)))) 1112 END DO 1113 DO J=NV+1,NQ 1114 QD(J)=CONTQ(J)+H*S*(CONTQ(J+NQ)+S*(CONTQ(J+NQ2)+ 1115 & S*(CONTQ(J+NQ3)+S*CONTQ(J+NQ4)))) 1116 END DO 1117 DO J=1,NU 1118 UD(J)=CONTU(J)+H*S*(CONTU(J+NU)+S*(CONTU(J+NU2)+ 1119 & S*(CONTU(J+NU3)+S*CONTU(J+NU4)))) 1120 END DO 1121 CALL FPROB(8,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1122 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 1123 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 1124 & G,GQ0,XD(1),XD(NV+1),GT,FL,QDOT3,UDOT3,B) 1125 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1126 & IPAR(3),IPAR(4),INDGD(1),INDGD(LDG+1),INDFL(1), 1127 & INDFL(LDF+1),TSH,QD,VD,UD,XL, 1128 & G,GD,XD(1),XD(NV+1),GTD,FL,QDOT3,UDOT3,B) 1129 FAC = -1.D0/(H*BS) 1130 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDGD, 1131 & GD,FAC,VD,GTD,XD(NV+1)) 1132 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 1133 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG1,INDGD, 1134 & INDFL,INDA,IP,IUMF,XUMF,B,GQ1,GD,FL,AVALUE,IER) 1135 IF (IER.NE.0) GOTO 176 1136 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 1137 & NZA,AVALUE,XUMF,B,XD,XD,XL) 1138 CALL FPROB(0,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1139 & IPAR(3),IPAR(4),INDG0(1),INDG0(LDG+1),INDFL(1), 1140 & INDFL(LDF+1),TCH,Q3,V3,U3,XL, 1141 & G,GQ0,XD(1),XD(NV+1),GT,FL,QDOT3,UDOT3,B) 1142 DO J=1,NU 1143 UD(J)=UD(J)+H*BS*UDOT3(J) 1144 END DO 1145 DO J=1,NV 1146 VD(J)=VD(J)+H*BS*XD(J) 1147 END DO 1148C-- STATISTICS 1149 ISTAT(4)=ISTAT(4)+1 1150 ISTAT(5)=ISTAT(5)+1 1151 ISTAT(6)=ISTAT(6)+1 1152 ISTAT(8)=ISTAT(8)+1 1153C --- PREPARES DOWK 1154 DO I=1,NV 1155 IP2=I+2 1156 DOWK(IP2)=Q1(I) 1157 DOWK(IP2+NQVU)=QDOT1(I) 1158 DOWK(IP2+NQVU2)=QDOT(I) 1159 DOWK(IP2+NQVU3)=QD(I)-Q1(I) 1160 DOWK(IP2+NQVU4)=Q(I)-Q1(I) 1161 IP3=IP2+NQ 1162 DOWK(IP3)=V1(I) 1163 DOWK(IP3+NQVU)=VP1(I) 1164 DOWK(IP3+NQVU2)=VP(I) 1165 DOWK(IP3+NQVU3)=VD(I)-V1(I) 1166 DOWK(IP3+NQVU4)=V(I)-V1(I) 1167 END DO 1168 DO I=NV+1,NQ 1169 IP2=I+2 1170 DOWK(IP2)=Q1(I) 1171 DOWK(IP2+NQVU)=QDOT1(I) 1172 DOWK(IP2+NQVU2)=QDOT(I) 1173 DOWK(IP2+NQVU3)=QD(I)-Q1(I) 1174 DOWK(IP2+NQVU4)=Q(I)-Q1(I) 1175 END DO 1176 NNN=NQ+NV+2 1177 DO I=1,NU 1178 IP4=NNN+I 1179 DOWK(IP4)=U1(I) 1180 DOWK(IP4+NQVU)=UDOT1(I) 1181 DOWK(IP4+NQVU2)=UDOT(I) 1182 DOWK(IP4+NQVU3)=UD(I)-U1(I) 1183 DOWK(IP4+NQVU4)=U(I)-U1(I) 1184 END DO 1185 DOWK(1)=T 1186 DOWK(2)=H 1187 CALL SOLOUT(MODE,ISTAT(2)+1,NQ,NV,NU,NL,LDG,LDF,LDA,LRDO, 1188 & LIDO,FPROB,Q,V,U,DOWK,IDOWK) 1189 END IF 1190 IF(IOUT.EQ.3) THEN 1191 DOWK(1)=T 1192 DOWK(2)=H 1193 CALL SOLOUT(MODE,ISTAT(2)+1,NQ,NV,NU,NL,LDG,LDF,LDA, 1194 & LRDO,LIDO,FPROB,Q,V,U,DOWK,IDOWK) 1195 END IF 1196 IF (IRTRN.LT.0) GOTO 176 1197 IF(ABS(HNEW).GT.HMAX)HNEW=POSNEG*HMAX 1198 IF(REJECT)HNEW=POSNEG*MIN(ABS(HNEW),ABS(H)) 1199 REJECT=.FALSE. 1200 TOLD=T 1201 T=TH 1202 IF (LAST) THEN 1203 IDID=1 1204 RETURN 1205 END IF 1206 DO I=1,NV 1207 Q1(I) = Q(I) 1208 QDOT1(I) = QDOT(I) 1209 V1(I) = V(I) 1210 VP(I) = VP1(I) 1211 END DO 1212 DO I=NV+1,NQ 1213 Q1(I) = Q(I) 1214 QDOT1(I) = QDOT(I) 1215 END DO 1216 DO I=1,NU 1217 U1(I) = U(I) 1218 END DO 1219 DO I=1,LDG 1220 DO J=1,NDIM 1221 GQ(I,J) = GQ0(I,J) 1222 END DO 1223 INDG(I)=INDG0(I) 1224 INDG(LDG+I)=INDG0(LDG+I) 1225 END DO 1226 ELSE 1227C -- STEP IS REJECTED AFTER A PROJECTION 1228 HNEW=0.7D0*H 1229 REJECT=.TRUE. 1230 IF (ISTAT(2).GE.1) ISTAT(3)=ISTAT(3)+1 1231 END IF 1232 ELSE 1233C --- STEP IS REJECTED WITHOUT PROJECTION 1234 HNEW=H/DMIN1(FACC1,FAC11/SAFE) 1235 REJECT=.TRUE. 1236 IF(ISTAT(2).GE.1)ISTAT(3)=ISTAT(3)+1 1237 END IF 1238 H = HNEW 1239 LAST=.FALSE. 1240 GOTO 909 1241C --- FAIL EXIT 1242 176 CONTINUE 1243 WRITE(6,979)T 1244 WRITE(6,*) ' MATRIX IN ASET IS SINGULAR, IER=',IER 1245 IDID=-4 1246 RETURN 1247 177 CONTINUE 1248 WRITE(6,979)T 1249 WRITE(6,*) ' STEP SIZE T0O SMALL, H=',H 1250 IDID=-3 1251 RETURN 1252 178 CONTINUE 1253 WRITE(6,979)T 1254 WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED' 1255 IDID=-2 1256 RETURN 1257C --- EXIT CAUSED BY SOLOUT 1258 179 CONTINUE 1259 WRITE(6,979)T 1260 WRITE(6,*) 'INITIAL PROJECTION: NO CONVERGENCE' 1261 IDID=-5 1262 RETURN 1263 979 FORMAT(' EXIT OF HEM5 AT X=',E18.4) 1264 END 1265C***************************************************************** 1266 SUBROUTINE GIINUM(MODE,NQ,NV,NU,NL,NDIM,MDIM,NMDIM,NM, 1267 & LDG,LDF,LDA,NBLK,NMRC,NPGP,NPFL,INDG,INDGD,INDFL,T, 1268 & UROUND,FPROB,Q,QD,QDOT,V,U,GQ,GQD,GT,GTD, 1269 & FL,B,GII,RES) 1270 IMPLICIT REAL*8(A-H,O-Z) 1271 EXTERNAL FPROB 1272 DIMENSION INDG(2*LDG),INDGD(2*LDG),Q(NQ),QD(NQ) 1273 DIMENSION QDOT(NQ),V(NV),U(NU),GQ(LDG,NDIM) 1274 DIMENSION GQD(LDG,NDIM),GT(NL),GTD(NL),FL(LDF,MDIM) 1275 DIMENSION B(LDA,NMDIM),GII(NL),RES(NL),INDFL(2*LDF) 1276 DA=DSQRT(UROUND) 1277 CALL FPROB(3,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1278 & NPGP,NPFL,INDGD(1),INDGD(LDG+1),INDFL(1), 1279 & INDFL(LDF+1),T,Q,V,U,RES, 1280 & RES,GQD,V,RES,GT,FL,QDOT,U,B) 1281 DO I=1,NQ 1282 QD(I)=Q(I)+DA*QDOT(I) 1283 END DO 1284 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1285 & NPGP,NPFL,INDGD(1),INDGD(LDG+1),INDFL(1), 1286 & INDFL(LDF+1),T,QD,V,U,RES, 1287 & RES,GQD,V,RES,GTD,FL,QDOT,U,B) 1288 CALL GPM2(MODE,NV,NL,NDIM,MDIM,LDG,INDG,INDGD, 1289 & UROUND,GQ,GQD,V,RES) 1290 DO I=1,NL 1291 GII(I)=RES(I)+(GTD(I)-GT(I))/DA 1292 END DO 1293 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1294 & NPGP,NPFL,INDGD(1),INDGD(LDG+1),INDFL(1), 1295 & INDFL(LDF+1),T+DA,Q,V,U, 1296 & RES,RES,GQD,V,RES,GTD,FL,QDOT,U,B) 1297 CALL GPM2(MODE,NV,NL,NDIM,MDIM,LDG,INDG,INDGD, 1298 & UROUND,GQ,GQD,V,RES) 1299 DO I=1,NL 1300 GII(I)=-(GII(I)+RES(I)+(GTD(I)-GT(I))/DA) 1301 END DO 1302 RETURN 1303 END 1304C********************************************************** 1305 SUBROUTINE GPM2(MODE,N,M,NDIM,MDIM,LDG,INDG,INDGD, 1306 & UROUND,GQ,GQD,V,RES) 1307 IMPLICIT REAL*8 (A-H,O-Z) 1308 DIMENSION INDG(2*LDG),INDGD(2*LDG),GQ(LDG,NDIM) 1309 DIMENSION GQD(LDG,NDIM),V(N),RES(M) 1310 DA=DSQRT(UROUND) 1311 IF ((MODE.GE.0).AND.(MODE.LE.3)) THEN 1312 GOTO 10 1313 ELSE IF ((MODE.EQ.4).OR.(MODE.EQ.5)) THEN 1314 GOTO 20 1315 ELSE 1316 STOP ' GPM2: INVALID MODE' 1317 ENDIF 1318 10 CONTINUE 1319 DO I=1,M 1320 SUM=0.D0 1321 DO J=1,N 1322 SUM=SUM+(GQD(I,J)-GQ(I,J))*V(J)/DA 1323 END DO 1324 RES(I)=SUM 1325 END DO 1326 RETURN 1327 20 CONTINUE 1328 DO I=1,M 1329 RES(I)=0.D0 1330 END DO 1331 DO K=1,LDG 1332 I=INDG(K) 1333 J=INDG(LDG+K) 1334 RES(I)=RES(I)+(GQD(K,1)-GQ(K,1))*V(J)/DA 1335 END DO 1336 RETURN 1337 END 1338C******************************************************************* 1339C* SUBROUTINE ASET CONSTRUCTS AND DECOMPOSES THE MATRIX 1340C* WITH : 1341C* | AM G0^t | 1342C* B = | | 1343C* | G1 0 | 1344C* 1345C******************************************************************* 1346 SUBROUTINE ASET(MODE,N,M,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA, 1347 & NZA,LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG0,INDG1, 1348 & INDFL,INDA,IP,IUMF,XUMF,B,G0,G1,FL,AVALUE,IER) 1349 IMPLICIT LOGICAL (A-Z) 1350 INTEGER MODE,N,M,NM,NDIM,MDIM,NMDIM,LDG,LDF,LDA,IP(NM), 1351 & LIPAR,LIUMF,LXUMF,IPAR(LIPAR),IUMF(LIUMF),INDG0(2*LDG), 1352 & INDG1(2*LDG),INDFL(2*LDF),NZA,INDA(2*NZA),IIRN,IICN,IIK, 1353 & IREFAC,I,J,LISTA,ISTAT(LISTA),IER,LICN,IIW,ERROR,IICNM1 1354 DOUBLE PRECISION B(LDA,NMDIM),G0(LDG,NDIM),G1(LDG,NDIM), 1355 & AVALUE(2*NZA),XUMF(LXUMF),FL(LDF,MDIM) 1356C For a description of the internal parameters of the package 1357C see: Harwell Subroutine Library Specification 1358C 1359C 1360 INTEGER LP, MP, IRNCP, ICNCP, MINIRN, MINICN, IRANK, IDISP(2) 1361 LOGICAL QBLOCK, QGROW, QABRT1, QABRT2 1362 DOUBLE PRECISION EPSMA, RMIN, RESID,THRSH1,THRSH2 1363C 1364 COMMON /MA28ED/ LP, MP, QBLOCK, QGROW 1365 COMMON /MA28FD/ EPSMA, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN, 1366 $ IRANK, QABRT1, QABRT2 1367 COMMON /MA28GD/ IDISP 1368C 1369C Description see Harwell Subroutine 1370C Library Specification 1371C 1372 SAVE /MA28ED/, /MA28FD/, /MA28GD/ 1373C 1374 IER=0 1375 IF (MODE.EQ.0) THEN 1376 GOTO 10 1377 ELSE IF (MODE.EQ.1) THEN 1378 GOTO 15 1379 ELSE IF (MODE.EQ.2) THEN 1380 GOTO 17 1381 ELSE IF (MODE.EQ.3) THEN 1382 GOTO 18 1383 ELSE IF ((MODE.EQ.4).OR.(MODE.EQ.5)) THEN 1384 GOTO 20 1385 ELSE 1386 STOP ' ASET: INVALID MODE' 1387 END IF 1388 10 CONTINUE 1389 DO I=1,N 1390 DO J=1,M 1391 B(I,N+J) = G0(J,I) 1392 END DO 1393 END DO 1394 DO I=1,M 1395 DO J=1,N 1396 B(N+I,J) = G1(I,J) 1397 END DO 1398 DO J=1,M 1399 B(N+I,N+J) = 0.D0 1400 END DO 1401 END DO 1402 CALL DEC(NM,NM,B,IP,IER) 1403 ISTAT(7)=ISTAT(7)+1 1404 RETURN 1405 15 CONTINUE 1406 DO I=1,N 1407 DO J=1,M 1408 B(I,N+J) = FL(I,J) 1409 END DO 1410 END DO 1411 DO I=1,M 1412 DO J=1,N 1413 B(N+I,J) = G1(I,J) 1414 END DO 1415 DO J=1,M 1416 B(N+I,N+J) = 0.D0 1417 END DO 1418 END DO 1419 CALL DEC(NM,NM,B,IP,IER) 1420 ISTAT(7)=ISTAT(7)+1 1421 RETURN 1422 17 CONTINUE 1423 DO I=1,N 1424 DO J=1,M 1425 B(I,N+J) = G0(J,I) 1426 END DO 1427 END DO 1428 DO I=1,M 1429 DO J=1,N 1430 B(N+I,J) = G1(I,J) 1431 END DO 1432 DO J=1,M 1433 B(N+I,N+J) = 0.D0 1434 END DO 1435 END DO 1436 CALL DGETRF(NM,NM,B,NM,IP,IER) 1437 ISTAT(7)=ISTAT(7)+1 1438 RETURN 1439 18 CONTINUE 1440 DO I=1,N 1441 DO J=1,M 1442 B(I,N+J) = FL(I,J) 1443 END DO 1444 END DO 1445 DO I=1,M 1446 DO J=1,N 1447 B(N+I,J) = G1(I,J) 1448 END DO 1449 DO J=1,M 1450 B(N+I,N+J) = 0.D0 1451 END DO 1452 END DO 1453 CALL DGETRF(NM,NM,B,NM,IP,IER) 1454 ISTAT(7)=ISTAT(7)+1 1455 RETURN 1456 20 CONTINUE 1457C 1458C 1459C 1460C RAPPEL INDG(1->LDG)=LIGNE, INDG(LDG+1->2*LDG)=COLONNE 1461C 1462C 1463 THRSH1 = 1.D-2 1464 THRSH2 = 1.D-6 1465 LICN=2*NZA 1466 IIRN=1 1467 IICN=IIRN+LICN 1468 IIK=IICN+LICN 1469 IIW=IIK+5*NM 1470 IICNM1=IICN-1 1471C ?? IFDEC = 0 1472 EPSMA =DMIN1 (1.D0, 10.D0*THRSH2) 1473C 1474 IF (IPAR(9).EQ.0) THEN 1475 IREFAC=IPAR(3)+IPAR(4) 1476 ELSE 1477 IPAR(9)=0 1478 IREFAC=1 1479 ENDIF 1480 100 CALL ACOPY(MODE,IPAR(1),IPAR(2),LDG,LDF,LDA,NZA,N,NMDIM, 1481 & INDG0,INDG1,INDFL,INDA,AVALUE,B,G0,G1,FL) 1482 IF (IREFAC.GT.0) THEN 1483C WRITE(6,*)'MA28AD' 1484 DO I=1,NZA 1485 IUMF(I)=INDA(I) 1486 IUMF(IICNM1+I)=INDA(NZA+I) 1487 END DO 1488 CALL MA28AD(NM,NZA,AVALUE,LICN,IUMF(IIRN),LICN,IUMF(IICN), 1489 & THRSH1,IUMF(IIK),IUMF(IIW),XUMF,ERROR) 1490 ISTAT(7)=ISTAT(7)+1 1491 ELSE 1492 CALL MA28BD(NM,NZA,AVALUE,LICN,INDA(1),INDA(NZA+1),IUMF(IICN), 1493 & IUMF(IIK),IUMF(IIW),XUMF,ERROR ) 1494 ISTAT(7)=ISTAT(7)+1 1495 IF(ERROR.LT.0) THEN 1496 WRITE(6,*)'ERROR IN UMDREFAC' 1497 IREFAC=1 1498 GOTO 100 1499 END IF 1500 IF(RMIN.LT.THRSH2) THEN 1501 IREFAC=1 1502 WRITE(6,*)'RMIN TOO SMALL' 1503 GOTO 100 1504 END IF 1505 END IF 1506 RETURN 1507 END 1508C 1509C******************************************************************* 1510C* SUBROUTINE GPMULT COMPUTES FAC*(AG*VECT+GT) = RES 1511C****************************************************************** 1512 SUBROUTINE GPMULT(MODE,NDIM,N,M,LDG,INDG, 1513 & AG,FAC,VECT,GT,RES) 1514 IMPLICIT LOGICAL (A-Z) 1515 INTEGER MODE,NDIM,N,M,LDG,INDG(2*LDG), 1516 & I,J,K 1517 DOUBLE PRECISION AG(LDG,NDIM),VECT(N),GT(M), 1518 & RES(M),FAC,SUM 1519 IF ((MODE.GE.0).AND.(MODE.LE.3)) THEN 1520 GOTO 10 1521 ELSE IF ((MODE.EQ.4).OR.(MODE.EQ.5)) THEN 1522 GOTO 20 1523 ELSE 1524 STOP ' GPMULT: INVALID MODE' 1525 END IF 1526 10 CONTINUE 1527 DO I=1,M 1528 SUM=0.D0 1529 DO J=1,N 1530 SUM=SUM+AG(I,J)*VECT(J) 1531 END DO 1532 RES(I)=FAC*(SUM+GT(I)) 1533 END DO 1534 RETURN 1535 20 CONTINUE 1536 DO I=1,M 1537 RES(I)=FAC*GT(I) 1538 END DO 1539 DO K=1,LDG 1540 I=INDG(K) 1541 J=INDG(LDG+K) 1542 RES(I)=RES(I)+FAC*AG(K,1)*VECT(J) 1543 END DO 1544 RETURN 1545 END 1546C******************************************************************* 1547C* SUBROUTINE AMULT COMPUTES FAC*AM*VECT = RES 1548C****************************************************************** 1549 SUBROUTINE AMULT(MODE,N,NMDIM,LDA,NBLK,NMRC,FAC,AM,VECT,RES) 1550 IMPLICIT LOGICAL (A-Z) 1551 INTEGER MODE,N,NMDIM,NBLK,NMRC,I,J,K,KN,LDA 1552 DOUBLE PRECISION AM(LDA,NMDIM),VECT(N),FAC,SUM,RES(N) 1553 IF ((MODE.GE.0).AND.(MODE.LE.3)) THEN 1554 GOTO 10 1555 ELSE IF ((MODE.EQ.4).OR.(MODE.EQ.5)) THEN 1556 GOTO 20 1557 ELSE 1558 STOP ' AMULT: INVALID MODE' 1559 END IF 1560 10 CONTINUE 1561 DO I=1,N 1562 SUM=0.D0 1563 DO J=1,N 1564 SUM=SUM+AM(I,J)*VECT(J) 1565 END DO 1566 RES(I)=FAC*SUM 1567 END DO 1568 RETURN 1569 20 CONTINUE 1570 DO K=1,NBLK 1571 KN=(K-1)*NMRC 1572 DO I=1,NMRC 1573 SUM=0.D0 1574 DO J=1,NMRC 1575 SUM=SUM+AM(KN+I,J)*VECT(KN+J) 1576 END DO 1577 RES(KN+I)=FAC*SUM 1578 END DO 1579 END DO 1580 RETURN 1581 END 1582C******************************************************************* 1583C* SUBROUTINE ACOPY 1584C******************************************************************* 1585 SUBROUTINE ACOPY(MODE,NMRC,NBLK,LDG,LDF,LDA,NZA,N,NMDIM, 1586 & INDG1,INDG2,INDFL,INDA,AVALUE,AM,G1,G2,FL) 1587 IMPLICIT LOGICAL (A-Z) 1588 INTEGER MODE,NMRC,NBLK,LDG,LDF,LDA,NZA,N,NMDIM,NMNB,I,J,K, 1589 & L,KN,INDG1(2*LDG),INDG2(2*LDG),INDA(2*NZA),INDFL(2*LDF) 1590 DOUBLE PRECISION AVALUE(2*NZA),G1(LDG),G2(LDG), 1591 & AM(LDA,NMDIM),FL(LDF) 1592 IF (MODE.EQ.4) THEN 1593 GOTO 10 1594 ELSE IF (MODE.EQ.5) THEN 1595 GOTO 20 1596 ELSE 1597 STOP ' ACOPY: INVALID MODE' 1598 END IF 1599 10 CONTINUE 1600 NMNB=NMRC*NBLK 1601 L=0 1602 DO K=1,NBLK 1603 KN=(K-1)*NMRC 1604 DO i=1,NMRC 1605 DO j=1,NMRC 1606 L=L+1 1607 AVALUE(L)=AM(I+KN,J) 1608 INDA(L)=KN+I 1609 INDA(NZA+L)=KN+J 1610 END DO 1611 END DO 1612 END DO 1613 DO I=1,LDG 1614 L=L+1 1615 AVALUE(L)=G1(I) 1616 INDA(L)=INDG1(LDG+I) 1617 INDA(NZA+L)=NMNB+INDG1(I) 1618 L=L+1 1619 AVALUE(L)=G2(I) 1620 INDA(L)=NMNB+INDG2(I) 1621 INDA(NZA+L)=INDG2(LDG+I) 1622 END DO 1623 RETURN 1624 20 CONTINUE 1625 NMNB=NMRC*NBLK 1626 L=0 1627 DO K=1,NBLK 1628 KN=(K-1)*NMRC 1629 DO i=1,NMRC 1630 DO j=1,NMRC 1631 L=L+1 1632 AVALUE(L)=AM(I+KN,J) 1633 INDA(L)=KN+I 1634 INDA(NZA+L)=KN+J 1635 END DO 1636 END DO 1637 END DO 1638 DO I=1,LDF 1639 L=L+1 1640 AVALUE(L)=FL(I) 1641 INDA(L)=INDFL(I) 1642 INDA(NZA+L)=NMNB+INDFL(LDF+I) 1643 END DO 1644 DO I=1,LDG 1645 L=L+1 1646 AVALUE(L)=G2(I) 1647 INDA(L)=NMNB+INDG2(I) 1648 INDA(NZA+L)=INDG2(LDG+I) 1649 END DO 1650 RETURN 1651 END 1652C******************************************************************* 1653C* SUBROUTINE ASOL SOLVES THE SYSTEM B*X0=R 1654C* WITH THE MATRIX B DECOMPOSED IN ASET 1655C* AND MOVES THE RESULT X0(I=1,..,N) -> RES(I) 1656C******************************************************************* 1657 SUBROUTINE ASOL(MODE,NM,N,M,LIPAR,LIUMF,LXUMF,IPAR, 1658 & IUMF,IP,NZA,AVALUE,XUMF,B,X0,RES1,RES2) 1659 IMPLICIT LOGICAL (A-Z) 1660 INTEGER MODE,N,NM,NZA,LIPAR,LIUMF,LXUMF,IPAR(LIPAR), 1661 & IUMF(LIUMF),IP(NM),IIW,I,M,LICN,IICN,IIK,MTYPE 1662 DOUBLE PRECISION B(NM,NM),X0(NM),RES1(N),RES2(M),XUMF(LXUMF), 1663 & AVALUE(2*NZA) 1664C -- umfpack variables: 1665C INTEGER MTYPE 1666C 1667 IF ((MODE.EQ.0).OR.(MODE.EQ.1)) THEN 1668 GOTO 10 1669 ELSE IF ((MODE.EQ.2).OR.(MODE.EQ.3)) THEN 1670 GOTO 15 1671 ELSE IF ((MODE.EQ.4).OR.(MODE.EQ.5)) THEN 1672 GOTO 20 1673 ELSE 1674 STOP ' ASOL: INVALID MODE' 1675 END IF 1676 10 CONTINUE 1677 CALL SOL(NM,NM,B,X0,IP) 1678 GOTO 100 1679 15 CONTINUE 1680 CALL DGETRS('N',NM,1,B,NM,IP,X0,NM,IER) 1681 GOTO 100 1682 20 CONTINUE 1683 LICN=2*NZA 1684 IICN=1+LICN 1685 IIK=IICN+LICN 1686 IIW=IIK+5*NM 1687 MTYPE=1 1688 CALL MA28CD(NM,AVALUE,LICN,IUMF(IICN),IUMF(IIK),X0, 1689 & XUMF,MTYPE) 1690 100 CONTINUE 1691 DO I=1,N 1692 RES1(I) = X0(I) 1693 END DO 1694 DO I=1,M 1695 RES2(I) = X0(N+I) 1696 END DO 1697 RETURN 1698 END 1699C*********************************************************** 1700C* SUBROUTINE PRO PROJECTION ON THE MANIFOLD DEFINED 1701C* BY {g(q) = 0} 1702C* 1703C* SOLVES BY NEWTON ITERATIONS A NON LINEAR SYSTEM 1704C*********************************************************** 1705 SUBROUTINE APROJ(MODE,NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM, 1706 * NBLK,NMRC,LDG,LDF,LDA,NZA,LIPAR,LIUMF,LXUMF,LISTA, 1707 * ISTAT,IPAR,IUMF,INDA,INDG,INDFL,FL,XUMF,AVALUE,T, 1708 * FPROB,Q,Q1,Q2,QDOT,V,V1,V2,U,XL,UDOT, G,GT,GQ,AM, 1709 * B,X0,IP,ATOL,RTOL,ITOL,ACCEPT) 1710 IMPLICIT REAL*8 (A-H,O-Z) 1711 DIMENSION Q(NQ),Q1(NQ),Q2(NQ),QDOT(NQ),V(NV),V1(NV) 1712 DIMENSION AM(LDA,NMDIM),G(NL),B(LDA,NMDIM),V2(NV),U(NU) 1713 DIMENSION X0(NM),GQ(LDG,NDIM),GT(NL),UDOT(NU),XL(NL) 1714 DIMENSION IP(NM),ATOL(1),RTOL(1),INDG(2*LDG),ISTAT(LISTA) 1715 DIMENSION INDFL(2*LDF),FL(LDF,MDIM),IPAR(LIPAR),IUMF(LIUMF) 1716 DIMENSION INDA(2*NZA),XUMF(LXUMF),AVALUE(2*NZA) 1717 EXTERNAL FPROB 1718 LOGICAL ACCEPT 1719 CALL FPROB(9,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1720 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 1721 & INDFL(LDF+1),T,Q,V,U,XL, 1722 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,AM) 1723C --- INITIAL PREPARATION ------------------------------------ 1724 ERROLD=1.D20 1725 FAC = -1.D0 1726 DO I=1,NV 1727 Q2(I)=Q(I) 1728 V2(I)=0.D0 1729 END DO 1730 DO I=NV+1,NQ 1731 Q2(I)=Q(I) 1732 END DO 1733 K=0 1734 ITMAX = 5 1735 EPS = 1.d-2 1736C --- COMPUTE DEFECT IN G -------------------------------------- 1737 CALL FPROB(4,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1738 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 1739 & INDFL(LDF+1),T,Q,V,U,XL, 1740 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 1741 DO I=1,NV 1742 X0(I)=0.D0 1743 END DO 1744 DO I=1,NL 1745 X0(NV+I)=-G(I) 1746 END DO 1747 106 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 1748 & NZA,AVALUE,XUMF,B,X0,X0,XL) 1749 ISTAT(8)=ISTAT(8)+1 1750 DO I=1,NV 1751 V2(I)=V2(I)+X0(I) 1752 END DO 1753 CALL FPROB(2,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1754 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 1755 & INDFL(LDF+1),T,Q,X0,U,XL, 1756 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 1757 DO I=1,NQ 1758 Q2(I)=Q2(I)+QDOT(I) 1759 END DO 1760 ERR=0.d0 1761 DO I=1,NQ 1762 IF (ITOL.EQ.0) THEN 1763 SK=ATOL(1)+RTOL(1)*MAX(ABS(Q(I)),ABS(Q1(I))) 1764 ELSE 1765 SK=ATOL(I)+RTOL(I)*MAX(ABS(Q(I)),ABS(Q1(I))) 1766 ENDIF 1767 ERR = ERR+(QDOT(I)/SK)**2 1768 END DO 1769 ERR=SQRT(ERR/NQ) 1770C 1771C --- TEST FOR CONVERGENCE -------------------------------------- 1772C 1773 IF(ERR.LT.EPS) GOTO 200 1774 IF(ERR.GE.ERROLD*2.D0 .OR. K.GE.ITMAX) THEN 1775 WRITE(6,*)'NO CONVERGENCE IN PROJECTION' 1776C --- NO CONVERGENCE STEP IS REJECTED -------------------------------------------- 1777 ACCEPT = .FALSE. 1778 RETURN 1779 ENDIF 1780C 1781C --- NEXT ITERATION -------------------------------------------- 1782C 1783 ERROLD=ERR 1784 K=K+1 1785 CALL FPROB(4,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1786 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 1787 & INDFL(LDF+1),T,Q2,V,U,XL, 1788 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 1789 CALL AMULT(MODE,NV,NMDIM,LDA,NBLK,NMRC,FAC,AM,V2,X0) 1790 DO I=1,NL 1791 X0(NV+I)=-G(I) 1792 END DO 1793 GOTO 106 1794 200 CONTINUE 1795C 1796c --- PROJECTION OF VELOCITY ----------------------------------- 1797c 1798cc 1799cc 1800 DO I=1,NQ 1801 Q(I)=Q2(I) 1802 END DO 1803 CALL FPROB(11,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1804 & IPAR(3),IPAR(4),INDG(1),INDG(LDG+1),INDFL(1), 1805 & INDFL(LDF+1),T,Q,V,U,XL, 1806 & G,GQ,X0(1),X0(NV+1),GT,FL,QDOT,UDOT,B) 1807 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA, 1808 & NZA,LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG,INDG, 1809 & INDFL,INDA,IP,IUMF,XUMF,B,GQ,GQ,FL,AVALUE,IER) 1810 IF (IER.NE.0) THEN 1811 WRITE(6,*)'MATRIX SINGULAR IN APROJ' 1812 ACCEPT=.FALSE. 1813 RETURN 1814 END IF 1815 DO I=1,NV 1816 X0(I)=0.D0 1817 END DO 1818 FAC=-1.D0 1819 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDG, 1820 & GQ,FAC,V,GT,X0(NV+1)) 1821 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 1822 & NZA,AVALUE,XUMF,B,X0,X0,XL) 1823 ISTAT(8)=ISTAT(8)+1 1824 DO I=1,NV 1825 V(I)=V(I)+X0(I) 1826 END DO 1827C -- ERROR ESTIMATION 1828 ERR=0.d0 1829 DO I=1,NV 1830 IF (ITOL.EQ.0) THEN 1831 SK=ATOL(1)+RTOL(1)*MAX(ABS(V(I)),ABS(V1(I))) 1832 ELSE 1833 SK=ATOL(I+NQ)+RTOL(I+NQ)*MAX(ABS(V(I)),ABS(V1(I))) 1834 ENDIF 1835 ERR = ERR+(X0(I)/SK)**2 1836 END DO 1837 ERR=SQRT(ERR/NV) 1838 IF(ERR.GT.EPS)ACCEPT=.FALSE. 1839 IF(ERR.GT.EPS)write(6,*)'reject after proj',err 1840 RETURN 1841 END 1842C 1843C*************************************************************** 1844C 1845C DENSE OUTPUT 1846C 1847C*************************************************************** 1848 FUNCTION CTQ4(MODE,I,NQ,NV,NU,NL,LDG,LDF,LDA,LRDO, 1849 & LIDO,IDOWK,X,DOWK,FPROB) 1850C ---------------------------------------------------------- 1851C THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONECTION 1852C WITH THE OUTPUT-SUBROUTINE FOR HEM5. IT PROVIDES AN 1853C APPROXIMATION TO THE I-TH COMPONENT OF THE Q-SOLUTION AT X. 1854C ---------------------------------------------------------- 1855 IMPLICIT REAL*8 (A-H,O-Z) 1856 DIMENSION DOWK(LRDO),IDOWK(LIDO) 1857 EXTERNAL FPROB 1858 XOLD=DOWK(1) 1859 H=DOWK(2) 1860 N=NQ 1861 NN=7+N 1862 NN2=NN+N 1863 NN3=NN2+N 1864 NN4=NN3+N 1865 S=(X-XOLD)/H 1866 CTQ4=DOWK(I+7)+H*S*(DOWK(I+NN)+S*(DOWK(I+NN2)+ 1867 & S*(DOWK(I+NN3)+S*DOWK(I+NN4)))) 1868 RETURN 1869 END 1870C 1871 FUNCTION CT4(MODE,I,NQ,NV,NU,NL,LDG,LDF,LDA,LRDO, 1872 & LIDO,IDOWK,X,DOWK,FPROB) 1873C ---------------------------------------------------------- 1874C THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONECTION 1875C WITH THE OUTPUT-SUBROUTINE FOR HEM5. IT PROVIDES AN 1876C APPROXIMATION TO THE I-TH COMPONENT OF THE (Q,V)-SOLUTION AT X. 1877C SET (1<=I<=N FOR Q(I) AND N<I<=2N FOR V(I-N)) 1878C ---------------------------------------------------------- 1879 IMPLICIT REAL*8 (A-H,O-Z) 1880 DIMENSION DOWK(LRDO),IDOWK(LIDO) 1881 EXTERNAL FPROB 1882 NM=NV+NL 1883 LIPAR = IDOWK(1) 1884 IPA = NM+2 1885 NMRC=IDOWK(IPA) 1886 NBLK=IDOWK(IPA+1) 1887 NZA = NBLK*NMRC**2+LDF+LDG 1888 IF (MODE.LE.3) THEN 1889 NDIM=NV 1890 MDIM=NL 1891 NMDIM=NM 1892 ELSE 1893 NDIM=1 1894 MDIM=1 1895 NMDIM=NMRC 1896 END IF 1897 ICQ=8 1898 ICV=ICQ+5*NQ 1899 ICU=ICV+5*NV 1900 IQD=ICU+5*NU 1901 IVD=IQD+NQ 1902 IUD=IVD+NV 1903 IGQ1=IUD+NU 1904 IGD=IGQ1+LDG*NDIM 1905 IB=IGD+LDG*NDIM 1906 IXD=IB+LDA*NMDIM 1907 IGTD=IXD+NM 1908 IFL = IGTD+NL 1909 IAV = IFL+LDF*MDIM 1910 IXU = IAV+2*NZA 1911 IST = IPA+LIPAR 1912 ING1 = IST+9 1913 INGD = ING1+2*LDG 1914 INGF = INGD+2*LDG 1915 INA = INGF+2*LDF 1916 INU = INA+2*NZA 1917 LIUMF = 13*NM+4*NZA 1918 LXUMF = NM 1919 CALL CMAT4(MODE,I,NQ,NV,NU,NL,NDIM,MDIM,NMDIM,LDG,LDF,LDA, 1920 & NZA,LIUMF,LXUMF,LIPAR,IDOWK(2),IDOWK(IPA),IDOWK(IST), 1921 & IDOWK(ING1),IDOWK(INGD),IDOWK(INGF),IDOWK(INA),IDOWK(INU), 1922 & FPROB,X,DOWK(1),DOWK(2),DOWK(3),DOWK(4),DOWK(ICQ), 1923 & DOWK(ICV),DOWK(ICU),DOWK(IQD),DOWK(IVD),DOWK(IUD), 1924 & DOWK(IGQ1),DOWK(IGD),DOWK(IB),DOWK(IXD),DOWK(IGTD), 1925 & DOWK(IFL),DOWK(IAV),DOWK(IXU),RES) 1926 CT4=RES 1927 RETURN 1928 END 1929C 1930 SUBROUTINE CMAT4(MODE,I,NQ,NV,NU,NL,NDIM,MDIM,NMDIM,LDG, 1931 & LDF,LDA,NZA,LIUMF,LXUMF,LIPAR,IP,IPAR,ISTAT,INDG1,INDGD, 1932 & INDFL,INDA,IUMF,FPROB,X,XOLD,H,TCH,COEFF,CONTQ,CONTV, 1933 & CONTU,QD,VD,UD,GQ1,GD,B,XD,GTD,FL,AVALUE,XUMF,RES) 1934 IMPLICIT REAL*8 (A-H,O-Z) 1935 DIMENSION B(LDA,NMDIM),CONTQ(5*NQ),CONTV(5*NV),COEFF(4) 1936 DIMENSION GD(LDG,NDIM),GTD(NL),XD(NV+NL),QD(NQ),VD(NV) 1937 DIMENSION IP(NV+NL),ISTAT(9),GQ1(LDG,NDIM) 1938 DIMENSION FL(LDF,MDIM),IUMF(LIUMF),XUMF(LXUMF) 1939 DIMENSION INDG1(2*LDG),INDGD(2*LDG),INDFL(2*LDF) 1940 DIMENSION INDA(2*NZA),AVALUE(2*NZA),IPAR(LIPAR) 1941 DIMENSION UD(NU),CONTU(5*NU) 1942 EXTERNAL FPROB 1943 EPS = 1.D-12 1944 S=(X-XOLD)/H 1945 LISTA=9 1946 NQ2=NQ*2 1947 NQ3=NQ*3 1948 NQ4=NQ*4 1949 NV2=NV*2 1950 NV3=NV*3 1951 NV4=NV*4 1952 NU2=NU*2 1953 NU3=NU*3 1954 NU4=NU*4 1955 NM=NV+NL 1956 NMRC = IPAR(1) 1957 NBLK =IPAR(2) 1958 B10=S*(COEFF(1)+S*(COEFF(2)+S*(COEFF(3)+S*COEFF(4)))) 1959C -- CDUM 1960 IF (S.LT.EPS) GOTO 10 1961c -- qd et vd contiennent q3 et v3 1962 CALL FPROB(8,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1963 & IPAR(3),IPAR(4),INDGD(1),INDGD(LDG+1),INDFL(1), 1964 & INDFL(LDF+1),TCH,QD,VD,UD,XL, 1965 & G,GD,XD(1),XD(NV+1),GTD,FL,QD,UD,B) 1966 DO J=1,NV 1967 QD(J)=CONTQ(J)+H*S*(CONTQ(J+NQ)+S*(CONTQ(J+NQ2)+ 1968 & S*(CONTQ(J+NQ3)+S*CONTQ(J+NQ4)))) 1969 VD(J)=CONTV(J)+H*S*(CONTV(J+NV)+S*(CONTV(J+NV2)+ 1970 & S*(CONTV(J+NV3)+S*CONTV(J+NV4)))) 1971 XD(J)=H*B10*XD(J) 1972 END DO 1973 DO J=NV+1,NQ 1974 QD(J)=CONTQ(J)+H*S*(CONTQ(J+NQ)+S*(CONTQ(J+NQ2)+ 1975 & S*(CONTQ(J+NQ3)+S*CONTQ(J+NQ4)))) 1976 END DO 1977 DO J=1,NU 1978 UD(J)=CONTU(J)+H*S*(CONTU(J+NU)+S*(CONTU(J+NU2)+ 1979 & S*(CONTU(J+NU3)+S*CONTU(J+NU4)))) 1980 END DO 1981 CALL FPROB(6,NQ,NV,NU,NL,LDG,LDF,LDA,NBLK,NMRC, 1982 & IPAR(3),IPAR(4),INDGD(1),INDGD(LDG+1),INDFL(1), 1983 & INDFL(LDF+1),X,QD,VD,UD,XL, 1984 & G,GD,XD(1),XD(NV+1),GTD,FL,QD,UD,B) 1985 FAC = -1.D0 1986 CALL GPMULT(MODE,NDIM,NV,NL,LDG,INDGD, 1987 & GD,FAC,VD,GTD,XD(NV+1)) 1988 CALL ASET(MODE,NV,NL,NM,NMDIM,NDIM,MDIM,LDG,LDF,LDA,NZA, 1989 & LIPAR,LIUMF,LXUMF,LISTA,ISTAT,IPAR,INDG1,INDGD, 1990 & INDFL,INDA,IP,IUMF,XUMF,B,GQ1,GD,FL,AVALUE,IER) 1991 IF (IER.NE.0) THEN 1992 WRITE(6,*)'MATRIX SINGULAR IN CMAT4' 1993 RETURN 1994 END IF 1995 CALL ASOL(MODE,NM,NV,NL,LIPAR,LIUMF,LXUMF,IPAR,IUMF,IP, 1996 & NZA,AVALUE,XUMF,B,XD,XD,XL) 1997 DO J=1,NV 1998 VD(J)=VD(J)+XD(J) 1999 END DO 2000 10 CONTINUE 2001 IF (I.LE.NQ) THEN 2002 RES=QD(I) 2003 ELSE 2004 IF (I.LE.NQ+NV) THEN 2005 RES=VD(I-NQ) 2006 ELSE 2007 RES=UD(I-NQ-NV) 2008 END IF 2009 END IF 2010C -- STATISTICS 2011 ISTAT(4)=ISTAT(4)+1 2012 ISTAT(5)=ISTAT(5)+1 2013 ISTAT(6)=ISTAT(6)+1 2014 ISTAT(8)=ISTAT(8)+1 2015 RETURN 2016 END 2017C****************************************************************** 2018C --- FUNCTION FOR DENSE OUTPUT -- IT IS BASED ON HERMITE 2019C INTERPOLATION 2020 FUNCTION POL(MODE,I,NQ,NV,NU,NL,LDG,LDF,LDA,LRDO, 2021 & LIDO,IDOWK,X,DOWK,FPROB) 2022 IMPLICIT REAL*8 (A-H,O-Z) 2023 DIMENSION DOWK(LRDO),IDOWK(LIDO) 2024 XOLD = DOWK(1) 2025 H= DOWK(2) 2026 S=(X-XOLD)/H 2027 NQVU=NQ+NV+NU 2028 NQVU2=2*NQVU 2029 NQVU3=3*NQVU 2030 NQVU4=4*NQVU 2031 CALL COPOL(S,CP1,CP2,CP3,CP4) 2032 IP2=I+2 2033 POL=DOWK(IP2)+H*CP1*DOWK(IP2+NQVU)+H*CP2*DOWK(IP2+NQVU2)+ 2034 & CP3*DOWK(IP2+NQVU3)+CP4*DOWK(IP2+NQVU4) 2035 RETURN 2036 END 2037C****************************************************************** 2038C --- COEFFICIENTS FOR HERMITE POLYNOMIAL 2039 SUBROUTINE COPOL(X,CP1,CP2,CP3,CP4) 2040 IMPLICIT REAL*8 (A-H,O-Z) 2041 T=0.540179418D0 2042 XM1=(X-1.D0) 2043 TMX=(T-X) 2044 FAC=X*XM1 2045 DEN1=T 2046 DEN2 = 0.4598205820000000D0 2047 DEN3 = 6.1695413425555630D-02 2048 DEN4 = 0.2114349676308187D0 2049 P1 = -1.285336261107544D0 2050 P2 = 3.416412392738363D0 2051 P3 = -1.919641164000000D0 2052 CP1 =FAC*XM1*TMX/DEN1 2053 CP2 =-FAC*X*TMX/DEN2 2054 CP3 =FAC**2/DEN3 2055 CP4 =X**2*(P1+X*(P2+X*P3))/DEN4 2056 RETURN 2057 END 2058C 2059C******************************************************************** 2060C --- COMPACTS SEVERAL VECTORS AND MATRIX IN A BIG ONE: CDUM 2061 SUBROUTINE CPCT(MODE,NQ,NV,NU,NL,NM,NDIM,MDIM,NMDIM,LDG,LDF, 2062 & LDA,NZA,LIPAR,LXUMF,LIUMF,LR,LI,IUMF,IPAR,ISTAT,INDG1, 2063 & INDFL,IP,IDUM,H,TOLD,TCH,B10,CONTQ,CONTV,CONTU, 2064 & Q,V,U,GQ,FL,AVALUE,XUMF,CDUM) 2065 IMPLICIT REAL*8(A-H,O-Z) 2066 DIMENSION B10(4),CONTQ(5*NQ),CDUM(LR),CONTV(5*NV) 2067 DIMENSION Q(NQ),V(NV),FL(LDF,MDIM),AVALUE(2*NZA) 2068 DIMENSION IP(NM),INDG1(2*LDG),INDFL(2*LDF),IDUM(LI) 2069 DIMENSION IPAR(LIPAR),IUMF(LIUMF),U(NU),ISTAT(9) 2070 DIMENSION CONTU(5*NU),GQ(LDG,NDIM),XUMF(LXUMF) 2071C -- THE ENTRY POINTS FOR CDUM 2072 NQ2=2*NQ 2073 NQ3=3*NQ 2074 NQ4=4*NQ 2075 NV2=2*NV 2076 NV3=3*NV 2077 NV4=4*NV 2078 NU2=2*NU 2079 NU3=3*NU 2080 NU4=4*NU 2081 IC1=7 2082 IC2=7+NQ 2083 IC3=IC2+NQ 2084 IC4=IC3+NQ 2085 IC5=IC4+NQ 2086 IC6=IC5+NQ 2087 IC7=IC6+NV 2088 IC8=IC7+NV 2089 IC9=IC8+NV 2090 IC10=IC9+NV 2091 IC11=IC10+NV 2092 IC12=IC11+NU 2093 IC13=IC12+NU 2094 IC14 =IC13+NU 2095 IC15 =IC14+NU 2096 IC16 =IC15+NU 2097 IC17 =IC16+NQ 2098 IC18 =IC17+NV 2099 IC19 =IC18+NU 2100 IC20=IC19+2*LDG*NDIM+LDA*NMDIM+NM+NL 2101 IC21=IC20+LDF*MDIM 2102 IC22=IC21+2*NZA 2103C -- THE CONTAIN OF CDUM 2104 CDUM(1)=TOLD 2105 CDUM(2)=H 2106 CDUM(3)=TCH 2107 DO I=1,4 2108 CDUM(3+I)=B10(I) 2109 END DO 2110 IF ((MODE.GE.0).AND.(MODE.LE.3)) THEN 2111 GOTO 10 2112 ELSE IF ((MODE.EQ.4).OR.(MODE.EQ.5)) THEN 2113 GOTO 20 2114 ELSE 2115 STOP ' CPCT: INVALID MODE' 2116 END IF 2117 10 CONTINUE 2118 DO I=1,NV 2119 CDUM(IC1+I)=CONTQ(I) 2120 CDUM(IC2+I)=CONTQ(I+NQ) 2121 CDUM(IC3+I)=CONTQ(I+NQ2) 2122 CDUM(IC4+I)=CONTQ(I+NQ3) 2123 CDUM(IC5+I)=CONTQ(I+NQ4) 2124 CDUM(IC6+I)=CONTV(I) 2125 CDUM(IC7+I)=CONTV(I+NV) 2126 CDUM(IC8+I)=CONTV(I+NV2) 2127 CDUM(IC9+I)=CONTV(I+NV3) 2128 CDUM(IC10+I)=CONTV(I+NV4) 2129 CDUM(IC16+I)=Q(I) 2130 CDUM(IC17+I)=V(I) 2131 DO J=1,NL 2132 IMJ=(I-1)*NL+J 2133 JNI=(J-1)*NV+I 2134 CDUM(IC19+IMJ)=GQ(J,I) 2135 CDUM(IC20+JNI)=FL(I,J) 2136 END DO 2137 END DO 2138 GOTO 100 2139 20 CONTINUE 2140 DO I=1,NV 2141 CDUM(IC1+I)=CONTQ(I) 2142 CDUM(IC2+I)=CONTQ(I+NQ) 2143 CDUM(IC3+I)=CONTQ(I+NQ2) 2144 CDUM(IC4+I)=CONTQ(I+NQ3) 2145 CDUM(IC5+I)=CONTQ(I+NQ4) 2146 CDUM(IC6+I)=CONTV(I) 2147 CDUM(IC7+I)=CONTV(I+NV) 2148 CDUM(IC8+I)=CONTV(I+NV2) 2149 CDUM(IC9+I)=CONTV(I+NV3) 2150 CDUM(IC10+I)=CONTV(I+NV4) 2151 CDUM(IC16+I)=Q(I) 2152 CDUM(IC17+I)=V(I) 2153 END DO 2154 DO I=1,LDG 2155 CDUM(IC19+I)=GQ(I,1) 2156 END DO 2157 DO I=1,LDF 2158 CDUM(IC20+I)=FL(I,1) 2159 END DO 2160 100 CONTINUE 2161 DO I=NV+1,NQ 2162 CDUM(IC1+I)=CONTQ(I) 2163 CDUM(IC2+I)=CONTQ(I+NQ) 2164 CDUM(IC3+I)=CONTQ(I+NQ2) 2165 CDUM(IC4+I)=CONTQ(I+NQ3) 2166 CDUM(IC5+I)=CONTQ(I+NQ4) 2167 CDUM(IC16+I)=Q(I) 2168 END DO 2169 DO I=1,NU 2170 CDUM(IC11+I)=CONTU(I) 2171 CDUM(IC12+I)=CONTU(I+NU) 2172 CDUM(IC13+I)=CONTU(I+NU2) 2173 CDUM(IC14+I)=CONTU(I+NU3) 2174 CDUM(IC15+I)=CONTU(I+NU4) 2175 CDUM(IC18+I)=U(I) 2176 END DO 2177 DO I=1,LXUMF 2178 CDUM(IC22+I)=XUMF(I) 2179 END DO 2180C -- ENTRY POINTS FOR IDUM 2181 IPA = NM+1 2182 IST=IPA+LIPAR 2183 IG1=IST+9 2184 IFL=IG1+4*LDG 2185 IUM=IFL+2*LDF+2*NZA 2186C -- THE CONTAIN OF IDUM 2187 IDUM(1)=LIPAR 2188 DO I=1,NM 2189 IDUM(I+1)=IP(I) 2190 END DO 2191 DO I=1,LIPAR 2192 IDUM(IPA+I)=IPAR(I) 2193 END DO 2194 DO I=1,9 2195 IDUM(IST+I)=ISTAT(I) 2196 END DO 2197 DO I=1,2*LDG 2198 IDUM(IG1+I)=INDG1(I) 2199 END DO 2200 DO I=1,2*LDF 2201 IDUM(IFL+I)=INDFL(I) 2202 END DO 2203 DO I=1,LIUMF 2204 IDUM(IUM+I)=IUMF(I) 2205 END DO 2206 RETURN 2207 END 2208C*********************************************************************** 2209C --- COEFFICIENTS OF THE METHOD 2210 SUBROUTINE COEF(C2,C3,C4,C5,C6,C7,C10,A21,A31,A32, 2211 & A41,A42,A43,A51,A52,A53,A54,A61,A62,A63,A64,A65, 2212 & A71,A72,A73,A74,A75,A76,A81,A82,A83,A84,A85,A86, 2213 & A87,A106,A107,A109,B1,B2,B3,B4,B5,B6,B7,B8, 2214 & D14,D15,D16,D17,D18,D19,D24,D25,D26,D27,D28,D29, 2215 & D34,D35,D36,D37,D38,D39,D44,D45,D46,D47,D48,D49) 2216 IMPLICIT REAL*8 (A-H,O-Z) 2217 C2 = 0.1000000000000000D+00 2218 C3 = 0.1500000000000000D+00 2219 C4 = 0.3614787346573635D+00 2220 C5 = 0.8535584335486351D-01 2221 C6 = 0.5000000000000000D+00 2222 C7 = 0.8000000000000000D+00 2223 A21 = 0.1000000000000000D+00 2224 A31 = 0.3749999999999999D-01 2225 A32 = 0.1125000000000000D+00 2226 A41 = 0.3222169236216038D+00 2227 A42 = -0.1188883322987607D+01 2228 A43 = 0.1228145134023366D+01 2229 A51 = -0.3501123898129943D-01 2230 A52 = 0.3725420601086163D+00 2231 A53 = -0.2721053535582034D+00 2232 A54 = 0.1993037578575077D-01 2233 A61 = -0.5576547055042005D+00 2234 A62 = 0.1367307289645883D+01 2235 A63 = -0.1732236360460725D+01 2236 A64 = 0.4587772007467548D+00 2237 A65 = 0.9638065755722880D+00 2238 A71 = 0.8654517193566155D-01 2239 A72 = -0.8810082847945416D-01 2240 A73 = 0.1981275547329404D+00 2241 A74 = -0.4645422679331083D+00 2242 A75 = 0.1615170091109488D+00 2243 A76 = 0.9064533606330119D+00 2244 A81 = 0.0000000000000000D+00 2245 A82 = 0.0000000000000000D+00 2246 A83 = 0.0000000000000000D+00 2247 A84 = 0.3624477248753816D+01 2248 A85 = -0.4617724189181256D+00 2249 A86 = -0.3198024628164272D+01 2250 A87 = 0.1035319798328740D+01 2251 B1 = 0.0000000000000000D+00 2252 B2 = 0.0000000000000000D+00 2253 B3 = 0.0000000000000000D+00 2254 B4 = 0.2467760667636791D+00 2255 B5 = 0.2106594087489728D+00 2256 B6 = 0.1769218149125021D+00 2257 B7 = 0.3064446444147922D+00 2258 B8 = 0.5919806516005373D-01 2259C-- COEFF FOR DENSE OUTPUT 2260 C10 = 0.05D0 2261 A106 = 0.2519444444444446D0 2262 A107 = -0.3861111111111117D0 2263 A109 = 0.1841666666666670D0 2264 D14 = -0.2661262387118526D0 2265 D15 = -0.4784628694418738D0 2266 D16 = -8.4922471157927090D-02 2267 D17 = 0.2500588298423810D0 2268 D18 = -0.1047577768468651D0 2269 D19 = 1.684210526316138D0 2270 D24 = 2.963408265225625D0 2271 D25 = 5.405521925140019D0 2272 D26 = 0.9129165649477744D0 2273 D27 = -2.963932600778991D0 2274 D28 = 1.261033213890067D0 2275 D29 = -7.578947368424495D0 2276 D34 = -4.141333547260734D0 2277 D35 = -8.533017606960350D0 2278 D36 = -0.8633784567714418D0 2279 D37 = 6.403467289689607D0 2280 D38 = -2.971000836599169D0 2281 D39 = 10.10526315790209D0 2282 D44 = 1.690827587510642D0 2283 D45 = 3.816617960011145D0 2284 D46 = 0.2123061778941055D0 2285 D47 = -3.383148874338215D0 2286 D48 = 1.873923464716025D0 2287 D49 = -4.210526315793703D0 2288 RETURN 2289 END 2290