 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 00002 $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 00003 $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 00004 $ IWORK, INFO ) 00005 * 00006 * -- LAPACK auxiliary routine (version 3.3.0) -- 00007 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00008 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00009 * November 2010 00010 * 00011 * .. Scalar Arguments .. 00012 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 00013 $ NR, SQRE 00014 DOUBLE PRECISION ALPHA, BETA, C, S 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 00018 $ PERM( * ) 00019 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ), 00020 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 00021 $ VF( * ), VL( * ), WORK( * ), Z( * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * DLASD6 computes the SVD of an updated upper bidiagonal matrix B 00028 * obtained by merging two smaller ones by appending a row. This 00029 * routine is used only for the problem which requires all singular 00030 * values and optionally singular vector matrices in factored form. 00031 * B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. 00032 * A related subroutine, DLASD1, handles the case in which all singular 00033 * values and singular vectors of the bidiagonal matrix are desired. 00034 * 00035 * DLASD6 computes the SVD as follows: 00036 * 00037 * ( D1(in) 0 0 0 ) 00038 * B = U(in) * ( Z1**T a Z2**T b ) * VT(in) 00039 * ( 0 0 D2(in) 0 ) 00040 * 00041 * = U(out) * ( D(out) 0) * VT(out) 00042 * 00043 * where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M 00044 * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 00045 * elsewhere; and the entry b is empty if SQRE = 0. 00046 * 00047 * The singular values of B can be computed using D1, D2, the first 00048 * components of all the right singular vectors of the lower block, and 00049 * the last components of all the right singular vectors of the upper 00050 * block. These components are stored and updated in VF and VL, 00051 * respectively, in DLASD6. Hence U and VT are not explicitly 00052 * referenced. 00053 * 00054 * The singular values are stored in D. The algorithm consists of two 00055 * stages: 00056 * 00057 * The first stage consists of deflating the size of the problem 00058 * when there are multiple singular values or if there is a zero 00059 * in the Z vector. For each such occurence the dimension of the 00060 * secular equation problem is reduced by one. This stage is 00061 * performed by the routine DLASD7. 00062 * 00063 * The second stage consists of calculating the updated 00064 * singular values. This is done by finding the roots of the 00065 * secular equation via the routine DLASD4 (as called by DLASD8). 00066 * This routine also updates VF and VL and computes the distances 00067 * between the updated singular values and the old singular 00068 * values. 00069 * 00070 * DLASD6 is called from DLASDA. 00071 * 00072 * Arguments 00073 * ========= 00074 * 00075 * ICOMPQ (input) INTEGER 00076 * Specifies whether singular vectors are to be computed in 00077 * factored form: 00078 * = 0: Compute singular values only. 00079 * = 1: Compute singular vectors in factored form as well. 00080 * 00081 * NL (input) INTEGER 00082 * The row dimension of the upper block. NL >= 1. 00083 * 00084 * NR (input) INTEGER 00085 * The row dimension of the lower block. NR >= 1. 00086 * 00087 * SQRE (input) INTEGER 00088 * = 0: the lower block is an NR-by-NR square matrix. 00089 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00090 * 00091 * The bidiagonal matrix has row dimension N = NL + NR + 1, 00092 * and column dimension M = N + SQRE. 00093 * 00094 * D (input/output) DOUBLE PRECISION array, dimension ( NL+NR+1 ). 00095 * On entry D(1:NL,1:NL) contains the singular values of the 00096 * upper block, and D(NL+2:N) contains the singular values 00097 * of the lower block. On exit D(1:N) contains the singular 00098 * values of the modified matrix. 00099 * 00100 * VF (input/output) DOUBLE PRECISION array, dimension ( M ) 00101 * On entry, VF(1:NL+1) contains the first components of all 00102 * right singular vectors of the upper block; and VF(NL+2:M) 00103 * contains the first components of all right singular vectors 00104 * of the lower block. On exit, VF contains the first components 00105 * of all right singular vectors of the bidiagonal matrix. 00106 * 00107 * VL (input/output) DOUBLE PRECISION array, dimension ( M ) 00108 * On entry, VL(1:NL+1) contains the last components of all 00109 * right singular vectors of the upper block; and VL(NL+2:M) 00110 * contains the last components of all right singular vectors of 00111 * the lower block. On exit, VL contains the last components of 00112 * all right singular vectors of the bidiagonal matrix. 00113 * 00114 * ALPHA (input/output) DOUBLE PRECISION 00115 * Contains the diagonal element associated with the added row. 00116 * 00117 * BETA (input/output) DOUBLE PRECISION 00118 * Contains the off-diagonal element associated with the added 00119 * row. 00120 * 00121 * IDXQ (output) INTEGER array, dimension ( N ) 00122 * This contains the permutation which will reintegrate the 00123 * subproblem just solved back into sorted order, i.e. 00124 * D( IDXQ( I = 1, N ) ) will be in ascending order. 00125 * 00126 * PERM (output) INTEGER array, dimension ( N ) 00127 * The permutations (from deflation and sorting) to be applied 00128 * to each block. Not referenced if ICOMPQ = 0. 00129 * 00130 * GIVPTR (output) INTEGER 00131 * The number of Givens rotations which took place in this 00132 * subproblem. Not referenced if ICOMPQ = 0. 00133 * 00134 * GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) 00135 * Each pair of numbers indicates a pair of columns to take place 00136 * in a Givens rotation. Not referenced if ICOMPQ = 0. 00137 * 00138 * LDGCOL (input) INTEGER 00139 * leading dimension of GIVCOL, must be at least N. 00140 * 00141 * GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00142 * Each number indicates the C or S value to be used in the 00143 * corresponding Givens rotation. Not referenced if ICOMPQ = 0. 00144 * 00145 * LDGNUM (input) INTEGER 00146 * The leading dimension of GIVNUM and POLES, must be at least N. 00147 * 00148 * POLES (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00149 * On exit, POLES(1,*) is an array containing the new singular 00150 * values obtained from solving the secular equation, and 00151 * POLES(2,*) is an array containing the poles in the secular 00152 * equation. Not referenced if ICOMPQ = 0. 00153 * 00154 * DIFL (output) DOUBLE PRECISION array, dimension ( N ) 00155 * On exit, DIFL(I) is the distance between I-th updated 00156 * (undeflated) singular value and the I-th (undeflated) old 00157 * singular value. 00158 * 00159 * DIFR (output) DOUBLE PRECISION array, 00160 * dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and 00161 * dimension ( N ) if ICOMPQ = 0. 00162 * On exit, DIFR(I, 1) is the distance between I-th updated 00163 * (undeflated) singular value and the I+1-th (undeflated) old 00164 * singular value. 00165 * 00166 * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 00167 * normalizing factors for the right singular vector matrix. 00168 * 00169 * See DLASD8 for details on DIFL and DIFR. 00170 * 00171 * Z (output) DOUBLE PRECISION array, dimension ( M ) 00172 * The first elements of this array contain the components 00173 * of the deflation-adjusted updating row vector. 00174 * 00175 * K (output) INTEGER 00176 * Contains the dimension of the non-deflated matrix, 00177 * This is the order of the related secular equation. 1 <= K <=N. 00178 * 00179 * C (output) DOUBLE PRECISION 00180 * C contains garbage if SQRE =0 and the C-value of a Givens 00181 * rotation related to the right null space if SQRE = 1. 00182 * 00183 * S (output) DOUBLE PRECISION 00184 * S contains garbage if SQRE =0 and the S-value of a Givens 00185 * rotation related to the right null space if SQRE = 1. 00186 * 00187 * WORK (workspace) DOUBLE PRECISION array, dimension ( 4 * M ) 00188 * 00189 * IWORK (workspace) INTEGER array, dimension ( 3 * N ) 00190 * 00191 * INFO (output) INTEGER 00192 * = 0: successful exit. 00193 * < 0: if INFO = -i, the i-th argument had an illegal value. 00194 * > 0: if INFO = 1, a singular value did not converge 00195 * 00196 * Further Details 00197 * =============== 00198 * 00199 * Based on contributions by 00200 * Ming Gu and Huan Ren, Computer Science Division, University of 00201 * California at Berkeley, USA 00202 * 00203 * ===================================================================== 00204 * 00205 * .. Parameters .. 00206 DOUBLE PRECISION ONE, ZERO 00207 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00208 * .. 00209 * .. Local Scalars .. 00210 INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M, 00211 $ N, N1, N2 00212 DOUBLE PRECISION ORGNRM 00213 * .. 00214 * .. External Subroutines .. 00215 EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA 00216 * .. 00217 * .. Intrinsic Functions .. 00218 INTRINSIC ABS, MAX 00219 * .. 00220 * .. Executable Statements .. 00221 * 00222 * Test the input parameters. 00223 * 00224 INFO = 0 00225 N = NL + NR + 1 00226 M = N + SQRE 00227 * 00228 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00229 INFO = -1 00230 ELSE IF( NL.LT.1 ) THEN 00231 INFO = -2 00232 ELSE IF( NR.LT.1 ) THEN 00233 INFO = -3 00234 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00235 INFO = -4 00236 ELSE IF( LDGCOL.LT.N ) THEN 00237 INFO = -14 00238 ELSE IF( LDGNUM.LT.N ) THEN 00239 INFO = -16 00240 END IF 00241 IF( INFO.NE.0 ) THEN 00242 CALL XERBLA( 'DLASD6', -INFO ) 00243 RETURN 00244 END IF 00245 * 00246 * The following values are for bookkeeping purposes only. They are 00247 * integer pointers which indicate the portion of the workspace 00248 * used by a particular array in DLASD7 and DLASD8. 00249 * 00250 ISIGMA = 1 00251 IW = ISIGMA + N 00252 IVFW = IW + M 00253 IVLW = IVFW + M 00254 * 00255 IDX = 1 00256 IDXC = IDX + N 00257 IDXP = IDXC + N 00258 * 00259 * Scale. 00260 * 00261 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 00262 D( NL+1 ) = ZERO 00263 DO 10 I = 1, N 00264 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 00265 ORGNRM = ABS( D( I ) ) 00266 END IF 00267 10 CONTINUE 00268 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00269 ALPHA = ALPHA / ORGNRM 00270 BETA = BETA / ORGNRM 00271 * 00272 * Sort and Deflate singular values. 00273 * 00274 CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF, 00275 $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA, 00276 $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ, 00277 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 00278 $ INFO ) 00279 * 00280 * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL. 00281 * 00282 CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, 00283 $ WORK( ISIGMA ), WORK( IW ), INFO ) 00284 * 00285 * Handle error returned 00286 * 00287 IF( INFO.NE.0 ) THEN 00288 CALL XERBLA( 'DLASD8', -INFO ) 00289 RETURN 00290 END IF 00291 * 00292 * Save the poles if ICOMPQ = 1. 00293 * 00294 IF( ICOMPQ.EQ.1 ) THEN 00295 CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 ) 00296 CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 ) 00297 END IF 00298 * 00299 * Unscale. 00300 * 00301 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00302 * 00303 * Prepare the IDXQ sorting permutation. 00304 * 00305 N1 = K 00306 N2 = N - K 00307 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 00308 * 00309 RETURN 00310 * 00311 * End of DLASD6 00312 * 00313 END
1.7.3 
