LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cspt03.f
Go to the documentation of this file.
00001 *> \brief \b CSPT03
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 CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            LDW, N
00017 *       REAL               RCOND, RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       REAL               RWORK( * )
00021 *       COMPLEX            A( * ), AINV( * ), WORK( LDW, * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> CSPT03 computes the residual for a complex symmetric packed matrix
00031 *> times its 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 *>          complex 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 COMPLEX array, dimension (N*(N+1)/2)
00057 *>          The original complex symmetric matrix A, stored as a packed
00058 *>          triangular matrix.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] AINV
00062 *> \verbatim
00063 *>          AINV is COMPLEX 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 COMPLEX array, dimension (LDW,N)
00071 *> \endverbatim
00072 *>
00073 *> \param[in] LDW
00074 *> \verbatim
00075 *>          LDW is INTEGER
00076 *>          The leading dimension of the array WORK.  LDW >= 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 complex_lin
00108 *
00109 *  =====================================================================
00110       SUBROUTINE CSPT03( UPLO, N, A, AINV, WORK, LDW, 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            LDW, N
00121       REAL               RCOND, RESID
00122 *     ..
00123 *     .. Array Arguments ..
00124       REAL               RWORK( * )
00125       COMPLEX            A( * ), AINV( * ), WORK( LDW, * )
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, ICOL, J, JCOL, K, KCOL, NALL
00136       REAL               AINVNM, ANORM, EPS
00137       COMPLEX            T
00138 *     ..
00139 *     .. External Functions ..
00140       LOGICAL            LSAME
00141       REAL               CLANGE, CLANSP, SLAMCH
00142       COMPLEX            CDOTU
00143       EXTERNAL           LSAME, CLANGE, CLANSP, SLAMCH, CDOTU
00144 *     ..
00145 *     .. Intrinsic Functions ..
00146       INTRINSIC          REAL
00147 *     ..
00148 *     .. Executable Statements ..
00149 *
00150 *     Quick exit if N = 0.
00151 *
00152       IF( N.LE.0 ) THEN
00153          RCOND = ONE
00154          RESID = ZERO
00155          RETURN
00156       END IF
00157 *
00158 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00159 *
00160       EPS = SLAMCH( 'Epsilon' )
00161       ANORM = CLANSP( '1', UPLO, N, A, RWORK )
00162       AINVNM = CLANSP( '1', UPLO, N, AINV, RWORK )
00163       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00164          RCOND = ZERO
00165          RESID = ONE / EPS
00166          RETURN
00167       END IF
00168       RCOND = ( ONE/ANORM ) / AINVNM
00169 *
00170 *     Case where both A and AINV are upper triangular:
00171 *     Each element of - A * AINV is computed by taking the dot product
00172 *     of a row of A with a column of AINV.
00173 *
00174       IF( LSAME( UPLO, 'U' ) ) THEN
00175          DO 70 I = 1, N
00176             ICOL = ( ( I-1 )*I ) / 2 + 1
00177 *
00178 *           Code when J <= I
00179 *
00180             DO 30 J = 1, I
00181                JCOL = ( ( J-1 )*J ) / 2 + 1
00182                T = CDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 )
00183                JCOL = JCOL + 2*J - 1
00184                KCOL = ICOL - 1
00185                DO 10 K = J + 1, I
00186                   T = T + A( KCOL+K )*AINV( JCOL )
00187                   JCOL = JCOL + K
00188    10          CONTINUE
00189                KCOL = KCOL + 2*I
00190                DO 20 K = I + 1, N
00191                   T = T + A( KCOL )*AINV( JCOL )
00192                   KCOL = KCOL + K
00193                   JCOL = JCOL + K
00194    20          CONTINUE
00195                WORK( I, J ) = -T
00196    30       CONTINUE
00197 *
00198 *           Code when J > I
00199 *
00200             DO 60 J = I + 1, N
00201                JCOL = ( ( J-1 )*J ) / 2 + 1
00202                T = CDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 )
00203                JCOL = JCOL - 1
00204                KCOL = ICOL + 2*I - 1
00205                DO 40 K = I + 1, J
00206                   T = T + A( KCOL )*AINV( JCOL+K )
00207                   KCOL = KCOL + K
00208    40          CONTINUE
00209                JCOL = JCOL + 2*J
00210                DO 50 K = J + 1, N
00211                   T = T + A( KCOL )*AINV( JCOL )
00212                   KCOL = KCOL + K
00213                   JCOL = JCOL + K
00214    50          CONTINUE
00215                WORK( I, J ) = -T
00216    60       CONTINUE
00217    70    CONTINUE
00218       ELSE
00219 *
00220 *        Case where both A and AINV are lower triangular
00221 *
00222          NALL = ( N*( N+1 ) ) / 2
00223          DO 140 I = 1, N
00224 *
00225 *           Code when J <= I
00226 *
00227             ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1
00228             DO 100 J = 1, I
00229                JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I )
00230                T = CDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 )
00231                KCOL = I
00232                JCOL = J
00233                DO 80 K = 1, J - 1
00234                   T = T + A( KCOL )*AINV( JCOL )
00235                   JCOL = JCOL + N - K
00236                   KCOL = KCOL + N - K
00237    80          CONTINUE
00238                JCOL = JCOL - J
00239                DO 90 K = J, I - 1
00240                   T = T + A( KCOL )*AINV( JCOL+K )
00241                   KCOL = KCOL + N - K
00242    90          CONTINUE
00243                WORK( I, J ) = -T
00244   100       CONTINUE
00245 *
00246 *           Code when J > I
00247 *
00248             ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2
00249             DO 130 J = I + 1, N
00250                JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1
00251                T = CDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 )
00252                KCOL = I
00253                JCOL = J
00254                DO 110 K = 1, I - 1
00255                   T = T + A( KCOL )*AINV( JCOL )
00256                   JCOL = JCOL + N - K
00257                   KCOL = KCOL + N - K
00258   110          CONTINUE
00259                KCOL = KCOL - I
00260                DO 120 K = I, J - 1
00261                   T = T + A( KCOL+K )*AINV( JCOL )
00262                   JCOL = JCOL + N - K
00263   120          CONTINUE
00264                WORK( I, J ) = -T
00265   130       CONTINUE
00266   140    CONTINUE
00267       END IF
00268 *
00269 *     Add the identity matrix to WORK .
00270 *
00271       DO 150 I = 1, N
00272          WORK( I, I ) = WORK( I, I ) + ONE
00273   150 CONTINUE
00274 *
00275 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
00276 *
00277       RESID = CLANGE( '1', N, N, WORK, LDW, RWORK )
00278 *
00279       RESID = ( ( RESID*RCOND )/EPS ) / REAL( N )
00280 *
00281       RETURN
00282 *
00283 *     End of CSPT03
00284 *
00285       END
 All Files Functions