Télécharger chole3.eso

Retour à la liste

Numérotation des lignes :

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

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