dhfti – bspline-fortran

private subroutine dhfti(a, mda, m, n, b, mdb, nb, tau, krank, rnorm, h, g, ip)

Solve a least squares problem for banded matrices using sequential accumulation of rows of the data matrix. Exactly one right-hand side vector is permitted.

This subroutine solves a linear least squares problem or a set of linear least squares problems having the same matrix but different right-side vectors. The problem data consists of an M by N matrix A, an M by NB matrix B, and an absolute tolerance parameter TAU whose usage is described below. The NB column vectors of B represent right-side vectors for NB distinct linear least squares problems.

This set of problems can also be written as the matrix least squares problem

A = B,

where X is the N by NB solution matrix.

Note that if B is the M by M identity matrix, then X will be the pseudo-inverse of A.

This subroutine first transforms the augmented matrix (A B) to a matrix (R C) using premultiplying Householder transformations with column interchanges. All subdiagonal elements in the matrix R are zero and its diagonal elements satisfy

  abs(r(i,i))>=abs(r(i+1,i+1)),
  i = 1,...,l-1, where
  l = min(m,n).

The subroutine will compute an integer, KRANK, equal to the number of diagonal terms of R that exceed TAU in magnitude. Then a solution of minimum Euclidean length is computed using the first KRANK rows of (R C).

To be specific we suggest that the user consider an easily computable matrix norm, such as, the maximum of all column sums of magnitudes.

Now if the relative uncertainty of B is EPS, (norm of uncertainty/ norm of B), it is suggested that TAU be set approximately equal to EPS*(norm of A).

References

C. L. Lawson and R. J. Hanson, Solving Least Squares Problems, Prentice-Hall, Inc., 1974, Chapter 14.

Revision history

790101 DATE WRITTEN. Lawson, C. L., (JPL), Hanson, R. J., (SNLA)
890531 Changed all specific intrinsics to generic. (WRB)
891006 Cosmetic changes to prologue. (WRB)
891006 REVISION DATE from Version 3.2
891214 Prologue converted to Version 4.0 format. (BAB)
900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
901005 Replace usage of DDIFF with usage of D1MACH. (RWC)
920501 Reformatted the REFERENCES section. (WRB)

Arguments

Type	Intent		Name
real(kind=wp),	intent(inout)	::	a(mda,*)	`A(MDA,N)`. The array A(,) initially contains the M by N matrix A of the least squares problem AX = B. The first dimensioning parameter of the array A(,) is MDA, which must satisfy MDA>=M Either M>=N or M<N is permitted. There is no restriction on the rank of A. The condition MDA<M is considered an error. The contents of the array A(,) will be modified by the subroutine. These contents are not generally required by the user.
integer,	intent(in)	::	mda	actual leading dimension of `a`
integer,	intent(in)	::	m
integer,	intent(in)	::	n
real(kind=wp),	intent(inout)	::	b(mdb,*)	`(B(MDB,NB) or B(M))`. If NB = 0 the subroutine will perform the orthogonal decomposition but will make no references to the array B(). If NB>0 the array B() must initially contain the M by NB matrix B of the least squares problem AX = B. If NB>=2 the array B() must be doubly subscripted with first dimensioning parameter MDB>=MAX(M,N). If NB = 1 the array B() may be either doubly or singly subscripted. In the latter case the value of MDB is arbitrary but it should be set to some valid integer value such as MDB = M. The condition of NB>1.AND.MDB< MAX(M,N) is considered an error. On return the array B(*) will contain the N by NB solution matrix X.
integer,	intent(in)	::	mdb	actual leading dimension of `b`
integer,	intent(in)	::	nb
real(kind=wp),	intent(in)	::	tau	Absolute tolerance parameter provided by user for pseudorank determination.
integer,	intent(out)	::	krank	Set by the subroutine to indicate the pseudorank of A.
real(kind=wp),	intent(out)	::	rnorm(*)	`RNORM(NB)`. On return, RNORM(J) will contain the Euclidean norm of the residual vector for the problem defined by the J-th column vector of the array B(,) for J = 1,...,NB.
real(kind=wp)		::	h(*)	`H(N)`. Array of working space used by DHFTI. On return, contains elements of the pre-multiplying Householder transformations used to compute the minimum Euclidean length solution. not generally required by the user.
real(kind=wp)		::	g(*)	`G(N)`. Array of working space used by DHFTI. On return, contain elements of the post-multiplying Householder transformations used to compute the minimum Euclidean length solution. not generally required by the user.
integer		::	ip(*)	`IP(N)`. Array of working space used by DHFTI. Array in which the subroutine records indices describing the permutation of column vectors. not generally required by the user.

Calls

Help

Called by

Help

Source Code

   subroutine dhfti (a, mda, m, n, b, mdb, nb, tau, krank, rnorm, h, g, ip)

   integer,intent(in) :: mda !! actual leading dimension of `a`
   integer,intent(in) :: mdb !! actual leading dimension of `b`
   real(wp),intent(inout) :: a(mda,*) !! `A(MDA,N)`.
                                      !! The array A(*,*) initially contains the M by N
                                      !! matrix A of the least squares problem AX = B.
                                      !! The first dimensioning parameter of the array
                                      !! A(*,*) is MDA, which must satisfy MDA>=M
                                      !! Either M>=N or M<N is permitted.  There
                                      !! is no restriction on the rank of A.  The
                                      !! condition MDA<M is considered an error.
                                      !!
                                      !! The contents of the array A(*,*) will be
                                      !! modified by the subroutine. These contents
                                      !! are not generally required by the user.
   integer,intent(in) :: m
   integer,intent(in) :: n
   real(wp),intent(inout) :: b(mdb,*) !! `(B(MDB,NB) or B(M))`.
                                      !! If NB = 0 the subroutine will perform the
                                      !! orthogonal decomposition but will make no
                                      !! references to the array B(*).  If NB>0
                                      !! the array B(*) must initially contain the M by
                                      !! NB matrix B of the least squares problem AX =
                                      !! B.  If NB>=2 the array B(*) must be doubly
                                      !! subscripted with first dimensioning parameter
                                      !! MDB>=MAX(M,N).  If NB = 1 the array B(*) may
                                      !! be either doubly or singly subscripted.  In
                                      !! the latter case the value of MDB is arbitrary
                                      !! but it should be set to some valid integer
                                      !! value such as MDB = M.
                                      !!
                                      !! The condition of NB>1.AND.MDB< MAX(M,N)
                                      !! is considered an error.
                                      !!
                                      !! On return the array B(*) will contain the N by
                                      !! NB solution matrix X.
   integer,intent(in) :: nb
   real(wp),intent(in) :: tau !! Absolute tolerance parameter provided by user
                              !! for pseudorank determination.
   integer,intent(out) :: krank !! Set by the subroutine to indicate the
                                !! pseudorank of A.
   real(wp),intent(out)  :: rnorm(*) !! `RNORM(NB)`.
                                     !! On return, RNORM(J) will contain the Euclidean
                                     !! norm of the residual vector for the problem
                                     !! defined by the J-th column vector of the array
                                     !! B(*,*) for J = 1,...,NB.
   real(wp) :: h(*) !! `H(N)`. Array of working space used by DHFTI.
                    !! On return, contains
                    !! elements of the pre-multiplying
                    !! Householder transformations used to compute
                    !! the minimum Euclidean length solution.
                    !! not generally required by the user.
   real(wp) :: g(*) !! `G(N)`. Array of working space used by DHFTI.
                    !! On return, contain
                    !! elements of the post-multiplying
                    !! Householder transformations used to compute
                    !! the minimum Euclidean length solution.
                    !! not generally required by the user.
   integer :: ip(*) !! `IP(N)`. Array of working space used by DHFTI.
                    !! Array in which the subroutine records indices
                    !! describing the permutation of column vectors.
                    !! not generally required by the user.

   integer :: i, ii, iopt, ip1, j, jb, jj, k, kp1, l, ldiag, lmax, nerr
   real(wp) :: dzero, factor, hmax, sm, sm1, szero, tmp
   logical :: lmax_found

   szero = 0.0_wp
   dzero = 0.0_wp
   factor = 0.001_wp

   k = 0
   ldiag = min(m,n)
   if (ldiag > 0) then

      if (mda < m) then
         nerr = 1
         iopt = 2
         write(*,*) 'MDA<M, PROBABLE ERROR.'
         return
      end if
      if (nb > 1 .and. max(m,n) > mdb) then
         nerr = 2
         iopt = 2
         write(*,*) 'MDB<MAX(M,N).AND.NB>1. PROBABLE ERROR.'
         return
      end if

      do j = 1, ldiag
         lmax_found = .false.
         if (j /= 1) then
            ! UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
            lmax = j
            do l = j, n
               h(l) = h(l) - a(j-1,l)**2
               if (h(l) > h(lmax)) lmax = l
            end do
            lmax_found = (factor*h(lmax) > hmax*drelpr)
         end if
         if (.not. lmax_found) then
            ! COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
            lmax = j
            do l = j, n
               h(l) = 0.0_wp
               do i = j, m
                  h(l) = h(l) + a(i,l)**2
               end do
               if (h(l) > h(lmax)) lmax = l
            end do
            hmax = h(lmax)
         end if
         ! LMAX HAS BEEN DETERMINED
         ! DO COLUMN INTERCHANGES IF NEEDED.
         ip(j) = lmax
         if (ip(j) /= j) then
            do i = 1, m
               tmp = a(i,j)
               a(i,j) = a(i,lmax)
               a(i,lmax) = tmp
            end do
            h(lmax) = h(j)
         end if
         ! COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A
         ! AND B.
         call dh12(1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j)
         call dh12(2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb)
      end do

      ! DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.
      do j = 1, ldiag
         if (abs(a(j,j)) <= tau) then
            k = j - 1
            exit
         else
            if (j==ldiag) k = ldiag
         end if
      end do

         kp1 = k + 1
         ! COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
         if (nb >= 1) then
            do jb = 1, nb
               tmp = szero
               if (m >= kp1) then
                  do i = kp1, m
                     tmp = tmp + b(i,jb)**2
                  end do
               end if
               rnorm(jb) = sqrt(tmp)
            end do
         end if
         ! SPECIAL FOR PSEUDORANK = 0
         if (k > 0) then
            ! IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER
            ! DECOMPOSITION OF FIRST K ROWS.

            if (k /= n) then
               do ii = 1, k
                  i = kp1 - ii
                  call dh12(1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1)
               end do
            end if

            if (nb >= 1) then
               do jb = 1, nb
                  ! SOLVE THE K BY K TRIANGULAR SYSTEM.
                  do l = 1, k
                     sm = dzero
                     i = kp1 - l
                     ip1 = i + 1
                     if (k >= ip1) then
                        do j = ip1, k
                           sm = sm + a(i,j)*b(j,jb)
                        end do
                     end if
                     sm1 = sm
                     b(i,jb) = (b(i,jb) - sm1)/a(i,i)
                  end do

                  ! COMPLETE COMPUTATION OF SOLUTION VECTOR.
                  if (k /= n) then
                     do j = kp1, n
                        b(j,jb) = szero
                     end do
                     do i = 1, k
                        call dh12(2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1)
                     end do
                  end if

                  ! RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
                  ! COLUMN INTERCHANGES.

                  do jj = 1, ldiag
                     j = ldiag + 1 - jj
                     if (ip(j) /= j) then
                        l = ip(j)
                        tmp = b(l,jb)
                        b(l,jb) = b(j,jb)
                        b(j,jb) = tmp
                     end if
                  end do

               end do

            end if

         elseif (nb >= 1) then
           do jb = 1, nb
              do i = 1, n
                 b(i,jb) = szero
             end do
           end do
         end if

   end if

   ! THE SOLUTION VECTORS, X, ARE NOW
   ! IN THE FIRST N ROWS OF THE ARRAY B(,).

   krank = k

   end subroutine dhfti

dhfti Subroutine

Contents

Source Code

private subroutine dhfti(a, mda, m, n, b, mdb, nb, tau, krank, rnorm, h, g, ip)

References

Revision history

Arguments

Calls

Called by

Source Code