Télécharger la plaquette Cast3M


fr en

Télécharger hbmnewt.eso

Retour à la liste

Numérotation des lignes : 

hbmnewt

  1. C HBMNEWT   SOURCE    OF166741  26/05/11    21:15:13     12538          
  2.  
  3. *=======================================================================
  4. *
  5. *                     Methode de Newton-Raphson
  6. *
  7. * Convergence vers une solution des equations d'equilibre (frequentiel).
  8. * TYPC(1) indique le type de calcul a realiser:
  9. *      = | 'I' : solution initiale
  10. *        | 'C' : pas correcteur de la continuation
  11. * TYPC(2) indique le type de calcul a realiser:
  12. *      = | 'F' : reponse Forcee
  13. *        | 'A' : reponse Autonome
  14. *        | 'N' : reponse libre : Non-linear normal modes
  15. *
  16. *=======================================================================
  17.  
  18.       SUBROUTINE HBMNEWT(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,
  19.      &                   KPARNUM,KTLIAA,MTLIAB,KTFEX,MTPAS,LOCLFA,
  20.      &                   LOCLB1,CHECK,TYPC,ITER)
  21.  
  22. *----- Declarations ----------------------------------------------------
  23.  
  24.       IMPLICIT INTEGER(I-N)
  25.       IMPLICIT REAL*8(A-H,O-Z)
  26.  
  27. -INC PPARAM
  28. -INC CCOPTIO
  29.  
  30. -INC SMLREEL
  31.       POINTEUR fvec.MLREEL
  32.  
  33. -INC TMDYNC
  34.  
  35.       SEGMENT,mwrkl
  36.         REAL*8 AVVP(NT),VCTCS(7),LAMBD(NDDL)
  37.         REAL*8 J2(NT+1,NT+1), Ja2(NT+2,NT+2)
  38.         REAL*8 gx(NT)
  39.         REAL*8 AiDi(2),Di(2),bi(2)
  40.       ENDSEGMENT
  41.  
  42.       CHARACTER*2 TYPC
  43.       LOGICAL check
  44.  
  45.       LOGICAL NNM,JANAL2
  46.       INTEGER INFO
  47.  
  48.       REAL*8 fref,fref0
  49.       REAL*8 ZERO,TWO,ONE
  50.       PARAMETER (ZERO=0.D0, TWO=2.D0,ONE=1.D0)
  51.  
  52. C Fonctions BLAS/LAPACK
  53.       REAL*8    DDOT, DNRM2
  54.       EXTERNAL  DDOT, DNRM2
  55.  
  56. c      NNM=(TYPC.EQ.'IN').OR.(TYPC.EQ.'CN')
  57.       NNM=TYPC(2:2).EQ.'N'
  58. c*dbg      write(ioimp,*) 'HBMNEWT : TYPC ='//TYPC//'='
  59.  
  60. * Segments
  61.       MTQ   = KTQ
  62.       MTKAM = KTKAM
  63.       MTEMP = KTEMP
  64.       PARNUM = KPARNUM
  65.       MTLIAA= KTLIAA
  66. c      MTLIAB= KTLIAB
  67.       MTFEX = KTFEX
  68.       ITER=-1
  69. c      RELAX=1.d0
  70.       JANAL2=JANAL
  71.  
  72.       SEGINI,mwrkl
  73.       JG=NT
  74.       SEGINI,fvec
  75.  
  76.       fref  = ZERO
  77.       fref0 = ZERO
  78.  
  79. c *------ Test de convergence pour la solution ------------------
  80. c       test=0.D0
  81. c       DO i=1,NT
  82. c         IF(abs(fvec.prog(i)).GT.test) test=abs(fvec.prog(i))
  83. c       ENDDO
  84.  
  85. *------ Calcul de la raideur dynamique tangente des efforts lineaires --
  86.  
  87. *   Matrice liee a la base modale (avec ou sans terme dissipatif selon NNM)
  88.       CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.not.NNM)
  89.  
  90. *   cas Granger-Paidoussis : rajouter raideur de couplage
  91. c   aussi pour cas NNM ? a voir...
  92.       IF (.NOT. NNM) THEN
  93.         DO J = 1,IPALA(/1)
  94.           IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  95.             DO IJ = 1,7
  96.               VCTCS(IJ) = XPALA(J,IJ)
  97.             ENDDO
  98.             IF (TYPC.EQ.'CA') THEN
  99.               VCTCS(4) = VV
  100.               XPALA(J,4) = VV
  101.             ENDIF
  102.             JG = NDDL
  103.             MLREE1 = IPALA(J,7)
  104.             MLREE2 = IPALA(J,4)
  105.             MLREE3 = IPALA(J,8)
  106. c            SEGACT, MLREE1,MLREE2,MLREE3
  107.             DO IJ = 1,JG
  108.               LAMBD(IJ) = MLREE1.PROG(IJ)
  109.             ENDDO
  110. c *           Nombre de termes fixe, a generaliser si besoin.
  111. c              JG = 2
  112. *           tentative de generalisation
  113.             JG = MLREE2.PROG(/1)
  114.             DO IJ = 1,JG
  115.               AiDi(IJ) = MLREE2.PROG(IJ)
  116.               Di(IJ)   = MLREE3.PROG(IJ)
  117.             ENDDO
  118. c            SEGDES,MLREE1,MLREE2,MLREE3
  119.             CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)
  120.           ENDIF
  121.         ENDDO
  122.       ENDIF
  123.  
  124. *------ Calcul du residu initial ---------------------------------------
  125. *       (lie a la prediction par exemple)
  126.       its=0
  127.       CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  128.      &         MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,JANAL2,fvec,fref0,NNM)
  129.  
  130.       IF (IIMPI.GE.2) THEN
  131.         Q1NRM2=dnrm2(NT,Q1,1)
  132.         res=dnrm2(NT,fvec.prog,1)/fref0
  133.         WRITE(IOIMP,990) its,OMEG,Q1NRM2,res,fref0,DS
  134. c        IF (IIMPI.GE.333) THEN
  135. c          WRITE(IOIMP,111) 'Q',its,(Q1(iou),iou=1,NT)
  136. c          WRITE(IOIMP,111) 'R',its,(fvec.prog(iou),iou=1,NT)
  137. c        ENDIF
  138.       ENDIF
  139.  
  140. *============================= Iterations ==============================
  141.       DO its=1,ITERMAX
  142.  
  143. *    ----- DIFFICULTE A CONVERGER : ON RELAXE (test non concluant) ----
  144. c            IF (its .GT. 10) THEN
  145. c              RELAX=0.222d0
  146. c            ENDIF
  147.  
  148. *    ----- CALCUL DE LA JACOBIENNE : dR/dX = Z(OMEGA) - dFNLdX ----
  149.         DO J=1,NT
  150.           DO I=1,NT
  151.             RX(I,J) = ZZ(I,J)-JAC(I,J)
  152.           ENDDO
  153.         ENDDO
  154. c        IF (IIMPI.GE.333)  THEN
  155. c          DO iou=1,NT
  156. c            WRITE(IOIMP,111) 'RX_',iou,(RX(iou,jou),jou=1,NT)
  157. c          ENDDO
  158. c        ENDIF
  159.  
  160. *    ----- CAS 1: INITIALISATION (systeme force) ------------------
  161.         IF (TYPC.EQ.'IF') THEN
  162. *           Calcul de la correction dX, telle que: Rx*dX = R
  163. *           Appel a "DGESV" (LAPACK)
  164.           CALL DGESV(NT,1,RX,NT,IPIV,fvec.prog,NT,INFO)
  165.           DO M = 1,NT
  166.             Q1(M) = Q1(M)-fvec.prog(M)
  167.           ENDDO
  168. c          IF (IIMPI.GE.333)
  169. c     &      WRITE(IOIMP,111) 'dQ^',its,(fvec.prog(iou),iou=1,NT)
  170. cc           si on souhaite verifier la convergence en correction dX
  171. cc           res=dnrm2(NT,fvec.prog,1)/dnrm2(NT,Q1,1)
  172.         ENDIF
  173.  
  174. *    ----- CAS 2: CONTINUATION (systeme force) --------------------
  175.         IF (TYPC.EQ.'CF') THEN
  176. *       Calcul de la correction dY, telle que: Ja*dY = {R,0}
  177. *       avec Ja la Jacobienne augmentee: Ja = {RX,Rw;dX,dw}
  178.            DO JJ=1,NT
  179.              DO II=1,NT
  180.                Ja(II,JJ) = RX(II,JJ)
  181.               ENDDO
  182.           ENDDO
  183. *       Calcul de la derivee: R,w = dR/dw
  184.           CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  185. *       Cas Granger-Paidoussis : R,w = R,w - dFf/dw
  186.           DO J = 1,IPALA(/1)
  187.             IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  188.               VCTCS(4) = VV
  189.               CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)
  190.             ENDIF
  191.           ENDDO
  192. c       Cas d'un balourd : force en OMEG**2 a deriver aussi
  193.           IF (BAL.EQ.1) THEN
  194.             DO I=1,NT
  195.               Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  196.             ENDDO
  197.           ENDIF
  198.           DO KK = 1,NT
  199.             Ja(KK,NT+1) = Rw(KK)
  200.             Ja(NT+1,KK) = dX(KK)
  201. ctest -> nechange rien :  Ja(NT+1,KK) = DS*dX(KK)
  202. ctest                Ja(NT+1,KK) = (Q1(KK)-QOLD(KK))
  203.             RHS(KK) = fvec.prog(KK)
  204. cdebug
  205. c           if(OMEG.gt.0.153.and.OMEG.lt.0.170)
  206. c      &      write(*,123) KK,(Ja(KK,iou),iou=1,6),RHS(KK)
  207.           ENDDO
  208.           Ja(NT+1,NT+1) = dw
  209. ctest -> nechange rien :      Ja(NT+1,NT+1) = DS*dw
  210. ctest            Ja(NT+1,NT+1) = OMEG-OMEGOLD
  211.           RHS(NT+1) = ZERO
  212. c          if(OMEG.gt.0.153.and.OMEG.lt.0.170)
  213. c      &     write(*,123) NT+1,(Ja(NT+1,iou),iou=1,6),RHS(NT+1)
  214. *        Appel a "DGESV" (LAPACK)
  215.           CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO)
  216.           DO MM = 1,NT
  217. c            Q1(MM) = Q1(MM) - RELAX*RHS(MM)
  218.             Q1(MM) = Q1(MM) - RHS(MM)
  219.           ENDDO
  220. c           OMEG = OMEG - RELAX*RHS(NT+1)
  221.           OMEG = OMEG - RHS(NT+1)
  222.           CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  223. c       Si on souhaite verifier la convergence en correction dX
  224. c             res=dnrm2(NT,RHS,1)/dnrm2(NT,Q1,1)
  225. c*          write(*,*) 'W,Q1=',OMEG,(Q1(iou),iou=1,NT)
  226.         ENDIF
  227.  
  228. *    ----- CAS 3: INITIALISATION (systeme autonome) --------------
  229.         IF (TYPC.EQ.'IA') THEN
  230. *          Calcul de la correction dY, telle que: Ja*dY = -[R;g]
  231. *          avec: Ja = [Rx Rw; gx 0]
  232.           DO JJ = 1,NT
  233.             DO II=1,NT
  234.               Ja(II,JJ) = RX(II,JJ)
  235.             ENDDO
  236.           ENDDO
  237.           VV = VCTCS(4)
  238. *       Calcul de la derivee: Rw = dR/dw
  239.           CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  240. *       cas Granger-Paidoussis : R,w = R,w - dFf/dw
  241.           DO J = 1,IPALA(/1)
  242.             IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  243.               CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)
  244.             ENDIF
  245.           ENDDO
  246.           DO KK = 1,NT
  247.             Ja(KK,NT+1) = Rw(KK)
  248.             RHS(KK) = fvec.prog(KK)
  249.           ENDDO
  250. *           Condition de phase
  251. *           On impose une vitesse nulle au 1er DDL en t = 0:
  252. *           dx1(0) = 0 ==> sum_i{i*Q1_si} = 0
  253.           RHS(NT+1) = ZERO
  254.           DO I = 1,NT+1
  255.             Ja(NT+1,I) = ZERO
  256.           ENDDO
  257.           DO JJ = 2,2*NHBM,2
  258.             Ja(NT+1,JJ*NDDL+1) = JJ/TWO
  259. *            RHS(NT+1)  = RHS(NT+1) + (JJ/TWO)*Q1(JJ*NDDL+1)
  260.           ENDDO
  261. c           si on souhaite verifier la convergence en correction dX
  262. c             res=dnrm2(NT,RHS,1)
  263.           Ja(NT+1,NT+1) = ZERO
  264.           CALL COPYMAT(NT+1,Ja,J2)
  265. *           Appel a "DGESV" (LAPACK)
  266.           CALL DGESV(NT+1,1,J2,NT+1,IPIV2,RHS,NT+1,INFO)
  267.           DO M = 1,NT
  268.             Q1(M) = Q1(M)-RHS(M)
  269.           ENDDO
  270.           OMEG = OMEG-RHS(NT+1)
  271.           CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  272. *           cas Granger-Paidoussis : rajouter raideur de couplage
  273.           DO J = 1,IPALA(/1)
  274.             IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  275.               CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)
  276.             ENDIF
  277.           ENDDO
  278.         ENDIF
  279.  
  280. *    ----- CAS 4: CONTINUATION (systeme autonome) -----------------
  281.         IF (TYPC.EQ.'CA') THEN
  282. *    Calcul de la correction dY, telle que: Jaa*dY = -[R;g;0]
  283. *    avec Jaa la Jacobienne augmentee: Jaa ={RX,Rw,Ra;gx,0,0;dX,dw,da}
  284.           DO JJ = 1,NT
  285.             DO II=1,NT
  286.               Jaa(II,JJ) = RX(II,JJ)
  287.             ENDDO
  288.           ENDDO
  289. *       Calcul de la derivee: Rw = dR/dw
  290.           CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  291. *       cas Granger-Paidoussis : R,w = R,w - dFf/dw
  292.           DO J = 1,IPALA(/1)
  293.             IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  294.               VCTCS(4) = VV
  295.               CALL HBMRWF(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Rw)
  296.             ENDIF
  297.           ENDDO
  298.           DO I = 1,NT+2
  299.             Jaa(NT+2,I) = ZERO
  300.           ENDDO
  301. *           Calcul de la derivee : R,Vmoy = dR/dVmoy
  302.           CALL HBMRV(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,Ra)
  303.           DO KK = 1,NT
  304.             Jaa(KK,NT+1) = Rw(KK)
  305.             Jaa(KK,NT+2) = Ra(KK)
  306.             Jaa(NT+2,KK) = dX(KK)
  307.             Jaa(NT+1,KK) = ZERO
  308.             RHS2(KK) = fvec.prog(KK)
  309.           ENDDO
  310.           Jaa(NT+2,NT+1) = dw
  311.           Jaa(NT+2,NT+2) = dv
  312. *       Condition de phase
  313. *       On impose une vitesse nulle au 1er DDL en t = 0:
  314. *         dx1(0) = 0 ==> g(Q) = sum_i{i*Q1_si} = 0
  315. *           dg/dw
  316.           Jaa(NT+1,NT+1) = ZERO
  317. *           dg/dv
  318.           Jaa(NT+1,NT+2) = ZERO
  319.           RHS2(NT+1) = ZERO
  320.           DO JJ = 2,2*NHBM,2
  321.             Jaa(NT+1,JJ*NDDL+1) = JJ/2
  322.             RHS2(NT+1)  = RHS2(NT+1) + (JJ/2)*Q1(JJ*NDDL+1)
  323.           ENDDO
  324.           RHS2(NT+2) = ZERO
  325. c*          WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)
  326. *            STOP
  327.           CALL COPYMAT(NT+2,Jaa,Ja2)
  328. *           Appel a "DGESV" (LAPACK)
  329.           CALL DGESV(NT+2,1,Ja2,NT+2,IPIV3,RHS2,NT+2,INFO)
  330. c           si on souhaite verifier la convergence en correction dX
  331. c             res=dnrm2(NT+1,RHS2,1)
  332. c*          WRITE(*,*) 'res=',res
  333. c*          WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)
  334.           DO M = 1,NT
  335.             Q1(M) = Q1(M)-RHS2(M)
  336.           ENDDO
  337.           OMEG = OMEG-RHS2(NT+1)
  338.           VV   = VV  -RHS2(NT+2)
  339.           CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  340.           DO J = 1,IPALA(/1)
  341.             IF (IPALA(J,3).EQ.101) THEN
  342.               VCTCS(4) = VV
  343.               XPALA(J,4) = VV
  344.               CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)
  345.             ENDIF
  346.           ENDDO
  347.         ENDIF
  348.  
  349. *    ----- CAS 5: INITIALISATION (systeme Hamiltonien) --------------
  350.         IF (TYPC.EQ.'IN') THEN
  351. *      Calcul de la correction dY, telle que: Ja*dY = -[R;g]
  352. *          avec: Ja = [Rx Rw; gx 0]
  353.           DO JJ = 1,NT
  354.             DO II=1,NT
  355.               Ja(II,JJ) = RX(II,JJ)
  356.             ENDDO
  357.           ENDDO
  358. *       Calcul de la derivee: Rw = dR/dw
  359.           CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.false.)
  360.           DO KK = 1,NT
  361.             Ja(KK,NT+1) = Rw(KK)
  362.             RHS(KK) = fvec.prog(KK)
  363.           ENDDO
  364. *           Condition d'energie
  365. *           On initialise en imposant un certain niveau d'energie,
  366. *           i.e. on fixe la norme de Q1
  367.           RHS(NT+1) = 0.5D0*(DDOT(NT,Q1,1,Q1,1))-XPARA
  368. cbp    <=>  RHS(NT+1) = -0.5*(DDOT(NT,Q1,1,Q1,1))
  369. c           car XPARA=DDOT(Q1)=sqrt(Energie_initiale) dans hbmalo
  370. cbp : incoherence avec calcul energie realise dans hbmco3
  371. c           RHS(NT+1) = DDOT(NT,Q1,1,Q1,1)-XPARA**2
  372.           DO I = 1,NT
  373.             Ja(NT+1,I) = Q1(I)
  374.           ENDDO
  375. c           si on souhaite verifier la convergence en correction dX
  376.           res=dnrm2(NT+1,RHS,1)
  377.           Ja(NT+1,NT+1) = ZERO
  378. *            CALL COPYMAT(NT+1,Ja,J2)
  379. *           Appel a "DGESV" (LAPACK)
  380.           CALL DGESV(NT+1,1,Ja,NT+1,IPIV2,RHS,NT+1,INFO)
  381.           DO M = 1,NT
  382.             Q1(M) = Q1(M)-RHS(M)
  383.           ENDDO
  384.           OMEG = OMEG-RHS(NT+1)
  385.           CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.false.)
  386.         ENDIF
  387.  
  388. *    ----- CAS 6: CONTINUATION (systeme Hamiltonien) -----------------
  389.         IF (TYPC.EQ.'CN') THEN
  390. *          Calcul de la correction dY, telle que: Jaa*dY = -[R;g;0]
  391. *          avec Jaa la Jacobienne augmentee:
  392. *                       Jaa ={RX,Rw,Ra;gx,0,0;dX,dw,0}
  393. *          Ici, 'a' est un parametre artificiel qui ferme le probleme
  394. *          initialement sous-contraint.
  395.           DO JJ = 1,NT
  396.             DO II=1,NT
  397.               Jaa(II,JJ) = RX(II,JJ)
  398.             ENDDO
  399.           ENDDO
  400. *           Calcul de la derivee: Rw = dR/dw
  401.           CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.false.)
  402. *           Calcul de la derivee : Ra= (LoI)*X
  403.           CALL HBMDVEC(NT,NHBM,NDDL,Q1,1.D0,Ra)
  404.           DO KK = 1,NT
  405.             Jaa(KK,NT+1) = Rw(KK)
  406.             Jaa(KK,NT+2) = Ra(KK)
  407.             Jaa(NT+2,KK) = dX(KK)
  408.             RHS2(KK) = fvec.prog(KK)
  409.           ENDDO
  410.           Jaa(NT+2,NT+1) = dw
  411.           Jaa(NT+2,NT+2) = ZERO
  412. *           Condition de phase
  413. *           On minimise la variation de phase entre deux pas succesifs
  414. *           (condition integrale)
  415. *           dg/dw
  416.           Jaa(NT+1,NT+1) = ZERO
  417. *           dg/da
  418.           Jaa(NT+1,NT+2) = ZERO
  419.           CALL HBMDVEC(NT,NHBM,NDDL,QOLD,-1.D0,gx)
  420.           DO JJ = 1,NT
  421.             Jaa(NT+1,JJ) = gx(JJ)
  422.           ENDDO
  423.           RHS2(NT+1)  = DDOT(NT,gx,1,Q1,1)
  424.           RHS2(NT+2)  = ZERO
  425. *          RHS2(NT+1) = ZERO
  426. *          RHS2(NT+2) = ZERO
  427. *          WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)
  428. *          STOP
  429.           CALL COPYMAT(NT+2,Jaa,Ja2)
  430.           res=dnrm2(NT+2,RHS2,1)
  431. *       Appel a "DGESV" (LAPACK)
  432.           CALL DGESV(NT+2,1,Ja2,NT+2,IPIV3,RHS2,NT+2,INFO)
  433. c          si on souhaite verifier la convergence en correction dX
  434. *          res=dnrm2(NT+2,RHS2,1)
  435. c          WRITE(*,*) 'res=',res
  436. *          WRITE(*,*) 'RHS2=',(RHS2(IOU),IOU=1,NT+2)
  437.           DO M = 1,NT
  438.             Q1(M) = Q1(M)-RHS2(M)
  439.           ENDDO
  440.           OMEG = OMEG-RHS2(NT+1)
  441.           CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.false.)
  442.         ENDIF
  443.  
  444. *    ----- EVALUATION DU RESIDU -----------------------------------
  445.         CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  446.      &           MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,JANAL2,fvec,fref,NNM)
  447.  
  448. *          residu relatif (calcule autrement dans le cas NNM)
  449. c*        IF (TYPC(2:2).NE.'N') THEN
  450.         IF (.NOT. NNM) THEN
  451.           res=dnrm2(NT,fvec.prog,1)/fref0
  452.         ENDIF
  453.  
  454.         IF (IIMPI.GE.2) THEN
  455.           Q1NRM2=dnrm2(NT,Q1,1)
  456. c          IF (IIMPI.GE.333) THEN
  457. c           WRITE(IOIMP,111) 'Q^',its,(Q1(iou),iou=1,NT)
  458. c           WRITE(IOIMP,111) 'R^',its,(fvec.prog(iou),iou=1,NT)
  459. c          ENDIF
  460.           WRITE(IOIMP,999) its,OMEG,Q1NRM2,res
  461.         ENDIF
  462.  
  463. *    ----- TEST DE CONVERGENCE ------------------------------------
  464.         ITER = its
  465.         IF (res.LE.TOLMIN) THEN
  466.           check=.FALSE.
  467. c       Calcul de la Jacobienne pour la solution convergee
  468.           CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.not.NNM)
  469.           DO J = 1,IPALA(/1)
  470.             IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  471.               VCTCS(4) = VV
  472.               XPALA(J,4) = VV
  473.               CALL HBMZZF(NT,NDDL,OMEG,AiDi,Di,LAMBD,VCTCS,ZZ)
  474.             ENDIF
  475.           ENDDO
  476.           DO J=1,NT
  477.             DO I=1,NT
  478.               RX(I,J) = ZZ(I,J) - JAC(I,J)
  479.             ENDDO
  480.           ENDDO
  481.           KTQ = MTQ
  482.           KTEMP=MTEMP
  483.           RETURN
  484.         ELSE
  485.           check = .true.
  486.         ENDIF
  487.  
  488.       ENDDO
  489. *====== FIN DES ITERATIONS =============================================
  490.  
  491. * Menage :
  492.       SEGSUP,mwrkl,fvec
  493.  
  494. *=======================================================================
  495. *       Mise en forme des messages
  496. *=======================================================================
  497.  990  FORMAT(' + newt : ',I3,'  w=',F9.5,' |Q|=',F9.5,' |R|=',E10.3,
  498.      &         ' |ref|=',E10.3,' ds=',F9.5)
  499.  999  FORMAT(' + newt : ',I3,'  w=',F9.5,' |Q|=',F9.5,' |R|=',E10.3)
  500. cdebug
  501.  123  FORMAT('Ja(',I2,',1:6)=',6(E12.5,1X),'... F=',E12.5)
  502.  111  FORMAT(A,I2,'=',11(1X,F10.6))
  503.  
  504.       RETURN
  505.       END
  506.  
  507.  
  508.  

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