Télécharger chole3.eso

Retour à la liste

Numérotation des lignes :

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

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