![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLALSA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLALSA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsa.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsa.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsa.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, 00022 * LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, 00023 * GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, 00024 * IWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS, 00028 * $ SMLSIZ 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 00032 * $ K( * ), PERM( LDGCOL, * ) 00033 * DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ), 00034 * $ DIFL( LDU, * ), DIFR( LDU, * ), 00035 * $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ), 00036 * $ U( LDU, * ), VT( LDU, * ), WORK( * ), 00037 * $ Z( LDU, * ) 00038 * .. 00039 * 00040 * 00041 *> \par Purpose: 00042 * ============= 00043 *> 00044 *> \verbatim 00045 *> 00046 *> DLALSA is an itermediate step in solving the least squares problem 00047 *> by computing the SVD of the coefficient matrix in compact form (The 00048 *> singular vectors are computed as products of simple orthorgonal 00049 *> matrices.). 00050 *> 00051 *> If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector 00052 *> matrix of an upper bidiagonal matrix to the right hand side; and if 00053 *> ICOMPQ = 1, DLALSA applies the right singular vector matrix to the 00054 *> right hand side. The singular vector matrices were generated in 00055 *> compact form by DLALSA. 00056 *> \endverbatim 00057 * 00058 * Arguments: 00059 * ========== 00060 * 00061 *> \param[in] ICOMPQ 00062 *> \verbatim 00063 *> ICOMPQ is INTEGER 00064 *> Specifies whether the left or the right singular vector 00065 *> matrix is involved. 00066 *> = 0: Left singular vector matrix 00067 *> = 1: Right singular vector matrix 00068 *> \endverbatim 00069 *> 00070 *> \param[in] SMLSIZ 00071 *> \verbatim 00072 *> SMLSIZ is INTEGER 00073 *> The maximum size of the subproblems at the bottom of the 00074 *> computation tree. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] N 00078 *> \verbatim 00079 *> N is INTEGER 00080 *> The row and column dimensions of the upper bidiagonal matrix. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] NRHS 00084 *> \verbatim 00085 *> NRHS is INTEGER 00086 *> The number of columns of B and BX. NRHS must be at least 1. 00087 *> \endverbatim 00088 *> 00089 *> \param[in,out] B 00090 *> \verbatim 00091 *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS ) 00092 *> On input, B contains the right hand sides of the least 00093 *> squares problem in rows 1 through M. 00094 *> On output, B contains the solution X in rows 1 through N. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] LDB 00098 *> \verbatim 00099 *> LDB is INTEGER 00100 *> The leading dimension of B in the calling subprogram. 00101 *> LDB must be at least max(1,MAX( M, N ) ). 00102 *> \endverbatim 00103 *> 00104 *> \param[out] BX 00105 *> \verbatim 00106 *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS ) 00107 *> On exit, the result of applying the left or right singular 00108 *> vector matrix to B. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LDBX 00112 *> \verbatim 00113 *> LDBX is INTEGER 00114 *> The leading dimension of BX. 00115 *> \endverbatim 00116 *> 00117 *> \param[in] U 00118 *> \verbatim 00119 *> U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ). 00120 *> On entry, U contains the left singular vector matrices of all 00121 *> subproblems at the bottom level. 00122 *> \endverbatim 00123 *> 00124 *> \param[in] LDU 00125 *> \verbatim 00126 *> LDU is INTEGER, LDU = > N. 00127 *> The leading dimension of arrays U, VT, DIFL, DIFR, 00128 *> POLES, GIVNUM, and Z. 00129 *> \endverbatim 00130 *> 00131 *> \param[in] VT 00132 *> \verbatim 00133 *> VT is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ). 00134 *> On entry, VT**T contains the right singular vector matrices of 00135 *> all subproblems at the bottom level. 00136 *> \endverbatim 00137 *> 00138 *> \param[in] K 00139 *> \verbatim 00140 *> K is INTEGER array, dimension ( N ). 00141 *> \endverbatim 00142 *> 00143 *> \param[in] DIFL 00144 *> \verbatim 00145 *> DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ). 00146 *> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. 00147 *> \endverbatim 00148 *> 00149 *> \param[in] DIFR 00150 *> \verbatim 00151 *> DIFR is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ). 00152 *> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record 00153 *> distances between singular values on the I-th level and 00154 *> singular values on the (I -1)-th level, and DIFR(*, 2 * I) 00155 *> record the normalizing factors of the right singular vectors 00156 *> matrices of subproblems on I-th level. 00157 *> \endverbatim 00158 *> 00159 *> \param[in] Z 00160 *> \verbatim 00161 *> Z is DOUBLE PRECISION array, dimension ( LDU, NLVL ). 00162 *> On entry, Z(1, I) contains the components of the deflation- 00163 *> adjusted updating row vector for subproblems on the I-th 00164 *> level. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] POLES 00168 *> \verbatim 00169 *> POLES is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ). 00170 *> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old 00171 *> singular values involved in the secular equations on the I-th 00172 *> level. 00173 *> \endverbatim 00174 *> 00175 *> \param[in] GIVPTR 00176 *> \verbatim 00177 *> GIVPTR is INTEGER array, dimension ( N ). 00178 *> On entry, GIVPTR( I ) records the number of Givens 00179 *> rotations performed on the I-th problem on the computation 00180 *> tree. 00181 *> \endverbatim 00182 *> 00183 *> \param[in] GIVCOL 00184 *> \verbatim 00185 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ). 00186 *> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the 00187 *> locations of Givens rotations performed on the I-th level on 00188 *> the computation tree. 00189 *> \endverbatim 00190 *> 00191 *> \param[in] LDGCOL 00192 *> \verbatim 00193 *> LDGCOL is INTEGER, LDGCOL = > N. 00194 *> The leading dimension of arrays GIVCOL and PERM. 00195 *> \endverbatim 00196 *> 00197 *> \param[in] PERM 00198 *> \verbatim 00199 *> PERM is INTEGER array, dimension ( LDGCOL, NLVL ). 00200 *> On entry, PERM(*, I) records permutations done on the I-th 00201 *> level of the computation tree. 00202 *> \endverbatim 00203 *> 00204 *> \param[in] GIVNUM 00205 *> \verbatim 00206 *> GIVNUM is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ). 00207 *> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S- 00208 *> values of Givens rotations performed on the I-th level on the 00209 *> computation tree. 00210 *> \endverbatim 00211 *> 00212 *> \param[in] C 00213 *> \verbatim 00214 *> C is DOUBLE PRECISION array, dimension ( N ). 00215 *> On entry, if the I-th subproblem is not square, 00216 *> C( I ) contains the C-value of a Givens rotation related to 00217 *> the right null space of the I-th subproblem. 00218 *> \endverbatim 00219 *> 00220 *> \param[in] S 00221 *> \verbatim 00222 *> S is DOUBLE PRECISION array, dimension ( N ). 00223 *> On entry, if the I-th subproblem is not square, 00224 *> S( I ) contains the S-value of a Givens rotation related to 00225 *> the right null space of the I-th subproblem. 00226 *> \endverbatim 00227 *> 00228 *> \param[out] WORK 00229 *> \verbatim 00230 *> WORK is DOUBLE PRECISION array. 00231 *> The dimension must be at least N. 00232 *> \endverbatim 00233 *> 00234 *> \param[out] IWORK 00235 *> \verbatim 00236 *> IWORK is INTEGER array. 00237 *> The dimension must be at least 3 * N 00238 *> \endverbatim 00239 *> 00240 *> \param[out] INFO 00241 *> \verbatim 00242 *> INFO is INTEGER 00243 *> = 0: successful exit. 00244 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00245 *> \endverbatim 00246 * 00247 * Authors: 00248 * ======== 00249 * 00250 *> \author Univ. of Tennessee 00251 *> \author Univ. of California Berkeley 00252 *> \author Univ. of Colorado Denver 00253 *> \author NAG Ltd. 00254 * 00255 *> \date November 2011 00256 * 00257 *> \ingroup doubleOTHERcomputational 00258 * 00259 *> \par Contributors: 00260 * ================== 00261 *> 00262 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of 00263 *> California at Berkeley, USA \n 00264 *> Osni Marques, LBNL/NERSC, USA \n 00265 * 00266 * ===================================================================== 00267 SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, 00268 $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, 00269 $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, 00270 $ IWORK, INFO ) 00271 * 00272 * -- LAPACK computational routine (version 3.4.0) -- 00273 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00274 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00275 * November 2011 00276 * 00277 * .. Scalar Arguments .. 00278 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS, 00279 $ SMLSIZ 00280 * .. 00281 * .. Array Arguments .. 00282 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 00283 $ K( * ), PERM( LDGCOL, * ) 00284 DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ), 00285 $ DIFL( LDU, * ), DIFR( LDU, * ), 00286 $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ), 00287 $ U( LDU, * ), VT( LDU, * ), WORK( * ), 00288 $ Z( LDU, * ) 00289 * .. 00290 * 00291 * ===================================================================== 00292 * 00293 * .. Parameters .. 00294 DOUBLE PRECISION ZERO, ONE 00295 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00296 * .. 00297 * .. Local Scalars .. 00298 INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2, 00299 $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL, 00300 $ NR, NRF, NRP1, SQRE 00301 * .. 00302 * .. External Subroutines .. 00303 EXTERNAL DCOPY, DGEMM, DLALS0, DLASDT, XERBLA 00304 * .. 00305 * .. Executable Statements .. 00306 * 00307 * Test the input parameters. 00308 * 00309 INFO = 0 00310 * 00311 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00312 INFO = -1 00313 ELSE IF( SMLSIZ.LT.3 ) THEN 00314 INFO = -2 00315 ELSE IF( N.LT.SMLSIZ ) THEN 00316 INFO = -3 00317 ELSE IF( NRHS.LT.1 ) THEN 00318 INFO = -4 00319 ELSE IF( LDB.LT.N ) THEN 00320 INFO = -6 00321 ELSE IF( LDBX.LT.N ) THEN 00322 INFO = -8 00323 ELSE IF( LDU.LT.N ) THEN 00324 INFO = -10 00325 ELSE IF( LDGCOL.LT.N ) THEN 00326 INFO = -19 00327 END IF 00328 IF( INFO.NE.0 ) THEN 00329 CALL XERBLA( 'DLALSA', -INFO ) 00330 RETURN 00331 END IF 00332 * 00333 * Book-keeping and setting up the computation tree. 00334 * 00335 INODE = 1 00336 NDIML = INODE + N 00337 NDIMR = NDIML + N 00338 * 00339 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00340 $ IWORK( NDIMR ), SMLSIZ ) 00341 * 00342 * The following code applies back the left singular vector factors. 00343 * For applying back the right singular vector factors, go to 50. 00344 * 00345 IF( ICOMPQ.EQ.1 ) THEN 00346 GO TO 50 00347 END IF 00348 * 00349 * The nodes on the bottom level of the tree were solved 00350 * by DLASDQ. The corresponding left and right singular vector 00351 * matrices are in explicit form. First apply back the left 00352 * singular vector matrices. 00353 * 00354 NDB1 = ( ND+1 ) / 2 00355 DO 10 I = NDB1, ND 00356 * 00357 * IC : center row of each node 00358 * NL : number of rows of left subproblem 00359 * NR : number of rows of right subproblem 00360 * NLF: starting row of the left subproblem 00361 * NRF: starting row of the right subproblem 00362 * 00363 I1 = I - 1 00364 IC = IWORK( INODE+I1 ) 00365 NL = IWORK( NDIML+I1 ) 00366 NR = IWORK( NDIMR+I1 ) 00367 NLF = IC - NL 00368 NRF = IC + 1 00369 CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU, 00370 $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 00371 CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU, 00372 $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 00373 10 CONTINUE 00374 * 00375 * Next copy the rows of B that correspond to unchanged rows 00376 * in the bidiagonal matrix to BX. 00377 * 00378 DO 20 I = 1, ND 00379 IC = IWORK( INODE+I-1 ) 00380 CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX ) 00381 20 CONTINUE 00382 * 00383 * Finally go through the left singular vector matrices of all 00384 * the other subproblems bottom-up on the tree. 00385 * 00386 J = 2**NLVL 00387 SQRE = 0 00388 * 00389 DO 40 LVL = NLVL, 1, -1 00390 LVL2 = 2*LVL - 1 00391 * 00392 * find the first node LF and last node LL on 00393 * the current level LVL 00394 * 00395 IF( LVL.EQ.1 ) THEN 00396 LF = 1 00397 LL = 1 00398 ELSE 00399 LF = 2**( LVL-1 ) 00400 LL = 2*LF - 1 00401 END IF 00402 DO 30 I = LF, LL 00403 IM1 = I - 1 00404 IC = IWORK( INODE+IM1 ) 00405 NL = IWORK( NDIML+IM1 ) 00406 NR = IWORK( NDIMR+IM1 ) 00407 NLF = IC - NL 00408 NRF = IC + 1 00409 J = J - 1 00410 CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX, 00411 $ B( NLF, 1 ), LDB, PERM( NLF, LVL ), 00412 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00413 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 00414 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 00415 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK, 00416 $ INFO ) 00417 30 CONTINUE 00418 40 CONTINUE 00419 GO TO 90 00420 * 00421 * ICOMPQ = 1: applying back the right singular vector factors. 00422 * 00423 50 CONTINUE 00424 * 00425 * First now go through the right singular vector matrices of all 00426 * the tree nodes top-down. 00427 * 00428 J = 0 00429 DO 70 LVL = 1, NLVL 00430 LVL2 = 2*LVL - 1 00431 * 00432 * Find the first node LF and last node LL on 00433 * the current level LVL. 00434 * 00435 IF( LVL.EQ.1 ) THEN 00436 LF = 1 00437 LL = 1 00438 ELSE 00439 LF = 2**( LVL-1 ) 00440 LL = 2*LF - 1 00441 END IF 00442 DO 60 I = LL, LF, -1 00443 IM1 = I - 1 00444 IC = IWORK( INODE+IM1 ) 00445 NL = IWORK( NDIML+IM1 ) 00446 NR = IWORK( NDIMR+IM1 ) 00447 NLF = IC - NL 00448 NRF = IC + 1 00449 IF( I.EQ.LL ) THEN 00450 SQRE = 0 00451 ELSE 00452 SQRE = 1 00453 END IF 00454 J = J + 1 00455 CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB, 00456 $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ), 00457 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00458 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 00459 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 00460 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK, 00461 $ INFO ) 00462 60 CONTINUE 00463 70 CONTINUE 00464 * 00465 * The nodes on the bottom level of the tree were solved 00466 * by DLASDQ. The corresponding right singular vector 00467 * matrices are in explicit form. Apply them back. 00468 * 00469 NDB1 = ( ND+1 ) / 2 00470 DO 80 I = NDB1, ND 00471 I1 = I - 1 00472 IC = IWORK( INODE+I1 ) 00473 NL = IWORK( NDIML+I1 ) 00474 NR = IWORK( NDIMR+I1 ) 00475 NLP1 = NL + 1 00476 IF( I.EQ.ND ) THEN 00477 NRP1 = NR 00478 ELSE 00479 NRP1 = NR + 1 00480 END IF 00481 NLF = IC - NL 00482 NRF = IC + 1 00483 CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU, 00484 $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 00485 CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU, 00486 $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 00487 80 CONTINUE 00488 * 00489 90 CONTINUE 00490 * 00491 RETURN 00492 * 00493 * End of DLALSA 00494 * 00495 END