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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
subroutine fmc11a(a,n,z,sig,w,ir,mk,eps)
implicit double precision (a-h,o-z)
dimension a(*),z(n),w(n)
c update factors given in a by sig*z*ztranspose
if(n.gt.1)goto1
a(1)=a(1)+sig *z(1)**2
ir=1
if(a(1).gt.0.d0)return
a(1)=0.d0
ir=0
return
1 continue
np=n+1
if(sig.gt.0.d0)goto40
if(sig.eq.0.d0.or.ir.eq.0)return
ti=1.d0/sig
ij=1
if(mk.eq.0)goto10
do 7 i=1,n
if(a(ij).ne.0.d0)ti=ti+w(i)**2/a(ij)
7 ij=ij+np-i
goto20
10 continue
do 11 i=1,n
11 w(i)=z(i)
do 15 i=1,n
ip=i+1
v=w(i)
if(a(ij).gt.0.d0)goto12
w(i)=0.d0
ij=ij+np-i
goto15
12 continue
ti=ti+v**2/a(ij)
if(i.eq.n)goto14
do 13 j=ip,n
ij=ij+1
13 w(j)=w(j)-v*a(ij)
14 ij=ij+1
15 continue
20 continue
if(ir.le.0 )goto21
if(ti.gt.0.d0)goto22
if(mk-1)40,40,23
21 ti=0.d0
ir=-ir-1
goto23
22 ti=eps/sig
if(eps.eq.0.d0)ir=ir-1
23 continue
mm=1
tim=ti
do 30 i=1,n
j=np-i
ij=ij-i
if(a(ij).ne.0.d0)tim=ti-w(j)**2/a(ij)
w(j)=ti
30 ti=tim
goto41
40 continue
mm=0
tim=1.d0/sig
41 continue
ij=1
do 66 i=1,n
ip=i+1
v=z(i)
if(a(ij).gt.0.d0)goto53
if(ir.gt.0 .or.sig.lt.0.d0.or.v.eq.0.d0)goto52
ir=1-ir
a(ij)=v**2/tim
if(i.eq.n)return
do 51 j=ip,n
ij=ij+1
51 a(ij)=z(j)/v
return
52 continue
ti=tim
ij=ij+np-i
goto66
53 continue
al=v/a(ij)
if(mm)54,54,55
54 ti=tim+v*al
goto56
55 ti=w(i)
56 continue
r=ti/tim
a(ij)=a(ij)*r
if(r.eq.0.d0)goto70
if(i.eq.n)goto70
b=al/ti
if(r.gt.4.d0)goto62
do 61 j=ip,n
ij=ij+1
z(j)=z(j)-v*a(ij)
61 a(ij)=a(ij)+b*z(j)
goto64
62 gm=tim/ti
do 63 j=ip,n
ij=ij+1
y=a(ij)
a(ij)=b*z(j)+y*gm
63 z(j)=z(j)-v*y
64 continue
tim=ti
ij=ij+1
66 continue
70 continue
if(ir.lt.0)ir=-ir
return
end
|
