dbndac – bspline-fortran

private subroutine dbndac(g, Mdg, Nb, Ip, Ir, Mt, Jt)

These subroutines solve the least squares problem Ax = b for banded matrices A using sequential accumulation of rows of the data matrix. Exactly one right-hand side vector is permitted.

These subroutines are intended for the type of least squares systems that arise in applications such as curve or surface fitting of data. The least squares equations are accumulated and processed using only part of the data. This requires a certain user interaction during the solution of Ax = b.

Specifically, suppose the data matrix (A B) is row partitioned into Q submatrices. Let (E F) be the T-th one of these submatrices where E = (0 C 0). Here the dimension of E is MT by N and the dimension of C is MT by NB. The value of NB is the bandwidth of A. The dimensions of the leading block of zeros in E are MT by JT-1.

The user of the subroutine DBNDAC provides MT,JT,C and F for T=1,...,Q. Not all of this data must be supplied at once.

Following the processing of the various blocks (E F), the matrix (A B) has been transformed to the form (R D) where R is upper triangular and banded with bandwidth NB. The least squares system Rx = d is then easily solved using back substitution by executing the statement CALL DBNDSL(1,...). The sequence of values for JT must be nondecreasing. This may require some preliminary interchanges of rows and columns of the matrix A.

The primary reason for these subroutines is that the total processing can take place in a working array of dimension MU by NB+1. An acceptable value for MU is

                MU = MAX(MT + N + 1),

where N is the number of unknowns.

Here the maximum is taken over all values of MT for T=1,...,Q. Notice that MT can be taken to be a small as one, showing that MU can be as small as N+2. The subprogram DBNDAC processes the rows more efficiently if MU is large enough so that each new block (C F) has a distinct value of JT.

The four principle parts of these algorithms are obtained by the following call statements:

CALL [[DBNDAC]](...) Introduce new blocks of data.
CALL [[DBNDSL]](1,...) Compute solution vector and length of residual vector.
CALL [[DBNDSL]](2,...) Given any row vector H solve YR = H for the row vector Y.
CALL [[DBNDSL]](3,...) Given any column vector W solve RZ = W for the column vector Z.

Remarks

To obtain the upper triangular matrix and transformed right-hand side vector D so that the super diagonals of R form the columns of G(,), execute the following Fortran statements.

     nbp1=nb+1
     do j=1, nbp1
       g(ir,j) = 0.0
     end do
     mt=1
     jt=n+1
     call dbndac(g,mdg,nb,ip,ir,mt,jt)

References

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

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)
920501 Reformatted the REFERENCES section. (WRB)

Arguments

Type	Intent		Name
real(kind=wp),	intent(inout)	::	g(Mdg,*)	`G(MDG,NB+1)` Input The working array into which the user will place the `MT` by `NB+1` block `(C F)` in rows `IR` through `IR+MT-1`, columns 1 through `NB+1`. See descriptions of `IR` and `MT` below. Output The working array which will contain the processed rows of that part of the data matrix which has been passed to DBNDAC.
integer,	intent(in)	::	Mdg	The number of rows in the working array `G(,)`. The value of MDG should be `>= MU`. The value of `MU` is defined in the abstract of these subprograms.
integer,	intent(in)	::	Nb	The bandwidth of the data matrix `A`.
integer,	intent(inout)	::	Ip	Input Set by the user to the value 1 before the first call to DBNDAC. Its subsequent value is controlled by DBNDAC to set up for the next call to DBNDAC. Output The value of this argument is advanced by DBNDAC to be ready for storing and processing a new block of data in `G(,)`.
integer,	intent(inout)	::	Ir	Input Index of the row of `G(,)` where the user is to place the new block of data `(C F)`. Set by the user to the value 1 before the first call to DBNDAC. Its subsequent value is controlled by DBNDAC. A value of `IR > MDG` is considered an error. Output The value of this argument is advanced by DBNDAC to be ready for storing and processing a new block of data in `G(,)`.
integer,	intent(in)	::	Mt	Set by the user to indicate the number of new rows of data in the block
integer,	intent(in)	::	Jt	Set by the user to indicate the index of the first nonzero column in that set of rows `(E F) = (0 C 0 F)` being processed.

Calls

Help

Called by

Help

Source Code

   subroutine dbndac(g, Mdg, Nb, Ip, Ir, Mt, Jt)

      implicit none

      integer,intent(in) :: Mdg !! The number of rows in the working array
                                !! `G(*,*)`.  The value of MDG should be `>= MU`.
                                !! The value of `MU` is defined in the abstract
                                !! of these subprograms.
      real(wp),intent(inout) :: g(Mdg, *) !! `G(MDG,NB+1)`
                                          !!
                                          !! *Input*
                                          !! The working array into which the user will
                                          !! place the `MT` by `NB+1` block `(C F)` in rows `IR`
                                          !! through `IR+MT-1`, columns 1 through `NB+1`.
                                          !! See descriptions of `IR` and `MT` below.
                                          !!
                                          !! *Output*
                                          !! The working array which will contain the
                                          !! processed rows of that part of the data
                                          !! matrix which has been passed to [[DBNDAC]].
      integer,intent(in) :: Nb !! The bandwidth of the data matrix `A`.
      integer,intent(inout) :: Ip !! *Input*
                                  !! Set by the user to the value 1 before the
                                  !! first call to [[DBNDAC]].  Its subsequent value
                                  !! is controlled by [[DBNDAC]] to set up for the
                                  !! next call to [[DBNDAC]].
                                  !!
                                  !! *Output*
                                  !! The value of this argument is advanced by
                                  !! [[DBNDAC]] to be ready for storing and processing
                                  !! a new block of data in `G(*,*)`.
      integer,intent(inout) :: Ir !! *Input*
                                  !! Index of the row of `G(*,*)` where the user is
                                  !! to place the new block of data `(C F)`.  Set by
                                  !! the user to the value 1 before the first call
                                  !! to [[DBNDAC]].  Its subsequent value is controlled
                                  !! by [[DBNDAC]]. A value of `IR > MDG` is considered
                                  !! an error.
                                  !!
                                  !! *Output*
                                  !! The value of this argument is advanced by
                                  !! [[DBNDAC]] to be ready for storing and processing
                                  !! a new block of data in `G(*,*)`.
      integer,intent(in) :: Mt !! Set by the user to indicate the
                               !! number of new rows of data in the block
      integer,intent(in) :: Jt !! Set by the user to indicate
                               !! the index of the first nonzero column in that
                               !! set of rows `(E F) = (0 C 0 F)` being processed.

      real(wp) :: rho
      integer :: i, ie, ig, ig1, ig2, iopt, j, jg, &
                 k, kh, l, lp1, mh, mu, nbp1, nerr

      real(wp), parameter :: zero = 0.0_wp

      ! ALG. STEPS 1-4 ARE PERFORMED EXTERNAL TO THIS SUBROUTINE.

      nbp1 = Nb + 1
      if (Mt <= 0 .or. Nb <= 0) return
      if (.not. Mdg < Ir) then
         if (Jt /= Ip) then
            if (Jt > Ir) then
               do i = 1, Mt
                  ig1 = Jt + Mt - i
                  ig2 = Ir + Mt - i
                  do j = 1, nbp1
                     g(ig1, j) = g(ig2, j)
                  end do
               end do
               ie = Jt - Ir
               do i = 1, ie
                  ig = Ir + i - 1
                  do j = 1, nbp1
                     g(ig, j) = zero
                  end do
               end do
               Ir = Jt
            end if
            mu = min(Nb - 1, Ir - Ip - 1)
            if (mu /= 0) then
               do l = 1, mu
                  k = min(l, Jt - Ip)
                  lp1 = l + 1
                  ig = Ip + l
                  do i = lp1, Nb
                     jg = i - k
                     g(ig, jg) = g(ig, i)
                  end do
                  do i = 1, k
                     jg = nbp1 - i
                     g(ig, jg) = zero
                  end do
               end do
            end if
            Ip = Jt
         end if
         mh = Ir + Mt - Ip
         kh = min(nbp1, mh)
         do i = 1, kh
            call dh12(1, i, max(i + 1, Ir - Ip + 1), mh, g(Ip, i), 1, &
                      rho, g(Ip, i + 1), 1, Mdg, nbp1 - i)
         end do
         Ir = Ip + kh
         if (kh >= nbp1) then
            do i = 1, Nb
               g(Ir - 1, i) = zero
            end do
         end if
      else
         nerr = 1
         iopt = 2
         write (*, *) 'MDG<IR, Probable error.'
      end if

   end subroutine dbndac

dbndac Subroutine

Contents

Source Code

private subroutine dbndac(g, Mdg, Nb, Ip, Ir, Mt, Jt)

Remarks

References

Revision history

Arguments

Calls

Called by

Source Code