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 98... Line 98...
98 
c
 98 
c
99 
c     dsvdc uses the following functions and subprograms.
 99 
c     dsvdc uses the following functions and subprograms.
100 
c
 100 
c
101 
c     external drot
 101 
c     external drot
102 
c     blas daxpy,ddot,dscal,dswap,dnrm2,drotg
 102 
c     blas daxpy,ddot,dscal,dswap,dnrm2,drotg
103 
c     fortran dabs,dmax1,max0,min0,mod,dsqrt
 103 
c     fortran abs,max,min,mod,sqrt
104 
c
 104 
c
105 
      subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
 105 
      subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
106 
      integer ldx,n,p,ldu,ldv,job,info
 106 
      integer ldx,n,p,ldu,ldv,job,info
107 
      double precision x(ldx,p),s(min(n+1,p)),e(p),
 107 
      double precision x(ldx,p),s(min(n+1,p)),e(p),
108 
     + u(ldu,n),v(ldv,p),work(n)
 108 
     + u(ldu,n),v(ldv,p),work(n)
Line 130... Line 130...
130 
c
 130 
c
131 
      wantu = .false.
 131 
      wantu = .false.
132 
      wantv = .false.
 132 
      wantv = .false.
133 
      jobu = mod(job,100)/10
 133 
      jobu = mod(job,100)/10
134 
      ncu = n
 134 
      ncu = n
135 
      if (jobu .gt. 1) ncu = min0(n,p)
 135 
      if (jobu .gt. 1) ncu = min(n,p)
136 
      if (jobu .ne. 0) wantu = .true.
 136 
      if (jobu .ne. 0) wantu = .true.
137 
      if (mod(job,10) .ne. 0) wantv = .true.
 137 
      if (mod(job,10) .ne. 0) wantv = .true.
138 
c
 138 
c
139 
c     reduce x to bidiagonal form, storing the diagonal elements
 139 
c     reduce x to bidiagonal form, storing the diagonal elements
140 
c     in s and the super-diagonal elements in e.
 140 
c     in s and the super-diagonal elements in e.
141 
c
 141 
c
142 
      info = 0
 142 
      info = 0
143 
      nct = min0(n-1,p)
 143 
      nct = min(n-1,p)
144 
      nrt = max0(0,min0(p-2,n))
 144 
      nrt = max(0,min(p-2,n))
145 
      lu = max0(nct,nrt)
 145 
      lu = max(nct,nrt)
146 
      if (lu .lt. 1) go to 170
 146 
      if (lu .lt. 1) go to 170
147 
      do 160 l = 1, lu
 147 
      do 160 l = 1, lu
148 
         lp1 = l + 1
 148 
         lp1 = l + 1
149 
         if (l .gt. nct) go to 20
 149 
         if (l .gt. nct) go to 20
150 
c
 150 
c
Line 224... Line 224...
224 
  160 continue
 224 
  160 continue
225 
  170 continue
 225 
  170 continue
226 
c
 226 
c
227 
c     set up the final bidiagonal matrix or order m.
 227 
c     set up the final bidiagonal matrix or order m.
228 
c
 228 
c
229 
      m = min0(p,n+1)
 229 
      m = min(p,n+1)
230 
      nctp1 = nct + 1
 230 
      nctp1 = nct + 1
231 
      nrtp1 = nrt + 1
 231 
      nrtp1 = nrt + 1
232 
      if (nct .lt. p) s(nctp1) = x(nctp1,nctp1)
 232 
      if (nct .lt. p) s(nctp1) = x(nctp1,nctp1)
233 
      if (n .lt. m) s(m) = 0.0d0
 233 
      if (n .lt. m) s(m) = 0.0d0
234 
      if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m)
 234 
      if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m)
Line 327... Line 327...
327 
c
 327 
c
328 
         do 390 ll = 1, m
 328 
         do 390 ll = 1, m
329 
            l = m - ll
 329 
            l = m - ll
330 
c        ...exit
 330 
c        ...exit
331 
            if (l .eq. 0) go to 400
 331 
            if (l .eq. 0) go to 400
332 
            test = dabs(s(l)) + dabs(s(l+1))
 332 
            test = abs(s(l)) + abs(s(l+1))
333 
            ztest = test + dabs(e(l))
 333 
            ztest = test + abs(e(l))
334 
            acc = dabs(test - ztest)/(1.0d-100 + test)
 334 
            acc = abs(test - ztest)/(1.0d-100 + test)
335 
            if (acc .gt. 1.d-15) goto 380
 335 
            if (acc .gt. 1.d-15) goto 380
336 
c            if (ztest .ne. test) go to 380
 336 
c            if (ztest .ne. test) go to 380
337 
               e(l) = 0.0d0
 337 
               e(l) = 0.0d0
338 
c        ......exit
 338 
c        ......exit
339 
               go to 400
 339 
               go to 400
Line 349... Line 349...
349 
            do 430 lls = lp1, mp1
 349 
            do 430 lls = lp1, mp1
350 
               ls = m - lls + lp1
 350 
               ls = m - lls + lp1
351 
c           ...exit
 351 
c           ...exit
352 
               if (ls .eq. l) go to 440
 352 
               if (ls .eq. l) go to 440
353 
               test = 0.0d0
 353 
               test = 0.0d0
354 
               if (ls .ne. m) test = test + dabs(e(ls))
 354 
               if (ls .ne. m) test = test + abs(e(ls))
355 
               if (ls .ne. l + 1) test = test + dabs(e(ls-1))
 355 
               if (ls .ne. l + 1) test = test + abs(e(ls-1))
356 
               ztest = test + dabs(s(ls))
 356 
               ztest = test + abs(s(ls))
357 
c 1.0d-100 is to guard against a zero matrix, hence zero test
 357 
c 1.0d-100 is to guard against a zero matrix, hence zero test
358 
               acc = dabs(test - ztest)/(1.0d-100 + test)
 358 
               acc = abs(test - ztest)/(1.0d-100 + test)
359 
               if (acc .gt. 1.d-15) goto 420
 359 
               if (acc .gt. 1.d-15) goto 420
360 
c               if (ztest .ne. test) go to 420
 360 
c               if (ztest .ne. test) go to 420
361 
                  s(ls) = 0.0d0
 361 
                  s(ls) = 0.0d0
362 
c           ......exit
 362 
c           ......exit
363 
                  go to 440
 363 
                  go to 440
Line 430... Line 430...
430 
c
 430 
c
431 
  540    continue
 431 
  540    continue
432 
c
 432 
c
433 
c           calculate the shift.
 433 
c           calculate the shift.
434 
c
 434 
c
435 
            scale = dmax1(dabs(s(m)),dabs(s(m-1)),dabs(e(m-1)),
 435 
            scale = max(abs(s(m)),abs(s(m-1)),abs(e(m-1)),
436 
     *                    dabs(s(l)),dabs(e(l)))
 436 
     *                  abs(s(l)),dabs(e(l)))
437 
            sm = s(m)/scale
 437 
            sm = s(m)/scale
438 
            smm1 = s(m-1)/scale
 438 
            smm1 = s(m-1)/scale
439 
            emm1 = e(m-1)/scale
 439 
            emm1 = e(m-1)/scale
440 
            sl = s(l)/scale
 440 
            sl = s(l)/scale
441 
            el = e(l)/scale
 441 
            el = e(l)/scale
442 
            b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0
 442 
            b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0
443 
            c = (sm*emm1)**2
 443 
            c = (sm*emm1)**2
444 
            shift = 0.0d0
 444 
            shift = 0.0d0
445 
            if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550
 445 
            if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550
446 
               shift = dsqrt(b**2+c)
 446 
               shift = sqrt(b**2+c)
447 
               if (b .lt. 0.0d0) shift = -shift
 447 
               if (b .lt. 0.0d0) shift = -shift
448 
               shift = c/(b + shift)
 448 
               shift = c/(b + shift)
449 
  550       continue
 449 
  550       continue
450 
            f = (sl + sm)*(sl - sm) + shift
 450 
            f = (sl + sm)*(sl - sm) + shift
451 
            g = sl*el
 451 
            g = sl*el