LAPACK: TESTING/EIG/zbdt01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
00002      $                   RWORK, 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       INTEGER            KD, LDA, LDPT, LDQ, M, N
00010       DOUBLE PRECISION   RESID
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
00014       COMPLEX*16         A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
00015      $                   WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZBDT01 reconstructs a general matrix A from its bidiagonal form
00022 *     A = Q * B * P'
00023 *  where Q (m by min(m,n)) and P' (min(m,n) by n) are unitary
00024 *  matrices and B is bidiagonal.
00025 *
00026 *  The test ratio to test the reduction is
00027 *     RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS )
00028 *  where PT = P' and EPS is the machine precision.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  M       (input) INTEGER
00034 *          The number of rows of the matrices A and Q.
00035 *
00036 *  N       (input) INTEGER
00037 *          The number of columns of the matrices A and P'.
00038 *
00039 *  KD      (input) INTEGER
00040 *          If KD = 0, B is diagonal and the array E is not referenced.
00041 *          If KD = 1, the reduction was performed by xGEBRD; B is upper
00042 *          bidiagonal if M >= N, and lower bidiagonal if M < N.
00043 *          If KD = -1, the reduction was performed by xGBBRD; B is
00044 *          always upper bidiagonal.
00045 *
00046 *  A       (input) COMPLEX*16 array, dimension (LDA,N)
00047 *          The m by n matrix A.
00048 *
00049 *  LDA     (input) INTEGER
00050 *          The leading dimension of the array A.  LDA >= max(1,M).
00051 *
00052 *  Q       (input) COMPLEX*16 array, dimension (LDQ,N)
00053 *          The m by min(m,n) unitary matrix Q in the reduction
00054 *          A = Q * B * P'.
00055 *
00056 *  LDQ     (input) INTEGER
00057 *          The leading dimension of the array Q.  LDQ >= max(1,M).
00058 *
00059 *  D       (input) DOUBLE PRECISION array, dimension (min(M,N))
00060 *          The diagonal elements of the bidiagonal matrix B.
00061 *
00062 *  E       (input) DOUBLE PRECISION array, dimension (min(M,N)-1)
00063 *          The superdiagonal elements of the bidiagonal matrix B if
00064 *          m >= n, or the subdiagonal elements of B if m < n.
00065 *
00066 *  PT      (input) COMPLEX*16 array, dimension (LDPT,N)
00067 *          The min(m,n) by n unitary matrix P' in the reduction
00068 *          A = Q * B * P'.
00069 *
00070 *  LDPT    (input) INTEGER
00071 *          The leading dimension of the array PT.
00072 *          LDPT >= max(1,min(M,N)).
00073 *
00074 *  WORK    (workspace) COMPLEX*16 array, dimension (M+N)
00075 *
00076 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (M)
00077 *
00078 *  RESID   (output) DOUBLE PRECISION
00079 *          The test ratio:  norm(A - Q * B * P') / ( n * norm(A) * EPS )
00080 *
00081 *  =====================================================================
00082 *
00083 *     .. Parameters ..
00084       DOUBLE PRECISION   ZERO, ONE
00085       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00086 *     ..
00087 *     .. Local Scalars ..
00088       INTEGER            I, J
00089       DOUBLE PRECISION   ANORM, EPS
00090 *     ..
00091 *     .. External Functions ..
00092       DOUBLE PRECISION   DLAMCH, DZASUM, ZLANGE
00093       EXTERNAL           DLAMCH, DZASUM, ZLANGE
00094 *     ..
00095 *     .. External Subroutines ..
00096       EXTERNAL           ZCOPY, ZGEMV
00097 *     ..
00098 *     .. Intrinsic Functions ..
00099       INTRINSIC          DBLE, DCMPLX, MAX, MIN
00100 *     ..
00101 *     .. Executable Statements ..
00102 *
00103 *     Quick return if possible
00104 *
00105       IF( M.LE.0 .OR. N.LE.0 ) THEN
00106          RESID = ZERO
00107          RETURN
00108       END IF
00109 *
00110 *     Compute A - Q * B * P' one column at a time.
00111 *
00112       RESID = ZERO
00113       IF( KD.NE.0 ) THEN
00114 *
00115 *        B is bidiagonal.
00116 *
00117          IF( KD.NE.0 .AND. M.GE.N ) THEN
00118 *
00119 *           B is upper bidiagonal and M >= N.
00120 *
00121             DO 20 J = 1, N
00122                CALL ZCOPY( M, A( 1, J ), 1, WORK, 1 )
00123                DO 10 I = 1, N - 1
00124                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
00125    10          CONTINUE
00126                WORK( M+N ) = D( N )*PT( N, J )
00127                CALL ZGEMV( 'No transpose', M, N, -DCMPLX( ONE ), Q, LDQ,
00128      $                     WORK( M+1 ), 1, DCMPLX( ONE ), WORK, 1 )
00129                RESID = MAX( RESID, DZASUM( M, WORK, 1 ) )
00130    20       CONTINUE
00131          ELSE IF( KD.LT.0 ) THEN
00132 *
00133 *           B is upper bidiagonal and M < N.
00134 *
00135             DO 40 J = 1, N
00136                CALL ZCOPY( M, A( 1, J ), 1, WORK, 1 )
00137                DO 30 I = 1, M - 1
00138                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
00139    30          CONTINUE
00140                WORK( M+M ) = D( M )*PT( M, J )
00141                CALL ZGEMV( 'No transpose', M, M, -DCMPLX( ONE ), Q, LDQ,
00142      $                     WORK( M+1 ), 1, DCMPLX( ONE ), WORK, 1 )
00143                RESID = MAX( RESID, DZASUM( M, WORK, 1 ) )
00144    40       CONTINUE
00145          ELSE
00146 *
00147 *           B is lower bidiagonal.
00148 *
00149             DO 60 J = 1, N
00150                CALL ZCOPY( M, A( 1, J ), 1, WORK, 1 )
00151                WORK( M+1 ) = D( 1 )*PT( 1, J )
00152                DO 50 I = 2, M
00153                   WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
00154      $                          D( I )*PT( I, J )
00155    50          CONTINUE
00156                CALL ZGEMV( 'No transpose', M, M, -DCMPLX( ONE ), Q, LDQ,
00157      $                     WORK( M+1 ), 1, DCMPLX( ONE ), WORK, 1 )
00158                RESID = MAX( RESID, DZASUM( M, WORK, 1 ) )
00159    60       CONTINUE
00160          END IF
00161       ELSE
00162 *
00163 *        B is diagonal.
00164 *
00165          IF( M.GE.N ) THEN
00166             DO 80 J = 1, N
00167                CALL ZCOPY( M, A( 1, J ), 1, WORK, 1 )
00168                DO 70 I = 1, N
00169                   WORK( M+I ) = D( I )*PT( I, J )
00170    70          CONTINUE
00171                CALL ZGEMV( 'No transpose', M, N, -DCMPLX( ONE ), Q, LDQ,
00172      $                     WORK( M+1 ), 1, DCMPLX( ONE ), WORK, 1 )
00173                RESID = MAX( RESID, DZASUM( M, WORK, 1 ) )
00174    80       CONTINUE
00175          ELSE
00176             DO 100 J = 1, N
00177                CALL ZCOPY( M, A( 1, J ), 1, WORK, 1 )
00178                DO 90 I = 1, M
00179                   WORK( M+I ) = D( I )*PT( I, J )
00180    90          CONTINUE
00181                CALL ZGEMV( 'No transpose', M, M, -DCMPLX( ONE ), Q, LDQ,
00182      $                     WORK( M+1 ), 1, DCMPLX( ONE ), WORK, 1 )
00183                RESID = MAX( RESID, DZASUM( M, WORK, 1 ) )
00184   100       CONTINUE
00185          END IF
00186       END IF
00187 *
00188 *     Compute norm(A - Q * B * P') / ( n * norm(A) * EPS )
00189 *
00190       ANORM = ZLANGE( '1', M, N, A, LDA, RWORK )
00191       EPS = DLAMCH( 'Precision' )
00192 *
00193       IF( ANORM.LE.ZERO ) THEN
00194          IF( RESID.NE.ZERO )
00195      $      RESID = ONE / EPS
00196       ELSE
00197          IF( ANORM.GE.RESID ) THEN
00198             RESID = ( RESID / ANORM ) / ( DBLE( N )*EPS )
00199          ELSE
00200             IF( ANORM.LT.ONE ) THEN
00201                RESID = ( MIN( RESID, DBLE( N )*ANORM ) / ANORM ) /
00202      $                 ( DBLE( N )*EPS )
00203             ELSE
00204                RESID = MIN( RESID / ANORM, DBLE( N ) ) /
00205      $                 ( DBLE( N )*EPS )
00206             END IF
00207          END IF
00208       END IF
00209 *
00210       RETURN
00211 *
00212 *     End of ZBDT01
00213 *
00214       END