LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sppt03.f
Go to the documentation of this file.
00001 *> \brief \b SPPT03
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 SPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDWORK, N
00017 *       REAL               RCOND, RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               A( * ), AINV( * ), RWORK( * ),
00021 *      $                   WORK( LDWORK, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> SPPT03 computes the residual for a symmetric packed 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 (N*(N+1)/2)
00057 *>          The original symmetric matrix A, stored as a packed
00058 *>          triangular matrix.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] AINV
00062 *> \verbatim
00063 *>          AINV is REAL array, dimension (N*(N+1)/2)
00064 *>          The (symmetric) inverse of the matrix A, stored as a packed
00065 *>          triangular matrix.
00066 *> \endverbatim
00067 *>
00068 *> \param[out] WORK
00069 *> \verbatim
00070 *>          WORK is REAL array, dimension (LDWORK,N)
00071 *> \endverbatim
00072 *>
00073 *> \param[in] LDWORK
00074 *> \verbatim
00075 *>          LDWORK is INTEGER
00076 *>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00077 *> \endverbatim
00078 *>
00079 *> \param[out] RWORK
00080 *> \verbatim
00081 *>          RWORK is REAL array, dimension (N)
00082 *> \endverbatim
00083 *>
00084 *> \param[out] RCOND
00085 *> \verbatim
00086 *>          RCOND is REAL
00087 *>          The reciprocal of the condition number of A, computed as
00088 *>          ( 1/norm(A) ) / norm(AINV).
00089 *> \endverbatim
00090 *>
00091 *> \param[out] RESID
00092 *> \verbatim
00093 *>          RESID is REAL
00094 *>          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
00095 *> \endverbatim
00096 *
00097 *  Authors:
00098 *  ========
00099 *
00100 *> \author Univ. of Tennessee 
00101 *> \author Univ. of California Berkeley 
00102 *> \author Univ. of Colorado Denver 
00103 *> \author NAG Ltd. 
00104 *
00105 *> \date November 2011
00106 *
00107 *> \ingroup single_lin
00108 *
00109 *  =====================================================================
00110       SUBROUTINE SPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
00111      $                   RESID )
00112 *
00113 *  -- LAPACK test routine (version 3.4.0) --
00114 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00115 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00116 *     November 2011
00117 *
00118 *     .. Scalar Arguments ..
00119       CHARACTER          UPLO
00120       INTEGER            LDWORK, N
00121       REAL               RCOND, RESID
00122 *     ..
00123 *     .. Array Arguments ..
00124       REAL               A( * ), AINV( * ), RWORK( * ),
00125      $                   WORK( LDWORK, * )
00126 *     ..
00127 *
00128 *  =====================================================================
00129 *
00130 *     .. Parameters ..
00131       REAL               ZERO, ONE
00132       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00133 *     ..
00134 *     .. Local Scalars ..
00135       INTEGER            I, J, JJ
00136       REAL               AINVNM, ANORM, EPS
00137 *     ..
00138 *     .. External Functions ..
00139       LOGICAL            LSAME
00140       REAL               SLAMCH, SLANGE, SLANSP
00141       EXTERNAL           LSAME, SLAMCH, SLANGE, SLANSP
00142 *     ..
00143 *     .. Intrinsic Functions ..
00144       INTRINSIC          REAL
00145 *     ..
00146 *     .. External Subroutines ..
00147       EXTERNAL           SCOPY, SSPMV
00148 *     ..
00149 *     .. Executable Statements ..
00150 *
00151 *     Quick exit if N = 0.
00152 *
00153       IF( N.LE.0 ) THEN
00154          RCOND = ONE
00155          RESID = ZERO
00156          RETURN
00157       END IF
00158 *
00159 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00160 *
00161       EPS = SLAMCH( 'Epsilon' )
00162       ANORM = SLANSP( '1', UPLO, N, A, RWORK )
00163       AINVNM = SLANSP( '1', UPLO, N, AINV, RWORK )
00164       IF( ANORM.LE.ZERO .OR. AINVNM.EQ.ZERO ) THEN
00165          RCOND = ZERO
00166          RESID = ONE / EPS
00167          RETURN
00168       END IF
00169       RCOND = ( ONE / ANORM ) / AINVNM
00170 *
00171 *     UPLO = 'U':
00172 *     Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
00173 *     expand it to a full matrix, then multiply by A one column at a
00174 *     time, moving the result one column to the left.
00175 *
00176       IF( LSAME( UPLO, 'U' ) ) THEN
00177 *
00178 *        Copy AINV
00179 *
00180          JJ = 1
00181          DO 10 J = 1, N - 1
00182             CALL SCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
00183             CALL SCOPY( J-1, AINV( JJ ), 1, WORK( J, 2 ), LDWORK )
00184             JJ = JJ + J
00185    10    CONTINUE
00186          JJ = ( ( N-1 )*N ) / 2 + 1
00187          CALL SCOPY( N-1, AINV( JJ ), 1, WORK( N, 2 ), LDWORK )
00188 *
00189 *        Multiply by A
00190 *
00191          DO 20 J = 1, N - 1
00192             CALL SSPMV( 'Upper', N, -ONE, A, WORK( 1, J+1 ), 1, ZERO,
00193      $                  WORK( 1, J ), 1 )
00194    20    CONTINUE
00195          CALL SSPMV( 'Upper', N, -ONE, A, AINV( JJ ), 1, ZERO,
00196      $               WORK( 1, N ), 1 )
00197 *
00198 *     UPLO = 'L':
00199 *     Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
00200 *     and multiply by A, moving each column to the right.
00201 *
00202       ELSE
00203 *
00204 *        Copy AINV
00205 *
00206          CALL SCOPY( N-1, AINV( 2 ), 1, WORK( 1, 1 ), LDWORK )
00207          JJ = N + 1
00208          DO 30 J = 2, N
00209             CALL SCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
00210             CALL SCOPY( N-J, AINV( JJ+1 ), 1, WORK( J, J ), LDWORK )
00211             JJ = JJ + N - J + 1
00212    30    CONTINUE
00213 *
00214 *        Multiply by A
00215 *
00216          DO 40 J = N, 2, -1
00217             CALL SSPMV( 'Lower', N, -ONE, A, WORK( 1, J-1 ), 1, ZERO,
00218      $                  WORK( 1, J ), 1 )
00219    40    CONTINUE
00220          CALL SSPMV( 'Lower', N, -ONE, A, AINV( 1 ), 1, ZERO,
00221      $               WORK( 1, 1 ), 1 )
00222 *
00223       END IF
00224 *
00225 *     Add the identity matrix to WORK .
00226 *
00227       DO 50 I = 1, N
00228          WORK( I, I ) = WORK( I, I ) + ONE
00229    50 CONTINUE
00230 *
00231 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
00232 *
00233       RESID = SLANGE( '1', N, N, WORK, LDWORK, RWORK )
00234 *
00235       RESID = ( ( RESID*RCOND ) / EPS ) / REAL( N )
00236 *
00237       RETURN
00238 *
00239 *     End of SPPT03
00240 *
00241       END
 All Files Functions