LAPACK: TESTING/LIN/dqrt14.f Source File

Go to the documentation of this file.
00001       DOUBLE PRECISION FUNCTION DQRT14( TRANS, M, N, NRHS, A, LDA, X,
00002      $                 LDX, WORK, LWORK )
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          TRANS
00010       INTEGER            LDA, LDX, LWORK, M, N, NRHS
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), WORK( LWORK ), X( LDX, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DQRT14 checks whether X is in the row space of A or A'.  It does so
00020 *  by scaling both X and A such that their norms are in the range
00021 *  [sqrt(eps), 1/sqrt(eps)], then computing a QR factorization of [A,X]
00022 *  (if TRANS = 'T') or an LQ factorization of [A',X]' (if TRANS = 'N'),
00023 *  and returning the norm of the trailing triangle, scaled by
00024 *  MAX(M,N,NRHS)*eps.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  TRANS   (input) CHARACTER*1
00030 *          = 'N':  No transpose, check for X in the row space of A
00031 *          = 'T':  Transpose, check for X in the row space of A'.
00032 *
00033 *  M       (input) INTEGER
00034 *          The number of rows of the matrix A.
00035 *
00036 *  N       (input) INTEGER
00037 *          The number of columns of the matrix A.
00038 *
00039 *  NRHS    (input) INTEGER
00040 *          The number of right hand sides, i.e., the number of columns
00041 *          of X.
00042 *
00043 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00044 *          The M-by-N matrix A.
00045 *
00046 *  LDA     (input) INTEGER
00047 *          The leading dimension of the array A.
00048 *
00049 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
00050 *          If TRANS = 'N', the N-by-NRHS matrix X.
00051 *          IF TRANS = 'T', the M-by-NRHS matrix X.
00052 *
00053 *  LDX     (input) INTEGER
00054 *          The leading dimension of the array X.
00055 *
00056 *  WORK    (workspace) DOUBLE PRECISION array dimension (LWORK)
00057 *
00058 *  LWORK   (input) INTEGER
00059 *          length of workspace array required
00060 *          If TRANS = 'N', LWORK >= (M+NRHS)*(N+2);
00061 *          if TRANS = 'T', LWORK >= (N+NRHS)*(M+2).
00062 *
00063 *  =====================================================================
00064 *
00065 *     .. Parameters ..
00066       DOUBLE PRECISION   ZERO, ONE
00067       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00068 *     ..
00069 *     .. Local Scalars ..
00070       LOGICAL            TPSD
00071       INTEGER            I, INFO, J, LDWORK
00072       DOUBLE PRECISION   ANRM, ERR, XNRM
00073 *     ..
00074 *     .. Local Arrays ..
00075       DOUBLE PRECISION   RWORK( 1 )
00076 *     ..
00077 *     .. External Functions ..
00078       LOGICAL            LSAME
00079       DOUBLE PRECISION   DLAMCH, DLANGE
00080       EXTERNAL           LSAME, DLAMCH, DLANGE
00081 *     ..
00082 *     .. External Subroutines ..
00083       EXTERNAL           DGELQ2, DGEQR2, DLACPY, DLASCL, XERBLA
00084 *     ..
00085 *     .. Intrinsic Functions ..
00086       INTRINSIC          ABS, DBLE, MAX, MIN
00087 *     ..
00088 *     .. Executable Statements ..
00089 *
00090       DQRT14 = ZERO
00091       IF( LSAME( TRANS, 'N' ) ) THEN
00092          LDWORK = M + NRHS
00093          TPSD = .FALSE.
00094          IF( LWORK.LT.( M+NRHS )*( N+2 ) ) THEN
00095             CALL XERBLA( 'DQRT14', 10 )
00096             RETURN
00097          ELSE IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00098             RETURN
00099          END IF
00100       ELSE IF( LSAME( TRANS, 'T' ) ) THEN
00101          LDWORK = M
00102          TPSD = .TRUE.
00103          IF( LWORK.LT.( N+NRHS )*( M+2 ) ) THEN
00104             CALL XERBLA( 'DQRT14', 10 )
00105             RETURN
00106          ELSE IF( M.LE.0 .OR. NRHS.LE.0 ) THEN
00107             RETURN
00108          END IF
00109       ELSE
00110          CALL XERBLA( 'DQRT14', 1 )
00111          RETURN
00112       END IF
00113 *
00114 *     Copy and scale A
00115 *
00116       CALL DLACPY( 'All', M, N, A, LDA, WORK, LDWORK )
00117       ANRM = DLANGE( 'M', M, N, WORK, LDWORK, RWORK )
00118       IF( ANRM.NE.ZERO )
00119      $   CALL DLASCL( 'G', 0, 0, ANRM, ONE, M, N, WORK, LDWORK, INFO )
00120 *
00121 *     Copy X or X' into the right place and scale it
00122 *
00123       IF( TPSD ) THEN
00124 *
00125 *        Copy X into columns n+1:n+nrhs of work
00126 *
00127          CALL DLACPY( 'All', M, NRHS, X, LDX, WORK( N*LDWORK+1 ),
00128      $                LDWORK )
00129          XNRM = DLANGE( 'M', M, NRHS, WORK( N*LDWORK+1 ), LDWORK,
00130      $          RWORK )
00131          IF( XNRM.NE.ZERO )
00132      $      CALL DLASCL( 'G', 0, 0, XNRM, ONE, M, NRHS,
00133      $                   WORK( N*LDWORK+1 ), LDWORK, INFO )
00134          ANRM = DLANGE( 'One-norm', M, N+NRHS, WORK, LDWORK, RWORK )
00135 *
00136 *        Compute QR factorization of X
00137 *
00138          CALL DGEQR2( M, N+NRHS, WORK, LDWORK,
00139      $                WORK( LDWORK*( N+NRHS )+1 ),
00140      $                WORK( LDWORK*( N+NRHS )+MIN( M, N+NRHS )+1 ),
00141      $                INFO )
00142 *
00143 *        Compute largest entry in upper triangle of
00144 *        work(n+1:m,n+1:n+nrhs)
00145 *
00146          ERR = ZERO
00147          DO 20 J = N + 1, N + NRHS
00148             DO 10 I = N + 1, MIN( M, J )
00149                ERR = MAX( ERR, ABS( WORK( I+( J-1 )*M ) ) )
00150    10       CONTINUE
00151    20    CONTINUE
00152 *
00153       ELSE
00154 *
00155 *        Copy X' into rows m+1:m+nrhs of work
00156 *
00157          DO 40 I = 1, N
00158             DO 30 J = 1, NRHS
00159                WORK( M+J+( I-1 )*LDWORK ) = X( I, J )
00160    30       CONTINUE
00161    40    CONTINUE
00162 *
00163          XNRM = DLANGE( 'M', NRHS, N, WORK( M+1 ), LDWORK, RWORK )
00164          IF( XNRM.NE.ZERO )
00165      $      CALL DLASCL( 'G', 0, 0, XNRM, ONE, NRHS, N, WORK( M+1 ),
00166      $                   LDWORK, INFO )
00167 *
00168 *        Compute LQ factorization of work
00169 *
00170          CALL DGELQ2( LDWORK, N, WORK, LDWORK, WORK( LDWORK*N+1 ),
00171      $                WORK( LDWORK*( N+1 )+1 ), INFO )
00172 *
00173 *        Compute largest entry in lower triangle in
00174 *        work(m+1:m+nrhs,m+1:n)
00175 *
00176          ERR = ZERO
00177          DO 60 J = M + 1, N
00178             DO 50 I = J, LDWORK
00179                ERR = MAX( ERR, ABS( WORK( I+( J-1 )*LDWORK ) ) )
00180    50       CONTINUE
00181    60    CONTINUE
00182 *
00183       END IF
00184 *
00185       DQRT14 = ERR / ( DBLE( MAX( M, N, NRHS ) )*DLAMCH( 'Epsilon' ) )
00186 *
00187       RETURN
00188 *
00189 *     End of DQRT14
00190 *
00191       END