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