LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cpbt01.f
Go to the documentation of this file.
00001 *> \brief \b CPBT01
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            KD, LDA, LDAFAC, N
00017 *       REAL               RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               RWORK( * )
00021 *       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> CPBT01 reconstructs a Hermitian positive definite band matrix A from
00031 *> its L*L' or U'*U factorization and computes the residual
00032 *>    norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00033 *>    norm( U'*U - A ) / ( N * norm(A) * EPS ),
00034 *> where EPS is the machine epsilon, L' is the conjugate transpose of
00035 *> L, and U' is the conjugate transpose of U.
00036 *> \endverbatim
00037 *
00038 *  Arguments:
00039 *  ==========
00040 *
00041 *> \param[in] UPLO
00042 *> \verbatim
00043 *>          UPLO is CHARACTER*1
00044 *>          Specifies whether the upper or lower triangular part of the
00045 *>          Hermitian matrix A is stored:
00046 *>          = 'U':  Upper triangular
00047 *>          = 'L':  Lower triangular
00048 *> \endverbatim
00049 *>
00050 *> \param[in] N
00051 *> \verbatim
00052 *>          N is INTEGER
00053 *>          The number of rows and columns of the matrix A.  N >= 0.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] KD
00057 *> \verbatim
00058 *>          KD is INTEGER
00059 *>          The number of super-diagonals of the matrix A if UPLO = 'U',
00060 *>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] A
00064 *> \verbatim
00065 *>          A is COMPLEX array, dimension (LDA,N)
00066 *>          The original Hermitian band matrix A.  If UPLO = 'U', the
00067 *>          upper triangular part of A is stored as a band matrix; if
00068 *>          UPLO = 'L', the lower triangular part of A is stored.  The
00069 *>          columns of the appropriate triangle are stored in the columns
00070 *>          of A and the diagonals of the triangle are stored in the rows
00071 *>          of A.  See CPBTRF for further details.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] LDA
00075 *> \verbatim
00076 *>          LDA is INTEGER.
00077 *>          The leading dimension of the array A.  LDA >= max(1,KD+1).
00078 *> \endverbatim
00079 *>
00080 *> \param[in] AFAC
00081 *> \verbatim
00082 *>          AFAC is COMPLEX array, dimension (LDAFAC,N)
00083 *>          The factored form of the matrix A.  AFAC contains the factor
00084 *>          L or U from the L*L' or U'*U factorization in band storage
00085 *>          format, as computed by CPBTRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] LDAFAC
00089 *> \verbatim
00090 *>          LDAFAC is INTEGER
00091 *>          The leading dimension of the array AFAC.
00092 *>          LDAFAC >= max(1,KD+1).
00093 *> \endverbatim
00094 *>
00095 *> \param[out] RWORK
00096 *> \verbatim
00097 *>          RWORK is REAL array, dimension (N)
00098 *> \endverbatim
00099 *>
00100 *> \param[out] RESID
00101 *> \verbatim
00102 *>          RESID is REAL
00103 *>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00104 *>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00105 *> \endverbatim
00106 *
00107 *  Authors:
00108 *  ========
00109 *
00110 *> \author Univ. of Tennessee 
00111 *> \author Univ. of California Berkeley 
00112 *> \author Univ. of Colorado Denver 
00113 *> \author NAG Ltd. 
00114 *
00115 *> \date November 2011
00116 *
00117 *> \ingroup complex_lin
00118 *
00119 *  =====================================================================
00120       SUBROUTINE CPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
00121      $                   RESID )
00122 *
00123 *  -- LAPACK test routine (version 3.4.0) --
00124 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00125 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00126 *     November 2011
00127 *
00128 *     .. Scalar Arguments ..
00129       CHARACTER          UPLO
00130       INTEGER            KD, LDA, LDAFAC, N
00131       REAL               RESID
00132 *     ..
00133 *     .. Array Arguments ..
00134       REAL               RWORK( * )
00135       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * )
00136 *     ..
00137 *
00138 *  =====================================================================
00139 *
00140 *
00141 *     .. Parameters ..
00142       REAL               ZERO, ONE
00143       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00144 *     ..
00145 *     .. Local Scalars ..
00146       INTEGER            I, J, K, KC, KLEN, ML, MU
00147       REAL               AKK, ANORM, EPS
00148 *     ..
00149 *     .. External Functions ..
00150       LOGICAL            LSAME
00151       REAL               CLANHB, SLAMCH
00152       COMPLEX            CDOTC
00153       EXTERNAL           LSAME, CLANHB, SLAMCH, CDOTC
00154 *     ..
00155 *     .. External Subroutines ..
00156       EXTERNAL           CHER, CSSCAL, CTRMV
00157 *     ..
00158 *     .. Intrinsic Functions ..
00159       INTRINSIC          AIMAG, MAX, MIN, REAL
00160 *     ..
00161 *     .. Executable Statements ..
00162 *
00163 *     Quick exit if N = 0.
00164 *
00165       IF( N.LE.0 ) THEN
00166          RESID = ZERO
00167          RETURN
00168       END IF
00169 *
00170 *     Exit with RESID = 1/EPS if ANORM = 0.
00171 *
00172       EPS = SLAMCH( 'Epsilon' )
00173       ANORM = CLANHB( '1', UPLO, N, KD, A, LDA, RWORK )
00174       IF( ANORM.LE.ZERO ) THEN
00175          RESID = ONE / EPS
00176          RETURN
00177       END IF
00178 *
00179 *     Check the imaginary parts of the diagonal elements and return with
00180 *     an error code if any are nonzero.
00181 *
00182       IF( LSAME( UPLO, 'U' ) ) THEN
00183          DO 10 J = 1, N
00184             IF( AIMAG( AFAC( KD+1, J ) ).NE.ZERO ) THEN
00185                RESID = ONE / EPS
00186                RETURN
00187             END IF
00188    10    CONTINUE
00189       ELSE
00190          DO 20 J = 1, N
00191             IF( AIMAG( AFAC( 1, J ) ).NE.ZERO ) THEN
00192                RESID = ONE / EPS
00193                RETURN
00194             END IF
00195    20    CONTINUE
00196       END IF
00197 *
00198 *     Compute the product U'*U, overwriting U.
00199 *
00200       IF( LSAME( UPLO, 'U' ) ) THEN
00201          DO 30 K = N, 1, -1
00202             KC = MAX( 1, KD+2-K )
00203             KLEN = KD + 1 - KC
00204 *
00205 *           Compute the (K,K) element of the result.
00206 *
00207             AKK = CDOTC( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 )
00208             AFAC( KD+1, K ) = AKK
00209 *
00210 *           Compute the rest of column K.
00211 *
00212             IF( KLEN.GT.0 )
00213      $         CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', KLEN,
00214      $                     AFAC( KD+1, K-KLEN ), LDAFAC-1,
00215      $                     AFAC( KC, K ), 1 )
00216 *
00217    30    CONTINUE
00218 *
00219 *     UPLO = 'L':  Compute the product L*L', overwriting L.
00220 *
00221       ELSE
00222          DO 40 K = N, 1, -1
00223             KLEN = MIN( KD, N-K )
00224 *
00225 *           Add a multiple of column K of the factor L to each of
00226 *           columns K+1 through N.
00227 *
00228             IF( KLEN.GT.0 )
00229      $         CALL CHER( 'Lower', KLEN, ONE, AFAC( 2, K ), 1,
00230      $                    AFAC( 1, K+1 ), LDAFAC-1 )
00231 *
00232 *           Scale column K by the diagonal element.
00233 *
00234             AKK = AFAC( 1, K )
00235             CALL CSSCAL( KLEN+1, AKK, AFAC( 1, K ), 1 )
00236 *
00237    40    CONTINUE
00238       END IF
00239 *
00240 *     Compute the difference  L*L' - A  or  U'*U - A.
00241 *
00242       IF( LSAME( UPLO, 'U' ) ) THEN
00243          DO 60 J = 1, N
00244             MU = MAX( 1, KD+2-J )
00245             DO 50 I = MU, KD + 1
00246                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00247    50       CONTINUE
00248    60    CONTINUE
00249       ELSE
00250          DO 80 J = 1, N
00251             ML = MIN( KD+1, N-J+1 )
00252             DO 70 I = 1, ML
00253                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00254    70       CONTINUE
00255    80    CONTINUE
00256       END IF
00257 *
00258 *     Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
00259 *
00260       RESID = CLANHB( '1', UPLO, N, KD, AFAC, LDAFAC, RWORK )
00261 *
00262       RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00263 *
00264       RETURN
00265 *
00266 *     End of CPBT01
00267 *
00268       END
 All Files Functions