LAPACK: SRC/clalsa.f Source File

Go to the documentation of this file.
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