LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clalsa.f
Go to the documentation of this file.
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
 All Files Functions