LAPACK: TESTING/LIN/dgbt01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGBT01( 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       DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DGBT01 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 *          DGBTRF.  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 DGBTRF 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 DGBTRF.
00064 *
00065 *  WORK    (workspace) DOUBLE PRECISION 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, T
00079 *     ..
00080 *     .. External Functions ..
00081       DOUBLE PRECISION   DASUM, DLAMCH
00082       EXTERNAL           DASUM, DLAMCH
00083 *     ..
00084 *     .. External Subroutines ..
00085       EXTERNAL           DAXPY, DCOPY
00086 *     ..
00087 *     .. Intrinsic Functions ..
00088       INTRINSIC          DBLE, MAX, MIN
00089 *     ..
00090 *     .. Executable Statements ..
00091 *
00092 *     Quick exit if M = 0 or N = 0.
00093 *
00094       RESID = ZERO
00095       IF( M.LE.0 .OR. N.LE.0 )
00096      $   RETURN
00097 *
00098 *     Determine EPS and the norm of A.
00099 *
00100       EPS = DLAMCH( 'Epsilon' )
00101       KD = KU + 1
00102       ANORM = ZERO
00103       DO 10 J = 1, N
00104          I1 = MAX( KD+1-J, 1 )
00105          I2 = MIN( KD+M-J, KL+KD )
00106          IF( I2.GE.I1 )
00107      $      ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) )
00108    10 CONTINUE
00109 *
00110 *     Compute one column at a time of L*U - A.
00111 *
00112       KD = KL + KU + 1
00113       DO 40 J = 1, N
00114 *
00115 *        Copy the J-th column of U to WORK.
00116 *
00117          JU = MIN( KL+KU, J-1 )
00118          JL = MIN( KL, M-J )
00119          LENJ = MIN( M, J ) - J + JU + 1
00120          IF( LENJ.GT.0 ) THEN
00121             CALL DCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
00122             DO 20 I = LENJ + 1, JU + JL + 1
00123                WORK( I ) = ZERO
00124    20       CONTINUE
00125 *
00126 *           Multiply by the unit lower triangular matrix L.  Note that L
00127 *           is stored as a product of transformations and permutations.
00128 *
00129             DO 30 I = MIN( M-1, J ), J - JU, -1
00130                IL = MIN( KL, M-I )
00131                IF( IL.GT.0 ) THEN
00132                   IW = I - J + JU + 1
00133                   T = WORK( IW )
00134                   CALL DAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
00135      $                        1 )
00136                   IP = IPIV( I )
00137                   IF( I.NE.IP ) THEN
00138                      IP = IP - J + JU + 1
00139                      WORK( IW ) = WORK( IP )
00140                      WORK( IP ) = T
00141                   END IF
00142                END IF
00143    30       CONTINUE
00144 *
00145 *           Subtract the corresponding column of A.
00146 *
00147             JUA = MIN( JU, KU )
00148             IF( JUA+JL+1.GT.0 )
00149      $         CALL DAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1,
00150      $                     WORK( JU+1-JUA ), 1 )
00151 *
00152 *           Compute the 1-norm of the column.
00153 *
00154             RESID = MAX( RESID, DASUM( JU+JL+1, WORK, 1 ) )
00155          END IF
00156    40 CONTINUE
00157 *
00158 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00159 *
00160       IF( ANORM.LE.ZERO ) THEN
00161          IF( RESID.NE.ZERO )
00162      $      RESID = ONE / EPS
00163       ELSE
00164          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00165       END IF
00166 *
00167       RETURN
00168 *
00169 *     End of DGBT01
00170 *
00171       END