1: *> \brief \b DLALSA
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLALSA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsa.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsa.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsa.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
22: * LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
23: * GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
24: * IWORK, INFO )
25: *
26: * .. Scalar Arguments ..
27: * INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
28: * $ SMLSIZ
29: * ..
30: * .. Array Arguments ..
31: * INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
32: * $ K( * ), PERM( LDGCOL, * )
33: * DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ),
34: * $ DIFL( LDU, * ), DIFR( LDU, * ),
35: * $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
36: * $ U( LDU, * ), VT( LDU, * ), WORK( * ),
37: * $ Z( LDU, * )
38: * ..
39: *
40: *
41: *> \par Purpose:
42: * =============
43: *>
44: *> \verbatim
45: *>
46: *> DLALSA is an itermediate step in solving the least squares problem
47: *> by computing the SVD of the coefficient matrix in compact form (The
48: *> singular vectors are computed as products of simple orthorgonal
49: *> matrices.).
50: *>
51: *> If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
52: *> matrix of an upper bidiagonal matrix to the right hand side; and if
53: *> ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
54: *> right hand side. The singular vector matrices were generated in
55: *> compact form by DLALSA.
56: *> \endverbatim
57: *
58: * Arguments:
59: * ==========
60: *
61: *> \param[in] ICOMPQ
62: *> \verbatim
63: *> ICOMPQ is INTEGER
64: *> Specifies whether the left or the right singular vector
65: *> matrix is involved.
66: *> = 0: Left singular vector matrix
67: *> = 1: Right singular vector matrix
68: *> \endverbatim
69: *>
70: *> \param[in] SMLSIZ
71: *> \verbatim
72: *> SMLSIZ is INTEGER
73: *> The maximum size of the subproblems at the bottom of the
74: *> computation tree.
75: *> \endverbatim
76: *>
77: *> \param[in] N
78: *> \verbatim
79: *> N is INTEGER
80: *> The row and column dimensions of the upper bidiagonal matrix.
81: *> \endverbatim
82: *>
83: *> \param[in] NRHS
84: *> \verbatim
85: *> NRHS is INTEGER
86: *> The number of columns of B and BX. NRHS must be at least 1.
87: *> \endverbatim
88: *>
89: *> \param[in,out] B
90: *> \verbatim
91: *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS )
92: *> On input, B contains the right hand sides of the least
93: *> squares problem in rows 1 through M.
94: *> On output, B contains the solution X in rows 1 through N.
95: *> \endverbatim
96: *>
97: *> \param[in] LDB
98: *> \verbatim
99: *> LDB is INTEGER
100: *> The leading dimension of B in the calling subprogram.
101: *> LDB must be at least max(1,MAX( M, N ) ).
102: *> \endverbatim
103: *>
104: *> \param[out] BX
105: *> \verbatim
106: *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS )
107: *> On exit, the result of applying the left or right singular
108: *> vector matrix to B.
109: *> \endverbatim
110: *>
111: *> \param[in] LDBX
112: *> \verbatim
113: *> LDBX is INTEGER
114: *> The leading dimension of BX.
115: *> \endverbatim
116: *>
117: *> \param[in] U
118: *> \verbatim
119: *> U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
120: *> On entry, U contains the left singular vector matrices of all
121: *> subproblems at the bottom level.
122: *> \endverbatim
123: *>
124: *> \param[in] LDU
125: *> \verbatim
126: *> LDU is INTEGER, LDU = > N.
127: *> The leading dimension of arrays U, VT, DIFL, DIFR,
128: *> POLES, GIVNUM, and Z.
129: *> \endverbatim
130: *>
131: *> \param[in] VT
132: *> \verbatim
133: *> VT is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
134: *> On entry, VT**T contains the right singular vector matrices of
135: *> all subproblems at the bottom level.
136: *> \endverbatim
137: *>
138: *> \param[in] K
139: *> \verbatim
140: *> K is INTEGER array, dimension ( N ).
141: *> \endverbatim
142: *>
143: *> \param[in] DIFL
144: *> \verbatim
145: *> DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
146: *> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
147: *> \endverbatim
148: *>
149: *> \param[in] DIFR
150: *> \verbatim
151: *> DIFR is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
152: *> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
153: *> distances between singular values on the I-th level and
154: *> singular values on the (I -1)-th level, and DIFR(*, 2 * I)
155: *> record the normalizing factors of the right singular vectors
156: *> matrices of subproblems on I-th level.
157: *> \endverbatim
158: *>
159: *> \param[in] Z
160: *> \verbatim
161: *> Z is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
162: *> On entry, Z(1, I) contains the components of the deflation-
163: *> adjusted updating row vector for subproblems on the I-th
164: *> level.
165: *> \endverbatim
166: *>
167: *> \param[in] POLES
168: *> \verbatim
169: *> POLES is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
170: *> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
171: *> singular values involved in the secular equations on the I-th
172: *> level.
173: *> \endverbatim
174: *>
175: *> \param[in] GIVPTR
176: *> \verbatim
177: *> GIVPTR is INTEGER array, dimension ( N ).
178: *> On entry, GIVPTR( I ) records the number of Givens
179: *> rotations performed on the I-th problem on the computation
180: *> tree.
181: *> \endverbatim
182: *>
183: *> \param[in] GIVCOL
184: *> \verbatim
185: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
186: *> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
187: *> locations of Givens rotations performed on the I-th level on
188: *> the computation tree.
189: *> \endverbatim
190: *>
191: *> \param[in] LDGCOL
192: *> \verbatim
193: *> LDGCOL is INTEGER, LDGCOL = > N.
194: *> The leading dimension of arrays GIVCOL and PERM.
195: *> \endverbatim
196: *>
197: *> \param[in] PERM
198: *> \verbatim
199: *> PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
200: *> On entry, PERM(*, I) records permutations done on the I-th
201: *> level of the computation tree.
202: *> \endverbatim
203: *>
204: *> \param[in] GIVNUM
205: *> \verbatim
206: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
207: *> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
208: *> values of Givens rotations performed on the I-th level on the
209: *> computation tree.
210: *> \endverbatim
211: *>
212: *> \param[in] C
213: *> \verbatim
214: *> C is DOUBLE PRECISION array, dimension ( N ).
215: *> On entry, if the I-th subproblem is not square,
216: *> C( I ) contains the C-value of a Givens rotation related to
217: *> the right null space of the I-th subproblem.
218: *> \endverbatim
219: *>
220: *> \param[in] S
221: *> \verbatim
222: *> S is DOUBLE PRECISION array, dimension ( N ).
223: *> On entry, if the I-th subproblem is not square,
224: *> S( I ) contains the S-value of a Givens rotation related to
225: *> the right null space of the I-th subproblem.
226: *> \endverbatim
227: *>
228: *> \param[out] WORK
229: *> \verbatim
230: *> WORK is DOUBLE PRECISION array.
231: *> The dimension must be at least N.
232: *> \endverbatim
233: *>
234: *> \param[out] IWORK
235: *> \verbatim
236: *> IWORK is INTEGER array.
237: *> The dimension must be at least 3 * N
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: *> \date November 2011
256: *
257: *> \ingroup doubleOTHERcomputational
258: *
259: *> \par Contributors:
260: * ==================
261: *>
262: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
263: *> California at Berkeley, USA \n
264: *> Osni Marques, LBNL/NERSC, USA \n
265: *
266: * =====================================================================
267: SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
268: $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
269: $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
270: $ IWORK, INFO )
271: *
272: * -- LAPACK computational routine (version 3.4.0) --
273: * -- LAPACK is a software package provided by Univ. of Tennessee, --
274: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
275: * November 2011
276: *
277: * .. Scalar Arguments ..
278: INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
279: $ SMLSIZ
280: * ..
281: * .. Array Arguments ..
282: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
283: $ K( * ), PERM( LDGCOL, * )
284: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ),
285: $ DIFL( LDU, * ), DIFR( LDU, * ),
286: $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
287: $ U( LDU, * ), VT( LDU, * ), WORK( * ),
288: $ Z( LDU, * )
289: * ..
290: *
291: * =====================================================================
292: *
293: * .. Parameters ..
294: DOUBLE PRECISION ZERO, ONE
295: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
296: * ..
297: * .. Local Scalars ..
298: INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
299: $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
300: $ NR, NRF, NRP1, SQRE
301: * ..
302: * .. External Subroutines ..
303: EXTERNAL DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
304: * ..
305: * .. Executable Statements ..
306: *
307: * Test the input parameters.
308: *
309: INFO = 0
310: *
311: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
312: INFO = -1
313: ELSE IF( SMLSIZ.LT.3 ) THEN
314: INFO = -2
315: ELSE IF( N.LT.SMLSIZ ) THEN
316: INFO = -3
317: ELSE IF( NRHS.LT.1 ) THEN
318: INFO = -4
319: ELSE IF( LDB.LT.N ) THEN
320: INFO = -6
321: ELSE IF( LDBX.LT.N ) THEN
322: INFO = -8
323: ELSE IF( LDU.LT.N ) THEN
324: INFO = -10
325: ELSE IF( LDGCOL.LT.N ) THEN
326: INFO = -19
327: END IF
328: IF( INFO.NE.0 ) THEN
329: CALL XERBLA( 'DLALSA', -INFO )
330: RETURN
331: END IF
332: *
333: * Book-keeping and setting up the computation tree.
334: *
335: INODE = 1
336: NDIML = INODE + N
337: NDIMR = NDIML + N
338: *
339: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
340: $ IWORK( NDIMR ), SMLSIZ )
341: *
342: * The following code applies back the left singular vector factors.
343: * For applying back the right singular vector factors, go to 50.
344: *
345: IF( ICOMPQ.EQ.1 ) THEN
346: GO TO 50
347: END IF
348: *
349: * The nodes on the bottom level of the tree were solved
350: * by DLASDQ. The corresponding left and right singular vector
351: * matrices are in explicit form. First apply back the left
352: * singular vector matrices.
353: *
354: NDB1 = ( ND+1 ) / 2
355: DO 10 I = NDB1, ND
356: *
357: * IC : center row of each node
358: * NL : number of rows of left subproblem
359: * NR : number of rows of right subproblem
360: * NLF: starting row of the left subproblem
361: * NRF: starting row of the right subproblem
362: *
363: I1 = I - 1
364: IC = IWORK( INODE+I1 )
365: NL = IWORK( NDIML+I1 )
366: NR = IWORK( NDIMR+I1 )
367: NLF = IC - NL
368: NRF = IC + 1
369: CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
370: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
371: CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
372: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
373: 10 CONTINUE
374: *
375: * Next copy the rows of B that correspond to unchanged rows
376: * in the bidiagonal matrix to BX.
377: *
378: DO 20 I = 1, ND
379: IC = IWORK( INODE+I-1 )
380: CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
381: 20 CONTINUE
382: *
383: * Finally go through the left singular vector matrices of all
384: * the other subproblems bottom-up on the tree.
385: *
386: J = 2**NLVL
387: SQRE = 0
388: *
389: DO 40 LVL = NLVL, 1, -1
390: LVL2 = 2*LVL - 1
391: *
392: * find the first node LF and last node LL on
393: * the current level LVL
394: *
395: IF( LVL.EQ.1 ) THEN
396: LF = 1
397: LL = 1
398: ELSE
399: LF = 2**( LVL-1 )
400: LL = 2*LF - 1
401: END IF
402: DO 30 I = LF, LL
403: IM1 = I - 1
404: IC = IWORK( INODE+IM1 )
405: NL = IWORK( NDIML+IM1 )
406: NR = IWORK( NDIMR+IM1 )
407: NLF = IC - NL
408: NRF = IC + 1
409: J = J - 1
410: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
411: $ B( NLF, 1 ), LDB, PERM( NLF, LVL ),
412: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
413: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
414: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
415: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
416: $ INFO )
417: 30 CONTINUE
418: 40 CONTINUE
419: GO TO 90
420: *
421: * ICOMPQ = 1: applying back the right singular vector factors.
422: *
423: 50 CONTINUE
424: *
425: * First now go through the right singular vector matrices of all
426: * the tree nodes top-down.
427: *
428: J = 0
429: DO 70 LVL = 1, NLVL
430: LVL2 = 2*LVL - 1
431: *
432: * Find the first node LF and last node LL on
433: * the current level LVL.
434: *
435: IF( LVL.EQ.1 ) THEN
436: LF = 1
437: LL = 1
438: ELSE
439: LF = 2**( LVL-1 )
440: LL = 2*LF - 1
441: END IF
442: DO 60 I = LL, LF, -1
443: IM1 = I - 1
444: IC = IWORK( INODE+IM1 )
445: NL = IWORK( NDIML+IM1 )
446: NR = IWORK( NDIMR+IM1 )
447: NLF = IC - NL
448: NRF = IC + 1
449: IF( I.EQ.LL ) THEN
450: SQRE = 0
451: ELSE
452: SQRE = 1
453: END IF
454: J = J + 1
455: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
456: $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
457: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
458: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
459: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
460: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
461: $ INFO )
462: 60 CONTINUE
463: 70 CONTINUE
464: *
465: * The nodes on the bottom level of the tree were solved
466: * by DLASDQ. The corresponding right singular vector
467: * matrices are in explicit form. Apply them back.
468: *
469: NDB1 = ( ND+1 ) / 2
470: DO 80 I = NDB1, ND
471: I1 = I - 1
472: IC = IWORK( INODE+I1 )
473: NL = IWORK( NDIML+I1 )
474: NR = IWORK( NDIMR+I1 )
475: NLP1 = NL + 1
476: IF( I.EQ.ND ) THEN
477: NRP1 = NR
478: ELSE
479: NRP1 = NR + 1
480: END IF
481: NLF = IC - NL
482: NRF = IC + 1
483: CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
484: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
485: CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
486: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
487: 80 CONTINUE
488: *
489: 90 CONTINUE
490: *
491: RETURN
492: *
493: * End of DLALSA
494: *
495: END
CVSweb interface <joel.bertrand@systella.fr>