LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zget01.f
Go to the documentation of this file.
00001 *> \brief \b ZGET01
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 ZGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            LDA, LDAFAC, M, N
00016 *       DOUBLE PRECISION   RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IPIV( * )
00020 *       DOUBLE PRECISION   RWORK( * )
00021 *       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> ZGET01 reconstructs a matrix A from its L*U factorization and
00031 *> computes the residual
00032 *>    norm(L*U - A) / ( N * norm(A) * EPS ),
00033 *> where EPS is the machine epsilon.
00034 *> \endverbatim
00035 *
00036 *  Arguments:
00037 *  ==========
00038 *
00039 *> \param[in] M
00040 *> \verbatim
00041 *>          M is INTEGER
00042 *>          The number of rows of the matrix A.  M >= 0.
00043 *> \endverbatim
00044 *>
00045 *> \param[in] N
00046 *> \verbatim
00047 *>          N is INTEGER
00048 *>          The number of columns of the matrix A.  N >= 0.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] A
00052 *> \verbatim
00053 *>          A is COMPLEX*16 array, dimension (LDA,N)
00054 *>          The original M x N matrix A.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] LDA
00058 *> \verbatim
00059 *>          LDA is INTEGER
00060 *>          The leading dimension of the array A.  LDA >= max(1,M).
00061 *> \endverbatim
00062 *>
00063 *> \param[in,out] AFAC
00064 *> \verbatim
00065 *>          AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
00066 *>          The factored form of the matrix A.  AFAC contains the factors
00067 *>          L and U from the L*U factorization as computed by ZGETRF.
00068 *>          Overwritten with the reconstructed matrix, and then with the
00069 *>          difference L*U - A.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] LDAFAC
00073 *> \verbatim
00074 *>          LDAFAC is INTEGER
00075 *>          The leading dimension of the array AFAC.  LDAFAC >= max(1,M).
00076 *> \endverbatim
00077 *>
00078 *> \param[in] IPIV
00079 *> \verbatim
00080 *>          IPIV is INTEGER array, dimension (N)
00081 *>          The pivot indices from ZGETRF.
00082 *> \endverbatim
00083 *>
00084 *> \param[out] RWORK
00085 *> \verbatim
00086 *>          RWORK is DOUBLE PRECISION array, dimension (M)
00087 *> \endverbatim
00088 *>
00089 *> \param[out] RESID
00090 *> \verbatim
00091 *>          RESID is DOUBLE PRECISION
00092 *>          norm(L*U - A) / ( N * norm(A) * EPS )
00093 *> \endverbatim
00094 *
00095 *  Authors:
00096 *  ========
00097 *
00098 *> \author Univ. of Tennessee 
00099 *> \author Univ. of California Berkeley 
00100 *> \author Univ. of Colorado Denver 
00101 *> \author NAG Ltd. 
00102 *
00103 *> \date November 2011
00104 *
00105 *> \ingroup complex16_lin
00106 *
00107 *  =====================================================================
00108       SUBROUTINE ZGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
00109      $                   RESID )
00110 *
00111 *  -- LAPACK test routine (version 3.4.0) --
00112 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00113 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00114 *     November 2011
00115 *
00116 *     .. Scalar Arguments ..
00117       INTEGER            LDA, LDAFAC, M, N
00118       DOUBLE PRECISION   RESID
00119 *     ..
00120 *     .. Array Arguments ..
00121       INTEGER            IPIV( * )
00122       DOUBLE PRECISION   RWORK( * )
00123       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * )
00124 *     ..
00125 *
00126 *  =====================================================================
00127 *
00128 *     .. Parameters ..
00129       DOUBLE PRECISION   ZERO, ONE
00130       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00131       COMPLEX*16         CONE
00132       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00133 *     ..
00134 *     .. Local Scalars ..
00135       INTEGER            I, J, K
00136       DOUBLE PRECISION   ANORM, EPS
00137       COMPLEX*16         T
00138 *     ..
00139 *     .. External Functions ..
00140       DOUBLE PRECISION   DLAMCH, ZLANGE
00141       COMPLEX*16         ZDOTU
00142       EXTERNAL           DLAMCH, ZLANGE, ZDOTU
00143 *     ..
00144 *     .. External Subroutines ..
00145       EXTERNAL           ZGEMV, ZLASWP, ZSCAL, ZTRMV
00146 *     ..
00147 *     .. Intrinsic Functions ..
00148       INTRINSIC          DBLE, MIN
00149 *     ..
00150 *     .. Executable Statements ..
00151 *
00152 *     Quick exit if M = 0 or N = 0.
00153 *
00154       IF( M.LE.0 .OR. N.LE.0 ) THEN
00155          RESID = ZERO
00156          RETURN
00157       END IF
00158 *
00159 *     Determine EPS and the norm of A.
00160 *
00161       EPS = DLAMCH( 'Epsilon' )
00162       ANORM = ZLANGE( '1', M, N, A, LDA, RWORK )
00163 *
00164 *     Compute the product L*U and overwrite AFAC with the result.
00165 *     A column at a time of the product is obtained, starting with
00166 *     column N.
00167 *
00168       DO 10 K = N, 1, -1
00169          IF( K.GT.M ) THEN
00170             CALL ZTRMV( 'Lower', 'No transpose', 'Unit', M, AFAC,
00171      $                  LDAFAC, AFAC( 1, K ), 1 )
00172          ELSE
00173 *
00174 *           Compute elements (K+1:M,K)
00175 *
00176             T = AFAC( K, K )
00177             IF( K+1.LE.M ) THEN
00178                CALL ZSCAL( M-K, T, AFAC( K+1, K ), 1 )
00179                CALL ZGEMV( 'No transpose', M-K, K-1, CONE,
00180      $                     AFAC( K+1, 1 ), LDAFAC, AFAC( 1, K ), 1,
00181      $                     CONE, AFAC( K+1, K ), 1 )
00182             END IF
00183 *
00184 *           Compute the (K,K) element
00185 *
00186             AFAC( K, K ) = T + ZDOTU( K-1, AFAC( K, 1 ), LDAFAC,
00187      $                     AFAC( 1, K ), 1 )
00188 *
00189 *           Compute elements (1:K-1,K)
00190 *
00191             CALL ZTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC,
00192      $                  LDAFAC, AFAC( 1, K ), 1 )
00193          END IF
00194    10 CONTINUE
00195       CALL ZLASWP( N, AFAC, LDAFAC, 1, MIN( M, N ), IPIV, -1 )
00196 *
00197 *     Compute the difference  L*U - A  and store in AFAC.
00198 *
00199       DO 30 J = 1, N
00200          DO 20 I = 1, M
00201             AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00202    20    CONTINUE
00203    30 CONTINUE
00204 *
00205 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00206 *
00207       RESID = ZLANGE( '1', M, N, AFAC, LDAFAC, RWORK )
00208 *
00209       IF( ANORM.LE.ZERO ) THEN
00210          IF( RESID.NE.ZERO )
00211      $      RESID = ONE / EPS
00212       ELSE
00213          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00214       END IF
00215 *
00216       RETURN
00217 *
00218 *     End of ZGET01
00219 *
00220       END
 All Files Functions