LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cgbt01.f
Go to the documentation of this file.
00001 *> \brief \b CGBT01
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 CGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            KL, KU, LDA, LDAFAC, M, N
00016 *       REAL               RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IPIV( * )
00020 *       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> CGBT01 reconstructs a band matrix  A  from its L*U factorization and
00030 *> computes the residual:
00031 *>    norm(L*U - A) / ( N * norm(A) * EPS ),
00032 *> where EPS is the machine epsilon.
00033 *>
00034 *> The expression L*U - A is computed one column at a time, so A and
00035 *> AFAC are not modified.
00036 *> \endverbatim
00037 *
00038 *  Arguments:
00039 *  ==========
00040 *
00041 *> \param[in] M
00042 *> \verbatim
00043 *>          M is INTEGER
00044 *>          The number of rows of the matrix A.  M >= 0.
00045 *> \endverbatim
00046 *>
00047 *> \param[in] N
00048 *> \verbatim
00049 *>          N is INTEGER
00050 *>          The number of columns of the matrix A.  N >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] KL
00054 *> \verbatim
00055 *>          KL is INTEGER
00056 *>          The number of subdiagonals within the band of A.  KL >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] KU
00060 *> \verbatim
00061 *>          KU is INTEGER
00062 *>          The number of superdiagonals within the band of A.  KU >= 0.
00063 *> \endverbatim
00064 *>
00065 *> \param[in,out] A
00066 *> \verbatim
00067 *>          A is COMPLEX array, dimension (LDA,N)
00068 *>          The original matrix A in band storage, stored in rows 1 to
00069 *>          KL+KU+1.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] LDA
00073 *> \verbatim
00074 *>          LDA is INTEGER.
00075 *>          The leading dimension of the array A.  LDA >= max(1,KL+KU+1).
00076 *> \endverbatim
00077 *>
00078 *> \param[in] AFAC
00079 *> \verbatim
00080 *>          AFAC is COMPLEX array, dimension (LDAFAC,N)
00081 *>          The factored form of the matrix A.  AFAC contains the banded
00082 *>          factors L and U from the L*U factorization, as computed by
00083 *>          CGBTRF.  U is stored as an upper triangular band matrix with
00084 *>          KL+KU superdiagonals in rows 1 to KL+KU+1, and the
00085 *>          multipliers used during the factorization are stored in rows
00086 *>          KL+KU+2 to 2*KL+KU+1.  See CGBTRF for further details.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LDAFAC
00090 *> \verbatim
00091 *>          LDAFAC is INTEGER
00092 *>          The leading dimension of the array AFAC.
00093 *>          LDAFAC >= max(1,2*KL*KU+1).
00094 *> \endverbatim
00095 *>
00096 *> \param[in] IPIV
00097 *> \verbatim
00098 *>          IPIV is INTEGER array, dimension (min(M,N))
00099 *>          The pivot indices from CGBTRF.
00100 *> \endverbatim
00101 *>
00102 *> \param[out] WORK
00103 *> \verbatim
00104 *>          WORK is COMPLEX array, dimension (2*KL+KU+1)
00105 *> \endverbatim
00106 *>
00107 *> \param[out] RESID
00108 *> \verbatim
00109 *>          RESID is REAL
00110 *>          norm(L*U - A) / ( N * norm(A) * EPS )
00111 *> \endverbatim
00112 *
00113 *  Authors:
00114 *  ========
00115 *
00116 *> \author Univ. of Tennessee 
00117 *> \author Univ. of California Berkeley 
00118 *> \author Univ. of Colorado Denver 
00119 *> \author NAG Ltd. 
00120 *
00121 *> \date November 2011
00122 *
00123 *> \ingroup complex_lin
00124 *
00125 *  =====================================================================
00126       SUBROUTINE CGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
00127      $                   RESID )
00128 *
00129 *  -- LAPACK test routine (version 3.4.0) --
00130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00132 *     November 2011
00133 *
00134 *     .. Scalar Arguments ..
00135       INTEGER            KL, KU, LDA, LDAFAC, M, N
00136       REAL               RESID
00137 *     ..
00138 *     .. Array Arguments ..
00139       INTEGER            IPIV( * )
00140       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
00141 *     ..
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       REAL               ZERO, ONE
00147       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00148 *     ..
00149 *     .. Local Scalars ..
00150       INTEGER            I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
00151       REAL               ANORM, EPS
00152       COMPLEX            T
00153 *     ..
00154 *     .. External Functions ..
00155       REAL               SCASUM, SLAMCH
00156       EXTERNAL           SCASUM, SLAMCH
00157 *     ..
00158 *     .. External Subroutines ..
00159       EXTERNAL           CAXPY, CCOPY
00160 *     ..
00161 *     .. Intrinsic Functions ..
00162       INTRINSIC          CMPLX, MAX, MIN, REAL
00163 *     ..
00164 *     .. Executable Statements ..
00165 *
00166 *     Quick exit if M = 0 or N = 0.
00167 *
00168       RESID = ZERO
00169       IF( M.LE.0 .OR. N.LE.0 )
00170      $   RETURN
00171 *
00172 *     Determine EPS and the norm of A.
00173 *
00174       EPS = SLAMCH( 'Epsilon' )
00175       KD = KU + 1
00176       ANORM = ZERO
00177       DO 10 J = 1, N
00178          I1 = MAX( KD+1-J, 1 )
00179          I2 = MIN( KD+M-J, KL+KD )
00180          IF( I2.GE.I1 )
00181      $      ANORM = MAX( ANORM, SCASUM( I2-I1+1, A( I1, J ), 1 ) )
00182    10 CONTINUE
00183 *
00184 *     Compute one column at a time of L*U - A.
00185 *
00186       KD = KL + KU + 1
00187       DO 40 J = 1, N
00188 *
00189 *        Copy the J-th column of U to WORK.
00190 *
00191          JU = MIN( KL+KU, J-1 )
00192          JL = MIN( KL, M-J )
00193          LENJ = MIN( M, J ) - J + JU + 1
00194          IF( LENJ.GT.0 ) THEN
00195             CALL CCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
00196             DO 20 I = LENJ + 1, JU + JL + 1
00197                WORK( I ) = ZERO
00198    20       CONTINUE
00199 *
00200 *           Multiply by the unit lower triangular matrix L.  Note that L
00201 *           is stored as a product of transformations and permutations.
00202 *
00203             DO 30 I = MIN( M-1, J ), J - JU, -1
00204                IL = MIN( KL, M-I )
00205                IF( IL.GT.0 ) THEN
00206                   IW = I - J + JU + 1
00207                   T = WORK( IW )
00208                   CALL CAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
00209      $                        1 )
00210                   IP = IPIV( I )
00211                   IF( I.NE.IP ) THEN
00212                      IP = IP - J + JU + 1
00213                      WORK( IW ) = WORK( IP )
00214                      WORK( IP ) = T
00215                   END IF
00216                END IF
00217    30       CONTINUE
00218 *
00219 *           Subtract the corresponding column of A.
00220 *
00221             JUA = MIN( JU, KU )
00222             IF( JUA+JL+1.GT.0 )
00223      $         CALL CAXPY( JUA+JL+1, -CMPLX( ONE ), A( KU+1-JUA, J ), 1,
00224      $                     WORK( JU+1-JUA ), 1 )
00225 *
00226 *           Compute the 1-norm of the column.
00227 *
00228             RESID = MAX( RESID, SCASUM( JU+JL+1, WORK, 1 ) )
00229          END IF
00230    40 CONTINUE
00231 *
00232 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00233 *
00234       IF( ANORM.LE.ZERO ) THEN
00235          IF( RESID.NE.ZERO )
00236      $      RESID = ONE / EPS
00237       ELSE
00238          RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS
00239       END IF
00240 *
00241       RETURN
00242 *
00243 *     End of CGBT01
00244 *
00245       END
 All Files Functions