Sources Esope | Cast3M

triaks
C TRIAKS    SOURCE    FANDEUR   22/05/02    21:15:31     11359          
      SUBROUTINE TRIAKS(IZL)
C************************************************************************
C
C    TRIANGULATION D'UNE MATRICE SYMETRIQUE LIGNE A LIGNE
C
C    POINTEUR : EN ENTREE IZL  CONTIENT LA MATRICE TRIANGULEE
C
C    SUR LES CONSEILS DE PV, LORSQUE L'ON CALCULE LIJ*LKJ*DJ POUR
C    J VARIANT DE JMIN A J-1, ON UTILISE ZKJ=LKJ*DJ QUE L'ON DIVISE
C    ENSUITE PAR DJ POUR RATTRAPER LKJ. ON REMPLACE L'APPEL A SDT3
C    PAR UN APPEL A DDOT ET ON ECONOMISE 1 OPERATION SUR 3.
C
C    ON POURRAIT AMELIORER L'ALGORITHME EN CODANT PAR BLOCS. CETTE
C    PARTIE N'EST PAS CRUCIALE CAR ELLE N'EST EFFECTUEE QU'UNE SEULE
C    FOIS LORS D'UN CALCUL COMPLET.
C
C************************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
      REAL*8   SDT3,SOM
      EXTERNAL SDT3
      REAL*8   DDOT
      EXTERNAL DDOT
 
-INC CCREEL
*-
C-INC SMMATRAKANC
C*************************************************************************
C
C    REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES  puis assemblees
C
 
* LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange
* (points CENTRE ) pour chaque operateur de contrainte
* KGEOC SPG pour la totalite des points CENTRE.
* KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse)
* KLEMC Connectivites de l'ensemble des contraintes
* LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones
 
      SEGMENT MATRAK
      INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP)
      INTEGER LIZAFM(NBSOUS)
      INTEGER IKAM0 (NBSOUS)
      INTEGER IMEM  (NBELC)
      INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC
      ENDSEGMENT
 
      SEGMENT IZAFM
      REAL*8  AM(NNELP,NP,IESP),RPGI(NELAX)
      ENDSEGMENT
 
      POINTEUR IPMJ.IZAFM,IPMK.IZAFM
 
C*******************************************************************
-INC SMMATRK1
      POINTEUR IBLKK.IDBLK,IBLKI.IDBLK
 
      SEGMENT/IZAK/(AK(1)*D)
      SEGMENT/IZAI/(AI(1)*D)
      SEGMENT/IZD/(D(1)*D)
C
C
      GGIS=XPETIT
      SIG=XGRAND
 
      SEGACT IZL
      IDMAT=KZA1
      SEGACT IDMAT
      NBLK=IDESCR(/1)
C
C   CALCUL DE L(I,J) A LA LIGNE J
C    LA MATRICE EST SYMETRIQUE
 
      N=KZA(/1)
      IZD=IDIAG
      SEGACT IZD*MOD
      IF(N.NE.D(/1))CALL ARRET(0)
C
      DO 100 IBLK=1,NBLK
C NUMÉRO DES LIGNES DE DÉBUT ET FIN DE BLOC KJD,KJF
      KJD=NLDBLK(IBLK)
      KJF=NLDBLK(IBLK+1)-1
C ACTIVATION DE LA PORTION DE MATRICE
      IBLKK=IDESCR(IBLK)
      SEGACT IBLKK
      IZAK=IBLKK.IMAT
      SEGACT IZAK*MOD
      DO 1 K=KJD,KJF
      KLL=K-1
C LAK=LONGUEUR DE LA LIGNE K
      LAK=IBLKK.IDEBLK(K-KJD+2)-IBLKK.IDEBLK(K-KJD+1)
C DÉCALAGE DANS LE TABLEAU POUR ACCÉDER À LA LIGNE K
C - NUMÉRO DE LA LIGNE K DANS LE BLOC IBLKK=K-KJD+1
      IDECK=IBLKK.IDEBLK(K-KJD+1)-1
      JK=K-LAK
      L2=1-JK
      IF(LAK.EQ.0)GO TO 1
 
      DO 4 I=JK,KLL
C A QUEL BLOC APPARTIENT I ?
C ATTENTION AU CAS PARTICULIER OÙ I VAUT 1
      IF(I.EQ.1) THEN
         IZAI=IZD
         LAI=N
         IDECI=0
      ELSE
      IBLI=IBLK
      DO 9 IBL=1,NBLK
         IF(I.GE.NLDBLK(IBL).AND.I.LT.NLDBLK(IBL+1)) THEN
            IBLI=IBL
            GOTO 99
         ENDIF
 9    CONTINUE
 99   CONTINUE
      IJD=KJD
      IJF=KJF
      IZAI=IZAK
      IBLKI=IBLKK
C SI I APPARTIENT À UN BLOC DIFFÉRENT DE CELUI DE K
      IF(IBLI.NE.IBLK) THEN
         IJD=NLDBLK(IBLI)
         IJF=NLDBLK(IBLI+1)-1
         IBLKI=IDESCR(IBLI)
         SEGACT IBLKI
         IZAI=IBLKI.IMAT
         SEGACT IZAI*MOD
      ENDIF
C LAI=LONGUEUR DE LA LIGNE I
      LAI=IBLKI.IDEBLK(I-IJD+2)-IBLKI.IDEBLK(I-IJD+1)
C DÉCALAGE DANS LE TABLEAU POUR ACCÉDER À LA LIGNE I
C - NUMÉRO DE LA LIGNE I DANS LE BLOC IBLKI=I-IJD+1
      IDECI=IBLKI.IDEBLK(I-IJD+1)-1
      ENDIF
      JI=I-LAI
      L1=1-JI
      I1=I-1
      L=L2+I
      SOM=0.D0
      JMIN=MAX0(JI,JK)
      IF (JMIN.LE.I1) THEN
        JMI1=I1-JMIN+1
        SOM=DDOT(JMI1,AI(IDECI+L1+JMIN),1,AK(IDECK+L2+JMIN),1)
      ENDIF
      AK(IDECK+L)=AK(IDECK+L)-SOM
      IF(IZAI.NE.IZAK.AND.IZAI.NE.IZD) THEN
        SEGDES IZAI
        SEGDES IBLKI
      ENDIF
 4    CONTINUE
C IL FAUT DIVISER CHAQUE TERME DE LA LIGNE QUI VIENT D'ÊTRE CALCULÉE PAR
C LE TERME DIAGONAL CORRESPONDANT.
      DO 5 I=JK,KLL
         L=L2+I
         AK(IDECK+L)=AK(IDECK+L)/D(I)
 5    CONTINUE
C ON TRAITE LE TERME DIAGONAL, ON NE PEUT FAIRE LA RUSE DU TERME GÉNÉRAL
      SOM=SDT3(LAK,AK(IDECK+1),AK(IDECK+1),D(K-LAK))
      D(K)=D(K)-SOM
 
C! DEBUT
      ADK=ABS(D(K))
      IF(ADK.LE.GGIS)THEN
      WRITE(6,*)'SUB TRIAWS : DIAGONALE NULLE ADK=',ADK,' LIGNE K=',
     &NUNA(K)
      D(K)=SIG
      ENDIF
C! FIN
 
 1    CONTINUE
 
      SEGDES IZAK
      SEGDES IBLKK
 100  CONTINUE
 
      SEGDES IZD
      SEGDES IDMAT
      SEGDES IZL
 
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
      END