lsqrtest.f90

比较经典的求解线性方程的方法 原理是C.C. Paige and M.A. Sauders等你提出的

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     File lsqrtest.f90  (double precision)
!
!     Aprod    Aprod1   Aprod2   Hprod    lstp     test     MAIN
!
!     These routines define a class of least-squares test problems
!     for testing algorithms LSQR and CRAIG
!     (Paige and Saunders, ACM TOMS, 1982).
!
!     1982---1991:  Various versions implemented.
!     06 Feb 1992:  Test-problem generator lstp generalized to allow
!                   any m and n.  lstp is now the same as the generator
!                   for LSQR and CRAIG.
!     30 Nov 1993:  Modified lstp.
!                   For a while, damp = 0 implied r = damp*s = 0.
!                   This was a result of generating x and s.
!                   Reverted to generating x and r as in LSQR paper.
! 07 Sep 2007: Line by line translation for Fortran90 compilers
!              by Eric Badel <[e-mail]>.
! 16 Sep 2007: Polished a bit by Michael Saunders.
!
! Maintained by Michael Saunders <[e-mail]>.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Aprod (mode,m,n,x,y,leniw,lenrw,iw,rw)

      implicit          none
      integer           mode, m, n, leniw, lenrw
      integer           iw(leniw)
      double precision  x(n), y(m), rw(lenrw)

!------------------------------------------------------------------------
!     This is the matrix-vector product routine required by subroutines
!     LSQR and CRAIG for a test matrix of the form  A = HY*D*HZ.
!     The quantities defining D, HY, HZ are in the work array rw,
!     followed by a work array w.  These are passed to Aprod1 and Aprod2
!     in order to make the code readable.
!------------------------------------------------------------------------

      external          Aprod1, Aprod2
      integer           locd, lochy, lochz, locw, maxmn, minmn


      maxmn  = max(m,n)
      minmn  = min(m,n)
      locd   = 1
      lochy  = locd+minmn
      lochz  = lochy+m
      locw   = lochz+n

      if (mode == 1) then
         call Aprod1(m,n,maxmn,minmn,x,y,rw(locd),rw(lochy),rw(lochz),rw(locw))
      else
         call Aprod2(m,n,maxmn,minmn,x,y,rw(locd),rw(lochy),rw(lochz),rw(locw))
      end if

      end subroutine Aprod

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Aprod1(m,n,maxmn,minmn,x,y,d,hy,hz,w)

      implicit          none
      integer           m, n, maxmn, minmn
      double precision  x(n),y(m),d(minmn),hy(m),hz(n),w(maxmn)

!------------------------------------------------------------------------
!     Aprod1  computes  y = y + A*x  for subroutine Aprod,
!     where A is a test matrix of the form  A = HY*D*HZ,
!     and the latter matrices HY, D, HZ are represented by
!     input vectors with the same name.
!
! 17 Sep 2007: Fortran 90 version.
!------------------------------------------------------------------------

      external          Hprod
      double precision, parameter :: zero = 0.0d+0


      call Hprod (n,hz,x,w)
      w(1:minmn) = d(1:minmn)*w(1:minmn)
      w(n+1:m)   = zero
      call Hprod (m,hy,w,w)
      y = y + w(1:m)

      end subroutine Aprod1

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Aprod2(m,n,maxmn,minmn,x,y,d,hy,hz,w)

      integer           m, n, maxmn, minmn
      double precision  x(n),y(m),d(minmn),hy(m),hz(n),w(maxmn)

!-----------------------------------------------------------------------
!     Aprod2  computes  x = x + A(t)*y  for subroutine Aprod,
!     where  A  is a test matrix of the form  A = HY*D*HZ,
!     and the latter matrices  HY, D, HZ  are represented by
!     input vectors with the same name.
!
! 17 Sep 2007: Fortran 90 version.
!-----------------------------------------------------------------------

      external          Hprod
      double precision, parameter :: zero = 0.0d+0


      call Hprod (m,hy,y,w)
      w(1:minmn) = d(1:minmn)*w(1:minmn)
      w(m+1:n)   = zero
      call Hprod (n,hz,w,w)
      x = x + w(1:n)

      end subroutine Aprod2

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hprod (n,hz,x,y)

      integer           n
      double precision  hz(n), x(n), y(n)

!-----------------------------------------------------------------------
!     Hprod  applies a Householder transformation stored in  hz
!     to get  y = ( I - 2*hz*hz(transpose) ) * x.
!
! 17 Sep 2007: Fortran 90 version.
!-----------------------------------------------------------------------

      integer            i
      double precision   s
      double precision, parameter :: zero = 0.0d+0

      s = zero
      do i = 1,n
         s = hz(i)*x(i) + s
      end do

      s = s + s
      y = x - s*hz

      end subroutine Hprod

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lstp  (m,n,maxmn,minmn,nduplc,npower,damp,x, &
                        b,d,hy,hz,w,Acond,rnorm)

      implicit          none
      integer           m, n, maxmn, minmn, nduplc, npower
      double precision  damp, Acond, rnorm
      double precision  b(m), x(n), d(minmn), hy(m), hz(n), w(maxmn)

!-----------------------------------------------------------------------
!     lstp  generates a sparse least-squares test problem of the form
!                (   A    )*x = ( b ) 
!                ( damp*I )     ( 0 )
!     for solution by LSQR, or a sparse underdetermined system
!                   Ax + damp*s = b
!     for solution by CRAIG.  The matrix A is m by n and is
!     constructed in the form  A = HY*D*HZ,  where D is an m by n
!     diagonal matrix, and HY and HZ are Householder transformations.
!
!     m and n may contain any positive values.
!     The first 8 parameters are input to lstp.  The last 8 are output.
!     If m >= n  or  damp = 0, the true solution is x as given.
!     Otherwise, x is modified to contain the true solution.
!
!     Functions and subroutines
!
!     testprob           Aprod1, Hprod
!     blas               dnrm2 , dscal
!
! 17 Sep 2007: Fortran 90 version.
!-----------------------------------------------------------------------

      double precision   dnrm2
      external           dnrm2, dcopy, dscal, Aprod1, Hprod
      intrinsic          cos, min, sin, sqrt
      integer            i, j
      double precision   alfa,beta,dampsq,fourpi,t
      double precision, parameter :: zero = 0.0d+0, one = 1.0d+0

!     ------------------------------------------------------------------
!     Make two vectors of norm 1.0 for the Householder transformations.
!     fourpi  need not be exact.
!     ------------------------------------------------------------------
      minmn  = min(m,n)
      dampsq = damp**2
      fourpi = 4.0 * 3.141592
      alfa   = fourpi / m
      beta   = fourpi / n

      do i = 1,m
         hy(i) = sin( alfa*i )
      end do

      do i = 1, n
         hz(i) = cos( beta*i )
      end do

      alfa   = dnrm2 ( m, hy, 1 )
      beta   = dnrm2 ( n, hz, 1 )
      call dscal ( m, (- one / alfa), hy, 1 )
      call dscal ( n, (- one / beta), hz, 1 )
!
!     ------------------------------------------------------------------
!     Set the diagonal matrix  D.  These are the singular values of  A.
!     ------------------------------------------------------------------
      do i = 1, minmn
         j     = (i - 1 + nduplc) / nduplc
         t     =  one*j*nduplc
         t     =  t / minmn
         d(i)  =  t**npower
      end do

      Acond  = (d(minmn)**2 + dampsq) / (d(1)**2 + dampsq)
      Acond  = sqrt( Acond )

!     ------------------------------------------------------------------
!     Set the true solution   x.
!     It must be of the form  x = Z ( w )  for some  w.
!                                   ( 0 )
!     ------------------------------------------------------------------
      call Hprod (n,hz,x,w)
      w(m+1:n) = zero
      call Hprod (n,hz,w,x)

!     Solve D r1bar = dampsq x1bar 
!     where r1bar and x1bar are both in w.

      w(1:minmn) = dampsq * w(1:minmn)/d(1:minmn)

!     Set r2bar to be anything.  (It is empty if m <= n)