Sources Esope | Cast3M

chole1
C CHOLE1    SOURCE    PV090527  24/01/12    21:15:02     11821          
      FUNCTION CHOLE1(ILIGF,iprell,imasql,VALF,DAAG,IPKNO,IPPVF,KHG,
     >IVPOF,KIDEP,KI1,KQ,idep,prec,nbo)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC SMMATRI
-INC SMRIGID
-INC CCHOLE
-INC CCREEL
*   DAAG contient l'inverse de la diagonale de facon a faire des multiplications et non des divisions
      DIMENSION ILIGF(*),VALF(*),DAAG(*),IPKNO(*),IPPVF(*),IVPOF(*)
      dimension imasql(*)
*  constante pour la somme compensee utilisee pour l'accumulation. CF ddot2
      parameter (c=2.D0**26+1.D0)
      xmatri=matric
      nbnnma=nbnnmc
      IPPKHG=IPPVF(KHG)
      KBAS=IPKNO(KIDEP)
      KHAU=IPKNO(KI1)
      KDIAG=KI1+1
      DNORM=ABS(VALF(KDIAG))*PREC
      KPREM=IVPOF(KHG)-IPPKHG
      ILIG=IPRELL+KHG-NBNNMA-1
      IECAR=KQ-IPRELL+1
      DO 30 NNJ=MAX(1,KIDEP+IECAR),KI1+IECAR
       KK=NNJ-IECAR
       ICOL=KQ+KK-NBNNMA
       NNJJ=IPPVF(NNJ+1)
       NJ=NNJJ-IPPVF(NNJ)
       LLOL=MIN(NJ,KK)-1
       LLON=MIN(LLOL-KK+KPREM+1,LLOL-NNJJ+IVPOF(NNJ)+1)
       LLON=MIN(LLON,NBNNMA-KQ-KK+LLOL+1)
       IF (LLON.GT.0.and.kk.ge.1) THEN
        IEC1=KK-LLOL-1
        IEC2=NNJJ-IPPKHG-KK
        ideq=1+iec1+idep-1
        if (llon.gt.masdim+1) then
        p=ddotpw(llon,VALF(1+iec1),VALF(1+iec1+iec2),
     >     imasql(1),1+idep-1,nbo)
        else
         if (imasql(ideq/masdim+1).gt.0.or.
     >       imasql((ideq+llon-1)/masdim+1).gt.0) then
           p=ddotpv(llon,VALF(1+iec1),VALF(1+iec1+iec2))
           nbo=nbo+llon
         else
            p=0.d0
         endif
        endif
        VALF(KK)=VALF(KK)-P
        IF (ABS(VALF(KK)).GT.DNORM) then
         KPREM=KK
         imasql((kk+idep-1)/masdim+1) =1
*  si on remonte, on tombe au terme diagonal ou apres, mais ce n'est qu'un seul terme
         imasql(kk/masdim+1) =1
        ELSE
*  annuler le terme car on l'ignorera par la suite
         valf(kk)=0.d0
        ENDIF
       ENDIF
        if (ilig.ge.1.and.icol.ge.1) then
         RE(ILIG,ICOL,1)=VALF(KK)
         RE(ICOL,ILIG,1)=VALF(KK)
        endif
  30  CONTINUE
  50  CONTINUE
**    AUX1=0.D0
      s1=0.d0
      c1=0.d0
      kdeb=1
  43  continue
      kdebi=kdeb
  44  continue
      do 100 im=kdeb/masdim+1,kprem/masdim+1
         jm=imasql(im)
         if (jm.gt.0) goto 105
         if (jm.eq.0) goto 100
         jinio=-jm/masdim+1
        if (jinio.gt.im+jacc) then
*         write (6,*) 'saut kdeb jm ',kdeb,jm
          kdeb=-jm
          goto 44
        endif
 100  continue
 105  continue
      ideb=max((im-1)*masdim,kdebi)
      kdeb=ideb
 111  continue
      do 110 im=kdeb/masdim+1,kprem/masdim+1
         jm=imasql(im)
         if (jm.le.0) goto 115
         if (jm.eq.1) goto 110
         jfineo=jm/masdim+1
         if(jfineo.gt.im+jacc) then
           kdeb=jm
           goto 111        
         endif
 110  continue
 115  continue
      im=im-1
      ifin=min(im*masdim-1,kprem)
**    write (6,*) ' chole1 kdeb kprem ideb ifin ',kdeb,kprem,ideb,ifin
      DO 9 K=ideb,min(ifin,nbnnma-kq)
       AUX=VALF(K)
       if (aux.eq.0.d0) goto 9
       nbo=nbo+1
       VALFK=AUX*DAAG(K)
       VALF(K)=VALFK
**     AUX1=AUX1+AUX*VALF(K)
** twoprod
*pv    aux2=valf(k)
*pv    z=c*aux2
*pv    xx=z-(z-aux2)
*pv    xy=aux2-xx
*pv    z=c*aux
*pv    yx=z-(z-aux)
*pv    yy=aux-yx
*pv    s2=aux2*aux
*pv    c2=xy*yy-(((s2-xx*yx)-xx*yy)-xy*yx)
** twosum
       s2=aux*valfk
       xx=s1+s2
       z=xx-s1
       xy=(s1-(xx-z))+(s2-z)
       s1=xx
*pv    c1=c1+(xy+c2)
       c1=c1+xy
   9  CONTINUE
      if (ifin.lt.kprem) then
       kdeb=ifin+1
       goto 43
      endif
      aux1=s1+c1
      ivpof(khg)=kprem+ippkhg
      CHOLE1=-AUX1
      RETURN
      END