LAPACK: TESTING/EIG/dbdt01.f Source File

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