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 48... Line 48...
48 
c     cleve moler, university of new mexico, argonne national lab.
 48 
c     cleve moler, university of new mexico, argonne national lab.
49 
c
 49 
c
50 
c     subroutines and functions
 50 
c     subroutines and functions
51 
c
 51 
c
52 
c     blas daxpy,dscal,dasum
 52 
c     blas daxpy,dscal,dasum
53 
c     fortran dabs,dmax1,dsign
 53 
c     fortran abs,max,sign
54 
c
 54 
c
55 
c     internal variables
 55 
c     internal variables
56 
c
 56 
c
57 
      double precision w,wk,wkm,ek
 57 
      double precision w,wk,wkm,ek
58 
      double precision tnorm,ynorm,s,sm,dasum
 58 
      double precision tnorm,ynorm,s,sm,dasum
Line 67... Line 67...
67 
      do 10 j = 1, n
 67 
      do 10 j = 1, n
68 
         l = j
 68 
         l = j
69 
         if (lower) l = n + 1 - j
 69 
         if (lower) l = n + 1 - j
70 
         i1 = 1
 70 
         i1 = 1
71 
         if (lower) i1 = j
 71 
         if (lower) i1 = j
72 
         tnorm = dmax1(tnorm,dasum(l,t(i1,j),1))
 72 
         tnorm = max(tnorm,dasum(l,t(i1,j),1))
73 
   10 continue
 73 
   10 continue
74 
c
 74 
c
75 
c     rcond = 1/(norm(t)*(estimate of norm(inverse(t)))) .
 75 
c     rcond = 1/(norm(t)*(estimate of norm(inverse(t)))) .
76 
c     estimate = norm(z)/norm(y) where  t*z = y  and  trans(t)*y = e .
 76 
c     estimate = norm(z)/norm(y) where  t*z = y  and  trans(t)*y = e .
77 
c     trans(t)  is the transpose of t .
 77 
c     trans(t)  is the transpose of t .
Line 86... Line 86...
86 
         z(j) = 0.0d0
 86 
         z(j) = 0.0d0
87 
   20 continue
 87 
   20 continue
88 
      do 100 kk = 1, n
 88 
      do 100 kk = 1, n
89 
         k = kk
 89 
         k = kk
90 
         if (lower) k = n + 1 - kk
 90 
         if (lower) k = n + 1 - kk
91 
         if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
 91 
         if (z(k) .ne. 0.0d0) ek = sign(ek,-z(k))
92 
         if (dabs(ek-z(k)) .le. dabs(t(k,k))) go to 30
 92 
         if (abs(ek-z(k)) .le. abs(t(k,k))) go to 30
93 
            s = dabs(t(k,k))/dabs(ek-z(k))
 93 
            s = abs(t(k,k))/abs(ek-z(k))
94 
            call dscal(n,s,z,1)
 94 
            call dscal(n,s,z,1)
95 
            ek = s*ek
 95 
            ek = s*ek
96 
   30    continue
 96 
   30    continue
97 
         wk = ek - z(k)
 97 
         wk = ek - z(k)
98 
         wkm = -ek - z(k)
 98 
         wkm = -ek - z(k)
99 
         s = dabs(wk)
 99 
         s = abs(wk)
100 
         sm = dabs(wkm)
 100 
         sm = abs(wkm)
101 
         if (t(k,k) .eq. 0.0d0) go to 40
 101 
         if (t(k,k) .eq. 0.0d0) go to 40
102 
            wk = wk/t(k,k)
 102 
            wk = wk/t(k,k)
103 
            wkm = wkm/t(k,k)
 103 
            wkm = wkm/t(k,k)
104 
         go to 50
 104 
         go to 50
105 
   40    continue
 105 
   40    continue
Line 110... Line 110...
110 
            j1 = k + 1
 110 
            j1 = k + 1
111 
            if (lower) j1 = 1
 111 
            if (lower) j1 = 1
112 
            j2 = n
 112 
            j2 = n
113 
            if (lower) j2 = k - 1
 113 
            if (lower) j2 = k - 1
114 
            do 60 j = j1, j2
 114 
            do 60 j = j1, j2
115 
               sm = sm + dabs(z(j)+wkm*t(k,j))
 115 
               sm = sm + abs(z(j)+wkm*t(k,j))
116 
               z(j) = z(j) + wk*t(k,j)
 116 
               z(j) = z(j) + wk*t(k,j)
117 
               s = s + dabs(z(j))
 117 
               s = s + abs(z(j))
118 
   60       continue
 118 
   60       continue
119 
            if (s .ge. sm) go to 80
 119 
            if (s .ge. sm) go to 80
120 
               w = wkm - wk
 120 
               w = wkm - wk
121 
               wk = wkm
 121 
               wk = wkm
122 
               do 70 j = j1, j2
 122 
               do 70 j = j1, j2
Line 134... Line 134...
134 
c     solve t*z = y
 134 
c     solve t*z = y
135 
c
 135 
c
136 
      do 130 kk = 1, n
 136 
      do 130 kk = 1, n
137 
         k = n + 1 - kk
 137 
         k = n + 1 - kk
138 
         if (lower) k = kk
 138 
         if (lower) k = kk
139 
         if (dabs(z(k)) .le. dabs(t(k,k))) go to 110
 139 
         if (abs(z(k)) .le. abs(t(k,k))) go to 110
140 
            s = dabs(t(k,k))/dabs(z(k))
 140 
            s = abs(t(k,k))/abs(z(k))
141 
            call dscal(n,s,z,1)
 141 
            call dscal(n,s,z,1)
142 
            ynorm = s*ynorm
 142 
            ynorm = s*ynorm
143 
  110    continue
 143 
  110    continue
144 
         if (t(k,k) .ne. 0.0d0) z(k) = z(k)/t(k,k)
 144 
         if (t(k,k) .ne. 0.0d0) z(k) = z(k)/t(k,k)
145 
         if (t(k,k) .eq. 0.0d0) z(k) = 1.0d0
 145 
         if (t(k,k) .eq. 0.0d0) z(k) = 1.0d0