1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
|
subroutine fmulb1(n,h,x,hx,tabaux,nmisaj,prosca,izs,rzs,dzs)
implicit double precision (a-h,o-z)
external prosca
c Copyright INRIA
c
c parametres
double precision un , deux
parameter ( un=1.0d+0, deux=2.0d+0 )
c declarations
double precision h(*), x(n), hx(n), tabaux(n), dzs(*)
real rzs(*)
integer izs(*)
double precision uscalx, sscalx, nu, eta, gamma, mu, sigma
integer n, nmisaj, memsup, ptnu, compt, iu, is, k
c
c calcul de la longueur d'un bloc
memsup=2*n+2
c calcul de h0*x=x=x
do 1000 k=1,n
hx(k)=x(k)
1000 continue
c
if (nmisaj.eq.0) then
return
else
ptnu=1
compt=1
endif
c
2000 iu=ptnu+1
is=iu+n
do 3000 k=1,n
tabaux(k)=h(iu+k)
3000 continue
call prosca(n,tabaux,x,uscalx,izs,rzs,dzs)
do 4000 k=1,n
tabaux(k)=h(is+k)
4000 continue
call prosca(n,tabaux,x,sscalx,izs,rzs,dzs)
nu=h(ptnu)
eta=h(ptnu+1)
c calcul du nouveau terme et addition dans hx
if (compt.eq.1) then
gamma=eta / nu
do 5000 k=1,n
hx(k)=gamma * hx(k)
5000 continue
mu=sscalx / nu
sigma=-(deux * sscalx / eta)+(uscalx / nu)
else
mu=sscalx / eta
sigma=-(un + nu / eta)* mu + uscalx / eta
endif
c
do 6000 k=1,n
hx(k)=hx(k) - mu * h(iu+k) - sigma * h(is+k)
6000 continue
c
compt=compt+1
if (compt.le.nmisaj) then
ptnu=ptnu+memsup
goto 2000
else
return
endif
end
|
