TESTING/LIN/cget01.f

   1 *> \brief \b CGET01
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE CGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
  12 *                          RESID )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER            LDA, LDAFAC, M, N
  16 *       REAL               RESID
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       INTEGER            IPIV( * )
  20 *       REAL               RWORK( * )
  21 *       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * )
  22 *       ..
  23 *
  24 *
  25 *> \par Purpose:
  26 *  =============
  27 *>
  28 *> \verbatim
  29 *>
  30 *> CGET01 reconstructs a matrix A from its L*U factorization and
  31 *> computes the residual
  32 *>    norm(L*U - A) / ( N * norm(A) * EPS ),
  33 *> where EPS is the machine epsilon.
  34 *> \endverbatim
  35 *
  36 *  Arguments:
  37 *  ==========
  38 *
  39 *> \param[in] M
  40 *> \verbatim
  41 *>          M is INTEGER
  42 *>          The number of rows of the matrix A.  M >= 0.
  43 *> \endverbatim
  44 *>
  45 *> \param[in] N
  46 *> \verbatim
  47 *>          N is INTEGER
  48 *>          The number of columns of the matrix A.  N >= 0.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] A
  52 *> \verbatim
  53 *>          A is COMPLEX array, dimension (LDA,N)
  54 *>          The original M x N matrix A.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] LDA
  58 *> \verbatim
  59 *>          LDA is INTEGER
  60 *>          The leading dimension of the array A.  LDA >= max(1,M).
  61 *> \endverbatim
  62 *>
  63 *> \param[in,out] AFAC
  64 *> \verbatim
  65 *>          AFAC is COMPLEX array, dimension (LDAFAC,N)
  66 *>          The factored form of the matrix A.  AFAC contains the factors
  67 *>          L and U from the L*U factorization as computed by CGETRF.
  68 *>          Overwritten with the reconstructed matrix, and then with the
  69 *>          difference L*U - A.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] LDAFAC
  73 *> \verbatim
  74 *>          LDAFAC is INTEGER
  75 *>          The leading dimension of the array AFAC.  LDAFAC >= max(1,M).
  76 *> \endverbatim
  77 *>
  78 *> \param[in] IPIV
  79 *> \verbatim
  80 *>          IPIV is INTEGER array, dimension (N)
  81 *>          The pivot indices from CGETRF.
  82 *> \endverbatim
  83 *>
  84 *> \param[out] RWORK
  85 *> \verbatim
  86 *>          RWORK is REAL array, dimension (M)
  87 *> \endverbatim
  88 *>
  89 *> \param[out] RESID
  90 *> \verbatim
  91 *>          RESID is REAL
  92 *>          norm(L*U - A) / ( N * norm(A) * EPS )
  93 *> \endverbatim
  94 *
  95 *  Authors:
  96 *  ========
  97 *
  98 *> \author Univ. of Tennessee
  99 *> \author Univ. of California Berkeley
 100 *> \author Univ. of Colorado Denver
 101 *> \author NAG Ltd.
 102 *
 103 *> \date November 2011
 104 *
 105 *> \ingroup complex_lin
 106 *
 107 *  =====================================================================
 108       SUBROUTINE CGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
 109      $                   RESID )
 110 *
 111 *  -- LAPACK test routine (version 3.4.0) --
 112 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 113 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 114 *     November 2011
 115 *
 116 *     .. Scalar Arguments ..
 117       INTEGER            LDA, LDAFAC, M, N
 118       REAL               RESID
 119 *     ..
 120 *     .. Array Arguments ..
 121       INTEGER            IPIV( * )
 122       REAL               RWORK( * )
 123       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * )
 124 *     ..
 125 *
 126 *  =====================================================================
 127 *
 128 *     .. Parameters ..
 129       REAL               ONE, ZERO
 130       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 131       COMPLEX            CONE
 132       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 133 *     ..
 134 *     .. Local Scalars ..
 135       INTEGER            I, J, K
 136       REAL               ANORM, EPS
 137       COMPLEX            T
 138 *     ..
 139 *     .. External Functions ..
 140       REAL               CLANGE, SLAMCH
 141       COMPLEX            CDOTU
 142       EXTERNAL           CLANGE, SLAMCH, CDOTU
 143 *     ..
 144 *     .. External Subroutines ..
 145       EXTERNAL           CGEMV, CLASWP, CSCAL, CTRMV
 146 *     ..
 147 *     .. Intrinsic Functions ..
 148       INTRINSIC          MIN, REAL
 149 *     ..
 150 *     .. Executable Statements ..
 151 *
 152 *     Quick exit if M = 0 or N = 0.
 153 *
 154       IF( M.LE.0 .OR. N.LE.0 ) THEN
 155          RESID = ZERO
 156          RETURN
 157       END IF
 158 *
 159 *     Determine EPS and the norm of A.
 160 *
 161       EPS = SLAMCH( 'Epsilon' )
 162       ANORM = CLANGE( '1', M, N, A, LDA, RWORK )
 163 *
 164 *     Compute the product L*U and overwrite AFAC with the result.
 165 *     A column at a time of the product is obtained, starting with
 166 *     column N.
 167 *
 168       DO 10 K = N, 1, -1
 169          IF( K.GT.M ) THEN
 170             CALL CTRMV( 'Lower', 'No transpose', 'Unit', M, AFAC,
 171      $                  LDAFAC, AFAC( 1, K ), 1 )
 172          ELSE
 173 *
 174 *           Compute elements (K+1:M,K)
 175 *
 176             T = AFAC( K, K )
 177             IF( K+1.LE.M ) THEN
 178                CALL CSCAL( M-K, T, AFAC( K+1, K ), 1 )
 179                CALL CGEMV( 'No transpose', M-K, K-1, CONE,
 180      $                     AFAC( K+1, 1 ), LDAFAC, AFAC( 1, K ), 1,
 181      $                     CONE, AFAC( K+1, K ), 1 )
 182             END IF
 183 *
 184 *           Compute the (K,K) element
 185 *
 186             AFAC( K, K ) = T + CDOTU( K-1, AFAC( K, 1 ), LDAFAC,
 187      $                     AFAC( 1, K ), 1 )
 188 *
 189 *           Compute elements (1:K-1,K)
 190 *
 191             CALL CTRMV( 'Lower', 'No transpose', 'Unit', K-1, AFAC,
 192      $                  LDAFAC, AFAC( 1, K ), 1 )
 193          END IF
 194    10 CONTINUE
 195       CALL CLASWP( N, AFAC, LDAFAC, 1, MIN( M, N ), IPIV, -1 )
 196 *
 197 *     Compute the difference  L*U - A  and store in AFAC.
 198 *
 199       DO 30 J = 1, N
 200          DO 20 I = 1, M
 201             AFAC( I, J ) = AFAC( I, J ) - A( I, J )
 202    20    CONTINUE
 203    30 CONTINUE
 204 *
 205 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
 206 *
 207       RESID = CLANGE( '1', M, N, AFAC, LDAFAC, RWORK )
 208 *
 209       IF( ANORM.LE.ZERO ) THEN
 210          IF( RESID.NE.ZERO )
 211      $      RESID = ONE / EPS
 212       ELSE
 213          RESID = ( ( RESID/REAL( N ) )/ANORM ) / EPS
 214       END IF
 215 *
 216       RETURN
 217 *
 218 *     End of CGET01
 219 *
 220       END