SUBROUTINE SCHKBB( NSIZES, MVAL, NVAL, NWDTHS, KK, NTYPES, DOTYPE, $ NRHS, ISEED, THRESH, NOUNIT, A, LDA, AB, LDAB, $ BD, BE, Q, LDQ, P, LDP, C, LDC, CC, WORK, $ LWORK, RESULT, INFO ) * * -- LAPACK test routine (release 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDAB, LDC, LDP, LDQ, LWORK, NOUNIT, $ NRHS, NSIZES, NTYPES, NWDTHS REAL THRESH * .. * .. Array Arguments .. LOGICAL DOTYPE( * ) INTEGER ISEED( 4 ), KK( * ), MVAL( * ), NVAL( * ) REAL A( LDA, * ), AB( LDAB, * ), BD( * ), BE( * ), $ C( LDC, * ), CC( LDC, * ), P( LDP, * ), $ Q( LDQ, * ), RESULT( * ), WORK( * ) * .. * * Purpose * ======= * * SCHKBB tests the reduction of a general real rectangular band * matrix to bidiagonal form. * * SGBBRD factors a general band matrix A as Q B P* , where * means * transpose, B is upper bidiagonal, and Q and P are orthogonal; * SGBBRD can also overwrite a given matrix C with Q* C . * * For each pair of matrix dimensions (M,N) and each selected matrix * type, an M by N matrix A and an M by NRHS matrix C are generated. * The problem dimensions are as follows * A: M x N * Q: M x M * P: N x N * B: min(M,N) x min(M,N) * C: M x NRHS * * For each generated matrix, 4 tests are performed: * * (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P' * * (2) | I - Q' Q | / ( M ulp ) * * (3) | I - PT PT' | / ( N ulp ) * * (4) | Y - Q' C | / ( |Y| max(M,NRHS) ulp ), where Y = Q' C. * * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. * Currently, the list of possible types is: * * The possible matrix types are * * (1) The zero matrix. * (2) The identity matrix. * * (3) A diagonal matrix with evenly spaced entries * 1, ..., ULP and random signs. * (ULP = (first number larger than 1) - 1 ) * (4) A diagonal matrix with geometrically spaced entries * 1, ..., ULP and random signs. * (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP * and random signs. * * (6) Same as (3), but multiplied by SQRT( overflow threshold ) * (7) Same as (3), but multiplied by SQRT( underflow threshold ) * * (8) A matrix of the form U D V, where U and V are orthogonal and * D has evenly spaced entries 1, ..., ULP with random signs * on the diagonal. * * (9) A matrix of the form U D V, where U and V are orthogonal and * D has geometrically spaced entries 1, ..., ULP with random * signs on the diagonal. * * (10) A matrix of the form U D V, where U and V are orthogonal and * D has "clustered" entries 1, ULP,..., ULP with random * signs on the diagonal. * * (11) Same as (8), but multiplied by SQRT( overflow threshold ) * (12) Same as (8), but multiplied by SQRT( underflow threshold ) * * (13) Rectangular matrix with random entries chosen from (-1,1). * (14) Same as (13), but multiplied by SQRT( overflow threshold ) * (15) Same as (13), but multiplied by SQRT( underflow threshold ) * * Arguments * ========= * * NSIZES (input) INTEGER * The number of values of M and N contained in the vectors * MVAL and NVAL. The matrix sizes are used in pairs (M,N). * If NSIZES is zero, SCHKBB does nothing. NSIZES must be at * least zero. * * MVAL (input) INTEGER array, dimension (NSIZES) * The values of the matrix row dimension M. * * NVAL (input) INTEGER array, dimension (NSIZES) * The values of the matrix column dimension N. * * NWDTHS (input) INTEGER * The number of bandwidths to use. If it is zero, * SCHKBB does nothing. It must be at least zero. * * KK (input) INTEGER array, dimension (NWDTHS) * An array containing the bandwidths to be used for the band * matrices. The values must be at least zero. * * NTYPES (input) INTEGER * The number of elements in DOTYPE. If it is zero, SCHKBB * does nothing. It must be at least zero. If it is MAXTYP+1 * and NSIZES is 1, then an additional type, MAXTYP+1 is * defined, which is to use whatever matrix is in A. This * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and * DOTYPE(MAXTYP+1) is .TRUE. . * * DOTYPE (input) LOGICAL array, dimension (NTYPES) * If DOTYPE(j) is .TRUE., then for each size in NN a * matrix of that size and of type j will be generated. * If NTYPES is smaller than the maximum number of types * defined (PARAMETER MAXTYP), then types NTYPES+1 through * MAXTYP will not be generated. If NTYPES is larger * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) * will be ignored. * * NRHS (input) INTEGER * The number of columns in the "right-hand side" matrix C. * If NRHS = 0, then the operations on the right-hand side will * not be tested. NRHS must be at least 0. * * ISEED (input/output) INTEGER array, dimension (4) * On entry ISEED specifies the seed of the random number * generator. The array elements should be between 0 and 4095; * if not they will be reduced mod 4096. Also, ISEED(4) must * be odd. The random number generator uses a linear * congruential sequence limited to small integers, and so * should produce machine independent random numbers. The * values of ISEED are changed on exit, and can be used in the * next call to SCHKBB to continue the same random number * sequence. * * THRESH (input) REAL * A test will count as "failed" if the "error", computed as * described above, exceeds THRESH. Note that the error * is scaled to be O(1), so THRESH should be a reasonably * small multiple of 1, e.g., 10 or 100. In particular, * it should not depend on the precision (single vs. double) * or the size of the matrix. It must be at least zero. * * NOUNIT (input) INTEGER * The FORTRAN unit number for printing out error messages * (e.g., if a routine returns IINFO not equal to 0.) * * A (input/workspace) REAL array, dimension * (LDA, max(NN)) * Used to hold the matrix A. * * LDA (input) INTEGER * The leading dimension of A. It must be at least 1 * and at least max( NN ). * * AB (workspace) REAL array, dimension (LDAB, max(NN)) * Used to hold A in band storage format. * * LDAB (input) INTEGER * The leading dimension of AB. It must be at least 2 (not 1!) * and at least max( KK )+1. * * BD (workspace) REAL array, dimension (max(NN)) * Used to hold the diagonal of the bidiagonal matrix computed * by SGBBRD. * * BE (workspace) REAL array, dimension (max(NN)) * Used to hold the off-diagonal of the bidiagonal matrix * computed by SGBBRD. * * Q (workspace) REAL array, dimension (LDQ, max(NN)) * Used to hold the orthogonal matrix Q computed by SGBBRD. * * LDQ (input) INTEGER * The leading dimension of Q. It must be at least 1 * and at least max( NN ). * * P (workspace) REAL array, dimension (LDP, max(NN)) * Used to hold the orthogonal matrix P computed by SGBBRD. * * LDP (input) INTEGER * The leading dimension of P. It must be at least 1 * and at least max( NN ). * * C (workspace) REAL array, dimension (LDC, max(NN)) * Used to hold the matrix C updated by SGBBRD. * * LDC (input) INTEGER * The leading dimension of U. It must be at least 1 * and at least max( NN ). * * CC (workspace) REAL array, dimension (LDC, max(NN)) * Used to hold a copy of the matrix C. * * WORK (workspace) REAL array, dimension (LWORK) * * LWORK (input) INTEGER * The number of entries in WORK. This must be at least * max( LDA+1, max(NN)+1 )*max(NN). * * RESULT (output) REAL array, dimension (4) * The values computed by the tests described above. * The values are currently limited to 1/ulp, to avoid * overflow. * * INFO (output) INTEGER * If 0, then everything ran OK. * *----------------------------------------------------------------------- * * Some Local Variables and Parameters: * ---- ----- --------- --- ---------- * ZERO, ONE Real 0 and 1. * MAXTYP The number of types defined. * NTEST The number of tests performed, or which can * be performed so far, for the current matrix. * NTESTT The total number of tests performed so far. * NMAX Largest value in NN. * NMATS The number of matrices generated so far. * NERRS The number of tests which have exceeded THRESH * so far. * COND, IMODE Values to be passed to the matrix generators. * ANORM Norm of A; passed to matrix generators. * * OVFL, UNFL Overflow and underflow thresholds. * ULP, ULPINV Finest relative precision and its inverse. * RTOVFL, RTUNFL Square roots of the previous 2 values. * The following four arrays decode JTYPE: * KTYPE(j) The general type (1-10) for type "j". * KMODE(j) The MODE value to be passed to the matrix * generator for type "j". * KMAGN(j) The order of magnitude ( O(1), * O(overflow^(1/2) ), O(underflow^(1/2) ) * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) INTEGER MAXTYP PARAMETER ( MAXTYP = 15 ) * .. * .. Local Scalars .. LOGICAL BADMM, BADNN, BADNNB INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JR, JSIZE, $ JTYPE, JWIDTH, K, KL, KMAX, KU, M, MMAX, MNMAX, $ MNMIN, MTYPES, N, NERRS, NMATS, NMAX, NTEST, $ NTESTT REAL AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL, ULP, $ ULPINV, UNFL * .. * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ), $ KMODE( MAXTYP ), KTYPE( MAXTYP ) * .. * .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH * .. * .. External Subroutines .. EXTERNAL SBDT01, SBDT02, SGBBRD, SLACPY, SLAHD2, SLASET, $ SLASUM, SLATMR, SLATMS, SORT01, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, REAL, SQRT * .. * .. Data statements .. DATA KTYPE / 1, 2, 5*4, 5*6, 3*9 / DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3 / DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, $ 0, 0 / * .. * .. Executable Statements .. * * Check for errors * NTESTT = 0 INFO = 0 * * Important constants * BADMM = .FALSE. BADNN = .FALSE. MMAX = 1 NMAX = 1 MNMAX = 1 DO 10 J = 1, NSIZES MMAX = MAX( MMAX, MVAL( J ) ) IF( MVAL( J ).LT.0 ) $ BADMM = .TRUE. NMAX = MAX( NMAX, NVAL( J ) ) IF( NVAL( J ).LT.0 ) $ BADNN = .TRUE. MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) ) 10 CONTINUE * BADNNB = .FALSE. KMAX = 0 DO 20 J = 1, NWDTHS KMAX = MAX( KMAX, KK( J ) ) IF( KK( J ).LT.0 ) $ BADNNB = .TRUE. 20 CONTINUE * * Check for errors * IF( NSIZES.LT.0 ) THEN INFO = -1 ELSE IF( BADMM ) THEN INFO = -2 ELSE IF( BADNN ) THEN INFO = -3 ELSE IF( NWDTHS.LT.0 ) THEN INFO = -4 ELSE IF( BADNNB ) THEN INFO = -5 ELSE IF( NTYPES.LT.0 ) THEN INFO = -6 ELSE IF( NRHS.LT.0 ) THEN INFO = -8 ELSE IF( LDA.LT.NMAX ) THEN INFO = -13 ELSE IF( LDAB.LT.2*KMAX+1 ) THEN INFO = -15 ELSE IF( LDQ.LT.NMAX ) THEN INFO = -19 ELSE IF( LDP.LT.NMAX ) THEN INFO = -21 ELSE IF( LDC.LT.NMAX ) THEN INFO = -23 ELSE IF( ( MAX( LDA, NMAX )+1 )*NMAX.GT.LWORK ) THEN INFO = -26 END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'SCHKBB', -INFO ) RETURN END IF * * Quick return if possible * IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 .OR. NWDTHS.EQ.0 ) $ RETURN * * More Important constants * UNFL = SLAMCH( 'Safe minimum' ) OVFL = ONE / UNFL ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) ULPINV = ONE / ULP RTUNFL = SQRT( UNFL ) RTOVFL = SQRT( OVFL ) * * Loop over sizes, widths, types * NERRS = 0 NMATS = 0 * DO 160 JSIZE = 1, NSIZES M = MVAL( JSIZE ) N = NVAL( JSIZE ) MNMIN = MIN( M, N ) AMNINV = ONE / REAL( MAX( 1, M, N ) ) * DO 150 JWIDTH = 1, NWDTHS K = KK( JWIDTH ) IF( K.GE.M .AND. K.GE.N ) $ GO TO 150 KL = MAX( 0, MIN( M-1, K ) ) KU = MAX( 0, MIN( N-1, K ) ) * IF( NSIZES.NE.1 ) THEN MTYPES = MIN( MAXTYP, NTYPES ) ELSE MTYPES = MIN( MAXTYP+1, NTYPES ) END IF * DO 140 JTYPE = 1, MTYPES IF( .NOT.DOTYPE( JTYPE ) ) $ GO TO 140 NMATS = NMATS + 1 NTEST = 0 * DO 30 J = 1, 4 IOLDSD( J ) = ISEED( J ) 30 CONTINUE * * Compute "A". * * Control parameters: * * KMAGN KMODE KTYPE * =1 O(1) clustered 1 zero * =2 large clustered 2 identity * =3 small exponential (none) * =4 arithmetic diagonal, (w/ singular values) * =5 random log (none) * =6 random nonhermitian, w/ singular values * =7 (none) * =8 (none) * =9 random nonhermitian * IF( MTYPES.GT.MAXTYP ) $ GO TO 90 * ITYPE = KTYPE( JTYPE ) IMODE = KMODE( JTYPE ) * * Compute norm * GO TO ( 40, 50, 60 )KMAGN( JTYPE ) * 40 CONTINUE ANORM = ONE GO TO 70 * 50 CONTINUE ANORM = ( RTOVFL*ULP )*AMNINV GO TO 70 * 60 CONTINUE ANORM = RTUNFL*MAX( M, N )*ULPINV GO TO 70 * 70 CONTINUE * CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) CALL SLASET( 'Full', LDAB, N, ZERO, ZERO, AB, LDAB ) IINFO = 0 COND = ULPINV * * Special Matrices -- Identity & Jordan block * * Zero * IF( ITYPE.EQ.1 ) THEN IINFO = 0 * ELSE IF( ITYPE.EQ.2 ) THEN * * Identity * DO 80 JCOL = 1, N A( JCOL, JCOL ) = ANORM 80 CONTINUE * ELSE IF( ITYPE.EQ.4 ) THEN * * Diagonal Matrix, singular values specified * CALL SLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, $ ANORM, 0, 0, 'N', A, LDA, WORK( M+1 ), $ IINFO ) * ELSE IF( ITYPE.EQ.6 ) THEN * * Nonhermitian, singular values specified * CALL SLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, $ ANORM, KL, KU, 'N', A, LDA, WORK( M+1 ), $ IINFO ) * ELSE IF( ITYPE.EQ.9 ) THEN * * Nonhermitian, random entries * CALL SLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, $ 'T', 'N', WORK( N+1 ), 1, ONE, $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, KL, $ KU, ZERO, ANORM, 'N', A, LDA, IDUMMA, $ IINFO ) * ELSE * IINFO = 1 END IF * * Generate Right-Hand Side * CALL SLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE, ONE, $ 'T', 'N', WORK( M+1 ), 1, ONE, $ WORK( 2*M+1 ), 1, ONE, 'N', IDUMMA, M, NRHS, $ ZERO, ONE, 'NO', C, LDC, IDUMMA, IINFO ) * IF( IINFO.NE.0 ) THEN WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, $ JTYPE, IOLDSD INFO = ABS( IINFO ) RETURN END IF * 90 CONTINUE * * Copy A to band storage. * DO 110 J = 1, N DO 100 I = MAX( 1, J-KU ), MIN( M, J+KL ) AB( KU+1+I-J, J ) = A( I, J ) 100 CONTINUE 110 CONTINUE * * Copy C * CALL SLACPY( 'Full', M, NRHS, C, LDC, CC, LDC ) * * Call SGBBRD to compute B, Q and P, and to update C. * CALL SGBBRD( 'B', M, N, NRHS, KL, KU, AB, LDAB, BD, BE, $ Q, LDQ, P, LDP, CC, LDC, WORK, IINFO ) * IF( IINFO.NE.0 ) THEN WRITE( NOUNIT, FMT = 9999 )'SGBBRD', IINFO, N, JTYPE, $ IOLDSD INFO = ABS( IINFO ) IF( IINFO.LT.0 ) THEN RETURN ELSE RESULT( 1 ) = ULPINV GO TO 120 END IF END IF * * Test 1: Check the decomposition A := Q * B * P' * 2: Check the orthogonality of Q * 3: Check the orthogonality of P * 4: Check the computation of Q' * C * CALL SBDT01( M, N, -1, A, LDA, Q, LDQ, BD, BE, P, LDP, $ WORK, RESULT( 1 ) ) CALL SORT01( 'Columns', M, M, Q, LDQ, WORK, LWORK, $ RESULT( 2 ) ) CALL SORT01( 'Rows', N, N, P, LDP, WORK, LWORK, $ RESULT( 3 ) ) CALL SBDT02( M, NRHS, C, LDC, CC, LDC, Q, LDQ, WORK, $ RESULT( 4 ) ) * * End of Loop -- Check for RESULT(j) > THRESH * NTEST = 4 120 CONTINUE NTESTT = NTESTT + NTEST * * Print out tests which fail. * DO 130 JR = 1, NTEST IF( RESULT( JR ).GE.THRESH ) THEN IF( NERRS.EQ.0 ) $ CALL SLAHD2( NOUNIT, 'SBB' ) NERRS = NERRS + 1 WRITE( NOUNIT, FMT = 9998 )M, N, K, IOLDSD, JTYPE, $ JR, RESULT( JR ) END IF 130 CONTINUE * 140 CONTINUE 150 CONTINUE 160 CONTINUE * * Summary * CALL SLASUM( 'SBB', NOUNIT, NERRS, NTESTT ) RETURN * 9999 FORMAT( ' SCHKBB: ', A, ' returned INFO=', I5, '.', / 9X, 'M=', $ I5, ' N=', I5, ' K=', I5, ', JTYPE=', I5, ', ISEED=(', $ 3( I5, ',' ), I5, ')' ) 9998 FORMAT( ' M =', I4, ' N=', I4, ', K=', I3, ', seed=', $ 4( I4, ',' ), ' type ', I2, ', test(', I2, ')=', G10.3 ) * * End of SCHKBB * END