LAPACK: TESTING/LIN/cspt03.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
00002      $                   RESID )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            LDW, N
00011       REAL               RCOND, RESID
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               RWORK( * )
00015       COMPLEX            A( * ), AINV( * ), WORK( LDW, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CSPT03 computes the residual for a complex symmetric packed matrix
00022 *  times its inverse:
00023 *     norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
00024 *  where EPS is the machine epsilon.
00025 *
00026 *  Arguments
00027 *  ==========
00028 *
00029 *  UPLO    (input) CHARACTER*1
00030 *          Specifies whether the upper or lower triangular part of the
00031 *          complex symmetric matrix A is stored:
00032 *          = 'U':  Upper triangular
00033 *          = 'L':  Lower triangular
00034 *
00035 *  N       (input) INTEGER
00036 *          The number of rows and columns of the matrix A.  N >= 0.
00037 *
00038 *  A       (input) COMPLEX array, dimension (N*(N+1)/2)
00039 *          The original complex symmetric matrix A, stored as a packed
00040 *          triangular matrix.
00041 *
00042 *  AINV    (input) COMPLEX array, dimension (N*(N+1)/2)
00043 *          The (symmetric) inverse of the matrix A, stored as a packed
00044 *          triangular matrix.
00045 *
00046 *  WORK    (workspace) COMPLEX array, dimension (LDWORK,N)
00047 *
00048 *  LDWORK  (input) INTEGER
00049 *          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00050 *
00051 *  RWORK   (workspace) REAL array, dimension (N)
00052 *
00053 *  RCOND   (output) REAL
00054 *          The reciprocal of the condition number of A, computed as
00055 *          ( 1/norm(A) ) / norm(AINV).
00056 *
00057 *  RESID   (output) REAL
00058 *          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
00059 *
00060 *  =====================================================================
00061 *
00062 *     .. Parameters ..
00063       REAL               ZERO, ONE
00064       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00065 *     ..
00066 *     .. Local Scalars ..
00067       INTEGER            I, ICOL, J, JCOL, K, KCOL, NALL
00068       REAL               AINVNM, ANORM, EPS
00069       COMPLEX            T
00070 *     ..
00071 *     .. External Functions ..
00072       LOGICAL            LSAME
00073       REAL               CLANGE, CLANSP, SLAMCH
00074       COMPLEX            CDOTU
00075       EXTERNAL           LSAME, CLANGE, CLANSP, SLAMCH, CDOTU
00076 *     ..
00077 *     .. Intrinsic Functions ..
00078       INTRINSIC          REAL
00079 *     ..
00080 *     .. Executable Statements ..
00081 *
00082 *     Quick exit if N = 0.
00083 *
00084       IF( N.LE.0 ) THEN
00085          RCOND = ONE
00086          RESID = ZERO
00087          RETURN
00088       END IF
00089 *
00090 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
00091 *
00092       EPS = SLAMCH( 'Epsilon' )
00093       ANORM = CLANSP( '1', UPLO, N, A, RWORK )
00094       AINVNM = CLANSP( '1', UPLO, N, AINV, RWORK )
00095       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00096          RCOND = ZERO
00097          RESID = ONE / EPS
00098          RETURN
00099       END IF
00100       RCOND = ( ONE/ANORM ) / AINVNM
00101 *
00102 *     Case where both A and AINV are upper triangular:
00103 *     Each element of - A * AINV is computed by taking the dot product
00104 *     of a row of A with a column of AINV.
00105 *
00106       IF( LSAME( UPLO, 'U' ) ) THEN
00107          DO 70 I = 1, N
00108             ICOL = ( ( I-1 )*I ) / 2 + 1
00109 *
00110 *           Code when J <= I
00111 *
00112             DO 30 J = 1, I
00113                JCOL = ( ( J-1 )*J ) / 2 + 1
00114                T = CDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 )
00115                JCOL = JCOL + 2*J - 1
00116                KCOL = ICOL - 1
00117                DO 10 K = J + 1, I
00118                   T = T + A( KCOL+K )*AINV( JCOL )
00119                   JCOL = JCOL + K
00120    10          CONTINUE
00121                KCOL = KCOL + 2*I
00122                DO 20 K = I + 1, N
00123                   T = T + A( KCOL )*AINV( JCOL )
00124                   KCOL = KCOL + K
00125                   JCOL = JCOL + K
00126    20          CONTINUE
00127                WORK( I, J ) = -T
00128    30       CONTINUE
00129 *
00130 *           Code when J > I
00131 *
00132             DO 60 J = I + 1, N
00133                JCOL = ( ( J-1 )*J ) / 2 + 1
00134                T = CDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 )
00135                JCOL = JCOL - 1
00136                KCOL = ICOL + 2*I - 1
00137                DO 40 K = I + 1, J
00138                   T = T + A( KCOL )*AINV( JCOL+K )
00139                   KCOL = KCOL + K
00140    40          CONTINUE
00141                JCOL = JCOL + 2*J
00142                DO 50 K = J + 1, N
00143                   T = T + A( KCOL )*AINV( JCOL )
00144                   KCOL = KCOL + K
00145                   JCOL = JCOL + K
00146    50          CONTINUE
00147                WORK( I, J ) = -T
00148    60       CONTINUE
00149    70    CONTINUE
00150       ELSE
00151 *
00152 *        Case where both A and AINV are lower triangular
00153 *
00154          NALL = ( N*( N+1 ) ) / 2
00155          DO 140 I = 1, N
00156 *
00157 *           Code when J <= I
00158 *
00159             ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1
00160             DO 100 J = 1, I
00161                JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I )
00162                T = CDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 )
00163                KCOL = I
00164                JCOL = J
00165                DO 80 K = 1, J - 1
00166                   T = T + A( KCOL )*AINV( JCOL )
00167                   JCOL = JCOL + N - K
00168                   KCOL = KCOL + N - K
00169    80          CONTINUE
00170                JCOL = JCOL - J
00171                DO 90 K = J, I - 1
00172                   T = T + A( KCOL )*AINV( JCOL+K )
00173                   KCOL = KCOL + N - K
00174    90          CONTINUE
00175                WORK( I, J ) = -T
00176   100       CONTINUE
00177 *
00178 *           Code when J > I
00179 *
00180             ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2
00181             DO 130 J = I + 1, N
00182                JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1
00183                T = CDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 )
00184                KCOL = I
00185                JCOL = J
00186                DO 110 K = 1, I - 1
00187                   T = T + A( KCOL )*AINV( JCOL )
00188                   JCOL = JCOL + N - K
00189                   KCOL = KCOL + N - K
00190   110          CONTINUE
00191                KCOL = KCOL - I
00192                DO 120 K = I, J - 1
00193                   T = T + A( KCOL+K )*AINV( JCOL )
00194                   JCOL = JCOL + N - K
00195   120          CONTINUE
00196                WORK( I, J ) = -T
00197   130       CONTINUE
00198   140    CONTINUE
00199       END IF
00200 *
00201 *     Add the identity matrix to WORK .
00202 *
00203       DO 150 I = 1, N
00204          WORK( I, I ) = WORK( I, I ) + ONE
00205   150 CONTINUE
00206 *
00207 *     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
00208 *
00209       RESID = CLANGE( '1', N, N, WORK, LDW, RWORK )
00210 *
00211       RESID = ( ( RESID*RCOND )/EPS ) / REAL( N )
00212 *
00213       RETURN
00214 *
00215 *     End of CSPT03
00216 *
00217       END