dchkbd.f

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

*
               IF( IINFO.NE.0 ) THEN
                  WRITE( NOUT, FMT = 9998 )'DGEBRD', IINFO, M, N,
     $               JTYPE, IOLDSD
                  INFO = ABS( IINFO )
                  RETURN
               END IF
*
               CALL DLACPY( ' ', 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 DORGBR( 'Q', M, MQ, N, Q, LDQ, WORK,
     $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
*
*              Check error code from DORGBR.
*
               IF( IINFO.NE.0 ) THEN
                  WRITE( NOUT, FMT = 9998 )'DORGBR(Q)', IINFO, M, N,
     $               JTYPE, IOLDSD
                  INFO = ABS( IINFO )
                  RETURN
               END IF
*
*              Generate P'
*
               CALL DORGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ),
     $                      WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
*
*              Check error code from DORGBR.
*
               IF( IINFO.NE.0 ) THEN
                  WRITE( NOUT, FMT = 9998 )'DORGBR(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 DGEMM( 'Transpose', 'No transpose', M, NRHS, M, ONE,
     $                     Q, LDQ, X, LDX, ZERO, Y, LDX )
*
*              Test 1:  Check the decomposition A := Q * B * PT
*                   2:  Check the orthogonality of Q
*                   3:  Check the orthogonality of PT
*
               CALL DBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT,
     $                      WORK, RESULT( 1 ) )
               CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
     $                      RESULT( 2 ) )
               CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
     $                      RESULT( 3 ) )
            END IF
*
*           Use DBDSQR to form the SVD of the bidiagonal matrix B:
*           B := U * S1 * VT, and compute Z = U' * Y.
*
            CALL DCOPY( MNMIN, BD, 1, S1, 1 )
            IF( MNMIN.GT.0 )
     $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
            CALL DLACPY( ' ', M, NRHS, Y, LDX, Z, LDX )
            CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT )
            CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT )
*
            CALL DBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, WORK, VT,
     $                   LDPT, U, LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO )
*
*           Check error code from DBDSQR.
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'DBDSQR(vects)', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 ) THEN
                  RETURN
               ELSE
                  RESULT( 4 ) = ULPINV
                  GO TO 170
               END IF
            END IF
*
*           Use DBDSQR to compute only the singular values of the
*           bidiagonal matrix B;  U, VT, and Z should not be modified.
*
            CALL DCOPY( MNMIN, BD, 1, S2, 1 )
            IF( MNMIN.GT.0 )
     $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
*
            CALL DBDSQR( UPLO, MNMIN, 0, 0, 0, S2, WORK, VT, LDPT, U,
     $                   LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO )
*
*           Check error code from DBDSQR.
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'DBDSQR(values)', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 ) THEN
                  RETURN
               ELSE
                  RESULT( 9 ) = ULPINV
                  GO TO 170
               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 DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
     $                   WORK, RESULT( 4 ) )
            CALL DBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK,
     $                   RESULT( 5 ) )
            CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
     $                   RESULT( 6 ) )
            CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
     $                   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 DBDSQR 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 DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
               IF( IINFO.EQ.0 )
     $            GO TO 140
               TEMP1 = TEMP1*TWO
  130       CONTINUE
*
  140       CONTINUE
            RESULT( 10 ) = TEMP1
*
*           Use DBDSQR to form the decomposition A := (QU) S (VT PT)
*           from the bidiagonal form A := Q B PT.
*
            IF( .NOT.BIDIAG ) THEN
               CALL DCOPY( MNMIN, BD, 1, S2, 1 )
               IF( MNMIN.GT.0 )
     $            CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
*
               CALL DBDSQR( UPLO, MNMIN, N, M, NRHS, S2, WORK, PT, LDPT,
     $                      Q, LDQ, Y, LDX, WORK( 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 DBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT,
     $                      LDPT, WORK, RESULT( 11 ) )
               CALL DBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK,
     $                      RESULT( 12 ) )
               CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
     $                      RESULT( 13 ) )
               CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
     $                      RESULT( 14 ) )
            END IF
*
*           Use DBDSDC to form the SVD of the bidiagonal matrix B:
*           B := U * S1 * VT
*
            CALL DCOPY( MNMIN, BD, 1, S1, 1 )
            IF( MNMIN.GT.0 )
     $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
            CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT )
            CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT )
*
            CALL DBDSDC( UPLO, 'I', MNMIN, S1, WORK, U, LDPT, VT, LDPT,
     $                   DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO )
*
*           Check error code from DBDSDC.
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'DBDSDC(vects)', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 ) THEN
                  RETURN
               ELSE
                  RESULT( 15 ) = ULPINV
                  GO TO 170
               END IF
            END IF
*
*           Use DBDSDC to compute only the singular values of the
*           bidiagonal matrix B;  U and VT should not be modified.
*
            CALL DCOPY( MNMIN, BD, 1, S2, 1 )
            IF( MNMIN.GT.0 )
     $         CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 )
*
            CALL DBDSDC( UPLO, 'N', MNMIN, S2, WORK, DUM, 1, DUM, 1,
     $                   DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO )
*
*           Check error code from DBDSDC.
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUT, FMT = 9998 )'DBDSDC(values)', IINFO, M, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 ) THEN
                  RETURN
               ELSE
                  RESULT( 18 ) = ULPINV
                  GO TO 170
               END IF
            END IF
*
*           Test 15:  Check the decomposition B := U * S1 * VT
*                16:  Check the orthogonality of U
*                17:  Check the orthogonality of VT
*
            CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
     $                   WORK, RESULT( 15 ) )
            CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
     $                   RESULT( 16 ) )
            CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
     $                   RESULT( 17 ) )
*
*           Test 18:  Check that the singular values are sorted in
*                     non-increasing order and are non-negative
*
            RESULT( 18 ) = ZERO
            DO 150 I = 1, MNMIN - 1
               IF( S1( I ).LT.S1( I+1 ) )
     $            RESULT( 18 ) = ULPINV
               IF( S1( I ).LT.ZERO )
     $            RESULT( 18 ) = ULPINV
  150       CONTINUE
            IF( MNMIN.GE.1 ) THEN
               IF( S1( MNMIN ).LT.ZERO )
     $            RESULT( 18 ) = ULPINV
            END IF
*
*           Test 19:  Compare DBDSQR with and without singular vectors
*
            TEMP2 = ZERO
*
            DO 160 J = 1, MNMIN
               TEMP1 = ABS( S1( J )-S2( J ) ) /
     $                 MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
     $                 ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) )
               TEMP2 = MAX( TEMP1, TEMP2 )
  160       CONTINUE
*
            RESULT( 19 ) = TEMP2
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  170       CONTINUE
            DO 180 J = 1, 19
               IF( RESULT( J ).GE.THRESH ) THEN
                  IF( NFAIL.EQ.0 )
     $               CALL DLAHD2( NOUT, PATH )
                  WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J,
     $               RESULT( J )
                  NFAIL = NFAIL + 1
               END IF
  180       CONTINUE
            IF( .NOT.BIDIAG ) THEN
               NTEST = NTEST + 19
            ELSE
               NTEST = NTEST + 5
            END IF
*
  190    CONTINUE
  200 CONTINUE
*
*     Summary
*
      CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 )
*
      RETURN
*
*     End of DCHKBD
*
 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=',
     $      4( I4, ',' ), ' test(', I2, ')=', G11.4 )
 9998 FORMAT( ' DCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
     $      I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
     $      I5, ')' )
*
      END