1: SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
2: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
3: $ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
4: *
5: * -- LAPACK routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
12: $ LDGNUM, NL, NR, NRHS, SQRE
13: DOUBLE PRECISION C, S
14: * ..
15: * .. Array Arguments ..
16: INTEGER GIVCOL( LDGCOL, * ), PERM( * )
17: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ),
18: $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
19: $ POLES( LDGNUM, * ), WORK( * ), Z( * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * DLALS0 applies back the multiplying factors of either the left or the
26: * right singular vector matrix of a diagonal matrix appended by a row
27: * to the right hand side matrix B in solving the least squares problem
28: * using the divide-and-conquer SVD approach.
29: *
30: * For the left singular vector matrix, three types of orthogonal
31: * matrices are involved:
32: *
33: * (1L) Givens rotations: the number of such rotations is GIVPTR; the
34: * pairs of columns/rows they were applied to are stored in GIVCOL;
35: * and the C- and S-values of these rotations are stored in GIVNUM.
36: *
37: * (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
38: * row, and for J=2:N, PERM(J)-th row of B is to be moved to the
39: * J-th row.
40: *
41: * (3L) The left singular vector matrix of the remaining matrix.
42: *
43: * For the right singular vector matrix, four types of orthogonal
44: * matrices are involved:
45: *
46: * (1R) The right singular vector matrix of the remaining matrix.
47: *
48: * (2R) If SQRE = 1, one extra Givens rotation to generate the right
49: * null space.
50: *
51: * (3R) The inverse transformation of (2L).
52: *
53: * (4R) The inverse transformation of (1L).
54: *
55: * Arguments
56: * =========
57: *
58: * ICOMPQ (input) INTEGER
59: * Specifies whether singular vectors are to be computed in
60: * factored form:
61: * = 0: Left singular vector matrix.
62: * = 1: Right singular vector matrix.
63: *
64: * NL (input) INTEGER
65: * The row dimension of the upper block. NL >= 1.
66: *
67: * NR (input) INTEGER
68: * The row dimension of the lower block. NR >= 1.
69: *
70: * SQRE (input) INTEGER
71: * = 0: the lower block is an NR-by-NR square matrix.
72: * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
73: *
74: * The bidiagonal matrix has row dimension N = NL + NR + 1,
75: * and column dimension M = N + SQRE.
76: *
77: * NRHS (input) INTEGER
78: * The number of columns of B and BX. NRHS must be at least 1.
79: *
80: * B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
81: * On input, B contains the right hand sides of the least
82: * squares problem in rows 1 through M. On output, B contains
83: * the solution X in rows 1 through N.
84: *
85: * LDB (input) INTEGER
86: * The leading dimension of B. LDB must be at least
87: * max(1,MAX( M, N ) ).
88: *
89: * BX (workspace) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
90: *
91: * LDBX (input) INTEGER
92: * The leading dimension of BX.
93: *
94: * PERM (input) INTEGER array, dimension ( N )
95: * The permutations (from deflation and sorting) applied
96: * to the two blocks.
97: *
98: * GIVPTR (input) INTEGER
99: * The number of Givens rotations which took place in this
100: * subproblem.
101: *
102: * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
103: * Each pair of numbers indicates a pair of rows/columns
104: * involved in a Givens rotation.
105: *
106: * LDGCOL (input) INTEGER
107: * The leading dimension of GIVCOL, must be at least N.
108: *
109: * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
110: * Each number indicates the C or S value used in the
111: * corresponding Givens rotation.
112: *
113: * LDGNUM (input) INTEGER
114: * The leading dimension of arrays DIFR, POLES and
115: * GIVNUM, must be at least K.
116: *
117: * POLES (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
118: * On entry, POLES(1:K, 1) contains the new singular
119: * values obtained from solving the secular equation, and
120: * POLES(1:K, 2) is an array containing the poles in the secular
121: * equation.
122: *
123: * DIFL (input) DOUBLE PRECISION array, dimension ( K ).
124: * On entry, DIFL(I) is the distance between I-th updated
125: * (undeflated) singular value and the I-th (undeflated) old
126: * singular value.
127: *
128: * DIFR (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
129: * On entry, DIFR(I, 1) contains the distances between I-th
130: * updated (undeflated) singular value and the I+1-th
131: * (undeflated) old singular value. And DIFR(I, 2) is the
132: * normalizing factor for the I-th right singular vector.
133: *
134: * Z (input) DOUBLE PRECISION array, dimension ( K )
135: * Contain the components of the deflation-adjusted updating row
136: * vector.
137: *
138: * K (input) INTEGER
139: * Contains the dimension of the non-deflated matrix,
140: * This is the order of the related secular equation. 1 <= K <=N.
141: *
142: * C (input) DOUBLE PRECISION
143: * C contains garbage if SQRE =0 and the C-value of a Givens
144: * rotation related to the right null space if SQRE = 1.
145: *
146: * S (input) DOUBLE PRECISION
147: * S contains garbage if SQRE =0 and the S-value of a Givens
148: * rotation related to the right null space if SQRE = 1.
149: *
150: * WORK (workspace) DOUBLE PRECISION array, dimension ( K )
151: *
152: * INFO (output) INTEGER
153: * = 0: successful exit.
154: * < 0: if INFO = -i, the i-th argument had an illegal value.
155: *
156: * Further Details
157: * ===============
158: *
159: * Based on contributions by
160: * Ming Gu and Ren-Cang Li, Computer Science Division, University of
161: * California at Berkeley, USA
162: * Osni Marques, LBNL/NERSC, USA
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION ONE, ZERO, NEGONE
168: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
169: * ..
170: * .. Local Scalars ..
171: INTEGER I, J, M, N, NLP1
172: DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
173: * ..
174: * .. External Subroutines ..
175: EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DROT, DSCAL,
176: $ XERBLA
177: * ..
178: * .. External Functions ..
179: DOUBLE PRECISION DLAMC3, DNRM2
180: EXTERNAL DLAMC3, DNRM2
181: * ..
182: * .. Intrinsic Functions ..
183: INTRINSIC MAX
184: * ..
185: * .. Executable Statements ..
186: *
187: * Test the input parameters.
188: *
189: INFO = 0
190: *
191: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
192: INFO = -1
193: ELSE IF( NL.LT.1 ) THEN
194: INFO = -2
195: ELSE IF( NR.LT.1 ) THEN
196: INFO = -3
197: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
198: INFO = -4
199: END IF
200: *
201: N = NL + NR + 1
202: *
203: IF( NRHS.LT.1 ) THEN
204: INFO = -5
205: ELSE IF( LDB.LT.N ) THEN
206: INFO = -7
207: ELSE IF( LDBX.LT.N ) THEN
208: INFO = -9
209: ELSE IF( GIVPTR.LT.0 ) THEN
210: INFO = -11
211: ELSE IF( LDGCOL.LT.N ) THEN
212: INFO = -13
213: ELSE IF( LDGNUM.LT.N ) THEN
214: INFO = -15
215: ELSE IF( K.LT.1 ) THEN
216: INFO = -20
217: END IF
218: IF( INFO.NE.0 ) THEN
219: CALL XERBLA( 'DLALS0', -INFO )
220: RETURN
221: END IF
222: *
223: M = N + SQRE
224: NLP1 = NL + 1
225: *
226: IF( ICOMPQ.EQ.0 ) THEN
227: *
228: * Apply back orthogonal transformations from the left.
229: *
230: * Step (1L): apply back the Givens rotations performed.
231: *
232: DO 10 I = 1, GIVPTR
233: CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
234: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
235: $ GIVNUM( I, 1 ) )
236: 10 CONTINUE
237: *
238: * Step (2L): permute rows of B.
239: *
240: CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
241: DO 20 I = 2, N
242: CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
243: 20 CONTINUE
244: *
245: * Step (3L): apply the inverse of the left singular vector
246: * matrix to BX.
247: *
248: IF( K.EQ.1 ) THEN
249: CALL DCOPY( NRHS, BX, LDBX, B, LDB )
250: IF( Z( 1 ).LT.ZERO ) THEN
251: CALL DSCAL( NRHS, NEGONE, B, LDB )
252: END IF
253: ELSE
254: DO 50 J = 1, K
255: DIFLJ = DIFL( J )
256: DJ = POLES( J, 1 )
257: DSIGJ = -POLES( J, 2 )
258: IF( J.LT.K ) THEN
259: DIFRJ = -DIFR( J, 1 )
260: DSIGJP = -POLES( J+1, 2 )
261: END IF
262: IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
263: $ THEN
264: WORK( J ) = ZERO
265: ELSE
266: WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
267: $ ( POLES( J, 2 )+DJ )
268: END IF
269: DO 30 I = 1, J - 1
270: IF( ( Z( I ).EQ.ZERO ) .OR.
271: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
272: WORK( I ) = ZERO
273: ELSE
274: WORK( I ) = POLES( I, 2 )*Z( I ) /
275: $ ( DLAMC3( POLES( I, 2 ), DSIGJ )-
276: $ DIFLJ ) / ( POLES( I, 2 )+DJ )
277: END IF
278: 30 CONTINUE
279: DO 40 I = J + 1, K
280: IF( ( Z( I ).EQ.ZERO ) .OR.
281: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
282: WORK( I ) = ZERO
283: ELSE
284: WORK( I ) = POLES( I, 2 )*Z( I ) /
285: $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+
286: $ DIFRJ ) / ( POLES( I, 2 )+DJ )
287: END IF
288: 40 CONTINUE
289: WORK( 1 ) = NEGONE
290: TEMP = DNRM2( K, WORK, 1 )
291: CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
292: $ B( J, 1 ), LDB )
293: CALL DLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
294: $ LDB, INFO )
295: 50 CONTINUE
296: END IF
297: *
298: * Move the deflated rows of BX to B also.
299: *
300: IF( K.LT.MAX( M, N ) )
301: $ CALL DLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
302: $ B( K+1, 1 ), LDB )
303: ELSE
304: *
305: * Apply back the right orthogonal transformations.
306: *
307: * Step (1R): apply back the new right singular vector matrix
308: * to B.
309: *
310: IF( K.EQ.1 ) THEN
311: CALL DCOPY( NRHS, B, LDB, BX, LDBX )
312: ELSE
313: DO 80 J = 1, K
314: DSIGJ = POLES( J, 2 )
315: IF( Z( J ).EQ.ZERO ) THEN
316: WORK( J ) = ZERO
317: ELSE
318: WORK( J ) = -Z( J ) / DIFL( J ) /
319: $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
320: END IF
321: DO 60 I = 1, J - 1
322: IF( Z( J ).EQ.ZERO ) THEN
323: WORK( I ) = ZERO
324: ELSE
325: WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
326: $ 2 ) )-DIFR( I, 1 ) ) /
327: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
328: END IF
329: 60 CONTINUE
330: DO 70 I = J + 1, K
331: IF( Z( J ).EQ.ZERO ) THEN
332: WORK( I ) = ZERO
333: ELSE
334: WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
335: $ 2 ) )-DIFL( I ) ) /
336: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
337: END IF
338: 70 CONTINUE
339: CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
340: $ BX( J, 1 ), LDBX )
341: 80 CONTINUE
342: END IF
343: *
344: * Step (2R): if SQRE = 1, apply back the rotation that is
345: * related to the right null space of the subproblem.
346: *
347: IF( SQRE.EQ.1 ) THEN
348: CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
349: CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
350: END IF
351: IF( K.LT.MAX( M, N ) )
352: $ CALL DLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
353: $ LDBX )
354: *
355: * Step (3R): permute rows of B.
356: *
357: CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
358: IF( SQRE.EQ.1 ) THEN
359: CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
360: END IF
361: DO 90 I = 2, N
362: CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
363: 90 CONTINUE
364: *
365: * Step (4R): apply back the Givens rotations performed.
366: *
367: DO 100 I = GIVPTR, 1, -1
368: CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
369: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
370: $ -GIVNUM( I, 1 ) )
371: 100 CONTINUE
372: END IF
373: *
374: RETURN
375: *
376: * End of DLALS0
377: *
378: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlals0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack