![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASD6 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASD6 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd6.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd6.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd6.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 00022 * IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 00023 * LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 00024 * IWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 00028 * $ NR, SQRE 00029 * REAL ALPHA, BETA, C, S 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 00033 * $ PERM( * ) 00034 * REAL D( * ), DIFL( * ), DIFR( * ), 00035 * $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 00036 * $ VF( * ), VL( * ), WORK( * ), Z( * ) 00037 * .. 00038 * 00039 * 00040 *> \par Purpose: 00041 * ============= 00042 *> 00043 *> \verbatim 00044 *> 00045 *> SLASD6 computes the SVD of an updated upper bidiagonal matrix B 00046 *> obtained by merging two smaller ones by appending a row. This 00047 *> routine is used only for the problem which requires all singular 00048 *> values and optionally singular vector matrices in factored form. 00049 *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. 00050 *> A related subroutine, SLASD1, handles the case in which all singular 00051 *> values and singular vectors of the bidiagonal matrix are desired. 00052 *> 00053 *> SLASD6 computes the SVD as follows: 00054 *> 00055 *> ( D1(in) 0 0 0 ) 00056 *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in) 00057 *> ( 0 0 D2(in) 0 ) 00058 *> 00059 *> = U(out) * ( D(out) 0) * VT(out) 00060 *> 00061 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M 00062 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 00063 *> elsewhere; and the entry b is empty if SQRE = 0. 00064 *> 00065 *> The singular values of B can be computed using D1, D2, the first 00066 *> components of all the right singular vectors of the lower block, and 00067 *> the last components of all the right singular vectors of the upper 00068 *> block. These components are stored and updated in VF and VL, 00069 *> respectively, in SLASD6. Hence U and VT are not explicitly 00070 *> referenced. 00071 *> 00072 *> The singular values are stored in D. The algorithm consists of two 00073 *> stages: 00074 *> 00075 *> The first stage consists of deflating the size of the problem 00076 *> when there are multiple singular values or if there is a zero 00077 *> in the Z vector. For each such occurence the dimension of the 00078 *> secular equation problem is reduced by one. This stage is 00079 *> performed by the routine SLASD7. 00080 *> 00081 *> The second stage consists of calculating the updated 00082 *> singular values. This is done by finding the roots of the 00083 *> secular equation via the routine SLASD4 (as called by SLASD8). 00084 *> This routine also updates VF and VL and computes the distances 00085 *> between the updated singular values and the old singular 00086 *> values. 00087 *> 00088 *> SLASD6 is called from SLASDA. 00089 *> \endverbatim 00090 * 00091 * Arguments: 00092 * ========== 00093 * 00094 *> \param[in] ICOMPQ 00095 *> \verbatim 00096 *> ICOMPQ is INTEGER 00097 *> Specifies whether singular vectors are to be computed in 00098 *> factored form: 00099 *> = 0: Compute singular values only. 00100 *> = 1: Compute singular vectors in factored form as well. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] NL 00104 *> \verbatim 00105 *> NL is INTEGER 00106 *> The row dimension of the upper block. NL >= 1. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] NR 00110 *> \verbatim 00111 *> NR is INTEGER 00112 *> The row dimension of the lower block. NR >= 1. 00113 *> \endverbatim 00114 *> 00115 *> \param[in] SQRE 00116 *> \verbatim 00117 *> SQRE is INTEGER 00118 *> = 0: the lower block is an NR-by-NR square matrix. 00119 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00120 *> 00121 *> The bidiagonal matrix has row dimension N = NL + NR + 1, 00122 *> and column dimension M = N + SQRE. 00123 *> \endverbatim 00124 *> 00125 *> \param[in,out] D 00126 *> \verbatim 00127 *> D is REAL array, dimension (NL+NR+1). 00128 *> On entry D(1:NL,1:NL) contains the singular values of the 00129 *> upper block, and D(NL+2:N) contains the singular values 00130 *> of the lower block. On exit D(1:N) contains the singular 00131 *> values of the modified matrix. 00132 *> \endverbatim 00133 *> 00134 *> \param[in,out] VF 00135 *> \verbatim 00136 *> VF is REAL array, dimension (M) 00137 *> On entry, VF(1:NL+1) contains the first components of all 00138 *> right singular vectors of the upper block; and VF(NL+2:M) 00139 *> contains the first components of all right singular vectors 00140 *> of the lower block. On exit, VF contains the first components 00141 *> of all right singular vectors of the bidiagonal matrix. 00142 *> \endverbatim 00143 *> 00144 *> \param[in,out] VL 00145 *> \verbatim 00146 *> VL is REAL array, dimension (M) 00147 *> On entry, VL(1:NL+1) contains the last components of all 00148 *> right singular vectors of the upper block; and VL(NL+2:M) 00149 *> contains the last components of all right singular vectors of 00150 *> the lower block. On exit, VL contains the last components of 00151 *> all right singular vectors of the bidiagonal matrix. 00152 *> \endverbatim 00153 *> 00154 *> \param[in,out] ALPHA 00155 *> \verbatim 00156 *> ALPHA is REAL 00157 *> Contains the diagonal element associated with the added row. 00158 *> \endverbatim 00159 *> 00160 *> \param[in,out] BETA 00161 *> \verbatim 00162 *> BETA is REAL 00163 *> Contains the off-diagonal element associated with the added 00164 *> row. 00165 *> \endverbatim 00166 *> 00167 *> \param[out] IDXQ 00168 *> \verbatim 00169 *> IDXQ is INTEGER array, dimension (N) 00170 *> This contains the permutation which will reintegrate the 00171 *> subproblem just solved back into sorted order, i.e. 00172 *> D( IDXQ( I = 1, N ) ) will be in ascending order. 00173 *> \endverbatim 00174 *> 00175 *> \param[out] PERM 00176 *> \verbatim 00177 *> PERM is INTEGER array, dimension ( N ) 00178 *> The permutations (from deflation and sorting) to be applied 00179 *> to each block. Not referenced if ICOMPQ = 0. 00180 *> \endverbatim 00181 *> 00182 *> \param[out] GIVPTR 00183 *> \verbatim 00184 *> GIVPTR is INTEGER 00185 *> The number of Givens rotations which took place in this 00186 *> subproblem. Not referenced if ICOMPQ = 0. 00187 *> \endverbatim 00188 *> 00189 *> \param[out] GIVCOL 00190 *> \verbatim 00191 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 ) 00192 *> Each pair of numbers indicates a pair of columns to take place 00193 *> in a Givens rotation. Not referenced if ICOMPQ = 0. 00194 *> \endverbatim 00195 *> 00196 *> \param[in] LDGCOL 00197 *> \verbatim 00198 *> LDGCOL is INTEGER 00199 *> leading dimension of GIVCOL, must be at least N. 00200 *> \endverbatim 00201 *> 00202 *> \param[out] GIVNUM 00203 *> \verbatim 00204 *> GIVNUM is REAL array, dimension ( LDGNUM, 2 ) 00205 *> Each number indicates the C or S value to be used in the 00206 *> corresponding Givens rotation. Not referenced if ICOMPQ = 0. 00207 *> \endverbatim 00208 *> 00209 *> \param[in] LDGNUM 00210 *> \verbatim 00211 *> LDGNUM is INTEGER 00212 *> The leading dimension of GIVNUM and POLES, must be at least N. 00213 *> \endverbatim 00214 *> 00215 *> \param[out] POLES 00216 *> \verbatim 00217 *> POLES is REAL array, dimension ( LDGNUM, 2 ) 00218 *> On exit, POLES(1,*) is an array containing the new singular 00219 *> values obtained from solving the secular equation, and 00220 *> POLES(2,*) is an array containing the poles in the secular 00221 *> equation. Not referenced if ICOMPQ = 0. 00222 *> \endverbatim 00223 *> 00224 *> \param[out] DIFL 00225 *> \verbatim 00226 *> DIFL is REAL array, dimension ( N ) 00227 *> On exit, DIFL(I) is the distance between I-th updated 00228 *> (undeflated) singular value and the I-th (undeflated) old 00229 *> singular value. 00230 *> \endverbatim 00231 *> 00232 *> \param[out] DIFR 00233 *> \verbatim 00234 *> DIFR is REAL array, 00235 *> dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and 00236 *> dimension ( N ) if ICOMPQ = 0. 00237 *> On exit, DIFR(I, 1) is the distance between I-th updated 00238 *> (undeflated) singular value and the I+1-th (undeflated) old 00239 *> singular value. 00240 *> 00241 *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 00242 *> normalizing factors for the right singular vector matrix. 00243 *> 00244 *> See SLASD8 for details on DIFL and DIFR. 00245 *> \endverbatim 00246 *> 00247 *> \param[out] Z 00248 *> \verbatim 00249 *> Z is REAL array, dimension ( M ) 00250 *> The first elements of this array contain the components 00251 *> of the deflation-adjusted updating row vector. 00252 *> \endverbatim 00253 *> 00254 *> \param[out] K 00255 *> \verbatim 00256 *> K is INTEGER 00257 *> Contains the dimension of the non-deflated matrix, 00258 *> This is the order of the related secular equation. 1 <= K <=N. 00259 *> \endverbatim 00260 *> 00261 *> \param[out] C 00262 *> \verbatim 00263 *> C is REAL 00264 *> C contains garbage if SQRE =0 and the C-value of a Givens 00265 *> rotation related to the right null space if SQRE = 1. 00266 *> \endverbatim 00267 *> 00268 *> \param[out] S 00269 *> \verbatim 00270 *> S is REAL 00271 *> S contains garbage if SQRE =0 and the S-value of a Givens 00272 *> rotation related to the right null space if SQRE = 1. 00273 *> \endverbatim 00274 *> 00275 *> \param[out] WORK 00276 *> \verbatim 00277 *> WORK is REAL array, dimension ( 4 * M ) 00278 *> \endverbatim 00279 *> 00280 *> \param[out] IWORK 00281 *> \verbatim 00282 *> IWORK is INTEGER array, dimension ( 3 * N ) 00283 *> \endverbatim 00284 *> 00285 *> \param[out] INFO 00286 *> \verbatim 00287 *> INFO is INTEGER 00288 *> = 0: successful exit. 00289 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00290 *> > 0: if INFO = 1, a singular value did not converge 00291 *> \endverbatim 00292 * 00293 * Authors: 00294 * ======== 00295 * 00296 *> \author Univ. of Tennessee 00297 *> \author Univ. of California Berkeley 00298 *> \author Univ. of Colorado Denver 00299 *> \author NAG Ltd. 00300 * 00301 *> \date November 2011 00302 * 00303 *> \ingroup auxOTHERauxiliary 00304 * 00305 *> \par Contributors: 00306 * ================== 00307 *> 00308 *> Ming Gu and Huan Ren, Computer Science Division, University of 00309 *> California at Berkeley, USA 00310 *> 00311 * ===================================================================== 00312 SUBROUTINE SLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 00313 $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 00314 $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 00315 $ IWORK, INFO ) 00316 * 00317 * -- LAPACK auxiliary routine (version 3.4.0) -- 00318 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00319 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00320 * November 2011 00321 * 00322 * .. Scalar Arguments .. 00323 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 00324 $ NR, SQRE 00325 REAL ALPHA, BETA, C, S 00326 * .. 00327 * .. Array Arguments .. 00328 INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 00329 $ PERM( * ) 00330 REAL D( * ), DIFL( * ), DIFR( * ), 00331 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 00332 $ VF( * ), VL( * ), WORK( * ), Z( * ) 00333 * .. 00334 * 00335 * ===================================================================== 00336 * 00337 * .. Parameters .. 00338 REAL ONE, ZERO 00339 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00340 * .. 00341 * .. Local Scalars .. 00342 INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M, 00343 $ N, N1, N2 00344 REAL ORGNRM 00345 * .. 00346 * .. External Subroutines .. 00347 EXTERNAL SCOPY, SLAMRG, SLASCL, SLASD7, SLASD8, XERBLA 00348 * .. 00349 * .. Intrinsic Functions .. 00350 INTRINSIC ABS, MAX 00351 * .. 00352 * .. Executable Statements .. 00353 * 00354 * Test the input parameters. 00355 * 00356 INFO = 0 00357 N = NL + NR + 1 00358 M = N + SQRE 00359 * 00360 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00361 INFO = -1 00362 ELSE IF( NL.LT.1 ) THEN 00363 INFO = -2 00364 ELSE IF( NR.LT.1 ) THEN 00365 INFO = -3 00366 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00367 INFO = -4 00368 ELSE IF( LDGCOL.LT.N ) THEN 00369 INFO = -14 00370 ELSE IF( LDGNUM.LT.N ) THEN 00371 INFO = -16 00372 END IF 00373 IF( INFO.NE.0 ) THEN 00374 CALL XERBLA( 'SLASD6', -INFO ) 00375 RETURN 00376 END IF 00377 * 00378 * The following values are for bookkeeping purposes only. They are 00379 * integer pointers which indicate the portion of the workspace 00380 * used by a particular array in SLASD7 and SLASD8. 00381 * 00382 ISIGMA = 1 00383 IW = ISIGMA + N 00384 IVFW = IW + M 00385 IVLW = IVFW + M 00386 * 00387 IDX = 1 00388 IDXC = IDX + N 00389 IDXP = IDXC + N 00390 * 00391 * Scale. 00392 * 00393 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 00394 D( NL+1 ) = ZERO 00395 DO 10 I = 1, N 00396 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 00397 ORGNRM = ABS( D( I ) ) 00398 END IF 00399 10 CONTINUE 00400 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00401 ALPHA = ALPHA / ORGNRM 00402 BETA = BETA / ORGNRM 00403 * 00404 * Sort and Deflate singular values. 00405 * 00406 CALL SLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF, 00407 $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA, 00408 $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ, 00409 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 00410 $ INFO ) 00411 * 00412 * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL. 00413 * 00414 CALL SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, 00415 $ WORK( ISIGMA ), WORK( IW ), INFO ) 00416 * 00417 * Handle error returned 00418 * 00419 IF( INFO.NE.0 ) THEN 00420 CALL XERBLA( 'SLASD8', -INFO ) 00421 RETURN 00422 END IF 00423 * 00424 * Save the poles if ICOMPQ = 1. 00425 * 00426 IF( ICOMPQ.EQ.1 ) THEN 00427 CALL SCOPY( K, D, 1, POLES( 1, 1 ), 1 ) 00428 CALL SCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 ) 00429 END IF 00430 * 00431 * Unscale. 00432 * 00433 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00434 * 00435 * Prepare the IDXQ sorting permutation. 00436 * 00437 N1 = K 00438 N2 = N - K 00439 CALL SLAMRG( N1, N2, D, 1, -1, IDXQ ) 00440 * 00441 RETURN 00442 * 00443 * End of SLASD6 00444 * 00445 END