 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, 00002 $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, 00003 $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, 00004 $ IWORK, INFO ) 00005 * 00006 * -- LAPACK routine (version 3.2) -- 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 2006 00010 * 00011 * .. Scalar Arguments .. 00012 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS, 00013 $ SMLSIZ 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 00017 $ K( * ), PERM( LDGCOL, * ) 00018 REAL B( LDB, * ), BX( LDBX, * ), C( * ), 00019 $ DIFL( LDU, * ), DIFR( LDU, * ), 00020 $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ), 00021 $ U( LDU, * ), VT( LDU, * ), WORK( * ), 00022 $ Z( LDU, * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * SLALSA is an itermediate step in solving the least squares problem 00029 * by computing the SVD of the coefficient matrix in compact form (The 00030 * singular vectors are computed as products of simple orthorgonal 00031 * matrices.). 00032 * 00033 * If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector 00034 * matrix of an upper bidiagonal matrix to the right hand side; and if 00035 * ICOMPQ = 1, SLALSA applies the right singular vector matrix to the 00036 * right hand side. The singular vector matrices were generated in 00037 * compact form by SLALSA. 00038 * 00039 * Arguments 00040 * ========= 00041 * 00042 * 00043 * ICOMPQ (input) INTEGER 00044 * Specifies whether the left or the right singular vector 00045 * matrix is involved. 00046 * = 0: Left singular vector matrix 00047 * = 1: Right singular vector matrix 00048 * 00049 * SMLSIZ (input) INTEGER 00050 * The maximum size of the subproblems at the bottom of the 00051 * computation tree. 00052 * 00053 * N (input) INTEGER 00054 * The row and column dimensions of the upper bidiagonal matrix. 00055 * 00056 * NRHS (input) INTEGER 00057 * The number of columns of B and BX. NRHS must be at least 1. 00058 * 00059 * B (input/output) REAL array, dimension ( LDB, NRHS ) 00060 * On input, B contains the right hand sides of the least 00061 * squares problem in rows 1 through M. 00062 * On output, B contains the solution X in rows 1 through N. 00063 * 00064 * LDB (input) INTEGER 00065 * The leading dimension of B in the calling subprogram. 00066 * LDB must be at least max(1,MAX( M, N ) ). 00067 * 00068 * BX (output) REAL array, dimension ( LDBX, NRHS ) 00069 * On exit, the result of applying the left or right singular 00070 * vector matrix to B. 00071 * 00072 * LDBX (input) INTEGER 00073 * The leading dimension of BX. 00074 * 00075 * U (input) REAL array, dimension ( LDU, SMLSIZ ). 00076 * On entry, U contains the left singular vector matrices of all 00077 * subproblems at the bottom level. 00078 * 00079 * LDU (input) INTEGER, LDU = > N. 00080 * The leading dimension of arrays U, VT, DIFL, DIFR, 00081 * POLES, GIVNUM, and Z. 00082 * 00083 * VT (input) REAL array, dimension ( LDU, SMLSIZ+1 ). 00084 * On entry, VT**T contains the right singular vector matrices of 00085 * all subproblems at the bottom level. 00086 * 00087 * K (input) INTEGER array, dimension ( N ). 00088 * 00089 * DIFL (input) REAL array, dimension ( LDU, NLVL ). 00090 * where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. 00091 * 00092 * DIFR (input) REAL array, dimension ( LDU, 2 * NLVL ). 00093 * On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record 00094 * distances between singular values on the I-th level and 00095 * singular values on the (I -1)-th level, and DIFR(*, 2 * I) 00096 * record the normalizing factors of the right singular vectors 00097 * matrices of subproblems on I-th level. 00098 * 00099 * Z (input) REAL array, dimension ( LDU, NLVL ). 00100 * On entry, Z(1, I) contains the components of the deflation- 00101 * adjusted updating row vector for subproblems on the I-th 00102 * level. 00103 * 00104 * POLES (input) REAL array, dimension ( LDU, 2 * NLVL ). 00105 * On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old 00106 * singular values involved in the secular equations on the I-th 00107 * level. 00108 * 00109 * GIVPTR (input) INTEGER array, dimension ( N ). 00110 * On entry, GIVPTR( I ) records the number of Givens 00111 * rotations performed on the I-th problem on the computation 00112 * tree. 00113 * 00114 * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ). 00115 * On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the 00116 * locations of Givens rotations performed on the I-th level on 00117 * the computation tree. 00118 * 00119 * LDGCOL (input) INTEGER, LDGCOL = > N. 00120 * The leading dimension of arrays GIVCOL and PERM. 00121 * 00122 * PERM (input) INTEGER array, dimension ( LDGCOL, NLVL ). 00123 * On entry, PERM(*, I) records permutations done on the I-th 00124 * level of the computation tree. 00125 * 00126 * GIVNUM (input) REAL array, dimension ( LDU, 2 * NLVL ). 00127 * On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S- 00128 * values of Givens rotations performed on the I-th level on the 00129 * computation tree. 00130 * 00131 * C (input) REAL array, dimension ( N ). 00132 * On entry, if the I-th subproblem is not square, 00133 * C( I ) contains the C-value of a Givens rotation related to 00134 * the right null space of the I-th subproblem. 00135 * 00136 * S (input) REAL array, dimension ( N ). 00137 * On entry, if the I-th subproblem is not square, 00138 * S( I ) contains the S-value of a Givens rotation related to 00139 * the right null space of the I-th subproblem. 00140 * 00141 * WORK (workspace) REAL array. 00142 * The dimension must be at least N. 00143 * 00144 * IWORK (workspace) INTEGER array. 00145 * The dimension must be at least 3 * N 00146 * 00147 * INFO (output) INTEGER 00148 * = 0: successful exit. 00149 * < 0: if INFO = -i, the i-th argument had an illegal value. 00150 * 00151 * Further Details 00152 * =============== 00153 * 00154 * Based on contributions by 00155 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00156 * California at Berkeley, USA 00157 * Osni Marques, LBNL/NERSC, USA 00158 * 00159 * ===================================================================== 00160 * 00161 * .. Parameters .. 00162 REAL ZERO, ONE 00163 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00164 * .. 00165 * .. Local Scalars .. 00166 INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2, 00167 $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL, 00168 $ NR, NRF, NRP1, SQRE 00169 * .. 00170 * .. External Subroutines .. 00171 EXTERNAL SCOPY, SGEMM, SLALS0, SLASDT, XERBLA 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Test the input parameters. 00176 * 00177 INFO = 0 00178 * 00179 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00180 INFO = -1 00181 ELSE IF( SMLSIZ.LT.3 ) THEN 00182 INFO = -2 00183 ELSE IF( N.LT.SMLSIZ ) THEN 00184 INFO = -3 00185 ELSE IF( NRHS.LT.1 ) THEN 00186 INFO = -4 00187 ELSE IF( LDB.LT.N ) THEN 00188 INFO = -6 00189 ELSE IF( LDBX.LT.N ) THEN 00190 INFO = -8 00191 ELSE IF( LDU.LT.N ) THEN 00192 INFO = -10 00193 ELSE IF( LDGCOL.LT.N ) THEN 00194 INFO = -19 00195 END IF 00196 IF( INFO.NE.0 ) THEN 00197 CALL XERBLA( 'SLALSA', -INFO ) 00198 RETURN 00199 END IF 00200 * 00201 * Book-keeping and setting up the computation tree. 00202 * 00203 INODE = 1 00204 NDIML = INODE + N 00205 NDIMR = NDIML + N 00206 * 00207 CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00208 $ IWORK( NDIMR ), SMLSIZ ) 00209 * 00210 * The following code applies back the left singular vector factors. 00211 * For applying back the right singular vector factors, go to 50. 00212 * 00213 IF( ICOMPQ.EQ.1 ) THEN 00214 GO TO 50 00215 END IF 00216 * 00217 * The nodes on the bottom level of the tree were solved 00218 * by SLASDQ. The corresponding left and right singular vector 00219 * matrices are in explicit form. First apply back the left 00220 * singular vector matrices. 00221 * 00222 NDB1 = ( ND+1 ) / 2 00223 DO 10 I = NDB1, ND 00224 * 00225 * IC : center row of each node 00226 * NL : number of rows of left subproblem 00227 * NR : number of rows of right subproblem 00228 * NLF: starting row of the left subproblem 00229 * NRF: starting row of the right subproblem 00230 * 00231 I1 = I - 1 00232 IC = IWORK( INODE+I1 ) 00233 NL = IWORK( NDIML+I1 ) 00234 NR = IWORK( NDIMR+I1 ) 00235 NLF = IC - NL 00236 NRF = IC + 1 00237 CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU, 00238 $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 00239 CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU, 00240 $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 00241 10 CONTINUE 00242 * 00243 * Next copy the rows of B that correspond to unchanged rows 00244 * in the bidiagonal matrix to BX. 00245 * 00246 DO 20 I = 1, ND 00247 IC = IWORK( INODE+I-1 ) 00248 CALL SCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX ) 00249 20 CONTINUE 00250 * 00251 * Finally go through the left singular vector matrices of all 00252 * the other subproblems bottom-up on the tree. 00253 * 00254 J = 2**NLVL 00255 SQRE = 0 00256 * 00257 DO 40 LVL = NLVL, 1, -1 00258 LVL2 = 2*LVL - 1 00259 * 00260 * find the first node LF and last node LL on 00261 * the current level LVL 00262 * 00263 IF( LVL.EQ.1 ) THEN 00264 LF = 1 00265 LL = 1 00266 ELSE 00267 LF = 2**( LVL-1 ) 00268 LL = 2*LF - 1 00269 END IF 00270 DO 30 I = LF, LL 00271 IM1 = I - 1 00272 IC = IWORK( INODE+IM1 ) 00273 NL = IWORK( NDIML+IM1 ) 00274 NR = IWORK( NDIMR+IM1 ) 00275 NLF = IC - NL 00276 NRF = IC + 1 00277 J = J - 1 00278 CALL SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX, 00279 $ B( NLF, 1 ), LDB, PERM( NLF, LVL ), 00280 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00281 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 00282 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 00283 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK, 00284 $ INFO ) 00285 30 CONTINUE 00286 40 CONTINUE 00287 GO TO 90 00288 * 00289 * ICOMPQ = 1: applying back the right singular vector factors. 00290 * 00291 50 CONTINUE 00292 * 00293 * First now go through the right singular vector matrices of all 00294 * the tree nodes top-down. 00295 * 00296 J = 0 00297 DO 70 LVL = 1, NLVL 00298 LVL2 = 2*LVL - 1 00299 * 00300 * Find the first node LF and last node LL on 00301 * the current level LVL. 00302 * 00303 IF( LVL.EQ.1 ) THEN 00304 LF = 1 00305 LL = 1 00306 ELSE 00307 LF = 2**( LVL-1 ) 00308 LL = 2*LF - 1 00309 END IF 00310 DO 60 I = LL, LF, -1 00311 IM1 = I - 1 00312 IC = IWORK( INODE+IM1 ) 00313 NL = IWORK( NDIML+IM1 ) 00314 NR = IWORK( NDIMR+IM1 ) 00315 NLF = IC - NL 00316 NRF = IC + 1 00317 IF( I.EQ.LL ) THEN 00318 SQRE = 0 00319 ELSE 00320 SQRE = 1 00321 END IF 00322 J = J + 1 00323 CALL SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB, 00324 $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ), 00325 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00326 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 00327 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 00328 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK, 00329 $ INFO ) 00330 60 CONTINUE 00331 70 CONTINUE 00332 * 00333 * The nodes on the bottom level of the tree were solved 00334 * by SLASDQ. The corresponding right singular vector 00335 * matrices are in explicit form. Apply them back. 00336 * 00337 NDB1 = ( ND+1 ) / 2 00338 DO 80 I = NDB1, ND 00339 I1 = I - 1 00340 IC = IWORK( INODE+I1 ) 00341 NL = IWORK( NDIML+I1 ) 00342 NR = IWORK( NDIMR+I1 ) 00343 NLP1 = NL + 1 00344 IF( I.EQ.ND ) THEN 00345 NRP1 = NR 00346 ELSE 00347 NRP1 = NR + 1 00348 END IF 00349 NLF = IC - NL 00350 NRF = IC + 1 00351 CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU, 00352 $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 00353 CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU, 00354 $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 00355 80 CONTINUE 00356 * 00357 90 CONTINUE 00358 * 00359 RETURN 00360 * 00361 * End of SLALSA 00362 * 00363 END
1.7.3 
