← Back to branch summary

~daniele-bigoni/pyorthpol/pyorthpol-python3

« back to all changes in this revision

Viewing changes to ORTHPOLxx/src/dgauss.f

  • Committer: Daniele Bigoni
  • Date: 2014-12-16 10:43:03 UTC
  • Revision ID: dabi@dtu.dk-20141216104303-ba4utrt7ln5g39lq
started with launchpad

Show diffs side-by-side

added added

removed removed

Lines of Context:
ORTHPOLxx/src/dgauss.f
 
1
c
 
2
c
 
3
      subroutine dgauss(n,dalpha,dbeta,deps,dzero,dweigh,ierr,de)
 
4
c
 
5
c This is a double-precision version of the routine  gauss.
 
6
c
 
7
      double precision dalpha,dbeta,deps,dzero,dweigh,de,dp,dg,dr,
 
8
     *ds,dc,df,db
 
9
      dimension dalpha(n),dbeta(n),dzero(n),dweigh(n),de(n)
 
10
      if(n.lt.1) then
 
11
        ierr=-1
 
12
        return
 
13
      end if
 
14
      ierr=0
 
15
      dzero(1)=dalpha(1)
 
16
      if(dbeta(1).lt.0.d0) then
 
17
        ierr=-2
 
18
        return
 
19
      end if
 
20
      dweigh(1)=dbeta(1)
 
21
      if (n.eq.1) return
 
22
      dweigh(1)=1.d0
 
23
      de(n)=0.d0
 
24
      do 100 k=2,n
 
25
        dzero(k)=dalpha(k)
 
26
        if(dbeta(k).lt.0.d0) then
 
27
          ierr=-2
 
28
          return
 
29
        end if
 
30
        de(k-1)=dsqrt(dbeta(k))
 
31
        dweigh(k)=0.d0
 
32
  100 continue
 
33
      do 240 l=1,n
 
34
        j=0
 
35
  105   do 110 m=l,n
 
36
          if(m.eq.n) goto 120
 
37
          if(dabs(de(m)).le.deps*(dabs(dzero(m))+dabs(dzero(m+1)))) 
 
38
     *      goto 120
 
39
  110   continue
 
40
  120   dp=dzero(l)
 
41
        if(m.eq.l) goto 240
 
42
        if(j.eq.30) goto 400
 
43
        j=j+1
 
44
        dg=(dzero(l+1)-dp)/(2.d0*de(l))
 
45
        dr=dsqrt(dg*dg+1.d0)
 
46
        dg=dzero(m)-dp+de(l)/(dg+dsign(dr,dg))
 
47
        ds=1.d0
 
48
        dc=1.d0
 
49
        dp=0.d0
 
50
        mml=m-l
 
51
        do 200 ii=1,mml
 
52
          i=m-ii
 
53
          df=ds*de(i)
 
54
          db=dc*de(i)
 
55
          if(dabs(df).lt.dabs(dg)) goto 150
 
56
          dc=dg/df
 
57
          dr=dsqrt(dc*dc+1.d0)
 
58
          de(i+1)=df*dr
 
59
          ds=1.d0/dr
 
60
          dc=dc*ds
 
61
          goto 160
 
62
  150     ds=df/dg
 
63
          dr=dsqrt(ds*ds+1.d0)
 
64
          de(i+1)=dg*dr
 
65
          dc=1.d0/dr
 
66
          ds=ds*dc
 
67
  160     dg=dzero(i+1)-dp
 
68
          dr=(dzero(i)-dg)*ds+2.d0*dc*db
 
69
          dp=ds*dr
 
70
          dzero(i+1)=dg+dp
 
71
          dg=dc*dr-db
 
72
          df=dweigh(i+1)
 
73
          dweigh(i+1)=ds*dweigh(i)+dc*df
 
74
          dweigh(i)=dc*dweigh(i)-ds*df
 
75
  200   continue
 
76
        dzero(l)=dzero(l)-dp
 
77
        de(l)=dg
 
78
        de(m)=0.d0
 
79
        goto 105
 
80
  240 continue
 
81
      do 300 ii=2,n
 
82
        i=ii-1
 
83
        k=i
 
84
        dp=dzero(i)
 
85
        do 260 j=ii,n
 
86
          if(dzero(j).ge.dp) goto 260
 
87
          k=j
 
88
          dp=dzero(j)
 
89
  260   continue
 
90
        if(k.eq.i) goto 300
 
91
        dzero(k)=dzero(i)
 
92
        dzero(i)=dp
 
93
        dp=dweigh(i)
 
94
        dweigh(i)=dweigh(k)
 
95
        dweigh(k)=dp
 
96
  300 continue
 
97
      do 310 k=1,n
 
98
        dweigh(k)=dbeta(1)*dweigh(k)*dweigh(k)
 
99
  310 continue
 
100
      return
 
101
  400 ierr=l
 
102
      return
 
103
      end
 
104