ddrvrf4.f

Go to the documentation of this file.

*> \brief \b DDRVRF4
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DDRVRF4( NOUT, NN, NVAL, THRESH, C1, C2, LDC, CRF, A,
*      +                    LDA, D_WORK_DLANGE )
* 
*       .. Scalar Arguments ..
*       INTEGER            LDA, LDC, NN, NOUT
*       DOUBLE PRECISION   THRESH
*       ..
*       .. Array Arguments ..
*       INTEGER            NVAL( NN )
*       DOUBLE PRECISION   A( LDA, * ), C1( LDC, * ), C2( LDC, *),
*      +                   CRF( * ), D_WORK_DLANGE( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DDRVRF4 tests the LAPACK RFP routines:
*>     DSFRK
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] NOUT
*> \verbatim
*>          NOUT is INTEGER
*>                The unit number for output.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*>          NN is INTEGER
*>                The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*>          NVAL is INTEGER array, dimension (NN)
*>                The values of the matrix dimension N.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*>          THRESH is DOUBLE PRECISION
*>                The threshold value for the test ratios.  A result is
*>                included in the output file if RESULT >= THRESH.  To
*>                have every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[out] C1
*> \verbatim
*>          C1 is DOUBLE PRECISION array,
*>                dimension (LDC,NMAX)
*> \endverbatim
*>
*> \param[out] C2
*> \verbatim
*>          C2 is DOUBLE PRECISION array,
*>                dimension (LDC,NMAX)
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>                The leading dimension of the array A.
*>                LDA >= max(1,NMAX).
*> \endverbatim
*>
*> \param[out] CRF
*> \verbatim
*>          CRF is DOUBLE PRECISION array,
*>                dimension ((NMAX*(NMAX+1))/2).
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*>          A is DOUBLE PRECISION array,
*>                dimension (LDA,NMAX)
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>                The leading dimension of the array A.  LDA >= max(1,NMAX).
*> \endverbatim
*>
*> \param[out] D_WORK_DLANGE
*> \verbatim
*>          D_WORK_DLANGE is DOUBLE PRECISION array, dimension (NMAX)
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup double_lin
*
*  =====================================================================
      SUBROUTINE ddrvrf4( NOUT, NN, NVAL, THRESH, C1, C2, LDC, CRF, A,
     +                    lda, d_work_dlange )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            LDA, LDC, NN, NOUT
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            NVAL( nn )
      DOUBLE PRECISION   A( lda, * ), C1( ldc, * ), C2( ldc, *),
     +                   crf( * ), d_work_dlange( * )
*     ..
*
*  =====================================================================
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter                ( zero = 0.0d+0, one  = 1.0d+0 )
      INTEGER            NTESTS
      parameter                ( ntests = 1 )
*     ..
*     .. Local Scalars ..
      CHARACTER          UPLO, CFORM, TRANS
      INTEGER            I, IFORM, IIK, IIN, INFO, IUPLO, J, K, N,
     +                   nfail, nrun, ialpha, itrans
      DOUBLE PRECISION   ALPHA, BETA, EPS, NORMA, NORMC
*     ..
*     .. Local Arrays ..
      CHARACTER          UPLOS( 2 ), FORMS( 2 ), TRANSS( 2 )
      INTEGER            ISEED( 4 ), ISEEDY( 4 )
      DOUBLE PRECISION   RESULT( ntests )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLARND, DLANGE
      EXTERNAL           dlamch, dlarnd, dlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           dsyrk, dsfrk, dtfttr, dtrttf
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Scalars in Common ..
      CHARACTER*32       SRNAMT
*     ..
*     .. Common blocks ..
      COMMON             / srnamc / srnamt
*     ..
*     .. Data statements ..
      DATA               iseedy / 1988, 1989, 1990, 1991 /
      DATA               uplos  / 'U', 'L' /
      DATA               forms  / 'N', 'T' /
      DATA               transs / 'N', 'T' /
*     ..
*     .. Executable Statements ..
*
*     Initialize constants and the random number seed.
*
      nrun = 0
      nfail = 0
      info = 0
      DO 10 i = 1, 4
         iseed( i ) = iseedy( i )
   10 CONTINUE
      eps = dlamch( 'Precision' )
*
      DO 150 iin = 1, nn
*
         n = nval( iin )
*
         DO 140 iik = 1, nn
*
            k = nval( iin )
*
            DO 130 iform = 1, 2
*
               cform = forms( iform )
*
               DO 120 iuplo = 1, 2
*
                  uplo = uplos( iuplo )
*
                  DO 110 itrans = 1, 2
*
                     trans = transs( itrans )
*
                     DO 100 ialpha = 1, 4
*
                        IF ( ialpha.EQ. 1) THEN
                           alpha = zero
                           beta = zero
                        ELSE IF ( ialpha.EQ. 2) THEN
                           alpha = one
                           beta = zero
                        ELSE IF ( ialpha.EQ. 3) THEN
                           alpha = zero
                           beta = one
                        ELSE
                           alpha = dlarnd( 2, iseed )
                           beta = dlarnd( 2, iseed )
                        END IF
*
*                       All the parameters are set:
*                          CFORM, UPLO, TRANS, M, N,
*                          ALPHA, and BETA
*                       READY TO TEST!
*
                        nrun = nrun + 1
*
                        IF ( itrans.EQ.1 ) THEN
*
*                          In this case we are NOTRANS, so A is N-by-K
*
                           DO j = 1, k
                              DO i = 1, n
                                 a( i, j) = dlarnd( 2, iseed )
                              END DO
                           END DO
*
                           norma = dlange( 'I', n, k, a, lda,
     +                                      d_work_dlange )
*
 
                        ELSE
*
*                          In this case we are TRANS, so A is K-by-N
*
                           DO j = 1,n 
                              DO i = 1, k
                                 a( i, j) = dlarnd( 2, iseed )
                              END DO
                           END DO
*
                           norma = dlange( 'I', k, n, a, lda,
     +                                      d_work_dlange )
*
                        END IF
*
*                       Generate C1 our N--by--N symmetric matrix. 
*                       Make sure C2 has the same upper/lower part,
*                       (the one that we do not touch), so
*                       copy the initial C1 in C2 in it.
*
                        DO j = 1, n
                           DO i = 1, n
                              c1( i, j) = dlarnd( 2, iseed )
                              c2(i,j) = c1(i,j)
                           END DO
                        END DO
*
*                       (See comment later on for why we use DLANGE and
*                       not DLANSY for C1.)
*
                        normc = dlange( 'I', n, n, c1, ldc,
     +                                      d_work_dlange )
*
                        srnamt = 'DTRTTF'
                        CALL dtrttf( cform, uplo, n, c1, ldc, crf,
     +                               info )
*
*                       call dsyrk the BLAS routine -> gives C1
*
                        srnamt = 'DSYRK '
                        CALL dsyrk( uplo, trans, n, k, alpha, a, lda,
     +                              beta, c1, ldc )
*
*                       call dsfrk the RFP routine -> gives CRF
*
                        srnamt = 'DSFRK '
                        CALL dsfrk( cform, uplo, trans, n, k, alpha, a,
     +                              lda, beta, crf )
*
*                       convert CRF in full format -> gives C2
*
                        srnamt = 'DTFTTR'
                        CALL dtfttr( cform, uplo, n, crf, c2, ldc,
     +                               info )
*
*                       compare C1 and C2
*
                        DO j = 1, n
                           DO i = 1, n
                              c1(i,j) = c1(i,j)-c2(i,j)
                           END DO
                        END DO
*
*                       Yes, C1 is symmetric so we could call DLANSY,
*                       but we want to check the upper part that is
*                       supposed to be unchanged and the diagonal that
*                       is supposed to be real -> DLANGE
*
                        result(1) = dlange( 'I', n, n, c1, ldc,
     +                                      d_work_dlange )
                        result(1) = result(1) 
     +                              / max( abs( alpha ) * norma
     +                                   + abs( beta ) , one )
     +                              / max( n , 1 ) / eps
*
                        IF( result(1).GE.thresh ) THEN
                           IF( nfail.EQ.0 ) THEN
                              WRITE( nout, * )
                              WRITE( nout, fmt = 9999 )
                           END IF
                           WRITE( nout, fmt = 9997 ) 'DSFRK', 
     +                        cform, uplo, trans, n, k, result(1)
                           nfail = nfail + 1
                        END IF
*
  100                CONTINUE
  110             CONTINUE
  120          CONTINUE
  130       CONTINUE
  140    CONTINUE
  150 CONTINUE
*
*     Print a summary of the results.
*
      IF ( nfail.EQ.0 ) THEN
         WRITE( nout, fmt = 9996 ) 'DSFRK', nrun
      ELSE
         WRITE( nout, fmt = 9995 ) 'DSFRK', nfail, nrun
      END IF
*
 9999 FORMAT( 1x, ' *** Error(s) or Failure(s) while testing DSFRK 
     +         ***')
 9997 FORMAT( 1x, '     Failure in ',a5,', CFORM=''',a1,''',',
     + ' UPLO=''',a1,''',',' TRANS=''',a1,''',', ' N=',i3,', K =', i3,
     + ', test=',g12.5)
 9996 FORMAT( 1x, 'All tests for ',a5,' auxiliary routine passed the ',
     +        'threshold ( ',i5,' tests run)')
 9995 FORMAT( 1x, a6, ' auxiliary routine: ',i5,' out of ',i5,
     +        ' tests failed to pass the threshold')
*
      RETURN
*
*     End of DDRVRF4
*
      END