1 : C OTTOB3 SOURCE CHAT 05/01/13 02:07:00 5004 2 : SUBROUTINE OTTOB3(A,B,LADIM,XX,KERRE) 3 : C====================================================================== 4 : C 5 : C ENTREES 6 : C ------- 7 : C A(3,3) = MATRICE SYMETRIQUE 8 : C B(3,3) = MATRICE SYMETRIQUE 9 : C LADIM = 2 OU 3 SI 2 ON NE S OCCUPE QUE DE A(2,2) 10 : C SI 3 DE A(3,3) 11 : C SORTIES 12 : C ------- 13 : C 14 : C ON SORT XX(1) < XX(2) < XX(3) 15 : C 16 : C=============================================================== 17 : IMPLICIT INTEGER(I-N) 18 : IMPLICIT REAL*8 (A-H,O-Z) 19 : -INC CCOPTIO 20 : DIMENSION A(3,3),B(3,3),XX(3) 21 : * 22 : KERRE=0 23 : * 24 : IF(LADIM.EQ.3) GOTO 700 25 : * 26 : * CAS DIMENSION 2 27 : * 28 : * A0 + A1*X + A2* X2 29 : * 30 : A0=A(1,1)*A(2,2)-A(1,2)*A(1,2) 31 : A1=A(1,1)*B(2,2)+A(2,2)*B(1,1)-2.D0*A(1,2)*B(1,2) 32 : A2=B(1,1)*B(2,2)-B(1,2)*B(1,2) 33 : 34 : * 35 : IF(IIMPI.EQ.42) THEN 36 : WRITE(IOIMP,70615) A0,A1,A2 37 : 70615 FORMAT(2X,'OTTOB3 - A0 A1 A2 ',3(1X,1PE12.5)/) 38 : ENDIF 39 : * 40 : IF(A2.EQ.0.D0) THEN 41 : IF(A1.EQ.0.D0) THEN 42 : WRITE(IOIMP,76000) 43 : 76000 FORMAT(2X,'OTTOB3 - A2 ET A1 SONT NULS ' ) 44 : KERRE=70 45 : RETURN 46 : ENDIF 47 : XX1=-A0/A1 48 : XX2=1.D4 49 : ELSE IF(A1.EQ.0.D0) THEN 50 : DIS=-A0/A2 51 : IF(DIS.LT.0.D0) THEN 52 : WRITE(IOIMP,76001) DIS 53 : 76001 FORMAT(2X,'OTTOB3 - RAPPORT DIS NEGATIF = ',1PE12.5/) 54 : KERRE=70 55 : RETURN 56 : ENDIF 57 : XX1=-SQRT(DIS) 58 : XX2= SQRT(DIS) 59 : ELSE 60 : AUX=4.D0*A2*A0/(A1*A1) 61 : IF(ABS(AUX).LT.1.D-8) THEN 62 : XX1=-A0/A1 - A0*A0*A2/(A1**3) 63 : XX2=-A1/A2 -XX1 64 : ELSE 65 : DIS = A1*A1-4.D0*A2*A0 66 : * 67 : IF(DIS.LT.0.D0) THEN 68 : AUX=1.D-12*A1*A1 69 : IF(ABS(DIS).LT.AUX ) THEN 70 : DIS= MAX(DIS,0.D0) 71 : ELSE 72 : WRITE(IOIMP,76002) DIS , AUX 73 : 76002 FORMAT(2X,'OTTOB3 - DISCRIMINANT NEGATIF = ', 74 : & 1PE12.5,2X,' AUX=',1PE12.5/) 75 : KERRE=70 76 : RETURN 77 : ENDIF 78 : ENDIF 79 : DIS=SQRT(DIS) 80 : XX1=(-A1-DIS)/(2.D0*A2) 81 : XX2=(-A1+DIS)/(2.D0*A2) 82 : ENDIF 83 : ENDIF 84 : XX(1)=MIN(XX1,XX2) 85 : XX(2)=MAX(XX1,XX2) 86 : RETURN 87 : * 88 : * CAS DIMENSION 3 89 : * 90 : * A0 + A1*X + A2* X2 + A3* X3 91 : * 92 : 700 CONTINUE 93 : * 94 : A0=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(2,3)) 95 : & -A(1,2)*(A(3,3)*A(1,2)-A(1,3)*A(2,3)) 96 : & +A(1,3)*(A(1,2)*A(2,3)-A(2,2)*A(1,3)) 97 : * 98 : A1=A(1,1)*(A(3,3)*B(2,2)+A(2,2)*B(3,3)-2.D0*A(2,3)*B(2,3)) 99 : & -A(1,2)*(A(3,3)*B(1,2)+A(1,2)*B(3,3) 100 : & -A(1,3)*B(2,3)-A(2,3)*B(1,3)) 101 : & +A(1,3)*(A(2,3)*B(1,2)+A(1,2)*B(2,3) 102 : & -A(1,3)*B(2,2)-A(2,2)*B(1,3)) 103 : & +B(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(2,3)) 104 : & -B(1,2)*(A(3,3)*A(1,2)-A(2,3)*A(1,3)) 105 : & +B(1,3)*(A(1,2)*A(2,3)-A(2,2)*A(1,3)) 106 : * 107 : A2=B(1,1)*(A(3,3)*B(2,2)+A(2,2)*B(3,3)-2.D0*A(2,3)*B(2,3)) 108 : & -B(1,2)*(A(3,3)*B(1,2)+A(1,2)*B(3,3) 109 : & -A(1,3)*B(2,3)-A(2,3)*B(1,3)) 110 : & +B(1,3)*(A(2,3)*B(1,2)+A(1,2)*B(2,3) 111 : & -A(1,3)*B(2,2)-A(2,2)*B(1,3)) 112 : & +A(1,1)*(B(2,2)*B(3,3)-B(2,3)*B(2,3)) 113 : & -A(1,2)*(B(3,3)*B(1,2)-B(2,3)*B(1,3)) 114 : & +A(1,3)*(B(1,2)*B(2,3)-B(2,2)*B(1,3)) 115 : * 116 : A3=B(1,1)*(B(2,2)*B(3,3)-B(2,3)*B(2,3)) 117 : & -B(1,2)*(B(3,3)*B(1,2)-B(1,3)*B(2,3)) 118 : & +B(1,3)*(B(1,2)*B(2,3)-B(2,2)*B(1,3)) 119 : * 120 : IF(IIMPI.EQ.42) THEN 121 : WRITE(IOIMP,70616) A0,A1,A2,A3 122 : 70616 FORMAT(2X,'OTTOB3 - A0 A1 A2 A3 ',4(1X,1PE12.5)/) 123 : ENDIF 124 : ****** 125 : IF(ABS(A3).LT.1.D-6) THEN 126 : A3=0.D0 127 : ENDIF 128 : ******** 129 : 130 : * 131 : IF(A3.EQ.0.D0) THEN 132 : * 133 : IF(A2.EQ.0.D0) THEN 134 : IF(A1.EQ.0.D0) THEN 135 : WRITE(IOIMP,76003) 136 : 76003 FORMAT(2X,'OTTOB3 - A3, A2 ET A1 SONT NULS ' ) 137 : KERRE=70 138 : RETURN 139 : ENDIF 140 : XX(1)=-A0/A1 141 : XX(2)=1.D4 142 : XX(3)=1.D4 143 : ELSE 144 : DIS = A1*A1-4.D0*A2*A0 145 : *********** 146 : DIS = MAX ( DIS, 0.D0) 147 : *********** 148 : 149 : 150 : IF(DIS.LT.0.D0) THEN 151 : WRITE(IOIMP,76001) DIS 152 : KERRE=70 153 : RETURN 154 : ENDIF 155 : DIS=SQRT(DIS) 156 : XX(1)=(-A1-DIS)/(2.D0*A2) 157 : XX(2)=(-A1+DIS)/(2.D0*A2) 158 : XX(3)=1.D4 159 : ENDIF 160 : RETURN 161 : * 162 : ELSE 163 : * 164 : * MISE SOUS FORME POUR DEGRE3 165 : * 166 : AS0=A0/A3 167 : AS1=A1/A3 168 : AS2=A2/A3 169 : * 169 : CALL DEGRE3(AS0,AS1,AS2,XR1,XI1,XR2,XI2,XR3,XI3) 170 : * 171 : IF(IIMPI.EQ.42) THEN 172 : WRITE(IOIMP,76005) XR1,XI1,XR2,XI2,XR3,XI3 173 : 76005 FORMAT(2X,'OTTOB3 - XR1 = ',1PE12.5,2X,'XI1=',1PE12.5/ 174 : & 2X,' XR2 = ',1PE12.5,2X,'XI2=',1PE12.5/ 175 : & 2X,' XR3 = ',1PE12.5,2X,'XI3=',1PE12.5/) 176 : ENDIF 177 : * 178 : * xi1 est toujours nul 179 : * 180 : IF(XI2.NE.0.D0) THEN 181 : XX(1)=XR1 182 : XX(2)=1.D4 183 : XX(3)=1.D4 184 : ELSE 185 : IF (XR1.LE.XR2.AND.XR1.LE.XR3)THEN 186 : XX(1)=XR1 187 : XX(2)=XR2 188 : XX(3)=XR3 189 : IF (XR3.LE.XR2) THEN 190 : XX(2)=XR3 191 : XX(3)=XR2 192 : ENDIF 193 : * 194 : ELSE IF (XR2.LE.XR3.AND.XR2.LE.XR1)THEN 195 : XX(1)=XR2 196 : XX(2)=XR3 197 : XX(3)=XR1 198 : IF (XR1.LE.XR3) THEN 199 : XX(2)=XR1 200 : XX(3)=XR3 201 : ENDIF 202 : * 203 : ELSE IF (XR3.LE.XR1.AND.XR3.LE.XR2)THEN 204 : XX(1)=XR3 205 : XX(2)=XR1 206 : XX(3)=XR2 207 : IF (XR2.LE.XR1) THEN 208 : XX(2)=XR2 209 : XX(3)=XR1 210 : ENDIF 211 : ENDIF 212 : ENDIF 213 : ENDIF 214 : RETURN 215 : END 216 : 217 : 218 :
© Cast3M 2003 - All rights reserved.
Disclaimer