Télécharger la plaquette Cast3M


fr en

Télécharger hbmco2.eso

Retour à la liste

Numérotation des lignes : 

hbmco2

  1. C HBMCO2    SOURCE    OF166741  26/05/11    21:15:07     12538          
  2.  
  3.       SUBROUTINE HBMCO2(KTKAM,KTQ,KTFEX,KTPAS,KTLIAA,KTEMP,KTLIAB,
  4.      &                  KTPHI,KCPR,KOCLFA,KOCLB1,NHBM,NFFT,KPARNUM,
  5.      &                  KSORT,NSPAS,ITER)
  6.  
  7. *=======================================================================
  8. *                Continuation par pseudo longueur d'arc
  9. *                     Cas des systemes autonomes
  10. *=======================================================================
  11.  
  12.       IMPLICIT INTEGER(I-N)
  13.       IMPLICIT REAL*8 (A-H,O-Z)
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17.  
  18. -INC SMLREEL
  19.  
  20. -INC TMDYNC
  21.  
  22.       SEGMENT mwork
  23.         REAL*8 XT(NDDL,NFFT), LAMBD(NDDL)
  24.         REAL*8 FTEST(4), FTEST0(4)
  25.         REAL*8 XAUX(NFFT)
  26.         REAL*8 VCTCS(7),AiDi(2),Di(2),bi(2)
  27.       ENDSEGMENT
  28.  
  29.       INTEGER NSPAS, IREDU, METSTB
  30.       LOGICAL CHECK, CHECKBIF
  31.       LOGICAL ZINIT,ZSTAB
  32.       CHARACTER*8 FLAG
  33.  
  34.       REAL*8 ZERO,ONE,TWO,FOUR
  35.       PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, FOUR=4.D0)
  36.  
  37. C Fonctions BLAS/LAPACK
  38.       REAL*8    DDOT, DNRM2
  39.       EXTERNAL  DDOT, DNRM2
  40.  
  41. *      Variables generalisees
  42.        MTQ = KTQ
  43. *      Partie lineaire
  44.        MTKAM = KTKAM
  45. *      Parametres numeriques
  46.        PARNUM = KPARNUM
  47. *      Reste des segments
  48.        MTPHI  = KTPHI
  49.        MTFEX  = KTFEX
  50.        MTPAS  = KTPAS
  51.        MTLIAA = KTLIAA
  52.        MTLIAB = KTLIAB
  53.        LOCLFA = KOCLFA
  54.        LOCLB1 = KOCLB1
  55.        MTEMP  = KTEMP
  56. *      Tableau des resultats
  57.        PSORT  = KSORT
  58. *
  59.        NT = Q1(/1)
  60.        NDDL = NT/(2*NHBM+1)
  61.        PDT = 1./NFFT
  62.  
  63.       SEGINI,MWORK
  64.  
  65. *=======================================================================
  66. *===   INITIALISATION
  67. *=======================================================================
  68.  
  69.       II=1
  70. *  Recuperation des coefficients si modele GP
  71.       DO J = 1,IPALA(/1)
  72.          IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  73.           DO IJ = 1,7
  74.             VCTCS(IJ) = XPALA(J,IJ)
  75.           ENDDO
  76.           JG = NDDL
  77.           MLREE1 = IPALA(J,7)
  78.           SEGACT, MLREE1
  79.           DO IJ = 1,JG
  80.             LAMBD(IJ) = MLREE1.PROG(IJ)
  81.           ENDDO
  82. *         Nombre de termes fixe, a generaliser si besoin.
  83.              JG = 2
  84.           MLREE2 = IPALA(J,4)
  85.           MLREE3 = IPALA(J,8)
  86.           SEGACT, MLREE2,MLREE3
  87.           DO IJ = 1,JG
  88.             AiDi(IJ) = MLREE2.PROG(IJ)
  89.             Di(IJ) = MLREE3.PROG(IJ)
  90.           ENDDO
  91.           SEGDES,MLREE1,MLREE2,MLREE3
  92.          ENDIF
  93.       ENDDO
  94.       VV = VCTCS(4)
  95.       DO I=1,NT
  96.         QOLD(I) = Q1(I)
  97.       ENDDO
  98.       OMEGOLD = OMEG
  99.       VVOLD = VV
  100. *      Derivee du residu par rapport au parametre de continuation: dR/da
  101. *      ATTENTION:
  102. *      L'appel direct a calcRv est a remplacer par l'appel a une sous-
  103. *      routine qui distingue selon les differentes possibilites associees
  104. *      a chaque liaison.
  105.        CALL HBMRV(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Ra)
  106.        DO I = 1,NT
  107.           RHS(I) = Ra(I)
  108.        ENDDO
  109.        RHS(NT+1) = ZERO
  110. *      DO I = 1,NT+1
  111. *       WRITE(*,*) 'Ja(',I,',:)=',(Ja(I,IOU),IOU=1,NT+1)
  112. *      ENDDO
  113. *      CALL COPYMAT(NT+1,Ja,J2)
  114. *      WRITE(*,*) 'RHS=',(RHS(IOU),IOU=1,NT+1)
  115.        CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO)
  116. *      WRITE(*,*) 'INFO=',INFO
  117. *      WRITE(*,*) 'pasT=',(RHS(IOU),IOU=1,NT+1)
  118. *      STOP
  119.        DO I=1,NT
  120.           dX(I) = -RHS(I)
  121.        ENDDO
  122.        dw = -RHS(NT+1)
  123. *       WRITE(*,*) 'dX=',(dX(IOU),IOU=1,NT)
  124.        a = dnrm2(NT,dX,1)
  125.        dv = ONE/sqrt(a**2+dw**2+ONE)
  126. *       WRITE(*,*) 'dw=',dw
  127. *       WRITE(*,*) 'dv=',dv
  128.        IF (ISENS.LT.ZERO) THEN
  129.          dv = -dv
  130.        ENDIF
  131.        t02(NT+2)=dv
  132.        DO J=1,NT
  133.           dX(J)  = dX(J)*dv
  134.           t02(J) = dX(J)
  135.        ENDDO
  136.        dw = dw*dv
  137.        t02(NT+1) = dw
  138. *       WRITE(*,*) 't0 = ', (t02(IOU),IOU=1,NT+2)
  139. *       STOP
  140. *      Sauvegarde des donnees de sortie
  141.        DO I=1,NT
  142.          QSAVE(I,II)=Q1(I)
  143.        ENDDO
  144.        WSAVE(II)=OMEG
  145.        VSAVE(II)= VV
  146.        ZSAVE(II)= .TRUE.
  147.  
  148. *      Message relatif au pas #0 (solution initiale)
  149.        IF(IIMPI.GE.1) THEN
  150.          WRITE(IOIMP,*) '+------+------+---------+---------+------+'
  151.          WRITE(IOIMP,*) '| Pas  | Iter |    w    |    Q    | Stab |'
  152.          WRITE(IOIMP,*) '+------+------+---------+---------+------+'
  153.          Q1NRM2=dnrm2(NT,Q1,1)
  154.          WRITE(IOIMP,666) (II-1),ITER,VV,Q1NRM2,ZSAVE(II)
  155.        ENDIF
  156.  666   FORMAT(' | ',I4,' | ',I4,' | ',F7.3,' | ',F7.3,' | ',L3,'  |')
  157.  
  158. *=======================================================================
  159. *      1ere prediction
  160. *=======================================================================
  161.        OMEG = OMEG + DS*dw
  162.        DO J=1,NT
  163.          Q1(J) = Q1(J)+DS*dX(J)
  164.        ENDDO
  165.        VV    = VV + DS*dv
  166.        IREDU = 0
  167.        II = 2
  168.  
  169. *=======================================================================
  170. *=============== Boucle sur le parametre de continuation ===============
  171. *=======================================================================
  172.       DO WHILE (II .LE. NBPAS)
  173.  
  174. *       === Correction =================================================
  175.         CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
  176.      &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'CA',ITER)
  177.  
  178. *       -----------------------------
  179. *       ---- si NEWT a converge, ----
  180. *       -----------------------------
  181.         IF (.NOT.CHECK) THEN
  182.  
  183. *          Sauvegarde des donnees de sortie:
  184. *          Frequence, coeffs de Fourier, Norme des coeffs de Fourier
  185. *          et parametre de continuation
  186.            DO I=1,NT
  187.              QSAVE(I,II)=Q1(I)
  188.            ENDDO
  189.            WSAVE(II)=OMEG
  190.            VSAVE(II)=VV
  191.            ZSAVE(II)= .TRUE.
  192. *
  193. *          === Prediction ==============================================
  194. *          Pas tangent a la courbe
  195.            CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  196.            DO J = 1,IPALA(/1)
  197.               IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  198.                 VCTCS(4) = VV
  199.                 CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)
  200.               ENDIF
  201.            ENDDO
  202.            CALL HBMRV(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Ra)
  203.            DO KK= 1,NT
  204.              DO JJ = 1,NT
  205.                Ja(JJ,KK) = RX(JJ,KK)
  206.              ENDDO
  207.              Ja(KK,NT+1) = Rw(KK)
  208.              Ja(NT+1,KK) = ZERO
  209.            ENDDO
  210.            DO JJ = 2,2*NHBM,2
  211.              Ja(NT+1,JJ*NDDL+1) = JJ/TWO
  212.            ENDDO
  213.            DO I = 1,NT
  214.              RHS(I) = Ra(I)
  215.            ENDDO
  216.            RHS(NT+1) = ZERO
  217.            Ja(NT+1,NT+1) = ZERO
  218.            CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO)
  219.            DO J=1,NT
  220.               dX(J) = -RHS(J)
  221.               tp2(J) = dX(J)
  222.            ENDDO
  223.            dw = -RHS(NT+1)
  224.            tp2(NT+1) = dw
  225.            tp2(NT+2) = ONE
  226.            a=dnrm2(NT,dX,1)
  227.            dv = ONE/sqrt(ONE+a**2+dw**2)
  228.            IF (II.EQ.2) THEN
  229.              IF (ISENS.LT.ZERO) THEN
  230.                dv = -dv
  231.              ENDIF
  232.            ELSE
  233.              aux = DDOT(NT+2,t02,1,tp2,1)
  234.              dv = (SIGN(dv,aux))
  235.            ENDIF
  236.            DO K = 1,NT
  237.               dX(K) = dX(K)*dv
  238.               tp2(K) = dX(K)
  239.            ENDDO
  240.            dw = dw*dv
  241.            tp2(NT+1) = dw
  242.            tp2(NT+2) = dv
  243.  
  244. *          === Message relatif au pas #II ==========
  245.            IF(IIMPI.GE.1) THEN
  246.               Q1NRM2=dnrm2(NT,Q1,1)
  247.               WRITE(IOIMP,666) (II-1),ITER,VV,Q1NRM2,.true.
  248.            ENDIF
  249. *
  250.  
  251. *          === Tests d'arret ===========================================
  252. *          On verifie la condition d'arret par rapport a VV
  253. c            IF ((VV.GT.PARFIN).OR.(VV.LT.PARINI)) THEN
  254.            IF ((VV-PARFIN)*(PARFIN-PARINI).GE.0) THEN
  255.                WRITE(IOIMP,*) 'Fin de la continuation apres',II,' pas.'
  256.                IF (II.LT.NBPAS) THEN
  257. *                  WRITE(*,*) 'NPAS sera egal a:',II
  258.                   NT1 = QSAVE(/1)
  259.                   NA1 = LSAVE(/2)/2
  260. *                  NPAS = QSAVE(/2)
  261.                   NBIFU = QBIFU(/1)
  262. *                  write(*,*) NT1,NA1,NPAS,NBIFU
  263.                   NPAS = II
  264.                   SEGADJ, PSORT
  265.                ENDIF
  266.                KSORT = PSORT
  267.                RETURN
  268.            ENDIF
  269.  
  270. *          === pour le prochain pas, ===================================
  271. *          ajustement automatique de la longueur du pas
  272.            CALL ADPAS(DS,DSMIN,DSMAX,ITER,ITERMOY,ANGMIN,ANGMAX,NT+2,
  273.      &                t02,tp2)
  274. *          donnees utiles
  275.            DO I=1,NT+2
  276.               t02(I) = tp2(I)
  277.            ENDDO
  278.            DO KK = 1,NT
  279.               QOLD(KK) = Q1(KK)+DS*dX(KK)/FOUR
  280.               Q1(KK) = Q1(KK)+DS*dX(KK)
  281.            ENDDO
  282.            OMEGOLD = OMEG + DS*dw/FOUR
  283.            VVOLD   = VV + DS*dv/FOUR
  284.            OMEG = OMEG + DS*dw
  285.            VV   = VV   + DS*dv
  286.            IREDU = 0
  287.  
  288. *       -----------------------------------
  289. *       ---- si NEWT n'a pas converge, ----
  290. *       -----------------------------------
  291.         ELSE
  292.  
  293. *          Revenir en arriere, en diminuant le pas
  294.            IF (II.GT.1) THEN
  295.              DO I=1,NT
  296.                Q1(I) = QOLD(I)
  297.              ENDDO
  298.              OMEG = OMEGOLD
  299.              VV = VVOLD
  300.              DS = DS/FOUR
  301.              II = II-1
  302.            ENDIF
  303. c          apres reduction de 0.25**10 = 1.E-6 , echec !
  304. c                             0.25**5  = 1.E-3
  305. c            IF (IREDU.GT.10) THEN
  306.            IF (IREDU.GT.5) THEN
  307.               WRITE(IOIMP,*)'Non-convergence apres ',IREDU,' reductions'
  308.               IF (II.LT.NBPAS) THEN
  309.                   NT1 = QSAVE(/1)
  310.                   NA1 = LSAVE(/2)/2
  311. *                  NPAS = QSAVE(/2)
  312.                   NBIFU = QBIFU(/1)
  313.                   NPAS = II
  314.                   SEGADJ, PSORT
  315.                ENDIF
  316.                KSORT = PSORT
  317.                RETURN
  318.            ENDIF
  319.            IREDU = IREDU + 1
  320.            IF (IIMPI.GE.3) THEN
  321.              WRITE(IOIMP,*) 'Reduction du pas #',IREDU
  322.            ENDIF
  323.  
  324.         ENDIF
  325. *       ----------------------------------------------
  326. *       ---- fin distinction NEWT converge ou pas ----
  327. *       ----------------------------------------------
  328.  
  329. *        incrementation du compteur
  330.          II = II+1
  331.       ENDDO
  332. *=======================================================================
  333. *       fin de la boucle sur le parametre de continuation
  334. *=======================================================================
  335.  
  336. *  Message
  337.       WRITE(IOIMP,*) NBPAS, 'pas de continuation realises!'
  338.       KSORT = PSORT
  339.  
  340.       SEGSUP,MWORK
  341.  
  342.       RETURN
  343.       END
  344.  
  345.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales