計算ルーチン: 複素エルミート正定値帯行列の平衡化と条件数を小さくする目的で行と列のスケーリングを計算する

ZPBEQU Example Program Data
  4  1                                                       :Values of N and KD
 (  9.39, 0.00) (  1.08,-1.73)
                (  1.69, 0.00) ( -0.04E+10, 0.29E+10)
                               (  2.64E+20, 0.00    ) ( -0.33E+10, 2.24E+10)
                                                      (  2.17,     0.00    )
                                                             :End of matrix A

 ZPBEQU Example Program Results

 Matrix A
                          1                       2                       3
 1  (  9.39E+00,  0.00E+00) (  1.08E+00, -1.73E+00)
 2                          (  1.69E+00,  0.00E+00) ( -4.00E+08,  2.90E+09)
 3                                                  (  2.64E+20,  0.00E+00)
 4
 
                          4
 1
 2
 3  ( -3.30E+09,  2.24E+10)
 4  (  2.17E+00,  0.00E+00)

 SCOND = 8.0E-11, AMAX = 2.6E+20

 Diagonal scaling factors
     3.3E-01    7.7E-01    6.2E-11    6.8E-01

 Scaled matrix
                    1                 2                 3                 4
 1  ( 1.0000, 0.0000) ( 0.2711,-0.4343)
 2                    ( 1.0000, 0.0000) (-0.0189, 0.1373)
 3                                      ( 1.0000, 0.0000) (-0.1379, 0.9359)
 4                                                        ( 1.0000, 0.0000)

    Program zpbequ_example

!     ZPBEQU Example Program Text

!     Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com

!     .. Use Statements ..
      Use blas_interfaces, Only: zdscal
      Use lapack_example_aux, Only: nagf_blas_zddscl, &
        nagf_file_print_matrix_complex_band_comp
      Use lapack_interfaces, Only: zpbequ
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=dp), Parameter :: one = 1.0_dp
      Real (Kind=dp), Parameter :: thresh = 0.1_dp
      Integer, Parameter :: nin = 5, nout = 6
      Character (1), Parameter :: uplo = 'U'
!     .. Local Scalars ..
      Real (Kind=dp) :: amax, big, scond, small
      Integer :: i, i0, i1, ifail, ilen, info, j, kd, ldab, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: ab(:, :)
      Real (Kind=dp), Allocatable :: s(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: epsilon, max, min, radix, real, tiny
!     .. Executable Statements ..
      Write (nout, *) 'ZPBEQU Example Program Results'
      Write (nout, *)
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n, kd
      ldab = kd + 1
      Allocate (ab(ldab,n), s(n))

!     Read the upper or lower triangular part of the band matrix A
!     from data file

      If (uplo=='U') Then
        Do i = 1, n
          Read (nin, *)(ab(kd+1+i-j,j), j=i, min(n,i+kd))
        End Do
      Else If (uplo=='L') Then
        Do i = 1, n
          Read (nin, *)(ab(1+i-j,j), j=max(1,i-kd), i)
        End Do
      End If

!     Print the matrix A

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      If (uplo=='U') Then
        Call nagf_file_print_matrix_complex_band_comp(n, n, 0, kd, ab, ldab, &
          'Bracketed', '1P,E10.2', 'Matrix A', 'Integer', rlabs, 'Integer', &
          clabs, 80, 0, ifail)
      Else If (uplo=='L') Then
        Call nagf_file_print_matrix_complex_band_comp(n, n, kd, 0, ab, ldab, &
          'Bracketed', '1P,E10.2', 'Matrix A', 'Integer', rlabs, 'Integer', &
          clabs, 80, 0, ifail)
      End If

      Write (nout, *)

!     Compute diagonal scaling factors

      Call zpbequ(uplo, n, kd, ab, ldab, s, scond, amax, info)

      If (info>0) Then
        Write (nout, 100) 'Diagonal element', info, ' of A is non positive'
      Else

!       Print SCOND, AMAX and the scale factors

        Write (nout, 110) 'SCOND =', scond, ', AMAX =', amax
        Write (nout, *)
        Write (nout, *) 'Diagonal scaling factors'
        Write (nout, 120) s(1:n)
        Write (nout, *)
        Flush (nout)

!       Compute values close to underflow and overflow

        small = tiny(1.0E0_dp)/(epsilon(1.0E0_dp)*real(radix(1.0E0_dp),kind=dp &
          ))
        big = one/small
        If ((scond<thresh) .Or. (amax<small) .Or. (amax>big)) Then

!         Scale A

          If (uplo=='U') Then
            Do j = 1, n
              i0 = max(1, j-kd)
              i1 = 1 + i0 - (j-kd)
              ilen = j - i0 + 1
              Call zdscal(ilen, s(j), ab(i1,j), 1)
              Call nagf_blas_zddscl(ilen, s(i0), 1, ab(i1,j), 1)
            End Do

          Else If (uplo=='L') Then
            Do j = 1, n
              i1 = 1
              ilen = min(n, j+kd) - j + 1
              Call zdscal(ilen, s(j), ab(i1,j), 1)
              Call nagf_blas_zddscl(ilen, s(j), 1, ab(i1,j), 1)
            End Do
          End If

!         Print the scaled matrix

          ifail = 0
          If (uplo=='U') Then
            Call nagf_file_print_matrix_complex_band_comp(n, n, 0, kd, ab, &
              ldab, 'Bracketed', 'F7.4', 'Scaled matrix', 'Integer', rlabs, &
              'Integer', clabs, 80, 0, ifail)
          Else If (uplo=='L') Then
            Call nagf_file_print_matrix_complex_band_comp(n, n, kd, 0, ab, &
              ldab, 'Bracketed', 'F7.4', 'Scaled matrix', 'Integer', rlabs, &
              'Integer', clabs, 80, 0, ifail)
          End If
        End If
      End If

100   Format (1X, A, I4, A)
110   Format (1X, 2(A,1P,E8.1))
120   Format ((1X,1P,7E11.1))
    End Program