LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
spot03.f
Go to the documentation of this file.
00001 *> \brief \b SPOT03
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 SPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
00012 *                          RWORK, RCOND, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDA, LDAINV, LDWORK, N
00017 *       REAL               RCOND, RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
00021 *      $                   WORK( LDWORK, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> SPOT03 computes the residual for a symmetric matrix times its
00031 *> inverse:
00032 *>    norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
00033 *> where EPS is the machine epsilon.
00034 *> \endverbatim
00035 *
00036 *  Arguments:
00037 *  ==========
00038 *
00039 *> \param[in] UPLO
00040 *> \verbatim
00041 *>          UPLO is CHARACTER*1
00042 *>          Specifies whether the upper or lower triangular part of the
00043 *>          symmetric matrix A is stored:
00044 *>          = 'U':  Upper triangular
00045 *>          = 'L':  Lower triangular
00046 *> \endverbatim
00047 *>
00048 *> \param[in] N
00049 *> \verbatim
00050 *>          N is INTEGER
00051 *>          The number of rows and columns of the matrix A.  N >= 0.
00052 *> \endverbatim
00053 *>
00054 *> \param[in] A
00055 *> \verbatim
00056 *>          A is REAL array, dimension (LDA,N)
00057 *>          The original symmetric matrix A.
00058 *> \endverbatim
00059 *>
00060 *> \param[in] LDA
00061 *> \verbatim
00062 *>          LDA is INTEGER
00063 *>          The leading dimension of the array A.  LDA >= max(1,N)
00064 *> \endverbatim
00065 *>
00066 *> \param[in,out] AINV
00067 *> \verbatim
00068 *>          AINV is REAL array, dimension (LDAINV,N)
00069 *>          On entry, the inverse of the matrix A, stored as a symmetric
00070 *>          matrix in the same format as A.
00071 *>          In this version, AINV is expanded into a full matrix and
00072 *>          multiplied by A, so the opposing triangle of AINV will be
00073 *>          changed; i.e., if the upper triangular part of AINV is
00074 *>          stored, the lower triangular part will be used as work space.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] LDAINV
00078 *> \verbatim
00079 *>          LDAINV is INTEGER
00080 *>          The leading dimension of the array AINV.  LDAINV >= max(1,N).
00081 *> \endverbatim
00082 *>
00083 *> \param[out] WORK
00084 *> \verbatim
00085 *>          WORK is REAL array, dimension (LDWORK,N)
00086 *> \endverbatim
00087 *>
00088 *> \param[in] LDWORK
00089 *> \verbatim
00090 *>          LDWORK is INTEGER
00091 *>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00092 *> \endverbatim
00093 *>
00094 *> \param[out] RWORK
00095 *> \verbatim
00096 *>          RWORK is REAL array, dimension (N)
00097 *> \endverbatim
00098 *>
00099 *> \param[out] RCOND
00100 *> \verbatim
00101 *>          RCOND is REAL
00102 *>          The reciprocal of the condition number of A, computed as
00103 *>          ( 1/norm(A) ) / norm(AINV).
00104 *> \endverbatim
00105 *>
00106 *> \param[out] RESID
00107 *> \verbatim
00108 *>          RESID is REAL
00109 *>          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
00110 *> \endverbatim
00111 *
00112 *  Authors:
00113 *  ========
00114 *
00115 *> \author Univ. of Tennessee 
00116 *> \author Univ. of California Berkeley 
00117 *> \author Univ. of Colorado Denver 
00118 *> \author NAG Ltd. 
00119 *
00120 *> \date November 2011
00121 *
00122 *> \ingroup single_lin
00123 *
00124 *  =====================================================================
00125       SUBROUTINE SPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
00126      $                   RWORK, RCOND, RESID )
00127 *
00128 *  -- LAPACK test routine (version 3.4.0) --
00129 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00130 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00131 *     November 2011
00132 *
00133 *     .. Scalar Arguments ..
00134       CHARACTER          UPLO
00135       INTEGER            LDA, LDAINV, LDWORK, N
00136       REAL               RCOND, RESID
00137 *     ..
00138 *     .. Array Arguments ..
00139       REAL               A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
00140      $                   WORK( LDWORK, * )
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, J
00151       REAL               AINVNM, ANORM, EPS
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       REAL               SLAMCH, SLANGE, SLANSY
00156       EXTERNAL           LSAME, SLAMCH, SLANGE, SLANSY
00157 *     ..
00158 *     .. External Subroutines ..
00159       EXTERNAL           SSYMM
00160 *     ..
00161 *     .. Intrinsic Functions ..
00162       INTRINSIC          REAL
00163 *     ..
00164 *     .. Executable Statements ..
00165 *
00166 *     Quick exit if N = 0.
00167 *
00168       IF( N.LE.0 ) THEN
00169          RCOND = ONE
00170          RESID = ZERO
00171          RETURN
00172       END IF
00173 *
00174 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00175 *
00176       EPS = SLAMCH( 'Epsilon' )
00177       ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
00178       AINVNM = SLANSY( '1', UPLO, N, AINV, LDAINV, RWORK )
00179       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00180          RCOND = ZERO
00181          RESID = ONE / EPS
00182          RETURN
00183       END IF
00184       RCOND = ( ONE / ANORM ) / AINVNM
00185 *
00186 *     Expand AINV into a full matrix and call SSYMM to multiply
00187 *     AINV on the left by A.
00188 *
00189       IF( LSAME( UPLO, 'U' ) ) THEN
00190          DO 20 J = 1, N
00191             DO 10 I = 1, J - 1
00192                AINV( J, I ) = AINV( I, J )
00193    10       CONTINUE
00194    20    CONTINUE
00195       ELSE
00196          DO 40 J = 1, N
00197             DO 30 I = J + 1, N
00198                AINV( J, I ) = AINV( I, J )
00199    30       CONTINUE
00200    40    CONTINUE
00201       END IF
00202       CALL SSYMM( 'Left', UPLO, N, N, -ONE, A, LDA, AINV, LDAINV, ZERO,
00203      $            WORK, LDWORK )
00204 *
00205 *     Add the identity matrix to WORK .
00206 *
00207       DO 50 I = 1, N
00208          WORK( I, I ) = WORK( I, I ) + ONE
00209    50 CONTINUE
00210 *
00211 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
00212 *
00213       RESID = SLANGE( '1', N, N, WORK, LDWORK, RWORK )
00214 *
00215       RESID = ( ( RESID*RCOND ) / EPS ) / REAL( N )
00216 *
00217       RETURN
00218 *
00219 *     End of SPOT03
00220 *
00221       END
 All Files Functions