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 :
© Cast3M 2003 - All rights reserved.
Disclaimer