C HBMCO2 SOURCE OF166741 26/05/11 21:15:07 12538 SUBROUTINE HBMCO2(KTKAM,KTQ,KTFEX,KTPAS,KTLIAA,KTEMP,KTLIAB, & KTPHI,KCPR,KOCLFA,KOCLB1,NHBM,NFFT,KPARNUM, & KSORT,NSPAS,ITER) *======================================================================= * Continuation par pseudo longueur d'arc * Cas des systemes autonomes *======================================================================= IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMLREEL -INC TMDYNC SEGMENT mwork REAL*8 XT(NDDL,NFFT), LAMBD(NDDL) REAL*8 FTEST(4), FTEST0(4) REAL*8 XAUX(NFFT) REAL*8 VCTCS(7),AiDi(2),Di(2),bi(2) ENDSEGMENT INTEGER NSPAS, IREDU, METSTB LOGICAL CHECK, CHECKBIF LOGICAL ZINIT,ZSTAB CHARACTER*8 FLAG REAL*8 ZERO,ONE,TWO,FOUR PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, FOUR=4.D0) C Fonctions BLAS/LAPACK REAL*8 DDOT, DNRM2 EXTERNAL DDOT, DNRM2 * Variables generalisees MTQ = KTQ * Partie lineaire MTKAM = KTKAM * Parametres numeriques PARNUM = KPARNUM * Reste des segments MTPHI = KTPHI MTFEX = KTFEX MTPAS = KTPAS MTLIAA = KTLIAA MTLIAB = KTLIAB LOCLFA = KOCLFA LOCLB1 = KOCLB1 MTEMP = KTEMP * Tableau des resultats PSORT = KSORT * NT = Q1(/1) NDDL = NT/(2*NHBM+1) PDT = 1./NFFT SEGINI,MWORK *======================================================================= *=== INITIALISATION *======================================================================= II=1 * Recuperation des coefficients si modele GP DO J = 1,IPALA(/1) IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN DO IJ = 1,7 VCTCS(IJ) = XPALA(J,IJ) ENDDO JG = NDDL MLREE1 = IPALA(J,7) SEGACT, MLREE1 DO IJ = 1,JG LAMBD(IJ) = MLREE1.PROG(IJ) ENDDO * Nombre de termes fixe, a generaliser si besoin. JG = 2 MLREE2 = IPALA(J,4) MLREE3 = IPALA(J,8) SEGACT, MLREE2,MLREE3 DO IJ = 1,JG AiDi(IJ) = MLREE2.PROG(IJ) Di(IJ) = MLREE3.PROG(IJ) ENDDO SEGDES,MLREE1,MLREE2,MLREE3 ENDIF ENDDO VV = VCTCS(4) DO I=1,NT QOLD(I) = Q1(I) ENDDO OMEGOLD = OMEG VVOLD = VV * Derivee du residu par rapport au parametre de continuation: dR/da * ATTENTION: * L'appel direct a calcRv est a remplacer par l'appel a une sous- * routine qui distingue selon les differentes possibilites associees * a chaque liaison. CALL HBMRV(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Ra) DO I = 1,NT RHS(I) = Ra(I) ENDDO RHS(NT+1) = ZERO * DO I = 1,NT+1 * WRITE(*,*) 'Ja(',I,',:)=',(Ja(I,IOU),IOU=1,NT+1) * ENDDO * CALL COPYMAT(NT+1,Ja,J2) * WRITE(*,*) 'RHS=',(RHS(IOU),IOU=1,NT+1) CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO) * WRITE(*,*) 'INFO=',INFO * WRITE(*,*) 'pasT=',(RHS(IOU),IOU=1,NT+1) * STOP DO I=1,NT dX(I) = -RHS(I) ENDDO dw = -RHS(NT+1) * WRITE(*,*) 'dX=',(dX(IOU),IOU=1,NT) a = dnrm2(NT,dX,1) dv = ONE/sqrt(a**2+dw**2+ONE) * WRITE(*,*) 'dw=',dw * WRITE(*,*) 'dv=',dv IF (ISENS.LT.ZERO) THEN dv = -dv ENDIF t02(NT+2)=dv DO J=1,NT dX(J) = dX(J)*dv t02(J) = dX(J) ENDDO dw = dw*dv t02(NT+1) = dw * WRITE(*,*) 't0 = ', (t02(IOU),IOU=1,NT+2) * STOP * Sauvegarde des donnees de sortie DO I=1,NT QSAVE(I,II)=Q1(I) ENDDO WSAVE(II)=OMEG VSAVE(II)= VV ZSAVE(II)= .TRUE. * Message relatif au pas #0 (solution initiale) IF(IIMPI.GE.1) THEN WRITE(IOIMP,*) '+------+------+---------+---------+------+' WRITE(IOIMP,*) '| Pas | Iter | w | Q | Stab |' WRITE(IOIMP,*) '+------+------+---------+---------+------+' Q1NRM2=dnrm2(NT,Q1,1) WRITE(IOIMP,666) (II-1),ITER,VV,Q1NRM2,ZSAVE(II) ENDIF 666 FORMAT(' | ',I4,' | ',I4,' | ',F7.3,' | ',F7.3,' | ',L3,' |') *======================================================================= * 1ere prediction *======================================================================= OMEG = OMEG + DS*dw DO J=1,NT Q1(J) = Q1(J)+DS*dX(J) ENDDO VV = VV + DS*dv IREDU = 0 II = 2 *======================================================================= *=============== Boucle sur le parametre de continuation =============== *======================================================================= DO WHILE (II .LE. NBPAS) * === Correction ================================================= CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM, & MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'CA',ITER) * ----------------------------- * ---- si NEWT a converge, ---- * ----------------------------- IF (.NOT.CHECK) THEN * Sauvegarde des donnees de sortie: * Frequence, coeffs de Fourier, Norme des coeffs de Fourier * et parametre de continuation DO I=1,NT QSAVE(I,II)=Q1(I) ENDDO WSAVE(II)=OMEG VSAVE(II)=VV ZSAVE(II)= .TRUE. * * === Prediction ============================================== * Pas tangent a la courbe CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) DO J = 1,IPALA(/1) IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN VCTCS(4) = VV CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw) ENDIF ENDDO CALL HBMRV(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Ra) DO KK= 1,NT DO JJ = 1,NT Ja(JJ,KK) = RX(JJ,KK) ENDDO Ja(KK,NT+1) = Rw(KK) Ja(NT+1,KK) = ZERO ENDDO DO JJ = 2,2*NHBM,2 Ja(NT+1,JJ*NDDL+1) = JJ/TWO ENDDO DO I = 1,NT RHS(I) = Ra(I) ENDDO RHS(NT+1) = ZERO Ja(NT+1,NT+1) = ZERO CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO) DO J=1,NT dX(J) = -RHS(J) tp2(J) = dX(J) ENDDO dw = -RHS(NT+1) tp2(NT+1) = dw tp2(NT+2) = ONE a=dnrm2(NT,dX,1) dv = ONE/sqrt(ONE+a**2+dw**2) IF (II.EQ.2) THEN IF (ISENS.LT.ZERO) THEN dv = -dv ENDIF ELSE aux = DDOT(NT+2,t02,1,tp2,1) dv = (SIGN(dv,aux)) ENDIF DO K = 1,NT dX(K) = dX(K)*dv tp2(K) = dX(K) ENDDO dw = dw*dv tp2(NT+1) = dw tp2(NT+2) = dv * === Message relatif au pas #II ========== IF(IIMPI.GE.1) THEN Q1NRM2=dnrm2(NT,Q1,1) WRITE(IOIMP,666) (II-1),ITER,VV,Q1NRM2,.true. ENDIF * * === Tests d'arret =========================================== * On verifie la condition d'arret par rapport a VV c IF ((VV.GT.PARFIN).OR.(VV.LT.PARINI)) THEN IF ((VV-PARFIN)*(PARFIN-PARINI).GE.0) THEN WRITE(IOIMP,*) 'Fin de la continuation apres',II,' pas.' IF (II.LT.NBPAS) THEN * WRITE(*,*) 'NPAS sera egal a:',II NT1 = QSAVE(/1) NA1 = LSAVE(/2)/2 * NPAS = QSAVE(/2) NBIFU = QBIFU(/1) * write(*,*) NT1,NA1,NPAS,NBIFU NPAS = II SEGADJ, PSORT ENDIF KSORT = PSORT RETURN ENDIF * === pour le prochain pas, =================================== * ajustement automatique de la longueur du pas CALL ADPAS(DS,DSMIN,DSMAX,ITER,ITERMOY,ANGMIN,ANGMAX,NT+2, & t02,tp2) * donnees utiles DO I=1,NT+2 t02(I) = tp2(I) ENDDO DO KK = 1,NT QOLD(KK) = Q1(KK)+DS*dX(KK)/FOUR Q1(KK) = Q1(KK)+DS*dX(KK) ENDDO OMEGOLD = OMEG + DS*dw/FOUR VVOLD = VV + DS*dv/FOUR OMEG = OMEG + DS*dw VV = VV + DS*dv IREDU = 0 * ----------------------------------- * ---- si NEWT n'a pas converge, ---- * ----------------------------------- ELSE * Revenir en arriere, en diminuant le pas IF (II.GT.1) THEN DO I=1,NT Q1(I) = QOLD(I) ENDDO OMEG = OMEGOLD VV = VVOLD DS = DS/FOUR II = II-1 ENDIF c apres reduction de 0.25**10 = 1.E-6 , echec ! c 0.25**5 = 1.E-3 c IF (IREDU.GT.10) THEN IF (IREDU.GT.5) THEN WRITE(IOIMP,*)'Non-convergence apres ',IREDU,' reductions' IF (II.LT.NBPAS) THEN NT1 = QSAVE(/1) NA1 = LSAVE(/2)/2 * NPAS = QSAVE(/2) NBIFU = QBIFU(/1) NPAS = II SEGADJ, PSORT ENDIF KSORT = PSORT RETURN ENDIF IREDU = IREDU + 1 IF (IIMPI.GE.3) THEN WRITE(IOIMP,*) 'Reduction du pas #',IREDU ENDIF ENDIF * ---------------------------------------------- * ---- fin distinction NEWT converge ou pas ---- * ---------------------------------------------- * incrementation du compteur II = II+1 ENDDO *======================================================================= * fin de la boucle sur le parametre de continuation *======================================================================= * Message WRITE(IOIMP,*) NBPAS, 'pas de continuation realises!' KSORT = PSORT SEGSUP,MWORK RETURN END