Sources Esope | Cast3M

elfli1
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 PPARAM
-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