C POLA2 SOURCE CB215821 26/08/24 21:17:45 12622 SUBROUTINE POLA2(F,R,U,N) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO DIMENSION F(*),R(*),U(*) DIMENSION C(3,3),C2(3,3),UI(9) * JEBOUC=0 IIMPI0=IIMPI 2020 JEBOUC=JEBOUC+1 KERRE=0 * * CALCUL DE C = U2 = FT*F * DO 77887 I=1,N DO 77886 J=1,N C(I,J)=0.D0 DO 1 K=1,N KI= (K-1)*N+I KJ= (K-1)*N+J C(I,J)=C(I,J)+F(KI)*F(KJ) 1 CONTINUE 77886 CONTINUE 77887 CONTINUE * IF(IIMPI.EQ.199) THEN WRITE(6,77771) N 77771 FORMAT(2X,'POLA2 - N=',I3/) N2=N*N WRITE(6,77772) (F(I),I=1,N2) 77772 FORMAT(2X,'F '/(3(1X,1PE12.5))) ENDIF * IF(N.EQ.2) THEN * * CAS 2D * TRC=C(1,1)+C(2,2) RDTC=SQRT(C(1,1)*C(2,2)-C(1,2)*C(2,1)) DTU=RDTC TRU=SQRT(TRC+2.D0*RDTC) * IF(IIMPI.EQ.199) THEN WRITE(6,77773) TRU,DTU 77773 FORMAT(2X,'POLA2 - TRU= ',1PE12.5,2X,'DTU= ',1PE12.5/) ENDIF IF(TRU.EQ.0.D0.OR.DTU.EQ.0.D0) THEN WRITE(6,77883) TRU,DTU 77883 FORMAT(2X,'TRU=',1PE12.5,2X,'DTU=',1PE12.5/) KERRE=26 GO TO 2021 ENDIF * U(1)=(DTU+C(1,1))/TRU U(2)=C(1,2)/TRU U(3)=C(2,1)/TRU U(4)=(DTU+C(2,2))/TRU * UI(1)=(TRU-U(1))/DTU UI(2)=-U(2)/DTU UI(3)=-U(3)/DTU UI(4)=(TRU-U(4))/DTU * ELSE IF(N.EQ.3) THEN * * CAS 3D * * CAS FAUX 3D * IF(C(1,3).EQ.0.D0.AND.C(2,3).EQ.0.D0.AND. & C(3,1).EQ.0.D0.AND.C(3,2).EQ.0.D0) THEN TRC=C(1,1)+C(2,2) RDTC=SQRT(C(1,1)*C(2,2)-C(1,2)*C(2,1)) DTU=RDTC TRU=SQRT(TRC+2.D0*RDTC) * IF(TRU.EQ.0.D0.OR.DTU.EQ.0.D0.OR.F(9).EQ.0.D0) THEN KERRE=26 GO TO 2021 ENDIF * U(1)=(DTU+C(1,1))/TRU U(2)=C(1,2)/TRU U(3)=0.D0 U(4)=C(2,1)/TRU U(5)=(DTU+C(2,2))/TRU U(6)=0.D0 U(7)=0.D0 U(8)=0.D0 U(9)=F(9) * UI(1)=(TRU-U(1))/DTU UI(2)=-U(2)/DTU UI(3)=0.D0 UI(4)=-U(4)/DTU UI(5)=(TRU-U(5))/DTU UI(6)=0.D0 UI(7)=0.D0 UI(8)=0.D0 UI(9)=1.D0/F(9) * ELSE * * CAS VRAI 3D * DO 77889 I=1,N DO 77888 J=1,N C2(I,J)=0.D0 DO 3 K=1,N C2(I,J)=C2(I,J)+C(I,K)*C(K,J) 3 CONTINUE 77888 CONTINUE 77889 CONTINUE TRC = C(1,1)+C(2,2)+C(3,3) TRC2 = C2(1,1)+C2(2,2)+C2(3,3) AUX=TRC**2 P2C = (AUX - TRC2)/2.D0 DTC = C(1,1)*(C(2,2)*C(3,3)-C(2,3)*C(3,2)) . -C(1,2)*(C(2,1)*C(3,3)-C(2,3)*C(3,1)) . +C(1,3)*(C(2,1)*C(3,2)-C(2,2)*C(3,1)) XK=AUX-3.D0*P2C TOL=1.D-6 * IF(IIMPI.EQ.199) THEN WRITE(6,77764) TRC,TRC2,AUX,P2C 77764 FORMAT(2X,'POLA2 - TRC= ',1PE12.5,2X,'TRC2= ',1PE12.5/ . 2X,'AUX= ',1PE12.5,2X,'P2C=',1PE12.5/) WRITE(6,77774) XK,TOL 77774 FORMAT(2X,'POLA2 - XK= ',1PE12.5,2X,'TOL= ',1PE12.5/) ENDIF * IF(XK.LT.TOL) THEN XLAM=SQRT(TRC/3.D0) CALL ZERO(U,9,1) CALL ZERO(UI,9,1) U(1)=XLAM U(5)=XLAM U(9)=XLAM * UNXLAM=1.D0/XLAM UI(1)=UNXLAM UI(5)=UNXLAM UI(9)=UNXLAM * ELSE XL=TRC*(AUX-4.5D0*P2C)+ 13.5D0*DTC AUX2=XL/SQRT(XK**3) IF(IIMPI.EQ.199) THEN WRITE(6,77777) XL,TRC,AUX2,DTC 77777 FORMAT(2X,'POLA2 - XL =',1PE12.5,2X,'TRC=',1PE12.5, . 2X,'AUX2=',1PE12.5,2X,'DTC=',1PE12.5/) ENDIF IF(ABS(AUX2).GT.1.D0) THEN TOLEPS=1.D-10 IF((ABS(AUX2)-1.D0).GE.TOLEPS) THEN ZZZ = ABS(AUX2)-1.D0 WRITE(6,77884) ZZZ ,TOLEPS 77884 FORMAT(2X,'ZOB =',1PE12.5,2X,'TOLEPS=',1PE12.5/) KERRE=26 GO TO 2021 ELSE IF(AUX2.GT.0.D0) THEN AUX2 = AUX2 - TOLEPS ELSE AUX2 = AUX2 + TOLEPS ENDIF ENDIF ENDIF PHI=ACOS(AUX2) XLAM2=(TRC+2.D0*SQRT(XK)*COS(PHI/3.D0))/3.D0 XLAM=SQRT(XLAM2) DTU=SQRT(DTC) TRU=XLAM+SQRT(TRC-XLAM2+(2.D0*DTU/XLAM)) P2U=(TRU**2-TRC)*0.5D0 * FAC1=TRU*P2U-DTU FAC2=TRU*DTU FAC3=TRU**2-P2U * IF(IIMPI.EQ.199) THEN WRITE(6,77775) FAC1,DTU,TRU,P2U 77775 FORMAT(2X,'POLA2 - FAC1 =',1PE12.5,2X,'DTU=',1PE12.5, . 2X,'TRU=',1PE12.5,2X,'P2U=',1PE12.5/) ENDIF IF(FAC1.EQ.0.D0.OR.DTU.EQ.0.D0) THEN WRITE(6,77885) FAC1,DTU 77885 FORMAT(2X,'FAC1=',1PE12.5,2X,'DTU=',1PE12.5/) KERRE=26 GO TO 2021 ENDIF * DO 77890 I=1,N IN=(I-1)*N DO 4 J=1,N IJ= IN+J IF(J.EQ.I) THEN U(IJ)=(FAC2+FAC3*C(I,I)-C2(I,I))/FAC1 UI(IJ)=(P2U-TRU*U(IJ)+C(I,I))/DTU ELSE U(IJ)=(FAC3*C(I,J)-C2(I,J))/FAC1 UI(IJ)=(-TRU*U(IJ)+C(I,J))/DTU ENDIF 4 CONTINUE 77890 CONTINUE * ENDIF ENDIF * ELSE KERRE=19 GO TO 2021 ENDIF * DO 77892 I=1,N IN=(I-1)*N DO 77891 J=1,N IJ= IN+J R(IJ)=0.D0 DO 2 K=1,N IK= IN+K KJ= (K-1)*N+J R(IJ)=R(IJ)+F(IK)*UI(KJ) 2 CONTINUE 77891 CONTINUE 77892 CONTINUE * * PROOF * IF(IIMPI.EQ.199) THEN DO 77894 I=1,N IN=(I-1)*N DO 77893 J=1,N IJ= IN+J UI(IJ)=F(IJ) DO 7 K=1,N IK= IN+K KJ= (K-1)*N+J UI(IJ)=UI(IJ)-R(IK)*U(KJ) 7 CONTINUE 77893 CONTINUE 77894 CONTINUE * WRITE(6,77830) (UI(K),K=1,N2) 77830 FORMAT(2X,'POLA2- PROOF '/(6(1X,1PE12.5))) WRITE(6,77831) (R(K),K=1,N2) 77831 FORMAT(2X,'POLA2- R '/(3(1X,1PE12.5))) WRITE(6,77832) (U(K),K=1,N2) 77832 FORMAT(2X,'POLA2- U '/(3(1X,1PE12.5))) ENDIF * 2021 CONTINUE IF(KERRE.EQ.0) GO TO 9999 * IF(JEBOUC.EQ.1.AND.IIMPI.EQ.1199) THEN IIMPI=199 GO TO 2020 ENDIF * 9999 CONTINUE IIMPI=IIMPI0 * IF(KERRE.NE.0) CALL ERREUR(KERRE) * RETURN END