C ELFLI1 SOURCE CHAT 05/01/12 23:32:09 5004 SUBROUTINE ELFLI1(C1,C2,C3,C4,N,B,IS,ISING,EPS,NLBLI,NBLI) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C C ===================================================================== C SOUS-PROGRAMME FORTRAN APPELE PAR ELFLIA C CALCUL DES MATRICES C3 ET C4 DE DIMENSION N C C EN ENTREE : C C C1 MATRICE BANDE A1 C -1 C C2 MATRICE BANDE A2 C C3 MATRICE A3 C C4 MATRICE A4 C B,IS TABLEAU DE TRAVAIL POUR L'INVERSION C ISING = 1 SI MATRICE SINGULIERE 0 SINON C EPS PRECISION C NLBLI TABLEAU POUR L'UTILISATION DES MATRICES BANDES C C EN SORTIE : C -1 -1 C C3 MATRICE (A3-A4*A2 *A1) C -1 -1 -1 C C4 MATRICE (A3-A4*A2 *A1) *A4*A2 C C CREATION : 04/08/87 C PROGRAMMEUR : GUILBAUD C ===================================================================== C -INC CCOPTIO DIMENSION C1(*),C2(*),C3(N,1),C4(N,1),B(*),IS(1),NLBLI(2,1) C -1 C 1 - A4 = A4 * A2 C C WRITE(IOIMP,*) ' A4 * A2-1 ' DO 20 I=1,N NB=0 DO 10 NBL=1,NBLI NX1=NLBLI(1,NBL) NX2=NLBLI(2,NBL) IF(NX1.EQ.NX2) THEN C4(I,NX1+1)=C2(NB+1)*C4(I,NX1+1) C4(I,NX1+2)=C2(NB+2)*C4(I,NX1+2) C4(I,NX1+3)=C2(NB+3)*C4(I,NX1+3) C4(I,NX1+4)=C2(NB+4)*C4(I,NX1+4) C4(I,NX1+5)=C2(NB+5)*C4(I,NX1+5) C4(I,NX1+6)=C2(NB+6)*C4(I,NX1+6) NB=NB+6 ELSE B( 1)=C2(NB+ 1)*C4(I,NX1+1)+C2(NB+21)*C4(I,NX2+1) B( 2)=C2(NB+ 2)*C4(I,NX1+2)+C2(NB+ 9)*C4(I,NX1+6) * +C2(NB+22)*C4(I,NX2+2)+C2(NB+29)*C4(I,NX2+6) B( 3)=C2(NB+ 4)*C4(I,NX1+3)+C2(NB+ 7)*C4(I,NX1+5) * +C2(NB+24)*C4(I,NX2+3)+C2(NB+27)*C4(I,NX2+5) B( 4)=C2(NB+ 6)*C4(I,NX1+4)+C2(NB+26)*C4(I,NX2+4) B( 5)=C2(NB+ 5)*C4(I,NX1+3)+C2(NB+ 8)*C4(I,NX1+5) * +C2(NB+25)*C4(I,NX2+3)+C2(NB+28)*C4(I,NX2+5) B( 6)=C2(NB+ 3)*C4(I,NX1+2)+C2(NB+10)*C4(I,NX1+6) * +C2(NB+23)*C4(I,NX2+2)+C2(NB+30)*C4(I,NX2+6) B( 7)=C2(NB+11)*C4(I,NX1+1)+C2(NB+31)*C4(I,NX2+1) B( 8)=C2(NB+12)*C4(I,NX1+2)+C2(NB+19)*C4(I,NX1+6) * +C2(NB+32)*C4(I,NX2+2)+C2(NB+39)*C4(I,NX2+6) B( 9)=C2(NB+14)*C4(I,NX1+3)+C2(NB+17)*C4(I,NX1+5) * +C2(NB+34)*C4(I,NX2+3)+C2(NB+37)*C4(I,NX2+5) B(10)=C2(NB+16)*C4(I,NX1+4)+C2(NB+36)*C4(I,NX2+4) B(11)=C2(NB+15)*C4(I,NX1+3)+C2(NB+18)*C4(I,NX1+5) * +C2(NB+35)*C4(I,NX2+3)+C2(NB+38)*C4(I,NX2+5) B(12)=C2(NB+13)*C4(I,NX1+2)+C2(NB+20)*C4(I,NX1+6) * +C2(NB+33)*C4(I,NX2+2)+C2(NB+40)*C4(I,NX2+6) C C4(I,NX1+1)=B( 1) C4(I,NX1+2)=B( 2) C4(I,NX1+3)=B( 3) C4(I,NX1+4)=B( 4) C4(I,NX1+5)=B( 5) C4(I,NX1+6)=B( 6) C4(I,NX2+1)=B( 7) C4(I,NX2+2)=B( 8) C4(I,NX2+3)=B( 9) C4(I,NX2+4)=B(10) C4(I,NX2+5)=B(11) C4(I,NX2+6)=B(12) C WRITE(IOIMP,*) I C WRITE(IOIMP,1000) (C4(I,NX1+K),K=1,6),(C4(I,NX2+K),K=1,6) C1000 FORMAT(12(1X,1PE10.3)) NB=NB+40 ENDIF 10 CONTINUE 20 CONTINUE C -1 C 2 - A3 = A3 - A4 * A2 * A1 C C WRITE(IOIMP,*) ' A3-A4 * A2-1 *A1 ' DO 40 I=1,N NA=0 DO 30 NBL=1,NBLI NX1=NLBLI(1,NBL) NX2=NLBLI(2,NBL) NX=NX1 IF(NX1.EQ.NX2) THEN C3(I,NX+1)=C3(I,NX+1)-C1(NA+1)*C4(I,NX+1) C3(I,NX+2)=C3(I,NX+2)-C1(NA+2)*C4(I,NX+2)-C1(NA+9)*C4(I,NX+6) C3(I,NX+3)=C3(I,NX+3)-C1(NA+4)*C4(I,NX+3)-C1(NA+7)*C4(I,NX+5) C3(I,NX+4)=C3(I,NX+4)-C1(NA+6)*C4(I,NX+4) C3(I,NX+5)=C3(I,NX+5)-C1(NA+5)*C4(I,NX+3)-C1(NA+8)*C4(I,NX+5) C3(I,NX+6)=C3(I,NX+6)-C1(NA+3)*C4(I,NX+2)-C1(NA+10)*C4(I,NX+6) NA=NA+10 ELSE C3(I,NX+1)=C3(I,NX+1)-C1(NA+ 1)*C4(I,NX+1)-C1(NA+21)*C4(I,NX2+1) C3(I,NX+2)=C3(I,NX+2)-C1(NA+ 2)*C4(I,NX+2)-C1(NA+ 9)*C4(I,NX+6) * -C1(NA+22)*C4(I,NX2+2)-C1(NA+29)*C4(I,NX2+6) C3(I,NX+3)=C3(I,NX+3)-C1(NA+ 4)*C4(I,NX+3)-C1(NA+ 7)*C4(I,NX+5) * -C1(NA+24)*C4(I,NX2+3)-C1(NA+27)*C4(I,NX2+5) C3(I,NX+4)=C3(I,NX+4)-C1(NA+ 6)*C4(I,NX+4)-C1(NA+26)*C4(I,NX2+4) C3(I,NX+5)=C3(I,NX+5)-C1(NA+ 5)*C4(I,NX+3)-C1(NA+ 8)*C4(I,NX+5) * -C1(NA+25)*C4(I,NX2+3)-C1(NA+28)*C4(I,NX2+5) C3(I,NX+6)=C3(I,NX+6)-C1(NA+ 3)*C4(I,NX+2)-C1(NA+10)*C4(I,NX+6) * -C1(NA+23)*C4(I,NX2+2)-C1(NA+30)*C4(I,NX2+6) C3(I,NX2+1)=C3(I,NX2+1)-C1(NA+11)*C4(I,NX+1)-C1(NA+31)*C4(I,NX2+1) C3(I,NX2+2)=C3(I,NX2+2)-C1(NA+12)*C4(I,NX+2)-C1(NA+19)*C4(I,NX+6) * -C1(NA+32)*C4(I,NX2+2)-C1(NA+39)*C4(I,NX2+6) C3(I,NX2+3)=C3(I,NX2+3)-C1(NA+14)*C4(I,NX+3)-C1(NA+17)*C4(I,NX+5) * -C1(NA+34)*C4(I,NX2+3)-C1(NA+37)*C4(I,NX2+5) C3(I,NX2+4)=C3(I,NX2+4)-C1(NA+16)*C4(I,NX+4)-C1(NA+36)*C4(I,NX2+4) C3(I,NX2+5)=C3(I,NX2+5)-C1(NA+15)*C4(I,NX+3)-C1(NA+18)*C4(I,NX+5) * -C1(NA+35)*C4(I,NX2+3)-C1(NA+38)*C4(I,NX2+5) C3(I,NX2+6)=C3(I,NX2+6)-C1(NA+13)*C4(I,NX+2)-C1(NA+20)*C4(I,NX+6) * -C1(NA+33)*C4(I,NX2+2)-C1(NA+40)*C4(I,NX2+6) C WRITE(IOIMP,*) I C WRITE(IOIMP,1000) (C3(I,NX1+K),K=1,6),(C3(I,NX2+K),K=1,6) NA=NA+40 ENDIF 30 CONTINUE 40 CONTINUE C DO 35 I=1,N C WRITE(6,1002) (C3(I,J),J=1,N) C1002 FORMAT(12(1X,1PE10.3)) C 35 CONTINUE C C -1 C 3 - C3 = A3 C CALL INVER(C3,N,ISING,B,IS,EPS) IF(ISING.EQ.1) RETURN C -1 -1 C 2 - C4 = A3 * A4 * A2 C DO 80 J=1,N DO 60 I=1,N AA=0.D0 DO 50 K=1,N AA=AA+C3(I,K)*C4(K,J) 50 CONTINUE B(I)=AA 60 CONTINUE DO 70 I=1,N C4(I,J)=B(I) 70 CONTINUE 80 CONTINUE C RETURN END