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