Sources Esope | Cast3M

extrap
C EXTRAP    SOURCE    CHAT      05/01/12    23:53:57     5004
      SUBROUTINE EXTRAP(SHPTOT,NBPGAU,NBNN,NBNO)
C================================================================
C   CALCULE LES FONCTIONS D EXTRAPOLATIONS A PARTIR DES FONCTIONS
C        D INTERPOLATIONS
C ENTREES
C   SHPTOT(6,NBNO,NBPGAU)   = FONCTIONS D INTERPOLATIONS
C   NBPGAU                  = NOMBRE DE POINTS DE GAUSS
C   NBNN                    = NOMBRE DE NOEUDS
C   NBNO                    = NOMBRE DE FONCTIONS D'INTERPOLATION
C SORTIES
C   SHPTOT(6,NBNO,NBPGAU)   = FONCTIONS D EXTRAPOLATIONS STOKEES
C                           SUR LA 6 IEME LIGNE
C  EBERSOLT NOVEMBRE 86    PAS PLUS DE 30 NOEUDS
C================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(XZER=0.D0,UN=1.D0)
      DIMENSION SHPTOT(6,NBNO,*)
      DIMENSION XMAT(30,30),XVEC(30)
C
C  PROTECTION PROVISOIRE
       IF(NBNN.NE.NBNO) RETURN
C
C      UN SEUL POINT DE GAUSS
C
      IF(NBPGAU.EQ.1) THEN
      DO 50 IA=1,NBNO
      SHPTOT(6,IA,1)=UN
  50  CONTINUE
C
C     PLUS D UN POINT DE GAUSS
C
      ELSE IF(NBPGAU.GT.1) THEN
      CALL ZERO(XMAT,30,30)
C
C     TRANSPOSE( A ) * A
C
      DO 100 IA=1,NBNN
      DO 100 IB=1,NBNN
      CC = XZER
      DO 300 IC=1,NBPGAU
      CC = CC + SHPTOT(1,IA,IC)*SHPTOT(1,IB,IC)
  300 CONTINUE
      XMAT(IA,IB)=CC
  100 CONTINUE
C
C      NOMBRE DE POINTS DE GAUSS DIFFERENTS DU NOMBRE DE NOEUDS
C
      IF(NBPGAU.NE.NBNN) THEN
C
C     SI ON A MOINS DE POINTS DE GAUSS  QUE DE NOEUDS ON RAJOUTE
C      UN PEU DE PENALISATION  EMPECHANT D OSCILLER SUR LES NOEUDS
C
      IF(NBPGAU.LT.NBNN) THEN
      DO 705 IA=1,NBNN
      CALL ZERO(XVEC,30,1)
      XVEC(IA)=NBNN
      DO 706 IB=1,NBNN
      XVEC(IB)=XVEC(IB)-UN
  706 CONTINUE
C
C     ON LE DEFLATIONNE  DE SES COMPOSANTES PARRALELLES AUX H
C
      DO 710 IB=1,NBPGAU
      SCAL=XZER
      XXNORM=XZER
      DO 730 IC=1,NBNN
      SCAL=SCAL+SHPTOT(1,IC,IB)*XVEC(IC)
      XXNORM=XXNORM+SHPTOT(1,IC,IB)*SHPTOT(1,IC,IB)
 730  CONTINUE
      IF(XXNORM.LT.1.E-7) GOTO 700
      SCAL=SCAL/XXNORM
C
      DO 720 IC=1,NBNN
      XVEC(IC)=XVEC(IC)-SCAL*SHPTOT(1,IC,IB)
 720  CONTINUE
 710  CONTINUE
C
C     ON RAJOUTE CES VECTEURS DANS LA PENALISATION
C
      DO 750 IB=1,NBNN
      DO 750 IC=1,NBNN
      XMAT(IB,IC)= XVEC(IB)*XVEC(IC)+XMAT(IB,IC)
  750 CONTINUE
C
  700 CONTINUE
  705 CONTINUE
C
      ENDIF
C
C      ( T A P A ) ** -1    ALGO WILSON
C
      DO 400 IEQ=1,NBNN
      DD = UN / XMAT(IEQ,IEQ)
      DO 410 IA=1,NBNN
      XMAT(IEQ,IA)=-XMAT(IEQ,IA)*DD
  410 CONTINUE
C
      DO 420 IA=1,NBNN
      IF(IA.EQ.IEQ) GOTO 420
      DO 430 IB=1,NBNN
      IF(IB.EQ.IEQ) GOTO 430
      XMAT(IA,IB)=XMAT(IA,IB)+XMAT(IA,IEQ)*XMAT(IEQ,IB)
  430 CONTINUE
  420 CONTINUE
C
      DO 440 IA=1,NBNN
      XMAT(IA,IEQ)= XMAT(IA,IEQ)*DD
  440 CONTINUE
      XMAT(IEQ,IEQ)= DD
  400 CONTINUE
C
C     (( T A . A ) ** -1 )  * ( T . A )
C
      DO 500 IA=1,NBNN
      DO 510 IB=1,NBPGAU
      CC=XZER
      DO 520 IC=1,NBNN
      CC=CC+XMAT(IA,IC)*SHPTOT(1,IC,IB)
  520 CONTINUE
      SHPTOT(6,IA,IB)=CC
  510 CONTINUE
  500 CONTINUE
C
C      NOMBRE DE POINTS DE GAUSS EGAL AUX NOMBRE DE NOEUDS
C
      ELSE IF(NBNN.EQ.NBPGAU) THEN
      DO 600 IA=1,NBNN
      DO 600 IB=1,NBNN
      SHPTOT(6,IA,IB)=XMAT(IA,IB)
  600 CONTINUE
      ENDIF
      ENDIF
      RETURN
      END