 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 00002 $ IDXQ, IWORK, WORK, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDU, LDVT, NL, NR, SQRE 00011 DOUBLE PRECISION ALPHA, BETA 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IDXQ( * ), IWORK( * ) 00015 DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, 00022 * where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0. 00023 * 00024 * A related subroutine DLASD7 handles the case in which the singular 00025 * values (and the singular vectors in factored form) are desired. 00026 * 00027 * DLASD1 computes the SVD as follows: 00028 * 00029 * ( D1(in) 0 0 0 ) 00030 * B = U(in) * ( Z1**T a Z2**T b ) * VT(in) 00031 * ( 0 0 D2(in) 0 ) 00032 * 00033 * = U(out) * ( D(out) 0) * VT(out) 00034 * 00035 * where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M 00036 * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 00037 * elsewhere; and the entry b is empty if SQRE = 0. 00038 * 00039 * The left singular vectors of the original matrix are stored in U, and 00040 * the transpose of the right singular vectors are stored in VT, and the 00041 * singular values are in D. The algorithm consists of three stages: 00042 * 00043 * The first stage consists of deflating the size of the problem 00044 * when there are multiple singular values or when there are zeros in 00045 * the Z vector. For each such occurence the dimension of the 00046 * secular equation problem is reduced by one. This stage is 00047 * performed by the routine DLASD2. 00048 * 00049 * The second stage consists of calculating the updated 00050 * singular values. This is done by finding the square roots of the 00051 * roots of the secular equation via the routine DLASD4 (as called 00052 * by DLASD3). This routine also calculates the singular vectors of 00053 * the current problem. 00054 * 00055 * The final stage consists of computing the updated singular vectors 00056 * directly using the updated singular values. The singular vectors 00057 * for the current problem are multiplied with the singular vectors 00058 * from the overall problem. 00059 * 00060 * Arguments 00061 * ========= 00062 * 00063 * NL (input) INTEGER 00064 * The row dimension of the upper block. NL >= 1. 00065 * 00066 * NR (input) INTEGER 00067 * The row dimension of the lower block. NR >= 1. 00068 * 00069 * SQRE (input) INTEGER 00070 * = 0: the lower block is an NR-by-NR square matrix. 00071 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00072 * 00073 * The bidiagonal matrix has row dimension N = NL + NR + 1, 00074 * and column dimension M = N + SQRE. 00075 * 00076 * D (input/output) DOUBLE PRECISION array, 00077 * dimension (N = NL+NR+1). 00078 * On entry D(1:NL,1:NL) contains the singular values of the 00079 * upper block; and D(NL+2:N) contains the singular values of 00080 * the lower block. On exit D(1:N) contains the singular values 00081 * of the modified matrix. 00082 * 00083 * ALPHA (input/output) DOUBLE PRECISION 00084 * Contains the diagonal element associated with the added row. 00085 * 00086 * BETA (input/output) DOUBLE PRECISION 00087 * Contains the off-diagonal element associated with the added 00088 * row. 00089 * 00090 * U (input/output) DOUBLE PRECISION array, dimension(LDU,N) 00091 * On entry U(1:NL, 1:NL) contains the left singular vectors of 00092 * the upper block; U(NL+2:N, NL+2:N) contains the left singular 00093 * vectors of the lower block. On exit U contains the left 00094 * singular vectors of the bidiagonal matrix. 00095 * 00096 * LDU (input) INTEGER 00097 * The leading dimension of the array U. LDU >= max( 1, N ). 00098 * 00099 * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M) 00100 * where M = N + SQRE. 00101 * On entry VT(1:NL+1, 1:NL+1)**T contains the right singular 00102 * vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains 00103 * the right singular vectors of the lower block. On exit 00104 * VT**T contains the right singular vectors of the 00105 * bidiagonal matrix. 00106 * 00107 * LDVT (input) INTEGER 00108 * The leading dimension of the array VT. LDVT >= max( 1, M ). 00109 * 00110 * IDXQ (output) INTEGER array, dimension(N) 00111 * This contains the permutation which will reintegrate the 00112 * subproblem just solved back into sorted order, i.e. 00113 * D( IDXQ( I = 1, N ) ) will be in ascending order. 00114 * 00115 * IWORK (workspace) INTEGER array, dimension( 4 * N ) 00116 * 00117 * WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M ) 00118 * 00119 * INFO (output) INTEGER 00120 * = 0: successful exit. 00121 * < 0: if INFO = -i, the i-th argument had an illegal value. 00122 * > 0: if INFO = 1, a singular value did not converge 00123 * 00124 * Further Details 00125 * =============== 00126 * 00127 * Based on contributions by 00128 * Ming Gu and Huan Ren, Computer Science Division, University of 00129 * California at Berkeley, USA 00130 * 00131 * ===================================================================== 00132 * 00133 * .. Parameters .. 00134 * 00135 DOUBLE PRECISION ONE, ZERO 00136 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00137 * .. 00138 * .. Local Scalars .. 00139 INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2, 00140 $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2 00141 DOUBLE PRECISION ORGNRM 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC ABS, MAX 00148 * .. 00149 * .. Executable Statements .. 00150 * 00151 * Test the input parameters. 00152 * 00153 INFO = 0 00154 * 00155 IF( NL.LT.1 ) THEN 00156 INFO = -1 00157 ELSE IF( NR.LT.1 ) THEN 00158 INFO = -2 00159 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00160 INFO = -3 00161 END IF 00162 IF( INFO.NE.0 ) THEN 00163 CALL XERBLA( 'DLASD1', -INFO ) 00164 RETURN 00165 END IF 00166 * 00167 N = NL + NR + 1 00168 M = N + SQRE 00169 * 00170 * The following values are for bookkeeping purposes only. They are 00171 * integer pointers which indicate the portion of the workspace 00172 * used by a particular array in DLASD2 and DLASD3. 00173 * 00174 LDU2 = N 00175 LDVT2 = M 00176 * 00177 IZ = 1 00178 ISIGMA = IZ + M 00179 IU2 = ISIGMA + N 00180 IVT2 = IU2 + LDU2*N 00181 IQ = IVT2 + LDVT2*M 00182 * 00183 IDX = 1 00184 IDXC = IDX + N 00185 COLTYP = IDXC + N 00186 IDXP = COLTYP + N 00187 * 00188 * Scale. 00189 * 00190 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 00191 D( NL+1 ) = ZERO 00192 DO 10 I = 1, N 00193 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 00194 ORGNRM = ABS( D( I ) ) 00195 END IF 00196 10 CONTINUE 00197 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00198 ALPHA = ALPHA / ORGNRM 00199 BETA = BETA / ORGNRM 00200 * 00201 * Deflate singular values. 00202 * 00203 CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU, 00204 $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2, 00205 $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ), 00206 $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO ) 00207 * 00208 * Solve Secular Equation and update singular vectors. 00209 * 00210 LDQ = K 00211 CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ), 00212 $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ), 00213 $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ), 00214 $ INFO ) 00215 IF( INFO.NE.0 ) THEN 00216 RETURN 00217 END IF 00218 * 00219 * Unscale. 00220 * 00221 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00222 * 00223 * Prepare the IDXQ sorting permutation. 00224 * 00225 N1 = K 00226 N2 = N - K 00227 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 00228 * 00229 RETURN 00230 * 00231 * End of DLASD1 00232 * 00233 END
1.7.3 
