![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASD1 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASD1 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd1.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd1.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd1.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 00022 * IDXQ, IWORK, WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDU, LDVT, NL, NR, SQRE 00026 * REAL ALPHA, BETA 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IDXQ( * ), IWORK( * ) 00030 * REAL D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, 00040 *> where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0. 00041 *> 00042 *> A related subroutine SLASD7 handles the case in which the singular 00043 *> values (and the singular vectors in factored form) are desired. 00044 *> 00045 *> SLASD1 computes the SVD as follows: 00046 *> 00047 *> ( D1(in) 0 0 0 ) 00048 *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in) 00049 *> ( 0 0 D2(in) 0 ) 00050 *> 00051 *> = U(out) * ( D(out) 0) * VT(out) 00052 *> 00053 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M 00054 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 00055 *> elsewhere; and the entry b is empty if SQRE = 0. 00056 *> 00057 *> The left singular vectors of the original matrix are stored in U, and 00058 *> the transpose of the right singular vectors are stored in VT, and the 00059 *> singular values are in D. The algorithm consists of three stages: 00060 *> 00061 *> The first stage consists of deflating the size of the problem 00062 *> when there are multiple singular values or when there are zeros in 00063 *> the Z vector. For each such occurence the dimension of the 00064 *> secular equation problem is reduced by one. This stage is 00065 *> performed by the routine SLASD2. 00066 *> 00067 *> The second stage consists of calculating the updated 00068 *> singular values. This is done by finding the square roots of the 00069 *> roots of the secular equation via the routine SLASD4 (as called 00070 *> by SLASD3). This routine also calculates the singular vectors of 00071 *> the current problem. 00072 *> 00073 *> The final stage consists of computing the updated singular vectors 00074 *> directly using the updated singular values. The singular vectors 00075 *> for the current problem are multiplied with the singular vectors 00076 *> from the overall problem. 00077 *> \endverbatim 00078 * 00079 * Arguments: 00080 * ========== 00081 * 00082 *> \param[in] NL 00083 *> \verbatim 00084 *> NL is INTEGER 00085 *> The row dimension of the upper block. NL >= 1. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] NR 00089 *> \verbatim 00090 *> NR is INTEGER 00091 *> The row dimension of the lower block. NR >= 1. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] SQRE 00095 *> \verbatim 00096 *> SQRE is INTEGER 00097 *> = 0: the lower block is an NR-by-NR square matrix. 00098 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00099 *> 00100 *> The bidiagonal matrix has row dimension N = NL + NR + 1, 00101 *> and column dimension M = N + SQRE. 00102 *> \endverbatim 00103 *> 00104 *> \param[in,out] D 00105 *> \verbatim 00106 *> D is REAL array, dimension (NL+NR+1). 00107 *> N = NL+NR+1 00108 *> On entry D(1:NL,1:NL) contains the singular values of the 00109 *> upper block; and D(NL+2:N) contains the singular values of 00110 *> the lower block. On exit D(1:N) contains the singular values 00111 *> of the modified matrix. 00112 *> \endverbatim 00113 *> 00114 *> \param[in,out] ALPHA 00115 *> \verbatim 00116 *> ALPHA is REAL 00117 *> Contains the diagonal element associated with the added row. 00118 *> \endverbatim 00119 *> 00120 *> \param[in,out] BETA 00121 *> \verbatim 00122 *> BETA is REAL 00123 *> Contains the off-diagonal element associated with the added 00124 *> row. 00125 *> \endverbatim 00126 *> 00127 *> \param[in,out] U 00128 *> \verbatim 00129 *> U is REAL array, dimension (LDU,N) 00130 *> On entry U(1:NL, 1:NL) contains the left singular vectors of 00131 *> the upper block; U(NL+2:N, NL+2:N) contains the left singular 00132 *> vectors of the lower block. On exit U contains the left 00133 *> singular vectors of the bidiagonal matrix. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] LDU 00137 *> \verbatim 00138 *> LDU is INTEGER 00139 *> The leading dimension of the array U. LDU >= max( 1, N ). 00140 *> \endverbatim 00141 *> 00142 *> \param[in,out] VT 00143 *> \verbatim 00144 *> VT is REAL array, dimension (LDVT,M) 00145 *> where M = N + SQRE. 00146 *> On entry VT(1:NL+1, 1:NL+1)**T contains the right singular 00147 *> vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains 00148 *> the right singular vectors of the lower block. On exit 00149 *> VT**T contains the right singular vectors of the 00150 *> bidiagonal matrix. 00151 *> \endverbatim 00152 *> 00153 *> \param[in] LDVT 00154 *> \verbatim 00155 *> LDVT is INTEGER 00156 *> The leading dimension of the array VT. LDVT >= max( 1, M ). 00157 *> \endverbatim 00158 *> 00159 *> \param[out] IDXQ 00160 *> \verbatim 00161 *> IDXQ is INTEGER array, dimension (N) 00162 *> This contains the permutation which will reintegrate the 00163 *> subproblem just solved back into sorted order, i.e. 00164 *> D( IDXQ( I = 1, N ) ) will be in ascending order. 00165 *> \endverbatim 00166 *> 00167 *> \param[out] IWORK 00168 *> \verbatim 00169 *> IWORK is INTEGER array, dimension (4*N) 00170 *> \endverbatim 00171 *> 00172 *> \param[out] WORK 00173 *> \verbatim 00174 *> WORK is REAL array, dimension (3*M**2+2*M) 00175 *> \endverbatim 00176 *> 00177 *> \param[out] INFO 00178 *> \verbatim 00179 *> INFO is INTEGER 00180 *> = 0: successful exit. 00181 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00182 *> > 0: if INFO = 1, a singular value did not converge 00183 *> \endverbatim 00184 * 00185 * Authors: 00186 * ======== 00187 * 00188 *> \author Univ. of Tennessee 00189 *> \author Univ. of California Berkeley 00190 *> \author Univ. of Colorado Denver 00191 *> \author NAG Ltd. 00192 * 00193 *> \date November 2011 00194 * 00195 *> \ingroup auxOTHERauxiliary 00196 * 00197 *> \par Contributors: 00198 * ================== 00199 *> 00200 *> Ming Gu and Huan Ren, Computer Science Division, University of 00201 *> California at Berkeley, USA 00202 *> 00203 * ===================================================================== 00204 SUBROUTINE SLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 00205 $ IDXQ, IWORK, WORK, INFO ) 00206 * 00207 * -- LAPACK auxiliary routine (version 3.4.0) -- 00208 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00210 * November 2011 00211 * 00212 * .. Scalar Arguments .. 00213 INTEGER INFO, LDU, LDVT, NL, NR, SQRE 00214 REAL ALPHA, BETA 00215 * .. 00216 * .. Array Arguments .. 00217 INTEGER IDXQ( * ), IWORK( * ) 00218 REAL D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 00219 * .. 00220 * 00221 * ===================================================================== 00222 * 00223 * .. Parameters .. 00224 * 00225 REAL ONE, ZERO 00226 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00227 * .. 00228 * .. Local Scalars .. 00229 INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2, 00230 $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2 00231 REAL ORGNRM 00232 * .. 00233 * .. External Subroutines .. 00234 EXTERNAL SLAMRG, SLASCL, SLASD2, SLASD3, XERBLA 00235 * .. 00236 * .. Intrinsic Functions .. 00237 INTRINSIC ABS, MAX 00238 * .. 00239 * .. Executable Statements .. 00240 * 00241 * Test the input parameters. 00242 * 00243 INFO = 0 00244 * 00245 IF( NL.LT.1 ) THEN 00246 INFO = -1 00247 ELSE IF( NR.LT.1 ) THEN 00248 INFO = -2 00249 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00250 INFO = -3 00251 END IF 00252 IF( INFO.NE.0 ) THEN 00253 CALL XERBLA( 'SLASD1', -INFO ) 00254 RETURN 00255 END IF 00256 * 00257 N = NL + NR + 1 00258 M = N + SQRE 00259 * 00260 * The following values are for bookkeeping purposes only. They are 00261 * integer pointers which indicate the portion of the workspace 00262 * used by a particular array in SLASD2 and SLASD3. 00263 * 00264 LDU2 = N 00265 LDVT2 = M 00266 * 00267 IZ = 1 00268 ISIGMA = IZ + M 00269 IU2 = ISIGMA + N 00270 IVT2 = IU2 + LDU2*N 00271 IQ = IVT2 + LDVT2*M 00272 * 00273 IDX = 1 00274 IDXC = IDX + N 00275 COLTYP = IDXC + N 00276 IDXP = COLTYP + N 00277 * 00278 * Scale. 00279 * 00280 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 00281 D( NL+1 ) = ZERO 00282 DO 10 I = 1, N 00283 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 00284 ORGNRM = ABS( D( I ) ) 00285 END IF 00286 10 CONTINUE 00287 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00288 ALPHA = ALPHA / ORGNRM 00289 BETA = BETA / ORGNRM 00290 * 00291 * Deflate singular values. 00292 * 00293 CALL SLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU, 00294 $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2, 00295 $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ), 00296 $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO ) 00297 * 00298 * Solve Secular Equation and update singular vectors. 00299 * 00300 LDQ = K 00301 CALL SLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ), 00302 $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ), 00303 $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ), 00304 $ INFO ) 00305 IF( INFO.NE.0 ) THEN 00306 RETURN 00307 END IF 00308 * 00309 * Unscale. 00310 * 00311 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00312 * 00313 * Prepare the IDXQ sorting permutation. 00314 * 00315 N1 = K 00316 N2 = N - K 00317 CALL SLAMRG( N1, N2, D, 1, -1, IDXQ ) 00318 * 00319 RETURN 00320 * 00321 * End of SLASD1 00322 * 00323 END