Télécharger la plaquette Cast3M


fr en

Télécharger subsolve.eso

Retour à la liste

Numérotation des lignes : 

subsolve

  1. C SUBSOLVE  SOURCE    FD218221  25/09/03    21:15:06     12351          
  2.       SUBROUTINE SUBSOLVE(M,N,MN,EPSIMIN,LOW,UPP,ALFA,BETA,P0,Q0,P,Q,
  3.      &                    A0,A,B,C,D,
  4.      &                    XMMA,YMMA,ZMMA,LAMMA,XSIMMA,ETAMMA,MUMMA,
  5.      &                    ZETMMA,SMMA)
  6.  
  7. C     Typages implicites habituels
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10.  
  11. C ------------------------------------------------------------------------------------------------
  12. C Ce code est une adaptation en FORTRAN de l'algorithme de la MMA (Methode des Asymptotes Mobiles)
  13. C propose initialement en langage Matlab par K. Svarberg
  14. C Les sources sont disponibles ici : https://www.smoptit.se/
  15. C ------------------------------------------------------------------------------------------------
  16. C    Cette subroutine resoud le sous probleme convexe suivant :
  17. C
  18. C    Minimiser :   SUM[ p0_j/(upp_j-x_j) + q0_j/(x_j-low_j) ] + a_0*z +
  19. C                + SUM[ c_i*y_i + 0.5*d_i*(y_i)^2 ]
  20. C
  21. C     Soumis a :   SUM[ p_ij/(upp_j-x_j) + q_ij/(x_j-low_j) ] - a_i*z - y_i <= b_i,  i = 1,...,m
  22. C                  alfa_j <= x_j <= beta_j,  j = 1,...,n
  23. C                  z >= 0,   y_i >= 0,       i = 1,...,m
  24. C
  25. C  Entrees :
  26. C  ---------
  27. C   m     = Nombre de contraintes
  28. C   n     = Nombre de variables x_j
  29. C   mn    = Min (m,n) (pour dimensionner les matrices AA, BB et SOLUT)
  30. C epsimin = Petit nombre positif afin de garantir la stabilite
  31. C  low    = Vecteur colonne (n x 1) des bornes inferieures
  32. C  upp    = Vecteur colonne (n x 1) des bornes superieures
  33. C  alfa   = Vecteur colonne (n x 1) des asymptotes inferieures
  34. C  beta   = Vecteur colonne (n x 1) des asymptotes superieures
  35. C  p0     = Vecteur colonne (n x 1) des coefficients pour les bornes inferieures
  36. C  q0     = Vecteur colonne (n x 1) des coefficients pour les bornes superieures
  37. C  P      = Matrice (m x n) des coefficients pour les bornes inferieures des contraintes
  38. C  Q      = Matrice (m x n) des coefficients pour les bornes superieures des contraintes
  39. C  a0     = Constante a_0
  40. C  a      = Vecteur colonne (m x 1) des constantes a_i
  41. C  b      = Vecteur colonne (m x 1) des constantes b_i
  42. C  c      = Vecteur colonne (m x 1) des constantes c_i
  43. C  d      = Vecteur colonne (m x 1) des constantes d_i
  44. C
  45. C  Sorties :
  46. C  ---------
  47. C  xmma  = Vecteur colonne (n x 1) des variables x_j optimisees
  48. C  ymma  = Vecteur colonne (m x 1) des variables y_j optimisees
  49. C  zmma  = Valeur de la variable z optimisee
  50. C  lamma = Vecteur colonne (m x 1) des multiplicateurs de Lagrange
  51. C          des contraintes generales
  52. C xsimma = Vecteur colonne (n x 1) des multiplicateurs de Lagrange
  53. C          des contraintes alfa_j - x_j <= 0
  54. C etamma = Vecteur colonne (n x 1) des multiplicateurs de Lagrange
  55. C          des contraintes x_j - beta_j <= 0
  56. C  mumma = Vecteur colonne (m x 1) des multiplicateurs de Lagrange
  57. C          des contraintes de positivite -y_i <= 0
  58. C zetmma = Multiplicateurs de Lagrange de la contrainte de positivite -z <= 0
  59. C   smma = Vecteur colonne (m x 1) des variables d'ecart des contraintes generales
  60. C ------------------------------------------------------------------------------------------------
  61.  
  62. C     Arguments de sortie
  63.       REAL*8 XMMA(N,1),YMMA(M,1),ZMMA
  64.       REAL*8 LAMMA(M,1),XSIMMA(N,1),ETAMMA(N,1),MUMMA(M,1),ZETMMA
  65.       REAL*8 SMMA(M,1)
  66. C     Arguments d'entree
  67.       REAL*8 EPSIMIN,LOW(N,1),UPP(N,1),ALFA(N,1),BETA(N,1)
  68.       REAL*8 P0(N,1),Q0(N,1),P(M,N),Q(M,N)
  69.       REAL*8 A0,A(M,1),B(M,1),C(M,1),D(M,1)
  70. C     Arguments internes
  71.       INTEGER ITTT,ITTO
  72.       INTEGER IPIV(MN+1)
  73.       REAL*8 EPSI,RESIDUNORM,RESIDUMAX
  74.       REAL*8 EEN(N,1),EEM(M,1),EPSVECM(M,1),EPSVECN(N,1)
  75.       REAL*8 X(N,1),Y(M,1),Z,DX(N,1),DY(M,1),DZ
  76.       REAL*8 XOLD(N,1),YOLD(M,1),ZOLD,DELX(N,1),DELY(M,1),DELZ
  77.       REAL*8 LAM(M,1),MU(M,1),XSI(N,1),ETA(N,1),S(M,1)
  78.       REAL*8 UX1(N,1),UX2(N,1),UX3(N,1)
  79.       REAL*8 XL1(N,1),XL2(N,1),XL3(N,1)
  80.       REAL*8 UXINV1(N,1),XLINV1(N,1),UXINV2(N,1),XLINV2(N,1)
  81.       REAL*8 DELLAM(M,1),DIAGX(N,1),DIAGY(M,1)
  82.       REAL*8 DIAGXINV(N,1),DIAGYINV(M,1)
  83.       REAL*8 DIAGLAM(M,1),DIAGLAMYI(M,1),PLAM(N,1),QLAM(N,1),GVEC(M,1)
  84.       REAL*8 DPSIDX(N,1),RES(M,1),REX(N,1),REY(M,1),REZ
  85.       REAL*8 RELAM(M,1),REXSI(N,1),REETA(N,1),REMU(M,1),REZET,RESINEW
  86.       REAL*8 BLAM(M,1)
  87.       REAL*8 GG(M,N),BX(N,1),BZ,AXX(N,N),AXZ(N,1),AZZ
  88.       REAL*8 TEMP(M,N),TEMP2(N,M)
  89.       REAL*8 ALAM(M,M),DLAM(M,1),DIAGLAMYIINV(M,1),DELLAMYI(M,1)
  90.       REAL*8 DXSI(N,1),DETA(N,1),DMU(M,1),ZET,DZET,DS(M,1)
  91.       REAL*8 STMXX,STEPALFA(N,1),STMALFA,STEPBETA(N,1),STMBETA
  92.       REAL*8 STMALBE,STMALBEXX,STMINV,STEG
  93.       REAL*8 LAMOLD(M,1),XSIOLD(N,1),ETAOLD(N,1)
  94.       REAL*8 MUOLD(M,1),ZETOLD,SOLD(M,1)
  95.       REAL*8 AA(MN+1,MN+1),BB(MN+1,1),SOLUT(MN+1,1)
  96.       REAL*8 XX(4*M+2*N+2,1),DXX(4*M+2*N+2,1)
  97.       REAL*8 STEPXX(4*M+2*N+2,1)
  98.       REAL*8 RESIDU1(N+M+1),RESIDU2(3*M+2*N+1),RESIDU(4*M+3*N+2)
  99.  
  100. C     Initialisation
  101.       CALL UN(EEN,N,1)
  102.       CALL UN(EEM,M,1)
  103.       EPSI = 1.0D0
  104.       EPSVECN = EPSI*EEN
  105.       EPSVECM = EPSI*EEM
  106.       X = 0.5D0*(ALFA+BETA)
  107.       Y = EEM
  108.       Z = 1
  109.       LAM = EEM
  110.       XSI = EEN/(X-ALFA)
  111.       XSI = MAX(XSI,EEN)
  112.       ETA = EEN/(BETA-X)
  113.       ETA = MAX(ETA,EEN)
  114.       MU  = MAX(EEM,0.5D0*C)
  115.       ZET = 1
  116.       S = EEM
  117.       ITERA = 0
  118.  
  119. C     Boucle de stabilite numerique
  120.       DO WHILE (EPSI.GT.EPSIMIN)
  121.           EPSVECN = EPSI*EEN
  122.           EPSVECM = EPSI*EEM
  123.           UX1 = UPP-X
  124.           XL1 = X-LOW
  125.           UX2 = UX1*UX1
  126.           XL2 = XL1*XL1
  127.           UXINV1 = EEN/UX1
  128.           XLINV1 = EEN/XL1
  129.           PLAM = P0 + MATMUL(TRANSPOSE(P),LAM) 
  130.           QLAM = Q0 + MATMUL(TRANSPOSE(Q),LAM) 
  131.           GVEC = MATMUL(P,UXINV1) + MATMUL(Q,XLINV1)
  132.           DPSIDX = PLAM/UX2 - QLAM/XL2
  133.           REX = DPSIDX - XSI + ETA
  134.           REY = C + D*Y - MU - LAM
  135.           REZ = A0 - ZET - MAXVAL(MATMUL(TRANSPOSE(A),LAM))
  136.           RELAM = GVEC - A*Z - Y + S - B
  137.           REXSI = XSI*(X-ALFA) - EPSVECN
  138.           REETA = ETA*(BETA-X) - EPSVECN
  139.           REMU = MU*Y - EPSVECM
  140.           REZET = ZET*Z - EPSI
  141.           RES = LAM*S - EPSVECM
  142. C         RESIDU1=[REX(N), REY(M), REZ]
  143.           NR1=N+M+1
  144.           RESIDU1(1:N)=REX(1:N,1)
  145.           RESIDU1(N+1:N+M)=REY(1:M,1)
  146.           RESIDU1(NR1)=REZ
  147. C         RESIDU2=[RELAM(M), REXSI(N), REETA(N), REMU(M), REZET, RES(M)]
  148.           NR2=2*N+3*M+1
  149.           RESIDU2(1:M)=RELAM(1:M,1)
  150.           RESIDU2(M+1:M+N)=REXSI(1:N,1)
  151.           RESIDU2(M+N+1:M+2*N)=REETA(1:N,1)
  152.           RESIDU2(M+2*N+1:2*M+2*N)=REMU(1:M,1)
  153.           RESIDU2(2*M+2*N+1)=REZET
  154.           RESIDU2(2*M+2*N+2:NR2)=RES(1:M,1)
  155. C         RESIDU=[RESIDU1, RESIDU2]
  156.           NR=NR1+NR2
  157.           RESIDU(1:NR1)=RESIDU1
  158.           RESIDU(NR1+1:NR)=RESIDU2
  159. C         Normes 2 et infinie du residu
  160.           RESIDUNORM = SQRT(DDOT(NR,RESIDU,1,RESIDU,1))
  161.           RESIDUMAX = VECMAX(ABS(RESIDU),NR)
  162.           ITTT = 0
  163.  
  164. C         Boucle d'optimisation
  165.           DO WHILE ((RESIDUMAX.GT.0.9D0*EPSI).AND.(ITTT.LT.200))
  166.               ITTT = ITTT + 1
  167.               ITERA = ITERA + 1
  168.               UX1 = UPP - X
  169.               XL1 = X - LOW
  170.               UX2 = UX1*UX1
  171.               XL2 = XL1*XL1
  172.               UX3 = UX1*UX2
  173.               XL3 = XL1*XL2
  174.               UXINV1 = EEN/UX1
  175.               XLINV1 = EEN/XL1
  176.               UXINV2 = EEN/UX2
  177.               XLINV2 = EEN/XL2
  178.               PLAM = P0 + MATMUL(TRANSPOSE(P),LAM)
  179.               QLAM = Q0 + MATMUL(TRANSPOSE(Q),LAM)
  180.               GVEC = MATMUL(P,UXINV1) + MATMUL(Q,XLINV1)
  181.               CALL MULMV(TEMP,P,UXINV2,M,N)
  182.               GG = TEMP
  183.               CALL MULMV(TEMP,Q,XLINV2,M,N)
  184.               GG = GG - TEMP
  185.               DPSIDX = PLAM/UX2 - QLAM/XL2
  186.               DELX = DPSIDX - EPSVECN/(X-ALFA) + EPSVECN/(BETA-X)
  187.               DELY = C + D*Y - LAM - EPSVECM/Y
  188.               DELZ = A0 - MAXVAL(MATMUL(TRANSPOSE(A),LAM)) - EPSI/Z
  189.               DELLAM = GVEC - A*Z - Y - B + EPSVECM/LAM
  190.               DIAGX = PLAM/UX3 + QLAM/XL3
  191.               DIAGX = 2*DIAGX + XSI/(X-ALFA) + ETA/(BETA-X)
  192.               DIAGXINV = EEN/DIAGX
  193.               DIAGY = D + MU/Y
  194.               DIAGYINV = EEM/DIAGY
  195.               DIAGLAM = S/LAM
  196.               DIAGLAMYI = DIAGLAM + DIAGYINV
  197.  
  198. C             Resolution du systeme lineaire : AA.SOLUT = BB
  199.  
  200. C             Cas ou il y moins de contraintes que de variables primales
  201. C             AA est une matrice (m+1 x m+1)
  202. C             BB est un vecteur  (m+1 x 1)
  203.               IF (M.LT.N) THEN
  204. C                 Assemblage de la matrice AA = [ALAM    A
  205. C                                                A'   -ZET/Z]
  206.                   CALL MULMV(TEMP,GG,DIAGXINV,M,N)
  207.                   ALAM=MATMUL(TEMP,TRANSPOSE(GG))
  208.                   DO I=1,M
  209.                       ALAM(I,I)=ALAM(I,I)+DIAGLAMYI(I,1)
  210.                   ENDDO
  211. C                 Ajout du bloc ALAM
  212.                   DO I=1,M
  213.                       DO J=1,M
  214.                           AA(I,J)=ALAM(I,J)
  215.                       ENDDO
  216.                   ENDDO
  217. C                 Ajout du bloc A
  218.                   DO I=1,M
  219.                       AA(I,M+1)=A(I,1)
  220.                   ENDDO
  221. C                 Ajout du bloc A'
  222.                   DO I=1,M
  223.                       AA(M+1,I)=A(I,1)
  224.                   ENDDO
  225. C                 Ajout du bloc -ZET/Z
  226.                   AA(M+1,M+1)=-ZET/Z
  227. C                 Assemblage du vecteur BB = [BLAM
  228. C                                             DELZ]
  229.                   BLAM=DELLAM + DELY/DIAGY - MATMUL(GG,(DELX/DIAGX))
  230.                   DO I=1,M
  231.                       BB(I,1)=BLAM(I,1)
  232.                   ENDDO
  233.                   BB(M+1,1)=DELZ
  234. C                 Appel au solveur de LAPACK
  235.                   CALL DGESV(M+1,1,AA,M+1,IPIV,BB,M+1,INFO)
  236. C                 Sortie en cas d'erreur
  237.                   IF (INFO.NE.0) THEN
  238.                       PRINT*,'ERROR IN DGESV SOLVER, INFO=',INFO
  239.                       PRINT*,'CHECK THE LINEAR SYSTEM'
  240.                       CALL ERREUR(26)
  241.                       RETURN
  242.                   ENDIF
  243.                   SOLUT=BB
  244.                   DLAM(1:M,1) = SOLUT(1:M,1)
  245.                   DZ = SOLUT(M+1,1)
  246.                   DX = -DELX/DIAGX - MATMUL(TRANSPOSE(GG),DLAM)/DIAGX
  247.  
  248. C             Cas ou il y plus de contraintes que de variables primales
  249. C             AA est une matrice (n+1 x n+1)
  250. C             BB est un vecteur  (n+1 x 1)
  251.               ELSE
  252. C                 Assemblage de la matrice AA = [AXX  AXZ
  253. C                                                AXZ' AZZ]
  254.                   DIAGLAMYIINV=EEM/DIAGLAMYI
  255.                   CALL MULMV(TEMP2,TRANSPOSE(GG),DIAGLAMYIINV,N,M)
  256.                   AXX=MATMUL(TEMP2,GG)
  257.                   DO I=1,N
  258.                       AXX(I,I)=AXX(I,I)+DIAGX(I,1)
  259.                   ENDDO
  260.                   AZZ=ZET/Z+MAXVAL(MATMUL(TRANSPOSE(A),(A/DIAGLAMYI)))
  261.                   AXZ=MATMUL(TRANSPOSE(-GG),(A/DIAGLAMYI))
  262. C                 Ajout du bloc AXX
  263.                   DO I=1,N
  264.                       DO J=1,N
  265.                           AA(I,J)=AXX(I,J)
  266.                       ENDDO
  267.                   ENDDO
  268. C                 Ajout du bloc AXZ
  269.                   DO I=1,N
  270.                       AA(I,N+1)=AXZ(I,1)
  271.                   ENDDO
  272. C                 Ajout du bloc AXZ'
  273.                   DO I=1,N
  274.                       AA(N+1,I)=AXZ(I,1)
  275.                   ENDDO
  276. C                 Ajout du bloc AZZ
  277.                   AA(N+1,N+1)=AZZ
  278. C                 Assemblage du vecteur BB = [-BX
  279. C                                             -BZ]
  280.                   DELLAMYI=DELLAM + DELY/DIAGY
  281.                   BX=DELX+MATMUL(TRANSPOSE(GG),(DELLAMYI/DIAGLAMYI))
  282.                   BZ=DELZ-MAXVAL(MATMUL(TRANSPOSE(A),
  283.      &                                  (DELLAMYI/DIAGLAMYI)))
  284.                   DO I=1,N
  285.                       BB(I,1)=-BX(I,1)
  286.                   ENDDO
  287.                   BB(N+1,1)=-BZ
  288. C                 Appel au solveur de LAPACK
  289.                   CALL DGESV(N+1,1,AA,N+1,IPIV,BB,N+1,INFO)
  290. C                 Sortie en cas d'erreur
  291.                   IF (INFO.NE.0) THEN
  292.                       PRINT*,'ERROR IN DGESV SOLVER, INFO=',INFO
  293.                       PRINT*,'CHECK THE LINEAR SYSTEM'
  294.                       CALL ERREUR(26)
  295.                       RETURN
  296.                   ENDIF
  297.                   SOLUT = BB
  298.                   DX(1:N,1) = SOLUT(1:N,1)
  299.                   DZ = SOLUT(N+1,1)
  300.                   DLAM =MATMUL(GG,DX)/DIAGLAMYI-DZ*(A/DIAGLAMYI)
  301.      &                                         +DELLAMYI/DIAGLAMYI
  302.               END IF
  303.  
  304.               DY = -DELY/DIAGY + DLAM/DIAGY
  305.               DXSI = -XSI + EPSVECN/(X-ALFA) - (XSI*DX)/(X-ALFA)
  306.               DETA = -ETA + EPSVECN/(BETA-X) + (ETA*DX)/(BETA-X)
  307.               DMU  = -MU  + EPSVECM/Y        - (MU*DY)/Y
  308.               DZET = -ZET + EPSI/Z           - ZET*DZ/Z
  309.               DS   = -S   + EPSVECM/LAM      - (S*DLAM)/LAM
  310.  
  311.               NX=4*M+2*N+2
  312. C             XX=[ Y, Z, LAM, XSI, ETA, MU, ZET, S]
  313.               XX(1:M,1)=Y(1:M,1)
  314.               XX(M+1,1)=Z
  315.               XX(M+2:2*M+1,1)=LAM(1:M,1)
  316.               XX(2*M+2:2*M+N+1,1)=XSI(1:N,1)
  317.               XX(2*M+N+2:2*M+2*N+1,1)=ETA(1:N,1)
  318.               XX(2*M+2*N+2:3*M+2*N+1,1)=MU(1:M,1)
  319.               XX(3*M+2*N+2,1)=ZET
  320.               XX(3*M+2*N+3:4*M+2*N+2,1)=S(1:M,1)
  321. C             DXX=[DY,DZ,DLAM,DXSI,DETA,DMU,DZET,DS]
  322.               DXX(1:M,1)=DY(1:M,1)
  323.               DXX(M+1,1)=DZ
  324.               DXX(M+2:2*M+1,1)=DLAM(1:M,1)
  325.               DXX(2*M+2:2*M+N+1,1)=DXSI(1:N,1)
  326.               DXX(2*M+N+2:2*M+2*N+1,1)=DETA(1:N,1)
  327.               DXX(2*M+2*N+2:3*M+2*N+1,1)=DMU(1:M,1)
  328.               DXX(3*M+2*N+2,1)=DZET
  329.               DXX(3*M+2*N+3:4*M+2*N+2,1)=DS(1:M,1)
  330.  
  331. C             Calcul de la longueur du pas
  332.               STEPXX = -1.01D0*DXX/XX
  333.               STMXX  = MAXVAL(STEPXX)
  334.               STEPALFA = -1.01D0*DX/(X-ALFA)
  335.               STMALFA = MAXVAL(STEPALFA)
  336.               STEPBETA = 1.01D0*DX/(BETA-X)
  337.               STMBETA = MAXVAL(STEPBETA)
  338.               STMALBE  = MAX(STMALFA,STMBETA)
  339.               STMALBEXX = MAX(STMALBE,STMXX)
  340.               STMINV = MAX(STMALBEXX,1.D0)
  341.               STEG = 1.D0/STMINV
  342.  
  343. C             Mise a jour des variables
  344.               XOLD = X
  345.               YOLD = Y
  346.               ZOLD = Z
  347.               LAMOLD = LAM
  348.               XSIOLD = XSI
  349.               ETAOLD = ETA
  350.               MUOLD  = MU
  351.               ZETOLD = ZET
  352.               SOLD = S
  353.  
  354.               ITTO = 0
  355.               RESINEW = 2*RESIDUNORM
  356.  
  357.               DO WHILE ((RESINEW.GT.RESIDUNORM).AND.(ITTO.LT.50))
  358.                   ITTO = ITTO + 1;
  359.                   X = XOLD + STEG*DX
  360.                   Y = YOLD + STEG*DY
  361.                   Z = ZOLD + STEG*DZ
  362.                   LAM = LAMOLD + STEG*DLAM
  363.                   XSI = XSIOLD + STEG*DXSI
  364.                   ETA = ETAOLD + STEG*DETA
  365.                   MU = MUOLD  + STEG*DMU
  366.                   ZET = ZETOLD + STEG*DZET
  367.                   S = SOLD + STEG*DS
  368.                   UX1 = UPP - X
  369.                   XL1 = X - LOW
  370.                   UX2 = UX1*UX1
  371.                   XL2 = XL1*XL1
  372.                   UXINV1 = EEN/UX1
  373.                   XLINV1 = EEN/XL1
  374.                   PLAM = P0 + MATMUL(TRANSPOSE(P),LAM) 
  375.                   QLAM = Q0 + MATMUL(TRANSPOSE(Q),LAM) 
  376.                   GVEC = MATMUL(P,UXINV1) + MATMUL(Q,XLINV1)
  377.                   DPSIDX = PLAM/UX2 - QLAM/XL2 
  378.                   REX = DPSIDX - XSI + ETA
  379.                   REY = C + D*Y - MU - LAM
  380.                   REZ = A0 - ZET - MAXVAL(MATMUL(TRANSPOSE(A),LAM))
  381.                   RELAM = GVEC - A*Z - Y + S - B
  382.                   REXSI = XSI*(X-ALFA) - EPSVECN
  383.                   REETA = ETA*(BETA-X) - EPSVECN
  384.                   REMU = MU*Y - EPSVECM
  385.                   REZET = ZET*Z - EPSI
  386.                   RES = LAM*S - EPSVECM
  387. C                 RESIDU1=[REX(N), REY(M), REZ]
  388.                   NR1=N+M+1
  389.                   RESIDU1(1:N)=REX(1:N,1)
  390.                   RESIDU1(N+1:N+M)=REY(1:M,1)
  391.                   RESIDU1(NR1)=REZ
  392. C                 RESIDU2=[RELAM(M), REXSI(N), REETA(N), REMU(M), REZET, RES(M)]
  393.                   NR2=2*N+3*M+1
  394.                   RESIDU2(1:M)=RELAM(1:M,1)
  395.                   RESIDU2(M+1:M+N)=REXSI(1:N,1)
  396.                   RESIDU2(M+N+1:M+2*N)=REETA(1:N,1)
  397.                   RESIDU2(M+2*N+1:2*M+2*N)=REMU(1:M,1)
  398.                   RESIDU2(2*M+2*N+1)=REZET
  399.                   RESIDU2(2*M+2*N+2:NR2)=RES(1:M,1)
  400. C                 RESIDU=[RESIDU1, RESIDU2]
  401.                   NR=NR1+NR2
  402.                   RESIDU(1:NR1)=RESIDU1
  403.                   RESIDU(NR1+1:NR)=RESIDU2
  404. C                 Normes 2 du residu
  405.                   NR=4*M+3*N+2
  406.                   RESINEW = SQRT(DDOT(NR,RESIDU,1,RESIDU,1))
  407.                   STEG = STEG/2.0D0
  408.               END DO
  409.               RESIDUNORM = RESINEW
  410. C             Norme infini du residu
  411.               NR=4*M+3*N+2
  412.               RESIDUMAX = VECMAX(ABS(RESIDU),NR)
  413.               STEG = 2.0D0*STEG
  414.           END DO
  415.           EPSI = 0.1D0*EPSI
  416.       END DO
  417.  
  418. C     Fin du programme
  419.       XMMA = X
  420.       YMMA = Y
  421.       ZMMA = Z
  422.       LAMMA = LAM
  423.       XSIMMA = XSI
  424.       ETAMMA = ETA
  425.       MUMMA = MU
  426.       ZETMMA = ZET
  427.       SMMA = S
  428.       RETURN
  429.       END
  430.  
  431.  
  432.  

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