LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dbdt04()

subroutine dbdt04 ( character  UPLO,
integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  S,
integer  NS,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldvt, * )  VT,
integer  LDVT,
double precision, dimension( * )  WORK,
double precision  RESID 
)

DBDT04

Purpose:
 DBDT04 reconstructs a bidiagonal matrix B from its (partial) SVD:
    S = U' * B * V
 where U and V are orthogonal matrices and S is diagonal.

 The test ratio to test the singular value decomposition is
    RESID = norm( S - U' * B * V ) / ( n * norm(B) * EPS )
 where VT = V' and EPS is the machine precision.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the matrix B is upper or lower bidiagonal.
          = 'U':  Upper bidiagonal
          = 'L':  Lower bidiagonal
[in]N
          N is INTEGER
          The order of the matrix B.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the bidiagonal matrix B.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) superdiagonal elements of the bidiagonal matrix B
          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
          UPLO = 'L'.
[in]S
          S is DOUBLE PRECISION array, dimension (NS)
          The singular values from the (partial) SVD of B, sorted in
          decreasing order.
[in]NS
          NS is INTEGER
          The number of singular values/vectors from the (partial)
          SVD of B.
[in]U
          U is DOUBLE PRECISION array, dimension (LDU,NS)
          The n by ns orthogonal matrix U in S = U' * B * V.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= max(1,N)
[in]VT
          VT is DOUBLE PRECISION array, dimension (LDVT,N)
          The n by ns orthogonal matrix V in S = U' * B * V.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
[out]RESID
          RESID is DOUBLE PRECISION
          The test ratio:  norm(S - U' * B * V) / ( n * norm(B) * EPS )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 133 of file dbdt04.f.

133 *
134 * -- LAPACK test routine (version 3.7.0) --
135 * -- LAPACK is a software package provided by Univ. of Tennessee, --
136 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137 * December 2016
138 *
139 * .. Scalar Arguments ..
140  CHARACTER uplo
141  INTEGER ldu, ldvt, n, ns
142  DOUBLE PRECISION resid
143 * ..
144 * .. Array Arguments ..
145  DOUBLE PRECISION d( * ), e( * ), s( * ), u( ldu, * ),
146  $ vt( ldvt, * ), work( * )
147 * ..
148 *
149 * ======================================================================
150 *
151 * .. Parameters ..
152  DOUBLE PRECISION zero, one
153  parameter( zero = 0.0d+0, one = 1.0d+0 )
154 * ..
155 * .. Local Scalars ..
156  INTEGER i, j, k
157  DOUBLE PRECISION bnorm, eps
158 * ..
159 * .. External Functions ..
160  LOGICAL lsame
161  INTEGER idamax
162  DOUBLE PRECISION dasum, dlamch
163  EXTERNAL lsame, idamax, dasum, dlamch
164 * ..
165 * .. External Subroutines ..
166  EXTERNAL dgemm
167 * ..
168 * .. Intrinsic Functions ..
169  INTRINSIC abs, dble, max, min
170 * ..
171 * .. Executable Statements ..
172 *
173 * Quick return if possible.
174 *
175  resid = zero
176  IF( n.LE.0 .OR. ns.LE.0 )
177  $ RETURN
178 *
179  eps = dlamch( 'Precision' )
180 *
181 * Compute S - U' * B * V.
182 *
183  bnorm = zero
184 *
185  IF( lsame( uplo, 'U' ) ) THEN
186 *
187 * B is upper bidiagonal.
188 *
189  k = 0
190  DO 20 i = 1, ns
191  DO 10 j = 1, n-1
192  k = k + 1
193  work( k ) = d( j )*vt( i, j ) + e( j )*vt( i, j+1 )
194  10 CONTINUE
195  k = k + 1
196  work( k ) = d( n )*vt( i, n )
197  20 CONTINUE
198  bnorm = abs( d( 1 ) )
199  DO 30 i = 2, n
200  bnorm = max( bnorm, abs( d( i ) )+abs( e( i-1 ) ) )
201  30 CONTINUE
202  ELSE
203 *
204 * B is lower bidiagonal.
205 *
206  k = 0
207  DO 50 i = 1, ns
208  k = k + 1
209  work( k ) = d( 1 )*vt( i, 1 )
210  DO 40 j = 1, n-1
211  k = k + 1
212  work( k ) = e( j )*vt( i, j ) + d( j+1 )*vt( i, j+1 )
213  40 CONTINUE
214  50 CONTINUE
215  bnorm = abs( d( n ) )
216  DO 60 i = 1, n-1
217  bnorm = max( bnorm, abs( d( i ) )+abs( e( i ) ) )
218  60 CONTINUE
219  END IF
220 *
221  CALL dgemm( 'T', 'N', ns, ns, n, -one, u, ldu, work( 1 ),
222  $ n, zero, work( 1+n*ns ), ns )
223 *
224 * norm(S - U' * B * V)
225 *
226  k = n*ns
227  DO 70 i = 1, ns
228  work( k+i ) = work( k+i ) + s( i )
229  resid = max( resid, dasum( ns, work( k+1 ), 1 ) )
230  k = k + ns
231  70 CONTINUE
232 *
233  IF( bnorm.LE.zero ) THEN
234  IF( resid.NE.zero )
235  $ resid = one / eps
236  ELSE
237  IF( bnorm.GE.resid ) THEN
238  resid = ( resid / bnorm ) / ( dble( n )*eps )
239  ELSE
240  IF( bnorm.LT.one ) THEN
241  resid = ( min( resid, dble( n )*bnorm ) / bnorm ) /
242  $ ( dble( n )*eps )
243  ELSE
244  resid = min( resid / bnorm, dble( n ) ) /
245  $ ( dble( n )*eps )
246  END IF
247  END IF
248  END IF
249 *
250  RETURN
251 *
252 * End of DBDT04
253 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:73
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:73
Here is the call graph for this function:
Here is the caller graph for this function: