lsqrtest.f90

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

!     Then form r = Y rbar (again in w).

      w(minmn+1:m) = one
      call Hprod (m,hy,w,w)

!     Compute the rhs    b = r  +  Ax.

      rnorm = dnrm2 (m,w,1)
      call dcopy ( m, w, 1, b, 1 )
      call Aprod1( m, n, maxmn, minmn, x, b, d, hy, hz, w )

      end subroutine lstp

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

      subroutine test  (m,n,nduplc,npower,damp)

      implicit          none
      integer           m, n, nduplc, npower
      double precision  damp

!------------------------------------------------------------------------
!     This is an example driver routine for running LSQR.
!     It generates a test problem, solves it, and examines the results.
!     Note that subroutine Aprod must be declared external
!     if it is used only in the call to LSQR (and acheck).
!
!     1982---1991:  Various versions implemented.
!     04 Sep 1991:  "wantse" added to argument list of LSQR,
!                   making standard errors optional.
!     10 Feb 1992:  Revised for use with lsqrcheck fortran.
!
! 17 Sep 2007: Fortran 90 version.
!------------------------------------------------------------------------

      external          acheck, Aprod, xcheck, dnrm2
      double precision                         dnrm2

      integer, parameter :: maxm=2000,  maxn =2000,  mxmn=2000
      double precision  b(maxm), u(maxm), v(maxn), w(mxmn),   &
                        x(maxn),se(maxn), xtrue(maxn), y(mxmn)
      double precision  eps,atol, btol, conlim, Anorm, Acond, &
                        rnorm, Arnorm, enorm, etol,           &
                        test1, test2, test3, wnorm, xnorm

      integer, parameter :: leniw = 1,  lenrw = 10000
      integer           iw(leniw)
      double precision  rw(lenrw)

      integer           inform, istop, itn, itnlim, j, &
                        maxmn, minmn, nout, nprint
      integer           locd, lochy, lochz, locw, ltotal
      logical           wantse

      double precision, parameter :: zero = 0.0d+0, one = 1.0d+0
      character(len=*), parameter :: line = '----------------------------------'


      if (m > maxm  .or.  n > maxn) go to 800

!     Set the output unit and the machine precision.

      nout   = 6      ! file unit number.  E.g. open(6,'toto.data')
      eps    = 2.22d-16

!     Set the desired solution xtrue.
!     For least-squares problems, this is it.
!     For underdetermined systems, lstp may alter it.

      do j = 1,n
!        xtrue(j) = one
         xtrue(j) = 0.1d+0*j
      end do

!     Generate the specified test problem.
!     The workspace array  iw  is not needed in this application.
!     The workspace array  rw  is used for the following vectors:
!        d(minmn), hy(m), hz(n), w(maxmn).
!     The vectors  d, hy, hz  will define the test matrix A.
!     w is needed for workspace in Aprod1 and Aprod2.

      maxmn  = max(m,n)
      minmn  = min(m,n)
      locd   = 1
      lochy  = locd  + minmn
      lochz  = lochy + m
      locw   = lochz + n
      ltotal = locw  + maxmn - 1
      if (ltotal > lenrw) go to 900

      call lstp  (m,n,maxmn,minmn,nduplc,npower,damp,xtrue, &
                  b,rw(locd),rw(lochy),rw(lochz),rw(locw),Acond,rnorm)

      write(nout,1000) line,line,m,n,nduplc,npower,damp,Acond,rnorm, &
                       line,line

!     Check that Aprod generates y + Ax and x + A'y consistently.

      call Acheck(m,n,nout,Aprod,eps,leniw,lenrw,iw,rw,v,w,x,y,inform)
          
      if (inform > 0) then
         write(nout,'(a)') ' Check eps and power in subroutine Acheck'
         stop
      end if

!     Solve the problem defined by Aprod, damp and b.
!     Copy the rhs vector b into u  (LSQR will overwrite u)
!     and set the other input parameters for LSQR.
!     We ask for standard errors only if they are well-defined.

      call dcopy (m,b,1,u,1)
!---  wantse = m .gt. n  .or.  damp .gt. zero
      wantse = .false.
      atol   = eps**0.99d+0  ! This asks for essentially maximum accuracy
      btol   = atol
      conlim = 1000.d+0 * Acond
      itnlim = 4*(m + n + 50)

      call LSQR  (m,n,Aprod,damp,wantse,         &
                  leniw,lenrw,iw,rw,u,v,w,x,se,  &
                  atol,btol,conlim,itnlim,nout,  &
                  istop,itn,Anorm,Acond,rnorm,Arnorm,xnorm )

      call xcheck(m,n,nout,Aprod,Anorm,damp,eps, &
                  leniw,lenrw,iw,rw,b,u,v,w,x,   &
                  inform,test1,test2,test3)

!     Print the solution and standard error estimates from  LSQR.

      nprint = min(m,n,8)
      write(nout,2500)    (j,x(j),j=1,nprint)
      if ( wantse ) then
         write(nout,2600) (j,se(j),j=1,nprint)
      end if

!     Print a clue about whether the solution looks OK.

      w(1:n) = x(1:n) - xtrue(1:n)
      wnorm  = dnrm2 (n,w,1)
      xnorm  = dnrm2 (n,xtrue,1)
      enorm  = wnorm/(one + xnorm)
      etol   = 1.0d-3
      if (enorm <= etol) then
         write(nout,3000) enorm
      else
         write(nout,3100) enorm
      end if

      return

!     m or n too large.

  800 write(nout,8000)
      return

!     Not enough workspace.

  900 write(nout,9000) ltotal
      return

 1000 format(1p &
      // 1x, 2a &
      /  ' Least-Squares Test Problem      P(', 4i5, e12.2, ' )'     &
      /  ' Condition no. =', e12.4,  '     Residual function =', e17.9 &
      /  1x, 2a)
 2000 format(1p &
      // ' We are solving    min norm(Ax - b)    with no damping.'   &
      // ' Estimates from LSQR:'                                     &
      /  '    norm(x)         =', e10.3, ' = xnorm'                  &
      /  '    norm(r)         =', e10.3, ' = rnorm'                  &
      /  '    norm(A''r)       =', e10.3, ' = Arnorm')
 2100 format(1p &
      // ' We are solving    min norm(Ax - b)    with damp =', e10.3 &
      /  '                           (damp*x)'                       &
      // ' Estimates from LSQR:'                                     &
      /  '    norm(x)         =', e10.3, ' = xnorm'                  &
      /  '    norm(rbar)      =', e10.3, ' = rnorm'                  &
      /  '    norm(Abar''rbar) =', e10.3, ' = Arnorm')
 2500 format(//' Solution  x:' / 4(i6, g14.6))
 2600 format(/ ' Standard errors  se:' / 4(i6, g14.6))
 3000 format(1p &
      // ' LSQR  appears to be successful.'                          &
      /' Relative error in  x  =', e10.2)
 3100 format(1p &
      /' LSQR  appears to have failed.  '                            &
      /' Relative error in  x  =', e10.2)
 8000 format(/ ' XXX  m or n is too large.')
 9000 format(/ ' XXX  Insufficient workspace.',                      &
      '  The length of  rw  should be at least', i6)

      end subroutine test

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

!     MAIN PROGRAM

      program           main

      implicit          none
      external          test
      integer           m,n,nbar,ndamp,nduplc,npower
      double precision  damp

      open(6,file='LSQR.txt',status='unknown')

      nbar   = 1000
      nduplc = 40

      m = 2*nbar
      n = nbar
      do ndamp = 2,7
         npower = ndamp
         damp   = 10.d+0**(-ndamp)
         call test(m,n,nduplc,npower,damp)
      end do

      m = nbar
      n = nbar
      do ndamp = 2, 7
         npower = ndamp
         damp   = 10.d+0**(-ndamp-6)
         call test(m,n,nduplc,npower,damp )
      end do

      m = nbar
      n = 2*nbar
      do ndamp = 2, 7
         npower = ndamp
         damp   = 10.d+0**(-ndamp-6)
         call test(m,n,nduplc,npower,damp )
      end do

      end program main