The R Project SVN R

Rev

Rev 76639 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 76639 Rev 86632
Line 56... Line 56...
56 
c
 56 
c
57 
c     subroutines and functions
 57 
c     subroutines and functions
58 
c
 58 
c
59 
c     linpack dpofa
 59 
c     linpack dpofa
60 
c     blas daxpy,ddot,dscal,dasum
 60 
c     blas daxpy,ddot,dscal,dasum
61 
c     fortran dabs,dmax1,dreal,dsign
 61 
c     fortran abs,max,sign
62 
c
 62 
c
63 
      subroutine dpoco(a,lda,n,rcond,z,info)
 63 
      subroutine dpoco(a,lda,n,rcond,z,info)
64 
      integer lda,n,info
 64 
      integer lda,n,info
65 
      double precision a(lda,n),z(n)
 65 
      double precision a(lda,n),z(n)
66 
      double precision rcond
 66 
      double precision rcond
Line 77... Line 77...
77 
      do 30 j = 1, n
 77 
      do 30 j = 1, n
78 
         z(j) = dasum(j,a(1,j),1)
 78 
         z(j) = dasum(j,a(1,j),1)
79 
         jm1 = j - 1
 79 
         jm1 = j - 1
80 
         if (jm1 .lt. 1) go to 20
 80 
         if (jm1 .lt. 1) go to 20
81 
         do 10 i = 1, jm1
 81 
         do 10 i = 1, jm1
82 
            z(i) = z(i) + dabs(a(i,j))
 82 
            z(i) = z(i) + abs(a(i,j))
83 
   10    continue
 83 
   10    continue
84 
   20    continue
 84 
   20    continue
85 
   30 continue
 85 
   30 continue
86 
      anorm = 0.0d0
 86 
      anorm = 0.0d0
87 
      do 40 j = 1, n
 87 
      do 40 j = 1, n
88 
         anorm = dmax1(anorm,z(j))
 88 
         anorm = max(anorm,z(j))
89 
   40 continue
 89 
   40 continue
90 
c
 90 
c
91 
c     factor
 91 
c     factor
92 
c
 92 
c
93 
      call dpofa(a,lda,n,info)
 93 
      call dpofa(a,lda,n,info)
Line 104... Line 104...
104 
         ek = 1.0d0
 104 
         ek = 1.0d0
105 
         do 50 j = 1, n
 105 
         do 50 j = 1, n
106 
            z(j) = 0.0d0
 106 
            z(j) = 0.0d0
107 
   50    continue
 107 
   50    continue
108 
         do 110 k = 1, n
 108 
         do 110 k = 1, n
109 
            if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
 109 
            if (z(k) .ne. 0.0d0) ek = sign(ek,-z(k))
110 
            if (dabs(ek-z(k)) .le. a(k,k)) go to 60
 110 
            if (abs(ek-z(k)) .le. a(k,k)) go to 60
111 
               s = a(k,k)/dabs(ek-z(k))
 111 
               s = a(k,k)/abs(ek-z(k))
112 
               call dscal(n,s,z,1)
 112 
               call dscal(n,s,z,1)
113 
               ek = s*ek
 113 
               ek = s*ek
114 
   60       continue
 114 
   60       continue
115 
            wk = ek - z(k)
 115 
            wk = ek - z(k)
116 
            wkm = -ek - z(k)
 116 
            wkm = -ek - z(k)
117 
            s = dabs(wk)
 117 
            s = abs(wk)
118 
            sm = dabs(wkm)
 118 
            sm = abs(wkm)
119 
            wk = wk/a(k,k)
 119 
            wk = wk/a(k,k)
120 
            wkm = wkm/a(k,k)
 120 
            wkm = wkm/a(k,k)
121 
            kp1 = k + 1
 121 
            kp1 = k + 1
122 
            if (kp1 .gt. n) go to 100
 122 
            if (kp1 .gt. n) go to 100
123 
               do 70 j = kp1, n
 123 
               do 70 j = kp1, n
124 
                  sm = sm + dabs(z(j)+wkm*a(k,j))
 124 
                  sm = sm + abs(z(j)+wkm*a(k,j))
125 
                  z(j) = z(j) + wk*a(k,j)
 125 
                  z(j) = z(j) + wk*a(k,j)
126 
                  s = s + dabs(z(j))
 126 
                  s = s + abs(z(j))
127 
   70          continue
 127 
   70          continue
128 
               if (s .ge. sm) go to 90
 128 
               if (s .ge. sm) go to 90
129 
                  t = wkm - wk
 129 
                  t = wkm - wk
130 
                  wk = wkm
 130 
                  wk = wkm
131 
                  do 80 j = kp1, n
 131 
                  do 80 j = kp1, n
Line 140... Line 140...
140 
c
 140 
c
141 
c        solve r*y = w
 141 
c        solve r*y = w
142 
c
 142 
c
143 
         do 130 kb = 1, n
 143 
         do 130 kb = 1, n
144 
            k = n + 1 - kb
 144 
            k = n + 1 - kb
145 
            if (dabs(z(k)) .le. a(k,k)) go to 120
 145 
            if (abs(z(k)) .le. a(k,k)) go to 120
146 
               s = a(k,k)/dabs(z(k))
 146 
               s = a(k,k)/abs(z(k))
147 
               call dscal(n,s,z,1)
 147 
               call dscal(n,s,z,1)
148 
  120       continue
 148 
  120       continue
149 
            z(k) = z(k)/a(k,k)
 149 
            z(k) = z(k)/a(k,k)
150 
            t = -z(k)
 150 
            t = -z(k)
151 
            call daxpy(k-1,t,a(1,k),1,z(1),1)
 151 
            call daxpy(k-1,t,a(1,k),1,z(1),1)
Line 157... Line 157...
157 
c
 157 
c
158 
c        solve trans(r)*v = y
 158 
c        solve trans(r)*v = y
159 
c
 159 
c
160 
         do 150 k = 1, n
 160 
         do 150 k = 1, n
161 
            z(k) = z(k) - ddot(k-1,a(1,k),1,z(1),1)
 161 
            z(k) = z(k) - ddot(k-1,a(1,k),1,z(1),1)
162 
            if (dabs(z(k)) .le. a(k,k)) go to 140
 162 
            if (abs(z(k)) .le. a(k,k)) go to 140
163 
               s = a(k,k)/dabs(z(k))
 163 
               s = a(k,k)/abs(z(k))
164 
               call dscal(n,s,z,1)
 164 
               call dscal(n,s,z,1)
165 
               ynorm = s*ynorm
 165 
               ynorm = s*ynorm
166 
  140       continue
 166 
  140       continue
167 
            z(k) = z(k)/a(k,k)
 167 
            z(k) = z(k)/a(k,k)
168 
  150    continue
 168 
  150    continue
Line 172... Line 172...
172 
c
 172 
c
173 
c        solve r*z = v
 173 
c        solve r*z = v
174 
c
 174 
c
175 
         do 170 kb = 1, n
 175 
         do 170 kb = 1, n
176 
            k = n + 1 - kb
 176 
            k = n + 1 - kb
177 
            if (dabs(z(k)) .le. a(k,k)) go to 160
 177 
            if (abs(z(k)) .le. a(k,k)) go to 160
178 
               s = a(k,k)/dabs(z(k))
 178 
               s = a(k,k)/abs(z(k))
179 
               call dscal(n,s,z,1)
 179 
               call dscal(n,s,z,1)
180 
               ynorm = s*ynorm
 180 
               ynorm = s*ynorm
181 
  160       continue
 181 
  160       continue
182 
            z(k) = z(k)/a(k,k)
 182 
            z(k) = z(k)/a(k,k)
183 
            t = -z(k)
 183 
            t = -z(k)