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