LAPACK: TESTING/LIN/zgbt01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
00002      $                   RESID )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            KL, KU, LDA, LDAFAC, M, N
00010       DOUBLE PRECISION   RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZGBT01 reconstructs a band matrix  A  from its L*U factorization and
00021 *  computes the residual:
00022 *     norm(L*U - A) / ( N * norm(A) * EPS ),
00023 *  where EPS is the machine epsilon.
00024 *
00025 *  The expression L*U - A is computed one column at a time, so A and
00026 *  AFAC are not modified.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  M       (input) INTEGER
00032 *          The number of rows of the matrix A.  M >= 0.
00033 *
00034 *  N       (input) INTEGER
00035 *          The number of columns of the matrix A.  N >= 0.
00036 *
00037 *  KL      (input) INTEGER
00038 *          The number of subdiagonals within the band of A.  KL >= 0.
00039 *
00040 *  KU      (input) INTEGER
00041 *          The number of superdiagonals within the band of A.  KU >= 0.
00042 *
00043 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00044 *          The original matrix A in band storage, stored in rows 1 to
00045 *          KL+KU+1.
00046 *
00047 *  LDA     (input) INTEGER.
00048 *          The leading dimension of the array A.  LDA >= max(1,KL+KU+1).
00049 *
00050 *  AFAC    (input) COMPLEX*16 array, dimension (LDAFAC,N)
00051 *          The factored form of the matrix A.  AFAC contains the banded
00052 *          factors L and U from the L*U factorization, as computed by
00053 *          ZGBTRF.  U is stored as an upper triangular band matrix with
00054 *          KL+KU superdiagonals in rows 1 to KL+KU+1, and the
00055 *          multipliers used during the factorization are stored in rows
00056 *          KL+KU+2 to 2*KL+KU+1.  See ZGBTRF for further details.
00057 *
00058 *  LDAFAC  (input) INTEGER
00059 *          The leading dimension of the array AFAC.
00060 *          LDAFAC >= max(1,2*KL*KU+1).
00061 *
00062 *  IPIV    (input) INTEGER array, dimension (min(M,N))
00063 *          The pivot indices from ZGBTRF.
00064 *
00065 *  WORK    (workspace) COMPLEX*16 array, dimension (2*KL+KU+1)
00066 *
00067 *  RESID   (output) DOUBLE PRECISION
00068 *          norm(L*U - A) / ( N * norm(A) * EPS )
00069 *
00070 *  =====================================================================
00071 *
00072 *     .. Parameters ..
00073       DOUBLE PRECISION   ZERO, ONE
00074       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00075 *     ..
00076 *     .. Local Scalars ..
00077       INTEGER            I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
00078       DOUBLE PRECISION   ANORM, EPS
00079       COMPLEX*16         T
00080 *     ..
00081 *     .. External Functions ..
00082       DOUBLE PRECISION   DLAMCH, DZASUM
00083       EXTERNAL           DLAMCH, DZASUM
00084 *     ..
00085 *     .. External Subroutines ..
00086       EXTERNAL           ZAXPY, ZCOPY
00087 *     ..
00088 *     .. Intrinsic Functions ..
00089       INTRINSIC          DBLE, DCMPLX, MAX, MIN
00090 *     ..
00091 *     .. Executable Statements ..
00092 *
00093 *     Quick exit if M = 0 or N = 0.
00094 *
00095       RESID = ZERO
00096       IF( M.LE.0 .OR. N.LE.0 )
00097      $   RETURN
00098 *
00099 *     Determine EPS and the norm of A.
00100 *
00101       EPS = DLAMCH( 'Epsilon' )
00102       KD = KU + 1
00103       ANORM = ZERO
00104       DO 10 J = 1, N
00105          I1 = MAX( KD+1-J, 1 )
00106          I2 = MIN( KD+M-J, KL+KD )
00107          IF( I2.GE.I1 )
00108      $      ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
00109    10 CONTINUE
00110 *
00111 *     Compute one column at a time of L*U - A.
00112 *
00113       KD = KL + KU + 1
00114       DO 40 J = 1, N
00115 *
00116 *        Copy the J-th column of U to WORK.
00117 *
00118          JU = MIN( KL+KU, J-1 )
00119          JL = MIN( KL, M-J )
00120          LENJ = MIN( M, J ) - J + JU + 1
00121          IF( LENJ.GT.0 ) THEN
00122             CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
00123             DO 20 I = LENJ + 1, JU + JL + 1
00124                WORK( I ) = ZERO
00125    20       CONTINUE
00126 *
00127 *           Multiply by the unit lower triangular matrix L.  Note that L
00128 *           is stored as a product of transformations and permutations.
00129 *
00130             DO 30 I = MIN( M-1, J ), J - JU, -1
00131                IL = MIN( KL, M-I )
00132                IF( IL.GT.0 ) THEN
00133                   IW = I - J + JU + 1
00134                   T = WORK( IW )
00135                   CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
00136      $                        1 )
00137                   IP = IPIV( I )
00138                   IF( I.NE.IP ) THEN
00139                      IP = IP - J + JU + 1
00140                      WORK( IW ) = WORK( IP )
00141                      WORK( IP ) = T
00142                   END IF
00143                END IF
00144    30       CONTINUE
00145 *
00146 *           Subtract the corresponding column of A.
00147 *
00148             JUA = MIN( JU, KU )
00149             IF( JUA+JL+1.GT.0 )
00150      $         CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ),
00151      $                     1, WORK( JU+1-JUA ), 1 )
00152 *
00153 *           Compute the 1-norm of the column.
00154 *
00155             RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) )
00156          END IF
00157    40 CONTINUE
00158 *
00159 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00160 *
00161       IF( ANORM.LE.ZERO ) THEN
00162          IF( RESID.NE.ZERO )
00163      $      RESID = ONE / EPS
00164       ELSE
00165          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00166       END IF
00167 *
00168       RETURN
00169 *
00170 *     End of ZGBT01
00171 *
00172       END