LAPACK 3.3.0

cpbt01.f

Go to the documentation of this file.
00001       SUBROUTINE CPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
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       CHARACTER          UPLO
00010       INTEGER            KD, LDA, LDAFAC, N
00011       REAL               RESID
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               RWORK( * )
00015       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CPBT01 reconstructs a Hermitian positive definite band matrix A from
00022 *  its L*L' or U'*U factorization and computes the residual
00023 *     norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00024 *     norm( U'*U - A ) / ( N * norm(A) * EPS ),
00025 *  where EPS is the machine epsilon, L' is the conjugate transpose of
00026 *  L, and U' is the conjugate transpose of U.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  UPLO    (input) CHARACTER*1
00032 *          Specifies whether the upper or lower triangular part of the
00033 *          Hermitian matrix A is stored:
00034 *          = 'U':  Upper triangular
00035 *          = 'L':  Lower triangular
00036 *
00037 *  N       (input) INTEGER
00038 *          The number of rows and columns of the matrix A.  N >= 0.
00039 *
00040 *  KD      (input) INTEGER
00041 *          The number of super-diagonals of the matrix A if UPLO = 'U',
00042 *          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
00043 *
00044 *  A       (input) COMPLEX array, dimension (LDA,N)
00045 *          The original Hermitian band matrix A.  If UPLO = 'U', the
00046 *          upper triangular part of A is stored as a band matrix; if
00047 *          UPLO = 'L', the lower triangular part of A is stored.  The
00048 *          columns of the appropriate triangle are stored in the columns
00049 *          of A and the diagonals of the triangle are stored in the rows
00050 *          of A.  See CPBTRF for further details.
00051 *
00052 *  LDA     (input) INTEGER.
00053 *          The leading dimension of the array A.  LDA >= max(1,KD+1).
00054 *
00055 *  AFAC    (input) COMPLEX array, dimension (LDAFAC,N)
00056 *          The factored form of the matrix A.  AFAC contains the factor
00057 *          L or U from the L*L' or U'*U factorization in band storage
00058 *          format, as computed by CPBTRF.
00059 *
00060 *  LDAFAC  (input) INTEGER
00061 *          The leading dimension of the array AFAC.
00062 *          LDAFAC >= max(1,KD+1).
00063 *
00064 *  RWORK   (workspace) REAL array, dimension (N)
00065 *
00066 *  RESID   (output) REAL
00067 *          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00068 *          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00069 *
00070 *  =====================================================================
00071 *
00072 *
00073 *     .. Parameters ..
00074       REAL               ZERO, ONE
00075       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00076 *     ..
00077 *     .. Local Scalars ..
00078       INTEGER            I, J, K, KC, KLEN, ML, MU
00079       REAL               AKK, ANORM, EPS
00080 *     ..
00081 *     .. External Functions ..
00082       LOGICAL            LSAME
00083       REAL               CLANHB, SLAMCH
00084       COMPLEX            CDOTC
00085       EXTERNAL           LSAME, CLANHB, SLAMCH, CDOTC
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           CHER, CSSCAL, CTRMV
00089 *     ..
00090 *     .. Intrinsic Functions ..
00091       INTRINSIC          AIMAG, MAX, MIN, REAL
00092 *     ..
00093 *     .. Executable Statements ..
00094 *
00095 *     Quick exit if N = 0.
00096 *
00097       IF( N.LE.0 ) THEN
00098          RESID = ZERO
00099          RETURN
00100       END IF
00101 *
00102 *     Exit with RESID = 1/EPS if ANORM = 0.
00103 *
00104       EPS = SLAMCH( 'Epsilon' )
00105       ANORM = CLANHB( '1', UPLO, N, KD, A, LDA, RWORK )
00106       IF( ANORM.LE.ZERO ) THEN
00107          RESID = ONE / EPS
00108          RETURN
00109       END IF
00110 *
00111 *     Check the imaginary parts of the diagonal elements and return with
00112 *     an error code if any are nonzero.
00113 *
00114       IF( LSAME( UPLO, 'U' ) ) THEN
00115          DO 10 J = 1, N
00116             IF( AIMAG( AFAC( KD+1, J ) ).NE.ZERO ) THEN
00117                RESID = ONE / EPS
00118                RETURN
00119             END IF
00120    10    CONTINUE
00121       ELSE
00122          DO 20 J = 1, N
00123             IF( AIMAG( AFAC( 1, J ) ).NE.ZERO ) THEN
00124                RESID = ONE / EPS
00125                RETURN
00126             END IF
00127    20    CONTINUE
00128       END IF
00129 *
00130 *     Compute the product U'*U, overwriting U.
00131 *
00132       IF( LSAME( UPLO, 'U' ) ) THEN
00133          DO 30 K = N, 1, -1
00134             KC = MAX( 1, KD+2-K )
00135             KLEN = KD + 1 - KC
00136 *
00137 *           Compute the (K,K) element of the result.
00138 *
00139             AKK = CDOTC( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 )
00140             AFAC( KD+1, K ) = AKK
00141 *
00142 *           Compute the rest of column K.
00143 *
00144             IF( KLEN.GT.0 )
00145      $         CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', KLEN,
00146      $                     AFAC( KD+1, K-KLEN ), LDAFAC-1,
00147      $                     AFAC( KC, K ), 1 )
00148 *
00149    30    CONTINUE
00150 *
00151 *     UPLO = 'L':  Compute the product L*L', overwriting L.
00152 *
00153       ELSE
00154          DO 40 K = N, 1, -1
00155             KLEN = MIN( KD, N-K )
00156 *
00157 *           Add a multiple of column K of the factor L to each of
00158 *           columns K+1 through N.
00159 *
00160             IF( KLEN.GT.0 )
00161      $         CALL CHER( 'Lower', KLEN, ONE, AFAC( 2, K ), 1,
00162      $                    AFAC( 1, K+1 ), LDAFAC-1 )
00163 *
00164 *           Scale column K by the diagonal element.
00165 *
00166             AKK = AFAC( 1, K )
00167             CALL CSSCAL( KLEN+1, AKK, AFAC( 1, K ), 1 )
00168 *
00169    40    CONTINUE
00170       END IF
00171 *
00172 *     Compute the difference  L*L' - A  or  U'*U - A.
00173 *
00174       IF( LSAME( UPLO, 'U' ) ) THEN
00175          DO 60 J = 1, N
00176             MU = MAX( 1, KD+2-J )
00177             DO 50 I = MU, KD + 1
00178                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00179    50       CONTINUE
00180    60    CONTINUE
00181       ELSE
00182          DO 80 J = 1, N
00183             ML = MIN( KD+1, N-J+1 )
00184             DO 70 I = 1, ML
00185                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00186    70       CONTINUE
00187    80    CONTINUE
00188       END IF
00189 *
00190 *     Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
00191 *
00192       RESID = CLANHB( '1', UPLO, N, KD, AFAC, LDAFAC, RWORK )
00193 *
00194       RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00195 *
00196       RETURN
00197 *
00198 *     End of CPBT01
00199 *
00200       END
 All Files Functions