1 : C INVER SOURCE CHAT 05/01/13 00:42:48 5004 2 : SUBROUTINE INVER(A,NZ,ICRIT,B,IS,EPS) 3 : C 4 : C ==================================================================== 5 : C SOUS-PROGRAMME FORTRAN APPELE PAR ELFINV 6 : C INVERSION DE MATRICE PAR RESOLUTION SUCCESSIVE DE NZ SYSTEMES 7 : C LINEAIRES 8 : C A MATRICE (NZ*NZ) A INVERSER EN ENTREE, MATRICE INVERSEE EN SORTIE 9 : C ICRIT=1 SI MATRICE SINGULIERE, 0 SINON 10 : C B TABLEAU DE REELS DE DIMENSION NZ*NZ 11 : C IS TABLEAU D'ENTIERS DE DIMENSION NZ 12 : C EPS PRECISION 13 : C ==================================================================== 14 : C 15 : IMPLICIT INTEGER(I-N) 16 : IMPLICIT REAL*8(A-H,O-Z) 17 : -INC CCOPTIO 18 : DIMENSION A(NZ,1),B(NZ,1) 19 : DIMENSION IS(1) 20 : C 21 : C INITIALISATIONS 22 : C IS SUITE REPRESENTANT L'INDICE J DE X(J) SOLU CORRESPONDANT 23 : C A LA I EME COLONNE DE LA MATRICE A TRIANGULARISEE 24 : C 25 : DO 20 I=1,NZ 26 : DO 10 J=1,NZ 27 : B(J,I)=0.D0 28 : 10 CONTINUE 29 : B(I,I)=1.D0 30 : IS(I)=I 31 : 20 CONTINUE 32 : ICRIT=0 33 : C 34 : C 1- TRIANGULARISATION 35 : C 36 : NZM1=NZ-1 37 : DO 100 NR=1,NZM1 38 : C 39 : C CHOIX DU PIVOT 40 : C 41 : PIVOT=0.D0 42 : DO 40 K=NR,NZ 43 : DO 40 L=NR,NZ 44 : ABSKL=ABS(A(K,L)) 45 : IF(ABSKL.GT.PIVOT) THEN 46 : I=K 47 : J=L 48 : PIVOT=ABSKL 49 : ENDIF 50 : 40 CONTINUE 51 : C 52 : C LE PIVOT EST-IL NUL? 53 : C 54 : IF(PIVOT.LE.EPS) THEN 55 : ICRIT=1 56 : RETURN 57 : ENDIF 58 : C 59 : C CHANGEMENT DE LIGNE : PLACE LE PIVOT EN NR EME LIGNE 60 : C 61 : DO 50 L=1,NZ 62 : D =A(NR,L) 63 : A(NR,L)=A(I,L) 64 : A(I,L)=D 65 : 50 CONTINUE 66 : DO 60 L=1,NZ 67 : E=B(NR,L) 68 : B(NR,L)=B(I,L) 69 : B(I,L)=E 70 : 60 CONTINUE 71 : C 72 : C CHANGEMENT DE COLONNE : PLACE LE PIVOT EN R EME COLONNE 73 : C 74 : DO 70 M=1,NZ 75 : CC =A(M,NR) 76 : A(M,NR)=A(M,J) 77 : A(M,J)=CC 78 : 70 CONTINUE 79 : C 80 : C INDICE DES VARIABLES CORRESPONDANT A LA J EME ET A LA R EME COLON 81 : C 82 : ISR=IS(NR) 83 : IS(NR)=IS(J) 84 : IS(J)=ISR 85 : C 86 : C CALCUL DE LA NOUVELLE MATRICE A 87 : C 88 : NRP1=NR+1 89 : DO 90 I=NRP1,NZ 90 : IF(A(I,NR).NE.0.D0)THEN 91 : G=A(I,NR)/A(NR,NR) 92 : DO 80 J=1,NZ 93 : A(I,J)=A(I,J)-G*A(NR,J) 94 : B(I,J)=B(I,J)-G*B(NR,J) 95 : 80 CONTINUE 96 : ENDIF 97 : 90 CONTINUE 98 : 100 CONTINUE 99 : C 100 : C 2- RESOLUTION DU SYSTEME TRIANGULARISE 101 : C 102 : DO 130 J=1,NZ 103 : B(NZ,J)=B(NZ,J)/A(NZ,NZ) 104 : DO 120 I= NZM1,1,-1 105 : F=0.D0 106 : I1=I+1 107 : DO 110 JJ=I1,NZ 108 : F=F-A(I,JJ)*B(JJ,J) 109 : 110 CONTINUE 110 : B(I,J)=(B(I,J)+F)/A(I,I) 111 : 120 CONTINUE 112 : 130 CONTINUE 113 : C 114 : DO 140 L=1,NZ 115 : IL=IS(L) 116 : DO 140 J=1,NZ 117 : A(IL,J)=B(L,J) 118 : 140 CONTINUE 119 : RETURN 120 : END 121 :
© Cast3M 2003 - All rights reserved.
Disclaimer