Télécharger la plaquette Cast3M


fr en

Télécharger chole3.eso

Retour à la liste

Numérotation des lignes : 

chole3

  1. C CHOLE3    SOURCE    PV090527  25/08/26    21:15:02     10699          
  2.       SUBROUTINE CHOLE3(IPREL,IDERL,lpl,IPPR,IDDR,IVPO,
  3.      # nbg1,VAL,VAL1,IVPO1              ,imasq,nbo,irondi,irondf,ivd)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6. -INC SMRIGID
  7. -INC CCHOLE
  8. -INC CCREEL
  9.       DIMENSION IVPO(*),VAL(*),VAL1(*),IVPO1(*)
  10.       DIMENSION imasq(*)
  11.       real*8 pt(36)
  12.  
  13. *  nombre max de lignes traitees simultanement
  14.       nbl=6
  15. *  recuperer en cas de super element dans cchole
  16.       xmatri=matric
  17.       nbnnma=nbnnmc
  18.  
  19.  
  20.  
  21. C  inconnues correspondant aux noeuds
  22.  
  23. **    write (6,*) ' chole3 irondi irondj ',irondi
  24.  
  25. C val1 en stockage compacte
  26. C val mis a jour et utilise imasq pour avoir les termes non nuls
  27. C  nombres de groupes (incluant la diagonale)
  28. **      nbg1=ippvv1(2)-1
  29.       nbg=1
  30. C  nb ligne
  31.       na=iderl-iprel+1
  32.       na1=iddr-ippr+1
  33. C  longueur de la premiere ligne incluant la diagonale
  34.       lpl1=ivpo1(2*(nbg1+1))-1              
  35. **    lpl1=ivpo1(2*(nbg1+1))-ivpo1(2*1)
  36. **    write(6,*) 'chole3 ivpo1 ',ivpo1(1),ivpo1(2)
  37. C  nb termes premiere ligne
  38.       nval1=ivpo1(2*(nbg1+1)-1)-1                 
  39. **    nval1=ivpo1(2*(nbg1+1)-1)-ivpo1(2*1-1)
  40. C  longueur de la premiere ligne de val
  41. **      lpl=ippvv(2)-ippvv(1)
  42. C  nb termes premiere ligne de val
  43.       nval=lpl
  44. C  position depart de val et val1
  45.       idepv=iprel-nval+1
  46.       idepv1=ippr-nval1+1
  47.       imb=idepv1-idepv
  48. C     write (6,*) 'chole3 idepv1 iml',idepv1,iml
  49. C  les groupes (hors groupe diagonal)
  50.       kidepg=ivpo(1)
  51.       do 121 im=2,na
  52.         kidepg=max(kidepg,ivpo(im))
  53.  121  continue
  54.  
  55.       do 10 ig1=1,nbg1-1
  56. C  il position dans la ligne compressee
  57. C  i position relative dans la ligne
  58.         ildeb1=ivpo1(2*ig1)
  59.         ilfin1=ivpo1(2*(ig1+1))-1
  60.         ideb1=ivpo1(2*ig1-1)
  61.         ifin1=ideb1+ilfin1-ildeb1
  62.         ifin1=min(ifin1,nbnnma-idepv1+1)
  63.  
  64.         ideb1n=max(1-imb,ideb1)
  65.         long=ifin1-ideb1n+1
  66.         lond=min(long,kidepg-imb-ideb1n+1)
  67. **      write(6,*) ' lond ',lond
  68.         ifin1=lond+ideb1n-1
  69.         ideb1n=max(irondi-imb,ideb1n)
  70. *  verif rapide avec imasq qu'il y a du travail a faire
  71.         jini=(ideb1n+imb)/masdim
  72.         if(-imasq(jini).gt.iddr) then
  73. **        write(6,*) ' chole3 sortie rapide ',-imasq(jini),ippr,iddr
  74.           goto 11
  75.         endif
  76.         if(-imasq(jini).gt.ifin1+imb) goto 10
  77.         ifin1 =min(irondf-imb,ifin1 )
  78.         lond=ifin1-ideb1n+1
  79.         if (ifin1.lt.ivd-imb) then
  80. **    fin d'operation avant le premier terme significatif de la rondelle
  81. **       write (6,*) ' chole3 saut',ivd,irondi,irondf
  82.          goto 10
  83.         endif
  84.         if (lond.le.0) goto 10
  85. C optimisation pour le cas na>1 na1>1
  86. C on decoupe les operations en groupes de 6x6 car au dela n'est pas programme
  87.        if (na.gt.1.or.na1.gt.1) then
  88.        do 301 ia=0,na-1,nbl
  89.         iposrb=imb+ia*lpl+(ia*(ia-1))/2
  90.        do 300 ia1=0,na1-1,nbl
  91.         ilpos1b=-ideb1+ildeb1+ia1*lpl1+(ia1*(ia1-1))/2
  92.            nboq=nbo
  93. **         write(6,*) 'appel mamupv ia na ia1 na1 ',ia,na,ia1,na1
  94.            npa=min(nbl,na-ia)
  95.            npa1=min(nbl,na1-ia1)
  96.            call mamupv(ideb1n,ifin1,val(1),iposrb,lpl+ia,val1(1),
  97.      >     ilpos1b,lpl1+ia1,imasq(1),imb,pt,
  98.      >     npa,npa1,nbo)
  99.            if (nbo.eq.nboq) then
  100. *  rien a mettre a jour
  101. ***        write (6,*) ' mamupv rien a mettre a jour '
  102.            goto 300
  103.            endif
  104. C  mise a jour val
  105.            do ic=1,npa
  106.              ivad=ippr-idepv+ia1+(ia+ic-1)*lpl+
  107.      >         ((ia+ic-1)*(ia+ic-2))/2
  108.              iaux=-ivad+(ic-1)*npa1
  109.            do iv=max(1,1+ivad),npa1+ivad
  110.              if(abs(pt(iv+iaux)).gt.xpetit)
  111.      >       val(iv)=val(iv)-pt(iv+iaux)
  112.            enddo
  113.            enddo
  114.   300  continue
  115.   301  continue
  116.        elseif (na.eq.1.and.na1.eq.1) then
  117.              ideqb= ideb1n+imb
  118.          im1=1
  119.              idebzc=ideb1n-ideb1+ildeb1+(im1-1)*lpl1+((im1-2)*(im1-1))/2
  120.          im=1
  121.            ideq= ideqb+(im-1)*lpl+((im-2)*(im-1))/2
  122.            nboq=nbo
  123.            p=ddotpw(lond,val(ideq),val1(idebzc),imasq(1),ideqb,nbo)
  124.            if (nbo.ne.nboq) then
  125. C  mise a jour val
  126.            ivad=ippr-idepv+im1+(im-1)*lpl+((im-2)*(im-1))/2
  127.            if (ivad.ge.1.and.abs(p).gt.xpetit) val(ivad)=val(ivad)-p
  128.            endif
  129.        endif
  130.   10  continue
  131.   11  continue
  132. 999   continue
  133.       if (irondi-imb.gt.nval1.or.irondf-imb.lt.nval1 ) return
  134.  
  135. C  le groupe diagonal
  136.       ig1=nbg1
  137. C les lignes de val1 on s'arrete avant le terme diagonal.
  138.           ivadb=ippr-idepv+1
  139. C*   >      ippr,iblcd
  140.           kidepb=ivpo(1)
  141.           do 210 im=1,na
  142.           icol=iprel +im-1-nbnnma
  143. **        write(6,*) 'super3 nbnnma icol ',nbnnma,icol
  144.            kidep=ivpo(im)
  145. C  comparaison au terme diagonal
  146.            dnorm = 
  147.      >       max(precc * abs(val(im*lpl+(im*(im-1))/2)),xpetit)
  148.           do 220 im1=1,na1
  149.           ilig=ippr+im1-1-nbnnma
  150.           p=0.d0
  151.  
  152.            ildeb1=ivpo1(2*ig1)
  153. C*           ilfin1=ildeb1+im1-2
  154.            ideb1=ivpo1(2*ig1-1)
  155.            ideb1n=max(1-imb,ideb1)
  156.            ifin1=ideb1+im1-2
  157.         ifin1=min(ifin1,nbnnma-idepv1+1)
  158. C le travail sur les colonnes de val
  159.             ilpos1=ideb1n-ideb1+ildeb1 +(im1-1)*lpl1+((im1-2)*(im1-1))/2
  160.             iposr=ideb1n+imb+(im-1)*lpl+((im-2)*(im-1))/2
  161.             ideqb=ideb1n+imb
  162.  
  163. **         p=ddotpw(ifin1-ideb1n+1,val(iposr),val1(ilpos1),imasq(1),
  164. **   >           ideqb,nbo)
  165.  
  166. * il y a trop peu de travail pour appeler ddotpw
  167.  
  168.            if (ifin1-ideb1n.ge.0) then
  169.            nbo=nbo+ifin1-ideb1n+1
  170.            do 200 ipos=ideb1n,ifin1
  171.             p=p+val(iposr)*val1(ilpos1)
  172.             ilpos1=ilpos1+1
  173.             iposr=iposr+1
  174.  200       continue
  175.            endif
  176. C  mise a jour val
  177.            ivad=ippr-idepv+im1+(im-1)*lpl+((im-2)*(im-1))/2
  178.            if (ivad.lt.1) goto 220
  179.            ivadb=ippr-idepv+im1
  180.            if (ivadb.lt.1) goto 220
  181.             if(abs(p).gt.xpetit) val(ivad)=val(ivad)-p
  182. **    write(6,*) 'super3 xmatri ilig icol re ',
  183. **   > xmatri,ilig,icol,val(ivad)
  184. *  cas du super element
  185.             if(xmatri.ne.0.and.ilig.ge.1.and.icol.ge.1) then
  186.             re(ilig,icol,1)=val(ivad)
  187.             re(icol,ilig,1)=val(ivad)
  188.             endif
  189.             if (abs(val(ivad)).gt.dnorm) then
  190.              imasq(ivad/masdim+1)=1
  191.              if (ivadb.le.lpl) imasq(ivadb/masdim+1)=1
  192. C  mise a jour imasq
  193.              do imt=kidep/masdim+1+1,ivad/masdim+1-1
  194.               imasq(imt)=-(ivad/masdim)*masdim+1
  195.              enddo
  196.              kidep=ivad
  197. C*           if (im.eq.1) then
  198. C*            kidepb=kidep
  199. C*           endif
  200.  
  201.              if (im.ne.1) then
  202.               do imt=kidepb/masdim+1+1,ivadb/masdim+1-1
  203.                imasq(imt)=-(ivadb/masdim)*masdim+1
  204.               enddo
  205.              endif
  206.              kidepb=max(ivadb,kidepb)
  207.            else
  208. C  si on note que la valeur est nulle, il faut qu'elle le soit vraiment
  209. **           val(ivad)=0.d0
  210. C*           write (6,*) ' chole3 ivad lpl ',ivad,lpl
  211.            endif
  212.  220      continue
  213.           ivpo(im)=kidep
  214.           ivpo(1)=kidepb
  215.  210      continue
  216. **       if(ivpo(1).ne.1) write(6,*) 'chole3 ivpo',ivpo(1)
  217.  
  218.       RETURN
  219.       END
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales