Télécharger inver1.eso

Retour à la liste

Numérotation des lignes :

  1. C INVER1 SOURCE CHAT 05/01/13 00:42:36 5004
  2. SUBROUTINE INVER1(TRAVAI,NZ,ICRIT,EPS)
  3. C
  4. C ====================================================================
  5. C SOUS-PROGRAMME APPELE PAR RAYE3 (RAYONNEMENT, oper. RAYN)
  6. C --> derive de inver.eso:
  7. C --> version esope
  8. C --> on travaille sur les matrices transposees
  9. C --> traitement des tableaux par colonne
  10. C RAPPEL:
  11. C INVERSION DE MATRICE PLEINE NON SYMETRIQUE PAR RESOLUTION
  12. C SUCCESSIVE DE NZ SYSTEMES LINEAIRES
  13. C A MATRICE (NZ*NZ) A INVERSER EN ENTREE, MATRICE INVERSEE EN SORTIE
  14. C ICRIT=1 SI MATRICE SINGULIERE, 0 SINON
  15. C BT TABLEAU DE REELS DE DIMENSION NZ*NZ
  16. C IS TABLEAU D'ENTIERS DE DIMENSION NZ
  17. C EPS PRECISION
  18. C
  19. C Le segment TRAVAI contenant A et BT est envoye actif par RAYE3
  20. C et retourné actif
  21. C
  22. C var locales : BTT transposee de BT dans segment DIVERS
  23. C ====================================================================
  24. C
  25. IMPLICIT INTEGER(I-N)
  26. IMPLICIT REAL*8(A-H,O-Z)
  27. -INC CCOPTIO
  28.  
  29. SEGMENT TRAVAI
  30. REAL *8 A(NEL,NEL), BT(NEL,NEL)
  31. INTEGER IS(NEL)
  32. ENDSEGMENT
  33.  
  34.  
  35. SEGMENT DIVERS
  36. REAL*8 BTT(NZ,NZ)
  37. ENDSEGMENT
  38. SEGINI DIVERS
  39. C
  40. C on peut stocker ABS(A) au lieu de le recalculer dans la boucle
  41. C 40 mais d'une part le gain en temps calcul n'est pas substantiel
  42. C et d'autre part cela implique une place memoire supplementaire
  43. C
  44. C DO 160 J=1,NZ
  45. C DO 160 I=1,NZ
  46. C ABA(I,J)=ABS(A(I,J))
  47. C 160 CONTINUE
  48. C
  49. C INITIALISATIONS
  50. C IS SUITE REPRESENTANT L'INDICE J DE X(J) SOLU CORRESPONDANT
  51. C A LA I EME COLONNE DE LA MATRICE A TRIANGULARISEE
  52. C
  53. DO 20 I=1,NZ
  54. DO 10 J=1,NZ
  55. BT(J,I)=0.D0
  56. 10 CONTINUE
  57. BT(I,I)=1.D0
  58. IS(I)=I
  59. 20 CONTINUE
  60. ICRIT=0
  61. C
  62. C 1- TRIANGULARISATION
  63. C
  64. NZM1=NZ-1
  65. DO 100 NR=1,NZM1
  66. C
  67. C CHOIX DU PIVOT
  68. C
  69. PIVOT=0.D0
  70. DO 40 K=NR,NZ
  71. DO 40 L=NR,NZ
  72. C ABSKL=ABA(L,K)
  73. ABSKL=ABS(A(L,K))
  74. IF(ABSKL.GT.PIVOT) THEN
  75. I=K
  76. J=L
  77. PIVOT=ABSKL
  78. ENDIF
  79. 40 CONTINUE
  80. C
  81. C LE PIVOT EST-IL NUL?
  82. C
  83. IF(PIVOT.LE.EPS) THEN
  84. ICRIT=1
  85. RETURN
  86. ENDIF
  87. C
  88. C CHANGEMENT DE LIGNE : PLACE LE PIVOT EN NR EME LIGNE
  89. C
  90. DO 50 L=1,NZ
  91. D=A(L,NR)
  92. A(L,NR)=A(L,I)
  93. A(L,I)=D
  94. 50 CONTINUE
  95. DO 60 L=1,NZ
  96. E=BT(L,NR)
  97. BT(L,NR)=BT(L,I)
  98. BT(L,I)=E
  99. 60 CONTINUE
  100. C
  101. C CHANGEMENT DE COLONNE : PLACE LE PIVOT EN R EME COLONNE
  102. C
  103. DO 70 M=1,NZ
  104. CC=A(NR,M)
  105. A(NR,M)=A(J,M)
  106. A(J,M)=CC
  107. 70 CONTINUE
  108. C
  109. C INDICE DES VARIABLES CORRESPONDANT A LA J EME ET A LA R EME COLON
  110. C
  111. ISR=IS(NR)
  112. IS(NR)=IS(J)
  113. IS(J)=ISR
  114. C
  115. C CALCUL DE LA NOUVELLE MATRICE A
  116. C
  117. NRP1=NR+1
  118. DO 90 I=NRP1,NZ
  119. IF(A(NR,I).NE.0.D0)THEN
  120. G=A(NR,I)/A(NR,NR)
  121. DO 80 J=1,NZ
  122. A(J,I)=A(J,I)-G*A(J,NR)
  123. BT(J,I)=BT(J,I)-G*BT(J,NR)
  124. 80 CONTINUE
  125. ENDIF
  126. 90 CONTINUE
  127. 100 CONTINUE
  128. C
  129. C 2- RESOLUTION DU SYSTEME TRIANGULARISE
  130. C
  131. * On introduit la transposée de B pour performance calcul fortran
  132. *
  133. DO 150 I=1,NZ
  134. DO 150 J=1,NZ
  135. BTT(J,I)=BT(I,J)
  136. 150 CONTINUE
  137.  
  138. DO 130 J=1,NZ
  139. BTT(NZ,J)=BTT(NZ,J)/A(NZ,NZ)
  140. DO 120 I= NZM1,1,-1
  141. F=0.D0
  142. I1=I+1
  143. DO 110 JJ=I1,NZ
  144. F=F-A(JJ,I)*BTT(JJ,J)
  145. 110 CONTINUE
  146. BTT(I,J)=(BTT(I,J)+F)/A(I,I)
  147. 120 CONTINUE
  148. 130 CONTINUE
  149. C
  150. DO 140 L=1,NZ
  151. IL=IS(L)
  152. DO 140 J=1,NZ
  153. A(J,IL)=BTT(L,J)
  154. 140 CONTINUE
  155.  
  156. SEGSUP DIVERS
  157. RETURN
  158. END
  159.  
  160.  
  161.  

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