cchkbd.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 864 行 · 第 1/3 页

               IF( M.GE.N ) THEN
                  UPLO = 'U'
               ELSE
                  UPLO = 'L'
               END IF
            ELSE
               IINFO = 1
            END IF
*
            IF( IINFO.EQ.0 ) THEN
*
*              Generate Right-Hand Side
*
               IF( BIDIAG ) THEN
                  CALL CLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6,
     $                         ONE, CONE, 'T', 'N', WORK( MNMIN+1 ), 1,
     $                         ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N',
     $                         IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y,
     $                         LDX, IWORK, IINFO )
               ELSE
                  CALL CLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE,
     $                         CONE, 'T', 'N', WORK( M+1 ), 1, ONE,
     $                         WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M,
     $                         NRHS, ZERO, ONE, 'NO', X, LDX, IWORK,
     $                         IINFO )
               END IF
            END IF
*
*           Error Exit
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               RETURN
            END IF
*
  100       CONTINUE
*
*           Call CGEBRD and CUNGBR to compute B, Q, and P, do tests.
*
            IF( .NOT.BIDIAG ) THEN
*
*              Compute transformations to reduce A to bidiagonal form:
*              B := Q' * A * P.
*
               CALL CLACPY( ' ', M, N, A, LDA, Q, LDQ )
               CALL CGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ),
     $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
*
*              Check error code from CGEBRD.
*
               IF( IINFO.NE.0 ) THEN
                  WRITE( NOUT, FMT = 9998 )'CGEBRD', IINFO, M, N,
     $               JTYPE, IOLDSD
                  INFO = ABS( IINFO )
                  RETURN
               END IF
*
               CALL CLACPY( ' ', M, N, Q, LDQ, PT, LDPT )
               IF( M.GE.N ) THEN
                  UPLO = 'U'
               ELSE
                  UPLO = 'L'
               END IF
*
*              Generate Q
*
               MQ = M
               IF( NRHS.LE.0 )
     $            MQ = MNMIN
               CALL CUNGBR( 'Q', M, MQ, N, Q, LDQ, WORK,
     $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
*
*              Check error code from CUNGBR.
*
               IF( IINFO.NE.0 ) THEN
                  WRITE( NOUT, FMT = 9998 )'CUNGBR(Q)', IINFO, M, N,
     $               JTYPE, IOLDSD
                  INFO = ABS( IINFO )
                  RETURN
               END IF
*
*              Generate P'
*
               CALL CUNGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ),
     $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
*
*              Check error code from CUNGBR.
*
               IF( IINFO.NE.0 ) THEN
                  WRITE( NOUT, FMT = 9998 )'CUNGBR(P)', IINFO, M, N,
     $               JTYPE, IOLDSD
                  INFO = ABS( IINFO )
                  RETURN
               END IF
*
*              Apply Q' to an M by NRHS matrix X:  Y := Q' * X.
*
               CALL CGEMM( 'Conjugate transpose', 'No transpose', M,
     $                     NRHS, M, CONE, Q, LDQ, X, LDX, CZERO, Y,
     $                     LDX )
*
*              Test 1:  Check the decomposition A := Q * B * PT
*                   2:  Check the orthogonality of Q
*                   3:  Check the orthogonality of PT
*
               CALL CBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT,
     $                      WORK, RWORK, RESULT( 1 ) )
               CALL CUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
     $                      RWORK, RESULT( 2 ) )
               CALL CUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
     $                      RWORK, RESULT( 3 ) )
            END IF
*
*           Use CBDSQR to form the SVD of the bidiagonal matrix B:
*           B := U * S1 * VT, and compute Z = U' * Y.
*
            CALL SCOPY( MNMIN, BD, 1, S1, 1 )
            IF( MNMIN.GT.0 )
     $         CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
            CALL CLACPY( ' ', M, NRHS, Y, LDX, Z, LDX )
            CALL CLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, U, LDPT )
            CALL CLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, VT, LDPT )
*
            CALL CBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, RWORK, VT,
     $                   LDPT, U, LDPT, Z, LDX, RWORK( MNMIN+1 ),
     $                   IINFO )
*
*           Check error code from CBDSQR.
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'CBDSQR(vects)', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 ) THEN
                  RETURN
               ELSE
                  RESULT( 4 ) = ULPINV
                  GO TO 150
               END IF
            END IF
*
*           Use CBDSQR to compute only the singular values of the
*           bidiagonal matrix B;  U, VT, and Z should not be modified.
*
            CALL SCOPY( MNMIN, BD, 1, S2, 1 )
            IF( MNMIN.GT.0 )
     $         CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
*
            CALL CBDSQR( UPLO, MNMIN, 0, 0, 0, S2, RWORK, VT, LDPT, U,
     $                   LDPT, Z, LDX, RWORK( MNMIN+1 ), IINFO )
*
*           Check error code from CBDSQR.
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'CBDSQR(values)', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 ) THEN
                  RETURN
               ELSE
                  RESULT( 9 ) = ULPINV
                  GO TO 150
               END IF
            END IF
*
*           Test 4:  Check the decomposition B := U * S1 * VT
*                5:  Check the computation Z := U' * Y
*                6:  Check the orthogonality of U
*                7:  Check the orthogonality of VT
*
            CALL CBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
     $                   WORK, RESULT( 4 ) )
            CALL CBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK,
     $                   RWORK, RESULT( 5 ) )
            CALL CUNT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
     $                   RWORK, RESULT( 6 ) )
            CALL CUNT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
     $                   RWORK, RESULT( 7 ) )
*
*           Test 8:  Check that the singular values are sorted in
*                    non-increasing order and are non-negative
*
            RESULT( 8 ) = ZERO
            DO 110 I = 1, MNMIN - 1
               IF( S1( I ).LT.S1( I+1 ) )
     $            RESULT( 8 ) = ULPINV
               IF( S1( I ).LT.ZERO )
     $            RESULT( 8 ) = ULPINV
  110       CONTINUE
            IF( MNMIN.GE.1 ) THEN
               IF( S1( MNMIN ).LT.ZERO )
     $            RESULT( 8 ) = ULPINV
            END IF
*
*           Test 9:  Compare CBDSQR with and without singular vectors
*
            TEMP2 = ZERO
*
            DO 120 J = 1, MNMIN
               TEMP1 = ABS( S1( J )-S2( J ) ) /
     $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
     $                 ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) )
               TEMP2 = MAX( TEMP1, TEMP2 )
  120       CONTINUE
*
            RESULT( 9 ) = TEMP2
*
*           Test 10:  Sturm sequence test of singular values
*                     Go up by factors of two until it succeeds
*
            TEMP1 = THRESH*( HALF-ULP )
*
            DO 130 J = 0, LOG2UI
               CALL SSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
               IF( IINFO.EQ.0 )
     $            GO TO 140
               TEMP1 = TEMP1*TWO
  130       CONTINUE
*
  140       CONTINUE
            RESULT( 10 ) = TEMP1
*
*           Use CBDSQR to form the decomposition A := (QU) S (VT PT)
*           from the bidiagonal form A := Q B PT.
*
            IF( .NOT.BIDIAG ) THEN
               CALL SCOPY( MNMIN, BD, 1, S2, 1 )
               IF( MNMIN.GT.0 )
     $            CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
*
               CALL CBDSQR( UPLO, MNMIN, N, M, NRHS, S2, RWORK, PT,
     $                      LDPT, Q, LDQ, Y, LDX, RWORK( MNMIN+1 ),
     $                      IINFO )
*
*              Test 11:  Check the decomposition A := Q*U * S2 * VT*PT
*                   12:  Check the computation Z := U' * Q' * X
*                   13:  Check the orthogonality of Q*U
*                   14:  Check the orthogonality of VT*PT
*
               CALL CBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT,
     $                      LDPT, WORK, RWORK, RESULT( 11 ) )
               CALL CBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK,
     $                      RWORK, RESULT( 12 ) )
               CALL CUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
     $                      RWORK, RESULT( 13 ) )
               CALL CUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
     $                      RWORK, RESULT( 14 ) )
            END IF
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  150       CONTINUE
            DO 160 J = 1, 14
               IF( RESULT( J ).GE.THRESH ) THEN
                  IF( NFAIL.EQ.0 )
     $               CALL SLAHD2( NOUT, PATH )
                  WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J,
     $               RESULT( J )
                  NFAIL = NFAIL + 1
               END IF
  160       CONTINUE
            IF( .NOT.BIDIAG ) THEN
               NTEST = NTEST + 14
            ELSE
               NTEST = NTEST + 5
            END IF
*
  170    CONTINUE
  180 CONTINUE
*
*     Summary
*
      CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 )
*
      RETURN
*
*     End of CCHKBD
*
 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=',
     $      4( I4, ',' ), ' test(', I2, ')=', G11.4 )
 9998 FORMAT( ' CCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
     $      I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
     $      I5, ')' )
*
      END