subroutine dormgqr_test ( )

#*****************************************************************************80
#
## dormgqr_test() tests dormgqr().
#
#  Discussion:
#
#    We want to solve the MxN linear system A*x=b using the QR approach:
#
#    Factor A:
#
#      A = Q * R                        (step 1)
#
#    Transform:
#
#               A * x =               b
#      ==>  Q * R * x =               b
#      ==>      R * x =          Q' * b  (step 2)
#      ==>          x = inv(R) * Q' * b. (step 3)
#
#    Step 1) DGEQRF computes the QR factorization of an M by N matrix A:
#    A(MxN) = Q(MxK) * R(KxN) where K = min ( M, N ).
#
#    Step 2) DORMQR can multiply Q' * b, putting the result back into b.
#
#    Step 3) We could call a LAPACK routine to solve the upper triangular
#    system R * x = Q' * b.  Instead, we will try this part ourselves.
#
#
#    LAPACK makes this process tricky because of two things it does
#    for efficiency:
#
#    *) LAPACK computes the Q and R factors in a
#       compressed and encoded form, overwriting the matrix A and
#       storing some extra information in a vector called TAU.
#
#    *) LAPACK defines K = min ( M, N ), and
#       does NOT compute the QR factorization as an MxM Q
#       times an MxN R.  Instead, it computes an MxK Q times
#       a KxN R.  This saves it worrying about zeroes, but it
#       means the programmer has to worry about proper storage
#       and correct dimensioning.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 January 2011
#
#  Author:
#
#    John Burkardt
#
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: m = 8
  integer, parameter :: n = 6

  integer, parameter :: k = min ( m, n )
  integer, parameter :: lwork = n

  real ( kind = rk8 ) a(m,n)
  real ( kind = rk8 ) b(m)
  integer i
  integer info
  integer j
  integer lda
  real ( kind = rk8 ) r(k,n)
  real ( kind = rk8 ) tau(k)
  real ( kind = rk8 ) work(lwork)
  real ( kind = rk8 ) x(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'dormgqr_test():'
  write ( *, '(a)' ) '  dormqr() computes Q'' * b.'
  write ( *, '(a)' ) '  after dgeqrf() computes the QR factorization:'
  write ( *, '(a)' ) '    A = Q * R'
  write ( *, '(a)' ) '  storing a double precision real matrix (D)'
  write ( *, '(a)' ) '  in general storage mode (GE).'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We use these routines to carry out a QR'
  write ( *, '(a)' ) '  solve of an M by N linear system A * x = b.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this case, our M x N matrix A has more rows'
  write ( *, '(a)' ) '  than columns:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  M = ', m
  write ( *, '(a,i8)' ) '  N = ', n
#
#  STEP 0: Set random A, simple x, and compute b.
#
  call random_number ( harvest = a(1:m,1:n) )

  do i = 1, n
    x(i) = real ( i, kind = rk8 )
  end do

  b(1:m) = matmul ( a(1:m,1:n), x(1:n) )
#
#  Wipe out X so we believe it when we get it back...
#
  x(1:n) =  0.0D+00

  call r8mat_print ( m, n, a, '  The matrix A:' )
#
#  STEP 1: Compute the QR factorization.
#
  lda = m

  call dgeqrf ( m, n, a, lda, tau, work, lwork, info )

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  DGEQRF returned nonzero INFO = ', info
    return
  end if
#
#  Extract and save the K by N matrix R.
#
  r(1:k,1:n) = 0.0D+00
  do i = 1, k
    r(i,i:n) = a(i,i:n)
  end do
#
#  STEP 2: Multiply Q' * b to get the N by 1 result in b.
#
  call dormqr ( 'L', 'T', m, 1, k, a, m, tau, b, m, work, lwork, info )

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  DORMQR returned nonzero INFO = ', info
    return
  end if
#
#  STEP 3: Compute inv(R) * Q' * b, or, equivalently,
#  solve R * x = Q' * b.
#
  do j = n, 1, -1
    x(j) = b(j) / r(j,j)
    do i = 1, n - 1
      b(i) = b(i) - r(i,j) * x(j)
    end do
  end do

  call r8vec_print ( n, x, '  The solution X:' )

  return
end


subroutine dpbtrf_test ( )

#*****************************************************************************80
#
## dpbtrf_test() tests dpbtrf().
#
#  Discussion:
#
#    We want to compute the lower triangular Cholesky factor L
#    of a symmetric positive definite (SPD) band matrix.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    12 September 2006
#
#  Author:
#
#    John Burkardt
#
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: n = 5
  integer, parameter :: nband = 1

  real ( kind = rk8 ) a(nband+1,n)
  integer i
  integer info
  integer j
  real ( kind = rk8 ) l(nband+1,n)
  real ( kind = rk8 ) l_row(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DPBTRF_TEST'
  write ( *, '(a)' ) '  DPBTRF computes'
  write ( *, '(a)' ) '    the lower Cholesky factor A = L*L'' or'
  write ( *, '(a)' ) '    the upper Cholesky factor A = U''*U;'
  write ( *, '(a)' ) '  For a double precision real matrix (D)'
  write ( *, '(a)' ) '  in positive definite band storage mode (PB):'
#
#  Zero out the matrix.
#
  a(1:nband+1,1:n) = 0.0D+00
#
#  Store the diagonal of a symmetric band matrix.
#
  a(1,1:n) = 2.0D+00
#
#  Store the subdiagonal of a symmetric band matrix.
#
  a(2,1:n-1) = -1.0D+00
#
#  Get the lower triangular Cholesky factor L:
#
  l(1:nband+1,1:n) = a(1:nband+1,1:n)

  call dpbtrf ( 'L', n, nband, l, nband+1, info )

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Factorization failed, INFO = ', info
    return
  end if
#
#  Print the relevant entries of L:
#
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The lower Cholesky factor L:'
  write ( *, '(a)' ) ' '

  do i = 1, n
    do j = 1, n

      if ( 0 <= i - j .and. i-j <= nband ) then
        l_row(j) = l(i-j+1,j)
      else
        l_row(j) = 0.0D+00
      end if

    end do

    write ( *, '(5f10.6)' ) l_row(1:n)

  end do

  return
end
subroutine dpbtrs_test ( )

#*****************************************************************************80
#
## dpbtrs_test() tests dpbtrs().
#
#  Discussion:
#
#    The problem is just an enlarged version of the
#    problem for n = 5, which is:
#
#    Matrix A is ( 2 -1  0  0  0)    right hand side b is  (1)
#                (-1  2 -1  0  0)                          (0)
#                ( 0 -1  2 -1  0)                          (0)
#                ( 0  0 -1  2 -1)                          (0)
#                ( 0  0  0 -1  2)                          (1)
#
#
#    Solution is   (1)
#                  (1)
#                  (1)
#                  (1)
#                  (1)
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    12 September 2006
#
#  Author:
#
#    John Burkardt
#
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: n = 25
  integer, parameter :: nband = 1

  integer, parameter :: lda = nband + 1

  real ( kind = rk8 ) a(lda,n)
  real ( kind = rk8 ) b(n)
  integer info
  integer nrhs

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DPBTRS_TEST'
  write ( *, '(a)' ) '  DPBTRS solves linear systems'
  write ( *, '(a)' ) '  for a positive definite symmetric band matrix,'
  write ( *, '(a)' ) '  stored as a double precision real matrix (D)'
  write ( *, '(a)' ) '  in positive definite band storage mode (PB):'
#
#  Zero out the matrix.
#
  a(1:lda,1:n) = 0.0D+00
#
#  Super (and sub) diagonal.
#
  a(1,2:n) = -1.0D+00
#
#  Diagonal.
#
  a(2,1:n) = 2.0D+00
#
#  Set the right hand side.
#
  b(1) = 1.0D+00
  b(2:n-1) = 0.0D+00
  b(n) = 1.0D+00
#
#  Factor the matrix.
#
  call dpbtrf ( 'u', n, nband, a, lda, info )

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Factorization failed, INFO = ', info
    return
  end if
#
#  Solve the linear system.
#
  nrhs = 1
  call dpbtrs ( 'u', n, nband, nrhs, a, lda, b, n, info )

  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Solution failed, INFO = ', info
  end if

  call r8vec_print_some ( n, b, 1, 5, '  Partial solution (all should be 1)' )

  return
end


subroutine dsbgvx_test ( )

#*****************************************************************************80
#
## dsbgvx_test() tests dsbgvx().
#
#  Discussion:
#
#    DSBGVX deals with the generalized eigenvalue problem:
#
#      A * x = lambda * B * x
#
#    where A and B are symmetric and banded (and stored in LAPACK symmetric
#    band storage mode).  B is additionally assumed to be positive definite.
#
#    This is an "expert" interface, and the user is requesting
#    only some of the eigenvalues and eigenvectors.  In this example,
#    only the largest and smallest (in magnitude) eigenvalues will
#    be requested.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    17 March 2003
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    real ( kind = rk8 ) AB(LDAB,N), contains, on input, the upper or lower
#    triangle of the symmetric band matrix A, stored in the first KA+1 rows
#    of the array AB.
#    If UPLO = 'U', then
#      AB(KA+1+I-J,J) = A(I,J) for max(1,J-KA) <= I <= J;
#    If UPLO = 'L', then
#      AB(1+I-J,J) = A(I,J) for J <= I <= min(N,J+KA).
#
#    real ( kind = rk8 ) ABSTOL, the absolute error tolerance for the
#    eigenvalues.  If the input value of ABSTOL is not positive, then an
#    appropriate value will be determined internally and used instead.
#
#    real ( kind = rk8 ) BB(LDBB,N), contains, on input, the upper or lower
#    triangle of the positive definite symmetric band matrix B, stored in
#    the first KB+1 rows of the array BB.
#    If UPLO = 'U', then
#      BB(KB+1+I-J,J) = B(I,J) for max(1,J-KB) <= I <= J;
#    If UPLO = 'L', then
#      BB(1+I-J,J) = B(I,J) for J <= I <= min(N,J+KB).
#
#    integer IFAIL(N), if JOBZ = 'V', then if INFO is 0, the first
#    M elements of IFAIL have been set to zero by DSBGVX, but if INFO
#    is nonzero, IFAIL contains the indices of the eigenvalues
#    for which the eigenvectors failed to converge.  If JOBZ = 'N',
#    then IFAIL is not referenced.
#
#    integer IL, IU, the indices of the first (smallest) and last
#    (largest) eigenvalues to be returned.  These values are only
#    used if RANGE = 'I'.  It must be the case that 1 <= IL <= IU <= N.
#
#    Integer INFO, is 0 for a successful computation,
#    negative if an input argument was illegal (the index of that
#    argument is the value of -INFO), or positive, in which case,
#    if 0 < INFO <= N, then INFO off diagonal elements of an
#    intermediate tridiagonal form did not converge to zero, or
#    if N < INFO, B is not positive definite and the computation
#    could not be completed.
#
#    integer IWORK(5*N), workspace.
#
#    character JOBZ, is 'N' if only eigenvalues are desired, or 'V'
#    if eigenvectors will also be required.
#
#    Integer KA, the number of superdiagonals (if UPLO = 'U') or
#    subdiagonals (if UPLO = 'L') of A that are nonzero.
#
#    integer KB, the number of superdiagonals (if UPLO = 'U') or
#    subdiagonals (if UPLO = 'L') of B that are nonzero.
#
#    integer LDAB, the leading dimension of the array AB, which
#    must be at least KA+1.
#
#    integer LDBB, the leading dimension of the array BB, which
#    must be at least KB+1.
#
#    integer LDQ, the leading dimension of the array Q.
#    If JOBZ = 'N', then Q is not used, and LDQ should be any
#    positive value such as 1.  If JOBZ = 'V', then LDQ must be
#    at least N.
#
#    integer LDZ, the leading dimension of the array Z.
#    If JOBZ = 'N', then Z is not used, and LDZ should be any
#    positive value such as 1.  If JOBZ = 'V', then LDZ must be
#    at least N.
#
#    integer M, the number of eigenvalues found by DSBGVX.
#
#    integer N, the order of the matrices A and B.
#
#    real ( kind = rk8 ) Q(LDQ,N), if JOBZ = 'V', the N by N matrix used to
#    reduce the problem to standard form: "C * x = lambda * x"
#    and then to reduce the matrix C to tridiagonal form.  But
#    if JOBZ is not 'V', Q is not referenced.
#
#    character RANGE, specifies which eigenvalues are desired.
#    'A' means all, 'V' means a real interval will be specified in which
#    eigenvalues are to be sought, 'I' means a range of indices will
#    be specified.
#
#    character UPLO, is 'U' if the upper triangles of A and B are stored,
#    'L' if the lower triangles are stored.
#
#    real ( kind = rk8 ) VL, VU, the lower and upper bounds of an interval to be
#    searched for eigenvalues.  In this case, VL must be less than VU.
#    These values are used only if RANGE = 'V'.
#
#    real ( kind = rk8 ) W(N), the requested eigenvalues, in ascending order.
#
#    real ( kind = rk8 ) WORK(7*N), workspace.
#
#    real ( kind = rk8 ) Z(LDZ,N), if JOBZ = 'V', the I-th column of Z contains
#    the eigenvector associated with the I-th eigenvalue W(I).
#
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: n = 11
  integer, parameter :: ka = 2
  integer, parameter :: kb = 1

  integer, parameter :: ldab = ka+1
  integer, parameter :: ldbb = kb+1
  integer, parameter :: ldq = 1
  integer, parameter :: ldz = 1

  real ( kind = rk8 ) ab(ldab,n)
  real ( kind = rk8 ), parameter :: abstol = 0.0D+00
  real ( kind = rk8 ) bb(ldbb,n)
  integer i
  integer ifail(n)
  integer il
  integer ilo
  integer info
  integer iu
  integer iwork(5*n)
  integer j
  character :: jobz = 'N'
  integer m
  real ( kind = rk8 ) q(ldq,n)
  character :: range = 'I'
  integer test
  character :: uplo = 'U'
  real ( kind = rk8 ) value
  real ( kind = rk8 ), parameter :: vl = 0.0D+00
  real ( kind = rk8 ), parameter :: vu = 1.0D+00
  real ( kind = rk8 ) w(n)
  real ( kind = rk8 ) work(7*n)
  real ( kind = rk8 ) z(ldz,n)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'dsbgvx_test():'
  write ( *, '(a)' ) '  dsbgvx() solves the generalized eigenvalue problem'
  write ( *, '(a)' ) '    A * X = LAMBDA * B * X'
  write ( *, '(a)' ) '  for a symmetric banded NxN matrix A, and a symmetric'
  write ( *, '(a)' ) '  banded positive definite NxN matrix B,'
  write ( *, '(a)' ) ' '

  do test = 1, 2
#
#  Set A.
#
    do j = 1, n
      ilo = max ( j - ka, 1 )
      do i = ilo, j

        if ( j == i-2 ) then
          value = -1.0D+00
        else if ( j == i-1 ) then
          value = -1.0D+00
        else if ( j == i ) then
          value = +4.0D+00
        else if ( j == i+1 ) then
          value = -1.0D+00
        else if ( j == i+2 ) then
          value = -1.0D+00
        else
          value = 0.0D+00
        end if

        ab(ka+1+i-j,j) = value

      end do
    end do
#
#  Set B.
#
    do j = 1, n
      ilo = max ( j - kb, 1 )
      do i = ilo, j

        if ( j == i-1 ) then
          value = -1.0D+00
        else if ( j == i ) then
          value = +2.0D+00
        else if ( j == i+1 ) then
          value = -1.0D+00
        else
          value = 0.0D+00
        end if

        bb(kb+1+i-j,j) = value

      end do
    end do
#
#  Request the value of the SMALLEST or LARGEST eigenvalue:
#
    if ( test == 1 ) then
      il = 1
      iu = 1
    else if ( test == 2 ) then
      il = n
      iu = n
    end if

    call dsbgvx ( jobz, range, uplo, n, ka, kb, ab, ldab, bb, ldbb, q, &
      ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info )

    if ( info < 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Illegal value for input argument ', -info
      return
    else if ( 0 < info ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The eigenvalue or eigenvector iterations'
      write ( *, '(a)' ) '  did not converge.'
      cycle
    end if

    call r8vec_print ( m, w, '  Computed eigenvalues' )

  end do

  return
end