LAPACK: TESTING/EIG/sckcsd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
00002      $                   MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
00003      $                   WORK, RWORK, NIN, NOUT, INFO )
00004       IMPLICIT NONE
00005 *
00006 *     Originally SCKGSV
00007 *  -- LAPACK test routine (version 3.1) --
00008 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00009 *     November 2006
00010 *
00011 *     Adapted to SCKCSD
00012 *     July 2010
00013 *
00014 *     .. Scalar Arguments ..
00015       INTEGER            INFO, NIN, NM, NMATS, MMAX, NOUT
00016       REAL               THRESH
00017 *     ..
00018 *     .. Array Arguments ..
00019       INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
00020      $                   QVAL( * )
00021       REAL               RWORK( * ), THETA( * )
00022       REAL               U1( * ), U2( * ), V1T( * ), V2T( * ),
00023      $                   WORK( * ), X( * ), XF( * )
00024 *     ..
00025 *
00026 *  Purpose
00027 *  =======
00028 *
00029 *  SCKCSD tests SORCSD:
00030 *         the CSD for an M-by-M orthogonal matrix X partitioned as
00031 *         [ X11 X12; X21 X22 ]. X11 is P-by-Q.
00032 *
00033 *  Arguments
00034 *  =========
00035 *
00036 *  NM      (input) INTEGER
00037 *          The number of values of M contained in the vector MVAL.
00038 *
00039 *  MVAL    (input) INTEGER array, dimension (NM)
00040 *          The values of the matrix row dimension M.
00041 *
00042 *  PVAL    (input) INTEGER array, dimension (NM)
00043 *          The values of the matrix row dimension P.
00044 *
00045 *  QVAL    (input) INTEGER array, dimension (NM)
00046 *          The values of the matrix column dimension Q.
00047 *
00048 *  NMATS   (input) INTEGER
00049 *          The number of matrix types to be tested for each combination
00050 *          of matrix dimensions.  If NMATS >= NTYPES (the maximum
00051 *          number of matrix types), then all the different types are
00052 *          generated for testing.  If NMATS < NTYPES, another input line
00053 *          is read to get the numbers of the matrix types to be used.
00054 *
00055 *  ISEED   (input/output) INTEGER array, dimension (4)
00056 *          On entry, the seed of the random number generator.  The array
00057 *          elements should be between 0 and 4095, otherwise they will be
00058 *          reduced mod 4096, and ISEED(4) must be odd.
00059 *          On exit, the next seed in the random number sequence after
00060 *          all the test matrices have been generated.
00061 *
00062 *  THRESH  (input) REAL
00063 *          The threshold value for the test ratios.  A result is
00064 *          included in the output file if RESULT >= THRESH.  To have
00065 *          every test ratio printed, use THRESH = 0.
00066 *
00067 *  MMAX    (input) INTEGER
00068 *          The maximum value permitted for M, used in dimensioning the
00069 *          work arrays.
00070 *
00071 *  X       (workspace) REAL array, dimension (MMAX*MMAX)
00072 *
00073 *  XF      (workspace) REAL array, dimension (MMAX*MMAX)
00074 *
00075 *  U1      (workspace) REAL array, dimension (MMAX*MMAX)
00076 *
00077 *  U2      (workspace) REAL array, dimension (MMAX*MMAX)
00078 *
00079 *  V1T     (workspace) REAL array, dimension (MMAX*MMAX)
00080 *
00081 *  V2T     (workspace) REAL array, dimension (MMAX*MMAX)
00082 *
00083 *  THETA   (workspace) REAL array, dimension (MMAX)
00084 *
00085 *  IWORK   (workspace) INTEGER array, dimension (MMAX)
00086 *
00087 *  WORK    (workspace) REAL array
00088 *
00089 *  RWORK   (workspace) REAL array
00090 *
00091 *  NIN     (input) INTEGER
00092 *          The unit number for input.
00093 *
00094 *  NOUT    (input) INTEGER
00095 *          The unit number for output.
00096 *
00097 *  INFO    (output) INTEGER
00098 *          = 0 :  successful exit
00099 *          > 0 :  If SLAROR returns an error code, the absolute value
00100 *                 of it is returned.
00101 *
00102 *  =====================================================================
00103 *
00104 *     .. Parameters ..
00105       INTEGER            NTESTS
00106       PARAMETER          ( NTESTS = 9 )
00107       INTEGER            NTYPES
00108       PARAMETER          ( NTYPES = 3 )
00109       REAL               GAPDIGIT, ORTH, PIOVER2, TEN
00110       PARAMETER          ( GAPDIGIT = 10.0E0, ORTH = 1.0E-4,
00111      $                     PIOVER2 = 1.57079632679489662E0,
00112      $                     TEN = 10.0D0 )
00113       REAL               ZERO, ONE
00114       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00115 *     ..
00116 *     .. Local Scalars ..
00117       LOGICAL            FIRSTT
00118       CHARACTER*3        PATH
00119       INTEGER            I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
00120      $                   LDV2T, LDX, LWORK, M, NFAIL, NRUN, NT, P, Q, R
00121 *     ..
00122 *     .. Local Arrays ..
00123       LOGICAL            DOTYPE( NTYPES )
00124       REAL               RESULT( NTESTS )
00125 *     ..
00126 *     .. External Subroutines ..
00127       EXTERNAL           ALAHDG, ALAREQ, ALASUM, SCSDTS, SLACSG, SLAROR,
00128      $                   SLASET
00129 *     ..
00130 *     .. Intrinsic Functions ..
00131       INTRINSIC          ABS, COS, MIN, SIN
00132 *     ..
00133 *     .. External Functions ..
00134       REAL               SLANGE, SLARND
00135       EXTERNAL           SLANGE, SLARND
00136 *     ..
00137 *     .. Executable Statements ..
00138 *
00139 *     Initialize constants and the random number seed.
00140 *
00141       PATH( 1: 3 ) = 'CSD'
00142       INFO = 0
00143       NRUN = 0
00144       NFAIL = 0
00145       FIRSTT = .TRUE.
00146       CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
00147       LDX = MMAX
00148       LDU1 = MMAX
00149       LDU2 = MMAX
00150       LDV1T = MMAX
00151       LDV2T = MMAX
00152       LWORK = MMAX*MMAX
00153 *
00154 *     Do for each value of M in MVAL.
00155 *
00156       DO 30 IM = 1, NM
00157          M = MVAL( IM )
00158          P = PVAL( IM )
00159          Q = QVAL( IM )
00160 *
00161          DO 20 IMAT = 1, NTYPES
00162 *
00163 *           Do the tests only if DOTYPE( IMAT ) is true.
00164 *
00165             IF( .NOT.DOTYPE( IMAT ) )
00166      $         GO TO 20
00167 *
00168 *           Generate X
00169 *
00170             IF( IMAT.EQ.1 ) THEN
00171                CALL SLAROR( 'L', 'I', M, M, X, LDX, ISEED, WORK, IINFO )
00172                IF( M .NE. 0 .AND. IINFO .NE. 0 ) THEN
00173                   WRITE( NOUT, FMT = 9999 ) M, IINFO
00174                   INFO = ABS( IINFO )
00175                   GO TO 20
00176                END IF
00177             ELSE IF( IMAT.EQ.2 ) THEN
00178                R = MIN( P, M-P, Q, M-Q )
00179                DO I = 1, R
00180                   THETA(I) = PIOVER2 * SLARND( 1, ISEED )
00181                END DO
00182                CALL SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
00183                DO I = 1, M
00184                   DO J = 1, M
00185                      X(I+(J-1)*LDX) = X(I+(J-1)*LDX) +
00186      $                                ORTH*SLARND(2,ISEED)
00187                   END DO
00188                END DO
00189             ELSE
00190                R = MIN( P, M-P, Q, M-Q )
00191                DO I = 1, R+1
00192                   THETA(I) = TEN**(-SLARND(1,ISEED)*GAPDIGIT)
00193                END DO
00194                DO I = 2, R+1
00195                   THETA(I) = THETA(I-1) + THETA(I)
00196                END DO
00197                DO I = 1, R
00198                   THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
00199                END DO
00200                CALL SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
00201             END IF
00202 *
00203             NT = 9
00204 *
00205             CALL SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
00206      $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
00207      $                   RWORK, RESULT )
00208 *
00209 *           Print information about the tests that did not
00210 *           pass the threshold.
00211 *
00212             DO 10 I = 1, NT
00213                IF( RESULT( I ).GE.THRESH ) THEN
00214                   IF( NFAIL.EQ.0 .AND. FIRSTT ) THEN
00215                      FIRSTT = .FALSE.
00216                      CALL ALAHDG( NOUT, PATH )
00217                   END IF
00218                   WRITE( NOUT, FMT = 9998 )M, P, Q, IMAT, I,
00219      $               RESULT( I )
00220                   NFAIL = NFAIL + 1
00221                END IF
00222    10       CONTINUE
00223             NRUN = NRUN + NT
00224    20    CONTINUE
00225    30 CONTINUE
00226 *
00227 *     Print a summary of the results.
00228 *
00229       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, 0 )
00230 *
00231  9999 FORMAT( ' SLAROR in SCKCSD: M = ', I5, ', INFO = ', I15 )
00232  9998 FORMAT( ' M=', I4, ' P=', I4, ', Q=', I4, ', type ', I2,
00233      $      ', test ', I2, ', ratio=', G13.6 )
00234       RETURN
00235 *
00236 *     End of SCKCSD
00237 *
00238       END
00239 *
00240 *
00241 *
00242       SUBROUTINE SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
00243       IMPLICIT NONE
00244 *
00245       INTEGER            LDX, M, P, Q
00246       INTEGER            ISEED( 4 )
00247       REAL               THETA( * )
00248       REAL               WORK( * ), X( LDX, * )
00249 *
00250       REAL               ONE, ZERO
00251       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
00252 *
00253       INTEGER            I, INFO, R
00254 *
00255       R = MIN( P, M-P, Q, M-Q )
00256 *
00257       CALL SLASET( 'Full', M, M, ZERO, ZERO, X, LDX )
00258 *
00259       DO I = 1, MIN(P,Q)-R
00260          X(I,I) = ONE
00261       END DO
00262       DO I = 1, R
00263          X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = COS(THETA(I))
00264       END DO
00265       DO I = 1, MIN(P,M-Q)-R
00266          X(P-I+1,M-I+1) = -ONE
00267       END DO
00268       DO I = 1, R
00269          X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
00270      $      -SIN(THETA(R-I+1))
00271       END DO
00272       DO I = 1, MIN(M-P,Q)-R
00273          X(M-I+1,Q-I+1) = ONE
00274       END DO
00275       DO I = 1, R
00276          X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
00277      $      SIN(THETA(R-I+1))
00278       END DO
00279       DO I = 1, MIN(M-P,M-Q)-R
00280          X(P+I,Q+I) = ONE
00281       END DO
00282       DO I = 1, R
00283          X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
00284      $      COS(THETA(I))
00285       END DO
00286       CALL SLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO )
00287       CALL SLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX,
00288      $             ISEED, WORK, INFO )
00289       CALL SLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED,
00290      $             WORK, INFO )
00291       CALL SLAROR( 'Right', 'No init', M, M-Q,
00292      $             X(1,Q+1), LDX, ISEED, WORK, INFO )
00293 *
00294       END
00295