001:       SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
002:      $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
003:      $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
004:      $                   IWORK, INFO )
005: *
006: *  -- LAPACK routine (version 3.2) --
007: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
012:      $                   SMLSIZ
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
016:      $                   K( * ), PERM( LDGCOL, * )
017:       DOUBLE PRECISION   B( LDB, * ), BX( LDBX, * ), C( * ),
018:      $                   DIFL( LDU, * ), DIFR( LDU, * ),
019:      $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
020:      $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
021:      $                   Z( LDU, * )
022: *     ..
023: *
024: *  Purpose
025: *  =======
026: *
027: *  DLALSA is an itermediate step in solving the least squares problem
028: *  by computing the SVD of the coefficient matrix in compact form (The
029: *  singular vectors are computed as products of simple orthorgonal
030: *  matrices.).
031: *
032: *  If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
033: *  matrix of an upper bidiagonal matrix to the right hand side; and if
034: *  ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
035: *  right hand side. The singular vector matrices were generated in
036: *  compact form by DLALSA.
037: *
038: *  Arguments
039: *  =========
040: *
041: *
042: *  ICOMPQ (input) INTEGER
043: *         Specifies whether the left or the right singular vector
044: *         matrix is involved.
045: *         = 0: Left singular vector matrix
046: *         = 1: Right singular vector matrix
047: *
048: *  SMLSIZ (input) INTEGER
049: *         The maximum size of the subproblems at the bottom of the
050: *         computation tree.
051: *
052: *  N      (input) INTEGER
053: *         The row and column dimensions of the upper bidiagonal matrix.
054: *
055: *  NRHS   (input) INTEGER
056: *         The number of columns of B and BX. NRHS must be at least 1.
057: *
058: *  B      (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
059: *         On input, B contains the right hand sides of the least
060: *         squares problem in rows 1 through M.
061: *         On output, B contains the solution X in rows 1 through N.
062: *
063: *  LDB    (input) INTEGER
064: *         The leading dimension of B in the calling subprogram.
065: *         LDB must be at least max(1,MAX( M, N ) ).
066: *
067: *  BX     (output) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
068: *         On exit, the result of applying the left or right singular
069: *         vector matrix to B.
070: *
071: *  LDBX   (input) INTEGER
072: *         The leading dimension of BX.
073: *
074: *  U      (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
075: *         On entry, U contains the left singular vector matrices of all
076: *         subproblems at the bottom level.
077: *
078: *  LDU    (input) INTEGER, LDU = > N.
079: *         The leading dimension of arrays U, VT, DIFL, DIFR,
080: *         POLES, GIVNUM, and Z.
081: *
082: *  VT     (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
083: *         On entry, VT' contains the right singular vector matrices of
084: *         all subproblems at the bottom level.
085: *
086: *  K      (input) INTEGER array, dimension ( N ).
087: *
088: *  DIFL   (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
089: *         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
090: *
091: *  DIFR   (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
092: *         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
093: *         distances between singular values on the I-th level and
094: *         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
095: *         record the normalizing factors of the right singular vectors
096: *         matrices of subproblems on I-th level.
097: *
098: *  Z      (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
099: *         On entry, Z(1, I) contains the components of the deflation-
100: *         adjusted updating row vector for subproblems on the I-th
101: *         level.
102: *
103: *  POLES  (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
104: *         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
105: *         singular values involved in the secular equations on the I-th
106: *         level.
107: *
108: *  GIVPTR (input) INTEGER array, dimension ( N ).
109: *         On entry, GIVPTR( I ) records the number of Givens
110: *         rotations performed on the I-th problem on the computation
111: *         tree.
112: *
113: *  GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
114: *         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
115: *         locations of Givens rotations performed on the I-th level on
116: *         the computation tree.
117: *
118: *  LDGCOL (input) INTEGER, LDGCOL = > N.
119: *         The leading dimension of arrays GIVCOL and PERM.
120: *
121: *  PERM   (input) INTEGER array, dimension ( LDGCOL, NLVL ).
122: *         On entry, PERM(*, I) records permutations done on the I-th
123: *         level of the computation tree.
124: *
125: *  GIVNUM (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
126: *         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
127: *         values of Givens rotations performed on the I-th level on the
128: *         computation tree.
129: *
130: *  C      (input) DOUBLE PRECISION array, dimension ( N ).
131: *         On entry, if the I-th subproblem is not square,
132: *         C( I ) contains the C-value of a Givens rotation related to
133: *         the right null space of the I-th subproblem.
134: *
135: *  S      (input) DOUBLE PRECISION array, dimension ( N ).
136: *         On entry, if the I-th subproblem is not square,
137: *         S( I ) contains the S-value of a Givens rotation related to
138: *         the right null space of the I-th subproblem.
139: *
140: *  WORK   (workspace) DOUBLE PRECISION array.
141: *         The dimension must be at least N.
142: *
143: *  IWORK  (workspace) INTEGER array.
144: *         The dimension must be at least 3 * N
145: *
146: *  INFO   (output) INTEGER
147: *          = 0:  successful exit.
148: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
149: *
150: *  Further Details
151: *  ===============
152: *
153: *  Based on contributions by
154: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
155: *       California at Berkeley, USA
156: *     Osni Marques, LBNL/NERSC, USA
157: *
158: *  =====================================================================
159: *
160: *     .. Parameters ..
161:       DOUBLE PRECISION   ZERO, ONE
162:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
163: *     ..
164: *     .. Local Scalars ..
165:       INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
166:      $                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
167:      $                   NR, NRF, NRP1, SQRE
168: *     ..
169: *     .. External Subroutines ..
170:       EXTERNAL           DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
171: *     ..
172: *     .. Executable Statements ..
173: *
174: *     Test the input parameters.
175: *
176:       INFO = 0
177: *
178:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
179:          INFO = -1
180:       ELSE IF( SMLSIZ.LT.3 ) THEN
181:          INFO = -2
182:       ELSE IF( N.LT.SMLSIZ ) THEN
183:          INFO = -3
184:       ELSE IF( NRHS.LT.1 ) THEN
185:          INFO = -4
186:       ELSE IF( LDB.LT.N ) THEN
187:          INFO = -6
188:       ELSE IF( LDBX.LT.N ) THEN
189:          INFO = -8
190:       ELSE IF( LDU.LT.N ) THEN
191:          INFO = -10
192:       ELSE IF( LDGCOL.LT.N ) THEN
193:          INFO = -19
194:       END IF
195:       IF( INFO.NE.0 ) THEN
196:          CALL XERBLA( 'DLALSA', -INFO )
197:          RETURN
198:       END IF
199: *
200: *     Book-keeping and  setting up the computation tree.
201: *
202:       INODE = 1
203:       NDIML = INODE + N
204:       NDIMR = NDIML + N
205: *
206:       CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
207:      $             IWORK( NDIMR ), SMLSIZ )
208: *
209: *     The following code applies back the left singular vector factors.
210: *     For applying back the right singular vector factors, go to 50.
211: *
212:       IF( ICOMPQ.EQ.1 ) THEN
213:          GO TO 50
214:       END IF
215: *
216: *     The nodes on the bottom level of the tree were solved
217: *     by DLASDQ. The corresponding left and right singular vector
218: *     matrices are in explicit form. First apply back the left
219: *     singular vector matrices.
220: *
221:       NDB1 = ( ND+1 ) / 2
222:       DO 10 I = NDB1, ND
223: *
224: *        IC : center row of each node
225: *        NL : number of rows of left  subproblem
226: *        NR : number of rows of right subproblem
227: *        NLF: starting row of the left   subproblem
228: *        NRF: starting row of the right  subproblem
229: *
230:          I1 = I - 1
231:          IC = IWORK( INODE+I1 )
232:          NL = IWORK( NDIML+I1 )
233:          NR = IWORK( NDIMR+I1 )
234:          NLF = IC - NL
235:          NRF = IC + 1
236:          CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
237:      $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
238:          CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
239:      $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
240:    10 CONTINUE
241: *
242: *     Next copy the rows of B that correspond to unchanged rows
243: *     in the bidiagonal matrix to BX.
244: *
245:       DO 20 I = 1, ND
246:          IC = IWORK( INODE+I-1 )
247:          CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
248:    20 CONTINUE
249: *
250: *     Finally go through the left singular vector matrices of all
251: *     the other subproblems bottom-up on the tree.
252: *
253:       J = 2**NLVL
254:       SQRE = 0
255: *
256:       DO 40 LVL = NLVL, 1, -1
257:          LVL2 = 2*LVL - 1
258: *
259: *        find the first node LF and last node LL on
260: *        the current level LVL
261: *
262:          IF( LVL.EQ.1 ) THEN
263:             LF = 1
264:             LL = 1
265:          ELSE
266:             LF = 2**( LVL-1 )
267:             LL = 2*LF - 1
268:          END IF
269:          DO 30 I = LF, LL
270:             IM1 = I - 1
271:             IC = IWORK( INODE+IM1 )
272:             NL = IWORK( NDIML+IM1 )
273:             NR = IWORK( NDIMR+IM1 )
274:             NLF = IC - NL
275:             NRF = IC + 1
276:             J = J - 1
277:             CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
278:      $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
279:      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
280:      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
281:      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
282:      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
283:      $                   INFO )
284:    30    CONTINUE
285:    40 CONTINUE
286:       GO TO 90
287: *
288: *     ICOMPQ = 1: applying back the right singular vector factors.
289: *
290:    50 CONTINUE
291: *
292: *     First now go through the right singular vector matrices of all
293: *     the tree nodes top-down.
294: *
295:       J = 0
296:       DO 70 LVL = 1, NLVL
297:          LVL2 = 2*LVL - 1
298: *
299: *        Find the first node LF and last node LL on
300: *        the current level LVL.
301: *
302:          IF( LVL.EQ.1 ) THEN
303:             LF = 1
304:             LL = 1
305:          ELSE
306:             LF = 2**( LVL-1 )
307:             LL = 2*LF - 1
308:          END IF
309:          DO 60 I = LL, LF, -1
310:             IM1 = I - 1
311:             IC = IWORK( INODE+IM1 )
312:             NL = IWORK( NDIML+IM1 )
313:             NR = IWORK( NDIMR+IM1 )
314:             NLF = IC - NL
315:             NRF = IC + 1
316:             IF( I.EQ.LL ) THEN
317:                SQRE = 0
318:             ELSE
319:                SQRE = 1
320:             END IF
321:             J = J + 1
322:             CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
323:      $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
324:      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
325:      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
326:      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
327:      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
328:      $                   INFO )
329:    60    CONTINUE
330:    70 CONTINUE
331: *
332: *     The nodes on the bottom level of the tree were solved
333: *     by DLASDQ. The corresponding right singular vector
334: *     matrices are in explicit form. Apply them back.
335: *
336:       NDB1 = ( ND+1 ) / 2
337:       DO 80 I = NDB1, ND
338:          I1 = I - 1
339:          IC = IWORK( INODE+I1 )
340:          NL = IWORK( NDIML+I1 )
341:          NR = IWORK( NDIMR+I1 )
342:          NLP1 = NL + 1
343:          IF( I.EQ.ND ) THEN
344:             NRP1 = NR
345:          ELSE
346:             NRP1 = NR + 1
347:          END IF
348:          NLF = IC - NL
349:          NRF = IC + 1
350:          CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
351:      $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
352:          CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
353:      $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
354:    80 CONTINUE
355: *
356:    90 CONTINUE
357: *
358:       RETURN
359: *
360: *     End of DLALSA
361: *
362:       END
363: