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

ZGBEQU Example Program Data
  4  1  2                                                :Values of N, KL and KU
 (-1.65, 2.26) (-2.05D-10,-8.50D-11) ( 9.70D-01,-2.84D+00)
 ( 0.00, 6.30) (-1.48D-10,-1.75D-10) (-2.99D+00, 3.01D+00) ( 0.59D+00,-0.48D+00)
               (-7.70D-01, 2.83D+00) (-1.06D+10, 1.94D+10) ( 3.33D+10,-1.04D+10)
                                     ( 4.48D+00,-1.09D+00) (-0.46D+00,-1.72D+00)
                                                                :End of matrix A

 ZGBEQU Example Program Results

 Matrix A
               1            2            3            4
 1   -1.6500E+00  -2.0500E-10   9.7000E-01
      2.2600E+00  -8.5000E-11  -2.8400E+00
 
 2    0.0000E+00  -1.4800E-10  -2.9900E+00   5.9000E-01
      6.3000E+00  -1.7500E-10   3.0100E+00  -4.8000E-01
 
 3                -7.7000E-01  -1.0600E+10   3.3300E+10
                   2.8300E+00   1.9400E+10  -1.0400E+10
 
 4                              4.4800E+00  -4.6000E-01
                               -1.0900E+00  -1.7200E+00

 ROWCND = 8.9E-11, COLCND = 8.2E-11, AMAX = 4.4E+10

 Row scale factors
    2.56E-01   1.59E-01   2.29E-11   1.80E-01

 Column scale factors
    1.00E+00   1.21E+10   1.00E+00   1.00E+00

 Scaled matrix
          1       2       3       4
 1  -0.4220 -0.6364  0.2481
     0.5780 -0.2639 -0.7263
 
 2   0.0000 -0.2852 -0.4746  0.0937
     1.0000 -0.3372  0.4778 -0.0762
 
 3          -0.2139 -0.2426  0.7620
             0.7861  0.4439 -0.2380
 
 4                   0.8043 -0.0826
                    -0.1957 -0.3088

    Program zgbequ_example

!     ZGBEQU 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
      Use lapack_interfaces, Only: zgbequ
      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
!     .. Local Scalars ..
      Real (Kind=dp) :: amax, big, colcnd, rowcnd, small
      Integer :: i, i0, i1, ifail, ilen, info, j, k, kl, ku, ldab, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: ab(:, :)
      Real (Kind=dp), Allocatable :: c(:), r(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: epsilon, max, min, radix, real, tiny
!     .. Executable Statements ..
      Write (nout, *) 'ZGBEQU Example Program Results'
      Write (nout, *)
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n, kl, ku
      ldab = kl + ku + 1
      Allocate (ab(ldab,n), c(n), r(n))

!     Read the band matrix A from data file

      k = ku + 1
      Read (nin, *)((ab(k+i-j,j),j=max(i-kl,1),min(i+ku,n)), i=1, n)

!     Print the matrix A

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call nagf_file_print_matrix_complex_band(n, n, kl, ku, ab, ldab, &
        'Matrix A', ifail)

      Write (nout, *)

!     Compute row and column scaling factors

      Call zgbequ(n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)

      If (info>0) Then
        If (info<=n) Then
          Write (nout, 100) 'Row ', info, ' of A is exactly zero'
        Else
          Write (nout, 100) 'Column ', info - n, ' of A is exactly zero'
        End If
      Else

!       Print ROWCND, COLCND, AMAX and the scale factors

        Write (nout, 110) 'ROWCND =', rowcnd, ', COLCND =', colcnd, &
          ', AMAX =', amax
        Write (nout, *)
        Write (nout, *) 'Row scale factors'
        Write (nout, 120) r(1:n)
        Write (nout, *)
        Write (nout, *) 'Column scale factors'
        Write (nout, 120) c(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 ((rowcnd>=thresh) .And. (amax>=small) .And. (amax<=big)) Then
          If (colcnd<thresh) Then

!           Just column scale A
            Do j = 1, n
              i1 = 1 + max(1, j-ku) - (j-ku)
              ilen = min(n, j+kl) - max(1, j-ku) + 1
              Call zdscal(ilen, c(j), ab(i1,j), 1)
            End Do

          End If
        Else If (colcnd>=thresh) Then

!         Just row scale A
          Do j = 1, n
            i0 = max(1, j-ku)
            i1 = 1 + i0 - (j-ku)
            ilen = min(n, j+kl) - i0 + 1
            Call nagf_blas_zddscl(ilen, r(i0), 1, ab(i1,j), 1)
          End Do

        Else

!         Row and column scale A
          Do j = 1, n
            i0 = max(1, j-ku)
            i1 = 1 + i0 - (j-ku)
            ilen = min(n, j+kl) - i0 + 1
            Call zdscal(ilen, c(j), ab(i1,j), 1)
            Call nagf_blas_zddscl(ilen, r(i0), 1, ab(i1,j), 1)
          End Do

        End If

!       Print the scaled matrix
        ifail = 0
        Call nagf_file_print_matrix_complex_band(n, n, kl, ku, ab, ldab, &
          'Scaled matrix', ifail)

      End If

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