modernize the Fortran

git-svn-id: https://svn.r-project.org/R/trunk@86631 00db46b3-68df-0310-9c12-caf00c1e9a41

src/appl/dqrdc2.f

@@ -1,14 +1,14 @@
 C     Dqrdc2 is a *modification* of Linpack's dqrdc ('DQRDC') for R
 c
-c     dqrdc2 uses householder transformations to compute the qr
-c     factorization of an n by p matrix x.  a limited column
+c     dqrdc2 uses Householder transformations to compute the qr
+c     factorization of an n by p matrix x.  A limited column
 c     pivoting strategy based on the 2-norms of the reduced columns
 c     moves columns with near-zero norm to the right-hand edge of
-c     the x matrix.  this strategy means that sequential one
+c     the x matrix.  This strategy means that sequential one
 c     degree-of-freedom effects can be computed in a natural way.
 c
-c     i am very nervous about modifying linpack code in this way.
-c     if you are a computational linear algebra guru and you really
+c     I am very nervous about modifying linpack code in this way.
+c     If you are a computational linear algebra guru and you really
 c     understand how to solve this problem please feel free to
 c     suggest improvements to this code.
 c
@@ -68,13 +68,14 @@
 c		 permutation of 1:p; it is called 'pivot' in R
 
 c
-c     original (dqrdc.f) linpack version dated 08/14/78 .
-c     g.w. stewart, university of maryland, argonne national lab.
+c     Original (dqrdc.f) linpack version dated 08/14/78 .
+c     G.W. Stewart, University of Maryland, Argonne National Lab.
 c
 C     This version dated 22 August 1995
 C     Ross Ihaka
 c
-c     bug fixes 29 September 1999 BDR (p > n case, inaccurate ranks)
+c     Bug fixes 29 September 1999 BDR (p > n case, inaccurate ranks)
+c     Fortran modernized to F77, 2024-05 BDR
 c
 c
 c     dqrdc2 uses the following functions and subprograms.
@@ -89,7 +90,7 @@ subroutine dqrdc2(x,ldx,n,p,tol,k,qraux,jpvt,work)
 c
 c     internal variables
 c
-      integer i,j,l,lp1,lup,k
+      integer i,j,l,lup,k
       double precision dnrm2,tt,ttt
       double precision ddot,nrmxl,t
 c
@@ -98,19 +99,19 @@ subroutine dqrdc2(x,ldx,n,p,tol,k,qraux,jpvt,work)
 c
       if (n .gt. 0) then
 c       avoid accessing element beyond the bound
-         do 70 j = 1, p
+         do j = 1, p
             qraux(j) = dnrm2(n,x(1,j),1)
             work(j,1) = qraux(j)
             work(j,2) = qraux(j)
             if(work(j,2) .eq. 0.0d0) work(j,2) = 1.0d0
-  70     continue
-      endif
+         end do
+      end if
 c
 c     perform the householder reduction of x.
 c
-      lup = min0(n,p)
+      lup = min(n,p)
       k = p + 1
-      do 200 l = 1, lup
+      do l = 1, lup
 c
 c     previous version only cycled l to lup
 c
@@ -120,81 +121,73 @@ subroutine dqrdc2(x,ldx,n,p,tol,k,qraux,jpvt,work)
 c     tol times its original norm.  the check for l .le. k
 c     avoids infinite cycling.
 c
-   80    continue
-         if (l .ge. k .or. qraux(l) .ge. work(l,2)*tol) go to 120
-            lp1 = l+1
-            do 100 i=1,n
+         do
+            if (l .ge. k .or. qraux(l) .ge. work(l,2)*tol) exit
+            do i=1,n
                t = x(i,l)
-               do 90 j=lp1,p
+               do j=l+1,p
                   x(i,j-1) = x(i,j)
-   90          continue
+               end do
                x(i,p) = t
-  100       continue
+            end do
             i = jpvt(l)
             t = qraux(l)
             tt = work(l,1)
             ttt = work(l,2)
-            do 110 j=lp1,p
+            do j=l+1,p
                jpvt(j-1) = jpvt(j)
                qraux(j-1) = qraux(j)
                work(j-1,1) = work(j,1)
                work(j-1,2) = work(j,2)
-  110       continue
+            end do
             jpvt(p) = i
             qraux(p) = t
             work(p,1) = tt
             work(p,2) = ttt
             k = k - 1
-            go to 80
-  120    continue
-         if (l .eq. n) go to 190
+         end do
+         if (l .ne. n) then
 c
 c           compute the householder transformation for column l.
 c
             nrmxl = dnrm2(n-l+1,x(l,l),1)
-            if (nrmxl .eq. 0.0d0) go to 180
-               if (x(l,l) .ne. 0.0d0) nrmxl = dsign(nrmxl,x(l,l))
+            if (nrmxl .ne. 0.0d0) then
+               if (x(l,l) .ne. 0.0d0) nrmxl = sign(nrmxl,x(l,l))
                call dscal(n-l+1,1.0d0/nrmxl,x(l,l),1)
                x(l,l) = 1.0d0 + x(l,l)
 c
 c              apply the transformation to the remaining columns,
 c              updating the norms.
 c
-               lp1 = l + 1
-               if (p .lt. lp1) go to 170
-               do 160 j = lp1, p
-                  t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
-                  call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
-                  if (qraux(j) .eq. 0.0d0) go to 150
-                     tt = 1.0d0 - (dabs(x(l,j))/qraux(j))**2
-                     tt = dmax1(tt,0.0d0)
-                     t = tt
-c
-c modified 9/99 by BDR. Re-compute norms if there is large reduction
+               if (p .ge. l+1) then
+                  do j = l+1, p
+                     t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
+                     call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
+                     if (qraux(j) .ne. 0.0d0) then
+                        tt = 1.0d0 - (abs(x(l,j))/qraux(j))**2
+                        tt = max(tt,0.0d0)
+                        t = tt
+c
+c Modified 9/99 by BDR. Re-compute norms if there is large reduction
 c The tolerance here is on the squared norm
 c In this version we need accurate norms, so re-compute often.
-c  work(j,1) is only updated in one case: looks like a bug -- no longer used
-c
-c                     tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j,1))**2
-c                     if (tt .eq. 1.0d0) go to 130
-                     if (dabs(t) .lt. 1d-6) go to 130
-                        qraux(j) = qraux(j)*dsqrt(t)
-                     go to 140
-  130                continue
-                        qraux(j) = dnrm2(n-l,x(l+1,j),1)
-                        work(j,1) = qraux(j)
-  140                continue
-  150             continue
-  160          continue
-  170          continue
-c
-c              save the transformation.
+                        if (dabs(t) .ge. 1d-6) then
+                           qraux(j) = qraux(j)*sqrt(t)
+                        else
+                           qraux(j) = dnrm2(n-l,x(l+1,j),1)
+                           work(j,1) = qraux(j)
+                        end if
+                     end if
+                  end do
+               end if
+c
+c              Save the transformation.
 c
                qraux(l) = x(l,l)
                x(l,l) = -nrmxl
-  180       continue
-  190    continue
-  200 continue
-      k = min0(k - 1, n)
+            end if
+         end if
+      end do
+      k = min(k - 1, n)
       return
       end