1: *> \brief \b DLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer SVD approach. Used by sgelsd.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLALS0 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlals0.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlals0.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlals0.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
22: * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
23: * POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
24: *
25: * .. Scalar Arguments ..
26: * INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
27: * $ LDGNUM, NL, NR, NRHS, SQRE
28: * DOUBLE PRECISION C, S
29: * ..
30: * .. Array Arguments ..
31: * INTEGER GIVCOL( LDGCOL, * ), PERM( * )
32: * DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ),
33: * $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
34: * $ POLES( LDGNUM, * ), WORK( * ), Z( * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> DLALS0 applies back the multiplying factors of either the left or the
44: *> right singular vector matrix of a diagonal matrix appended by a row
45: *> to the right hand side matrix B in solving the least squares problem
46: *> using the divide-and-conquer SVD approach.
47: *>
48: *> For the left singular vector matrix, three types of orthogonal
49: *> matrices are involved:
50: *>
51: *> (1L) Givens rotations: the number of such rotations is GIVPTR; the
52: *> pairs of columns/rows they were applied to are stored in GIVCOL;
53: *> and the C- and S-values of these rotations are stored in GIVNUM.
54: *>
55: *> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
56: *> row, and for J=2:N, PERM(J)-th row of B is to be moved to the
57: *> J-th row.
58: *>
59: *> (3L) The left singular vector matrix of the remaining matrix.
60: *>
61: *> For the right singular vector matrix, four types of orthogonal
62: *> matrices are involved:
63: *>
64: *> (1R) The right singular vector matrix of the remaining matrix.
65: *>
66: *> (2R) If SQRE = 1, one extra Givens rotation to generate the right
67: *> null space.
68: *>
69: *> (3R) The inverse transformation of (2L).
70: *>
71: *> (4R) The inverse transformation of (1L).
72: *> \endverbatim
73: *
74: * Arguments:
75: * ==========
76: *
77: *> \param[in] ICOMPQ
78: *> \verbatim
79: *> ICOMPQ is INTEGER
80: *> Specifies whether singular vectors are to be computed in
81: *> factored form:
82: *> = 0: Left singular vector matrix.
83: *> = 1: Right singular vector matrix.
84: *> \endverbatim
85: *>
86: *> \param[in] NL
87: *> \verbatim
88: *> NL is INTEGER
89: *> The row dimension of the upper block. NL >= 1.
90: *> \endverbatim
91: *>
92: *> \param[in] NR
93: *> \verbatim
94: *> NR is INTEGER
95: *> The row dimension of the lower block. NR >= 1.
96: *> \endverbatim
97: *>
98: *> \param[in] SQRE
99: *> \verbatim
100: *> SQRE is INTEGER
101: *> = 0: the lower block is an NR-by-NR square matrix.
102: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
103: *>
104: *> The bidiagonal matrix has row dimension N = NL + NR + 1,
105: *> and column dimension M = N + SQRE.
106: *> \endverbatim
107: *>
108: *> \param[in] NRHS
109: *> \verbatim
110: *> NRHS is INTEGER
111: *> The number of columns of B and BX. NRHS must be at least 1.
112: *> \endverbatim
113: *>
114: *> \param[in,out] B
115: *> \verbatim
116: *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS )
117: *> On input, B contains the right hand sides of the least
118: *> squares problem in rows 1 through M. On output, B contains
119: *> the solution X in rows 1 through N.
120: *> \endverbatim
121: *>
122: *> \param[in] LDB
123: *> \verbatim
124: *> LDB is INTEGER
125: *> The leading dimension of B. LDB must be at least
126: *> max(1,MAX( M, N ) ).
127: *> \endverbatim
128: *>
129: *> \param[out] BX
130: *> \verbatim
131: *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS )
132: *> \endverbatim
133: *>
134: *> \param[in] LDBX
135: *> \verbatim
136: *> LDBX is INTEGER
137: *> The leading dimension of BX.
138: *> \endverbatim
139: *>
140: *> \param[in] PERM
141: *> \verbatim
142: *> PERM is INTEGER array, dimension ( N )
143: *> The permutations (from deflation and sorting) applied
144: *> to the two blocks.
145: *> \endverbatim
146: *>
147: *> \param[in] GIVPTR
148: *> \verbatim
149: *> GIVPTR is INTEGER
150: *> The number of Givens rotations which took place in this
151: *> subproblem.
152: *> \endverbatim
153: *>
154: *> \param[in] GIVCOL
155: *> \verbatim
156: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
157: *> Each pair of numbers indicates a pair of rows/columns
158: *> involved in a Givens rotation.
159: *> \endverbatim
160: *>
161: *> \param[in] LDGCOL
162: *> \verbatim
163: *> LDGCOL is INTEGER
164: *> The leading dimension of GIVCOL, must be at least N.
165: *> \endverbatim
166: *>
167: *> \param[in] GIVNUM
168: *> \verbatim
169: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
170: *> Each number indicates the C or S value used in the
171: *> corresponding Givens rotation.
172: *> \endverbatim
173: *>
174: *> \param[in] LDGNUM
175: *> \verbatim
176: *> LDGNUM is INTEGER
177: *> The leading dimension of arrays DIFR, POLES and
178: *> GIVNUM, must be at least K.
179: *> \endverbatim
180: *>
181: *> \param[in] POLES
182: *> \verbatim
183: *> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
184: *> On entry, POLES(1:K, 1) contains the new singular
185: *> values obtained from solving the secular equation, and
186: *> POLES(1:K, 2) is an array containing the poles in the secular
187: *> equation.
188: *> \endverbatim
189: *>
190: *> \param[in] DIFL
191: *> \verbatim
192: *> DIFL is DOUBLE PRECISION array, dimension ( K ).
193: *> On entry, DIFL(I) is the distance between I-th updated
194: *> (undeflated) singular value and the I-th (undeflated) old
195: *> singular value.
196: *> \endverbatim
197: *>
198: *> \param[in] DIFR
199: *> \verbatim
200: *> DIFR is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
201: *> On entry, DIFR(I, 1) contains the distances between I-th
202: *> updated (undeflated) singular value and the I+1-th
203: *> (undeflated) old singular value. And DIFR(I, 2) is the
204: *> normalizing factor for the I-th right singular vector.
205: *> \endverbatim
206: *>
207: *> \param[in] Z
208: *> \verbatim
209: *> Z is DOUBLE PRECISION array, dimension ( K )
210: *> Contain the components of the deflation-adjusted updating row
211: *> vector.
212: *> \endverbatim
213: *>
214: *> \param[in] K
215: *> \verbatim
216: *> K is INTEGER
217: *> Contains the dimension of the non-deflated matrix,
218: *> This is the order of the related secular equation. 1 <= K <=N.
219: *> \endverbatim
220: *>
221: *> \param[in] C
222: *> \verbatim
223: *> C is DOUBLE PRECISION
224: *> C contains garbage if SQRE =0 and the C-value of a Givens
225: *> rotation related to the right null space if SQRE = 1.
226: *> \endverbatim
227: *>
228: *> \param[in] S
229: *> \verbatim
230: *> S is DOUBLE PRECISION
231: *> S contains garbage if SQRE =0 and the S-value of a Givens
232: *> rotation related to the right null space if SQRE = 1.
233: *> \endverbatim
234: *>
235: *> \param[out] WORK
236: *> \verbatim
237: *> WORK is DOUBLE PRECISION array, dimension ( K )
238: *> \endverbatim
239: *>
240: *> \param[out] INFO
241: *> \verbatim
242: *> INFO is INTEGER
243: *> = 0: successful exit.
244: *> < 0: if INFO = -i, the i-th argument had an illegal value.
245: *> \endverbatim
246: *
247: * Authors:
248: * ========
249: *
250: *> \author Univ. of Tennessee
251: *> \author Univ. of California Berkeley
252: *> \author Univ. of Colorado Denver
253: *> \author NAG Ltd.
254: *
255: *> \ingroup doubleOTHERcomputational
256: *
257: *> \par Contributors:
258: * ==================
259: *>
260: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
261: *> California at Berkeley, USA \n
262: *> Osni Marques, LBNL/NERSC, USA \n
263: *
264: * =====================================================================
265: SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
266: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
267: $ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
268: *
269: * -- LAPACK computational routine --
270: * -- LAPACK is a software package provided by Univ. of Tennessee, --
271: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
272: *
273: * .. Scalar Arguments ..
274: INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
275: $ LDGNUM, NL, NR, NRHS, SQRE
276: DOUBLE PRECISION C, S
277: * ..
278: * .. Array Arguments ..
279: INTEGER GIVCOL( LDGCOL, * ), PERM( * )
280: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ),
281: $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
282: $ POLES( LDGNUM, * ), WORK( * ), Z( * )
283: * ..
284: *
285: * =====================================================================
286: *
287: * .. Parameters ..
288: DOUBLE PRECISION ONE, ZERO, NEGONE
289: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
290: * ..
291: * .. Local Scalars ..
292: INTEGER I, J, M, N, NLP1
293: DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
294: * ..
295: * .. External Subroutines ..
296: EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DROT, DSCAL,
297: $ XERBLA
298: * ..
299: * .. External Functions ..
300: DOUBLE PRECISION DLAMC3, DNRM2
301: EXTERNAL DLAMC3, DNRM2
302: * ..
303: * .. Intrinsic Functions ..
304: INTRINSIC MAX
305: * ..
306: * .. Executable Statements ..
307: *
308: * Test the input parameters.
309: *
310: INFO = 0
311: N = NL + NR + 1
312: *
313: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
314: INFO = -1
315: ELSE IF( NL.LT.1 ) THEN
316: INFO = -2
317: ELSE IF( NR.LT.1 ) THEN
318: INFO = -3
319: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
320: INFO = -4
321: ELSE IF( NRHS.LT.1 ) THEN
322: INFO = -5
323: ELSE IF( LDB.LT.N ) THEN
324: INFO = -7
325: ELSE IF( LDBX.LT.N ) THEN
326: INFO = -9
327: ELSE IF( GIVPTR.LT.0 ) THEN
328: INFO = -11
329: ELSE IF( LDGCOL.LT.N ) THEN
330: INFO = -13
331: ELSE IF( LDGNUM.LT.N ) THEN
332: INFO = -15
333: ELSE IF( K.LT.1 ) THEN
334: INFO = -20
335: END IF
336: IF( INFO.NE.0 ) THEN
337: CALL XERBLA( 'DLALS0', -INFO )
338: RETURN
339: END IF
340: *
341: M = N + SQRE
342: NLP1 = NL + 1
343: *
344: IF( ICOMPQ.EQ.0 ) THEN
345: *
346: * Apply back orthogonal transformations from the left.
347: *
348: * Step (1L): apply back the Givens rotations performed.
349: *
350: DO 10 I = 1, GIVPTR
351: CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
352: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
353: $ GIVNUM( I, 1 ) )
354: 10 CONTINUE
355: *
356: * Step (2L): permute rows of B.
357: *
358: CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
359: DO 20 I = 2, N
360: CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
361: 20 CONTINUE
362: *
363: * Step (3L): apply the inverse of the left singular vector
364: * matrix to BX.
365: *
366: IF( K.EQ.1 ) THEN
367: CALL DCOPY( NRHS, BX, LDBX, B, LDB )
368: IF( Z( 1 ).LT.ZERO ) THEN
369: CALL DSCAL( NRHS, NEGONE, B, LDB )
370: END IF
371: ELSE
372: DO 50 J = 1, K
373: DIFLJ = DIFL( J )
374: DJ = POLES( J, 1 )
375: DSIGJ = -POLES( J, 2 )
376: IF( J.LT.K ) THEN
377: DIFRJ = -DIFR( J, 1 )
378: DSIGJP = -POLES( J+1, 2 )
379: END IF
380: IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
381: $ THEN
382: WORK( J ) = ZERO
383: ELSE
384: WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
385: $ ( POLES( J, 2 )+DJ )
386: END IF
387: DO 30 I = 1, J - 1
388: IF( ( Z( I ).EQ.ZERO ) .OR.
389: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
390: WORK( I ) = ZERO
391: ELSE
392: WORK( I ) = POLES( I, 2 )*Z( I ) /
393: $ ( DLAMC3( POLES( I, 2 ), DSIGJ )-
394: $ DIFLJ ) / ( POLES( I, 2 )+DJ )
395: END IF
396: 30 CONTINUE
397: DO 40 I = J + 1, K
398: IF( ( Z( I ).EQ.ZERO ) .OR.
399: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
400: WORK( I ) = ZERO
401: ELSE
402: WORK( I ) = POLES( I, 2 )*Z( I ) /
403: $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+
404: $ DIFRJ ) / ( POLES( I, 2 )+DJ )
405: END IF
406: 40 CONTINUE
407: WORK( 1 ) = NEGONE
408: TEMP = DNRM2( K, WORK, 1 )
409: CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
410: $ B( J, 1 ), LDB )
411: CALL DLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
412: $ LDB, INFO )
413: 50 CONTINUE
414: END IF
415: *
416: * Move the deflated rows of BX to B also.
417: *
418: IF( K.LT.MAX( M, N ) )
419: $ CALL DLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
420: $ B( K+1, 1 ), LDB )
421: ELSE
422: *
423: * Apply back the right orthogonal transformations.
424: *
425: * Step (1R): apply back the new right singular vector matrix
426: * to B.
427: *
428: IF( K.EQ.1 ) THEN
429: CALL DCOPY( NRHS, B, LDB, BX, LDBX )
430: ELSE
431: DO 80 J = 1, K
432: DSIGJ = POLES( J, 2 )
433: IF( Z( J ).EQ.ZERO ) THEN
434: WORK( J ) = ZERO
435: ELSE
436: WORK( J ) = -Z( J ) / DIFL( J ) /
437: $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
438: END IF
439: DO 60 I = 1, J - 1
440: IF( Z( J ).EQ.ZERO ) THEN
441: WORK( I ) = ZERO
442: ELSE
443: WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
444: $ 2 ) )-DIFR( I, 1 ) ) /
445: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
446: END IF
447: 60 CONTINUE
448: DO 70 I = J + 1, K
449: IF( Z( J ).EQ.ZERO ) THEN
450: WORK( I ) = ZERO
451: ELSE
452: WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
453: $ 2 ) )-DIFL( I ) ) /
454: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
455: END IF
456: 70 CONTINUE
457: CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
458: $ BX( J, 1 ), LDBX )
459: 80 CONTINUE
460: END IF
461: *
462: * Step (2R): if SQRE = 1, apply back the rotation that is
463: * related to the right null space of the subproblem.
464: *
465: IF( SQRE.EQ.1 ) THEN
466: CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
467: CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
468: END IF
469: IF( K.LT.MAX( M, N ) )
470: $ CALL DLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
471: $ LDBX )
472: *
473: * Step (3R): permute rows of B.
474: *
475: CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
476: IF( SQRE.EQ.1 ) THEN
477: CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
478: END IF
479: DO 90 I = 2, N
480: CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
481: 90 CONTINUE
482: *
483: * Step (4R): apply back the Givens rotations performed.
484: *
485: DO 100 I = GIVPTR, 1, -1
486: CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
487: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
488: $ -GIVNUM( I, 1 ) )
489: 100 CONTINUE
490: END IF
491: *
492: RETURN
493: *
494: * End of DLALS0
495: *
496: 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