![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLALSA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLALSA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clalsa.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clalsa.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clalsa.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLALSA( 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, RWORK, 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 * REAL C( * ), DIFL( LDU, * ), DIFR( LDU, * ), 00034 * $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ), 00035 * $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * ) 00036 * COMPLEX B( LDB, * ), BX( LDBX, * ) 00037 * .. 00038 * 00039 * 00040 *> \par Purpose: 00041 * ============= 00042 *> 00043 *> \verbatim 00044 *> 00045 *> CLALSA is an itermediate step in solving the least squares problem 00046 *> by computing the SVD of the coefficient matrix in compact form (The 00047 *> singular vectors are computed as products of simple orthorgonal 00048 *> matrices.). 00049 *> 00050 *> If ICOMPQ = 0, CLALSA applies the inverse of the left singular vector 00051 *> matrix of an upper bidiagonal matrix to the right hand side; and if 00052 *> ICOMPQ = 1, CLALSA applies the right singular vector matrix to the 00053 *> right hand side. The singular vector matrices were generated in 00054 *> compact form by CLALSA. 00055 *> \endverbatim 00056 * 00057 * Arguments: 00058 * ========== 00059 * 00060 *> \param[in] ICOMPQ 00061 *> \verbatim 00062 *> ICOMPQ is INTEGER 00063 *> Specifies whether the left or the right singular vector 00064 *> matrix is involved. 00065 *> = 0: Left singular vector matrix 00066 *> = 1: Right singular vector matrix 00067 *> \endverbatim 00068 *> 00069 *> \param[in] SMLSIZ 00070 *> \verbatim 00071 *> SMLSIZ is INTEGER 00072 *> The maximum size of the subproblems at the bottom of the 00073 *> computation tree. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] N 00077 *> \verbatim 00078 *> N is INTEGER 00079 *> The row and column dimensions of the upper bidiagonal matrix. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] NRHS 00083 *> \verbatim 00084 *> NRHS is INTEGER 00085 *> The number of columns of B and BX. NRHS must be at least 1. 00086 *> \endverbatim 00087 *> 00088 *> \param[in,out] B 00089 *> \verbatim 00090 *> B is COMPLEX array, dimension ( LDB, NRHS ) 00091 *> On input, B contains the right hand sides of the least 00092 *> squares problem in rows 1 through M. 00093 *> On output, B contains the solution X in rows 1 through N. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LDB 00097 *> \verbatim 00098 *> LDB is INTEGER 00099 *> The leading dimension of B in the calling subprogram. 00100 *> LDB must be at least max(1,MAX( M, N ) ). 00101 *> \endverbatim 00102 *> 00103 *> \param[out] BX 00104 *> \verbatim 00105 *> BX is COMPLEX array, dimension ( LDBX, NRHS ) 00106 *> On exit, the result of applying the left or right singular 00107 *> vector matrix to B. 00108 *> \endverbatim 00109 *> 00110 *> \param[in] LDBX 00111 *> \verbatim 00112 *> LDBX is INTEGER 00113 *> The leading dimension of BX. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] U 00117 *> \verbatim 00118 *> U is REAL array, dimension ( LDU, SMLSIZ ). 00119 *> On entry, U contains the left singular vector matrices of all 00120 *> subproblems at the bottom level. 00121 *> \endverbatim 00122 *> 00123 *> \param[in] LDU 00124 *> \verbatim 00125 *> LDU is INTEGER, LDU = > N. 00126 *> The leading dimension of arrays U, VT, DIFL, DIFR, 00127 *> POLES, GIVNUM, and Z. 00128 *> \endverbatim 00129 *> 00130 *> \param[in] VT 00131 *> \verbatim 00132 *> VT is REAL array, dimension ( LDU, SMLSIZ+1 ). 00133 *> On entry, VT**H contains the right singular vector matrices of 00134 *> all subproblems at the bottom level. 00135 *> \endverbatim 00136 *> 00137 *> \param[in] K 00138 *> \verbatim 00139 *> K is INTEGER array, dimension ( N ). 00140 *> \endverbatim 00141 *> 00142 *> \param[in] DIFL 00143 *> \verbatim 00144 *> DIFL is REAL array, dimension ( LDU, NLVL ). 00145 *> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. 00146 *> \endverbatim 00147 *> 00148 *> \param[in] DIFR 00149 *> \verbatim 00150 *> DIFR is REAL array, dimension ( LDU, 2 * NLVL ). 00151 *> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record 00152 *> distances between singular values on the I-th level and 00153 *> singular values on the (I -1)-th level, and DIFR(*, 2 * I) 00154 *> record the normalizing factors of the right singular vectors 00155 *> matrices of subproblems on I-th level. 00156 *> \endverbatim 00157 *> 00158 *> \param[in] Z 00159 *> \verbatim 00160 *> Z is REAL array, dimension ( LDU, NLVL ). 00161 *> On entry, Z(1, I) contains the components of the deflation- 00162 *> adjusted updating row vector for subproblems on the I-th 00163 *> level. 00164 *> \endverbatim 00165 *> 00166 *> \param[in] POLES 00167 *> \verbatim 00168 *> POLES is REAL array, dimension ( LDU, 2 * NLVL ). 00169 *> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old 00170 *> singular values involved in the secular equations on the I-th 00171 *> level. 00172 *> \endverbatim 00173 *> 00174 *> \param[in] GIVPTR 00175 *> \verbatim 00176 *> GIVPTR is INTEGER array, dimension ( N ). 00177 *> On entry, GIVPTR( I ) records the number of Givens 00178 *> rotations performed on the I-th problem on the computation 00179 *> tree. 00180 *> \endverbatim 00181 *> 00182 *> \param[in] GIVCOL 00183 *> \verbatim 00184 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ). 00185 *> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the 00186 *> locations of Givens rotations performed on the I-th level on 00187 *> the computation tree. 00188 *> \endverbatim 00189 *> 00190 *> \param[in] LDGCOL 00191 *> \verbatim 00192 *> LDGCOL is INTEGER, LDGCOL = > N. 00193 *> The leading dimension of arrays GIVCOL and PERM. 00194 *> \endverbatim 00195 *> 00196 *> \param[in] PERM 00197 *> \verbatim 00198 *> PERM is INTEGER array, dimension ( LDGCOL, NLVL ). 00199 *> On entry, PERM(*, I) records permutations done on the I-th 00200 *> level of the computation tree. 00201 *> \endverbatim 00202 *> 00203 *> \param[in] GIVNUM 00204 *> \verbatim 00205 *> GIVNUM is REAL array, dimension ( LDU, 2 * NLVL ). 00206 *> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S- 00207 *> values of Givens rotations performed on the I-th level on the 00208 *> computation tree. 00209 *> \endverbatim 00210 *> 00211 *> \param[in] C 00212 *> \verbatim 00213 *> C is REAL array, dimension ( N ). 00214 *> On entry, if the I-th subproblem is not square, 00215 *> C( I ) contains the C-value of a Givens rotation related to 00216 *> the right null space of the I-th subproblem. 00217 *> \endverbatim 00218 *> 00219 *> \param[in] S 00220 *> \verbatim 00221 *> S is REAL array, dimension ( N ). 00222 *> On entry, if the I-th subproblem is not square, 00223 *> S( I ) contains the S-value of a Givens rotation related to 00224 *> the right null space of the I-th subproblem. 00225 *> \endverbatim 00226 *> 00227 *> \param[out] RWORK 00228 *> \verbatim 00229 *> RWORK is REAL array, dimension at least 00230 *> MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ). 00231 *> \endverbatim 00232 *> 00233 *> \param[out] IWORK 00234 *> \verbatim 00235 *> IWORK is INTEGER array. 00236 *> The dimension must be at least 3 * N 00237 *> \endverbatim 00238 *> 00239 *> \param[out] INFO 00240 *> \verbatim 00241 *> INFO is INTEGER 00242 *> = 0: successful exit. 00243 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00244 *> \endverbatim 00245 * 00246 * Authors: 00247 * ======== 00248 * 00249 *> \author Univ. of Tennessee 00250 *> \author Univ. of California Berkeley 00251 *> \author Univ. of Colorado Denver 00252 *> \author NAG Ltd. 00253 * 00254 *> \date November 2011 00255 * 00256 *> \ingroup complexOTHERcomputational 00257 * 00258 *> \par Contributors: 00259 * ================== 00260 *> 00261 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of 00262 *> California at Berkeley, USA \n 00263 *> Osni Marques, LBNL/NERSC, USA \n 00264 * 00265 * ===================================================================== 00266 SUBROUTINE CLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, 00267 $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, 00268 $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK, 00269 $ IWORK, INFO ) 00270 * 00271 * -- LAPACK computational routine (version 3.4.0) -- 00272 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00274 * November 2011 00275 * 00276 * .. Scalar Arguments .. 00277 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS, 00278 $ SMLSIZ 00279 * .. 00280 * .. Array Arguments .. 00281 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 00282 $ K( * ), PERM( LDGCOL, * ) 00283 REAL C( * ), DIFL( LDU, * ), DIFR( LDU, * ), 00284 $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ), 00285 $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * ) 00286 COMPLEX B( LDB, * ), BX( LDBX, * ) 00287 * .. 00288 * 00289 * ===================================================================== 00290 * 00291 * .. Parameters .. 00292 REAL ZERO, ONE 00293 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00294 * .. 00295 * .. Local Scalars .. 00296 INTEGER I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL, 00297 $ JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML, 00298 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE 00299 * .. 00300 * .. External Subroutines .. 00301 EXTERNAL CCOPY, CLALS0, SGEMM, SLASDT, XERBLA 00302 * .. 00303 * .. Intrinsic Functions .. 00304 INTRINSIC AIMAG, CMPLX, REAL 00305 * .. 00306 * .. Executable Statements .. 00307 * 00308 * Test the input parameters. 00309 * 00310 INFO = 0 00311 * 00312 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00313 INFO = -1 00314 ELSE IF( SMLSIZ.LT.3 ) THEN 00315 INFO = -2 00316 ELSE IF( N.LT.SMLSIZ ) THEN 00317 INFO = -3 00318 ELSE IF( NRHS.LT.1 ) THEN 00319 INFO = -4 00320 ELSE IF( LDB.LT.N ) THEN 00321 INFO = -6 00322 ELSE IF( LDBX.LT.N ) THEN 00323 INFO = -8 00324 ELSE IF( LDU.LT.N ) THEN 00325 INFO = -10 00326 ELSE IF( LDGCOL.LT.N ) THEN 00327 INFO = -19 00328 END IF 00329 IF( INFO.NE.0 ) THEN 00330 CALL XERBLA( 'CLALSA', -INFO ) 00331 RETURN 00332 END IF 00333 * 00334 * Book-keeping and setting up the computation tree. 00335 * 00336 INODE = 1 00337 NDIML = INODE + N 00338 NDIMR = NDIML + N 00339 * 00340 CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00341 $ IWORK( NDIMR ), SMLSIZ ) 00342 * 00343 * The following code applies back the left singular vector factors. 00344 * For applying back the right singular vector factors, go to 170. 00345 * 00346 IF( ICOMPQ.EQ.1 ) THEN 00347 GO TO 170 00348 END IF 00349 * 00350 * The nodes on the bottom level of the tree were solved 00351 * by SLASDQ. The corresponding left and right singular vector 00352 * matrices are in explicit form. First apply back the left 00353 * singular vector matrices. 00354 * 00355 NDB1 = ( ND+1 ) / 2 00356 DO 130 I = NDB1, ND 00357 * 00358 * IC : center row of each node 00359 * NL : number of rows of left subproblem 00360 * NR : number of rows of right subproblem 00361 * NLF: starting row of the left subproblem 00362 * NRF: starting row of the right subproblem 00363 * 00364 I1 = I - 1 00365 IC = IWORK( INODE+I1 ) 00366 NL = IWORK( NDIML+I1 ) 00367 NR = IWORK( NDIMR+I1 ) 00368 NLF = IC - NL 00369 NRF = IC + 1 00370 * 00371 * Since B and BX are complex, the following call to SGEMM 00372 * is performed in two steps (real and imaginary parts). 00373 * 00374 * CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU, 00375 * $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 00376 * 00377 J = NL*NRHS*2 00378 DO 20 JCOL = 1, NRHS 00379 DO 10 JROW = NLF, NLF + NL - 1 00380 J = J + 1 00381 RWORK( J ) = REAL( B( JROW, JCOL ) ) 00382 10 CONTINUE 00383 20 CONTINUE 00384 CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU, 00385 $ RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1 ), NL ) 00386 J = NL*NRHS*2 00387 DO 40 JCOL = 1, NRHS 00388 DO 30 JROW = NLF, NLF + NL - 1 00389 J = J + 1 00390 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 00391 30 CONTINUE 00392 40 CONTINUE 00393 CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU, 00394 $ RWORK( 1+NL*NRHS*2 ), NL, ZERO, RWORK( 1+NL*NRHS ), 00395 $ NL ) 00396 JREAL = 0 00397 JIMAG = NL*NRHS 00398 DO 60 JCOL = 1, NRHS 00399 DO 50 JROW = NLF, NLF + NL - 1 00400 JREAL = JREAL + 1 00401 JIMAG = JIMAG + 1 00402 BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 00403 $ RWORK( JIMAG ) ) 00404 50 CONTINUE 00405 60 CONTINUE 00406 * 00407 * Since B and BX are complex, the following call to SGEMM 00408 * is performed in two steps (real and imaginary parts). 00409 * 00410 * CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU, 00411 * $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 00412 * 00413 J = NR*NRHS*2 00414 DO 80 JCOL = 1, NRHS 00415 DO 70 JROW = NRF, NRF + NR - 1 00416 J = J + 1 00417 RWORK( J ) = REAL( B( JROW, JCOL ) ) 00418 70 CONTINUE 00419 80 CONTINUE 00420 CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU, 00421 $ RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1 ), NR ) 00422 J = NR*NRHS*2 00423 DO 100 JCOL = 1, NRHS 00424 DO 90 JROW = NRF, NRF + NR - 1 00425 J = J + 1 00426 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 00427 90 CONTINUE 00428 100 CONTINUE 00429 CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU, 00430 $ RWORK( 1+NR*NRHS*2 ), NR, ZERO, RWORK( 1+NR*NRHS ), 00431 $ NR ) 00432 JREAL = 0 00433 JIMAG = NR*NRHS 00434 DO 120 JCOL = 1, NRHS 00435 DO 110 JROW = NRF, NRF + NR - 1 00436 JREAL = JREAL + 1 00437 JIMAG = JIMAG + 1 00438 BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 00439 $ RWORK( JIMAG ) ) 00440 110 CONTINUE 00441 120 CONTINUE 00442 * 00443 130 CONTINUE 00444 * 00445 * Next copy the rows of B that correspond to unchanged rows 00446 * in the bidiagonal matrix to BX. 00447 * 00448 DO 140 I = 1, ND 00449 IC = IWORK( INODE+I-1 ) 00450 CALL CCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX ) 00451 140 CONTINUE 00452 * 00453 * Finally go through the left singular vector matrices of all 00454 * the other subproblems bottom-up on the tree. 00455 * 00456 J = 2**NLVL 00457 SQRE = 0 00458 * 00459 DO 160 LVL = NLVL, 1, -1 00460 LVL2 = 2*LVL - 1 00461 * 00462 * find the first node LF and last node LL on 00463 * the current level LVL 00464 * 00465 IF( LVL.EQ.1 ) THEN 00466 LF = 1 00467 LL = 1 00468 ELSE 00469 LF = 2**( LVL-1 ) 00470 LL = 2*LF - 1 00471 END IF 00472 DO 150 I = LF, LL 00473 IM1 = I - 1 00474 IC = IWORK( INODE+IM1 ) 00475 NL = IWORK( NDIML+IM1 ) 00476 NR = IWORK( NDIMR+IM1 ) 00477 NLF = IC - NL 00478 NRF = IC + 1 00479 J = J - 1 00480 CALL CLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX, 00481 $ B( NLF, 1 ), LDB, PERM( NLF, LVL ), 00482 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00483 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 00484 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 00485 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK, 00486 $ INFO ) 00487 150 CONTINUE 00488 160 CONTINUE 00489 GO TO 330 00490 * 00491 * ICOMPQ = 1: applying back the right singular vector factors. 00492 * 00493 170 CONTINUE 00494 * 00495 * First now go through the right singular vector matrices of all 00496 * the tree nodes top-down. 00497 * 00498 J = 0 00499 DO 190 LVL = 1, NLVL 00500 LVL2 = 2*LVL - 1 00501 * 00502 * Find the first node LF and last node LL on 00503 * the current level LVL. 00504 * 00505 IF( LVL.EQ.1 ) THEN 00506 LF = 1 00507 LL = 1 00508 ELSE 00509 LF = 2**( LVL-1 ) 00510 LL = 2*LF - 1 00511 END IF 00512 DO 180 I = LL, LF, -1 00513 IM1 = I - 1 00514 IC = IWORK( INODE+IM1 ) 00515 NL = IWORK( NDIML+IM1 ) 00516 NR = IWORK( NDIMR+IM1 ) 00517 NLF = IC - NL 00518 NRF = IC + 1 00519 IF( I.EQ.LL ) THEN 00520 SQRE = 0 00521 ELSE 00522 SQRE = 1 00523 END IF 00524 J = J + 1 00525 CALL CLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB, 00526 $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ), 00527 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00528 $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ), 00529 $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ), 00530 $ Z( NLF, LVL ), K( J ), C( J ), S( J ), RWORK, 00531 $ INFO ) 00532 180 CONTINUE 00533 190 CONTINUE 00534 * 00535 * The nodes on the bottom level of the tree were solved 00536 * by SLASDQ. The corresponding right singular vector 00537 * matrices are in explicit form. Apply them back. 00538 * 00539 NDB1 = ( ND+1 ) / 2 00540 DO 320 I = NDB1, ND 00541 I1 = I - 1 00542 IC = IWORK( INODE+I1 ) 00543 NL = IWORK( NDIML+I1 ) 00544 NR = IWORK( NDIMR+I1 ) 00545 NLP1 = NL + 1 00546 IF( I.EQ.ND ) THEN 00547 NRP1 = NR 00548 ELSE 00549 NRP1 = NR + 1 00550 END IF 00551 NLF = IC - NL 00552 NRF = IC + 1 00553 * 00554 * Since B and BX are complex, the following call to SGEMM is 00555 * performed in two steps (real and imaginary parts). 00556 * 00557 * CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU, 00558 * $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) 00559 * 00560 J = NLP1*NRHS*2 00561 DO 210 JCOL = 1, NRHS 00562 DO 200 JROW = NLF, NLF + NLP1 - 1 00563 J = J + 1 00564 RWORK( J ) = REAL( B( JROW, JCOL ) ) 00565 200 CONTINUE 00566 210 CONTINUE 00567 CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU, 00568 $ RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO, RWORK( 1 ), 00569 $ NLP1 ) 00570 J = NLP1*NRHS*2 00571 DO 230 JCOL = 1, NRHS 00572 DO 220 JROW = NLF, NLF + NLP1 - 1 00573 J = J + 1 00574 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 00575 220 CONTINUE 00576 230 CONTINUE 00577 CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU, 00578 $ RWORK( 1+NLP1*NRHS*2 ), NLP1, ZERO, 00579 $ RWORK( 1+NLP1*NRHS ), NLP1 ) 00580 JREAL = 0 00581 JIMAG = NLP1*NRHS 00582 DO 250 JCOL = 1, NRHS 00583 DO 240 JROW = NLF, NLF + NLP1 - 1 00584 JREAL = JREAL + 1 00585 JIMAG = JIMAG + 1 00586 BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 00587 $ RWORK( JIMAG ) ) 00588 240 CONTINUE 00589 250 CONTINUE 00590 * 00591 * Since B and BX are complex, the following call to SGEMM is 00592 * performed in two steps (real and imaginary parts). 00593 * 00594 * CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU, 00595 * $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) 00596 * 00597 J = NRP1*NRHS*2 00598 DO 270 JCOL = 1, NRHS 00599 DO 260 JROW = NRF, NRF + NRP1 - 1 00600 J = J + 1 00601 RWORK( J ) = REAL( B( JROW, JCOL ) ) 00602 260 CONTINUE 00603 270 CONTINUE 00604 CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU, 00605 $ RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO, RWORK( 1 ), 00606 $ NRP1 ) 00607 J = NRP1*NRHS*2 00608 DO 290 JCOL = 1, NRHS 00609 DO 280 JROW = NRF, NRF + NRP1 - 1 00610 J = J + 1 00611 RWORK( J ) = AIMAG( B( JROW, JCOL ) ) 00612 280 CONTINUE 00613 290 CONTINUE 00614 CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU, 00615 $ RWORK( 1+NRP1*NRHS*2 ), NRP1, ZERO, 00616 $ RWORK( 1+NRP1*NRHS ), NRP1 ) 00617 JREAL = 0 00618 JIMAG = NRP1*NRHS 00619 DO 310 JCOL = 1, NRHS 00620 DO 300 JROW = NRF, NRF + NRP1 - 1 00621 JREAL = JREAL + 1 00622 JIMAG = JIMAG + 1 00623 BX( JROW, JCOL ) = CMPLX( RWORK( JREAL ), 00624 $ RWORK( JIMAG ) ) 00625 300 CONTINUE 00626 310 CONTINUE 00627 * 00628 320 CONTINUE 00629 * 00630 330 CONTINUE 00631 * 00632 RETURN 00633 * 00634 * End of CLALSA 00635 * 00636 END