1: *> \brief \b DLALSA computes the SVD of the coefficient matrix in compact form. 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 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, dimension (N)
231: *> \endverbatim
232: *>
233: *> \param[out] IWORK
234: *> \verbatim
235: *> IWORK is INTEGER array, dimension (3*N)
236: *> \endverbatim
237: *>
238: *> \param[out] INFO
239: *> \verbatim
240: *> INFO is INTEGER
241: *> = 0: successful exit.
242: *> < 0: if INFO = -i, the i-th argument had an illegal value.
243: *> \endverbatim
244: *
245: * Authors:
246: * ========
247: *
248: *> \author Univ. of Tennessee
249: *> \author Univ. of California Berkeley
250: *> \author Univ. of Colorado Denver
251: *> \author NAG Ltd.
252: *
253: *> \ingroup doubleOTHERcomputational
254: *
255: *> \par Contributors:
256: * ==================
257: *>
258: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
259: *> California at Berkeley, USA \n
260: *> Osni Marques, LBNL/NERSC, USA \n
261: *
262: * =====================================================================
263: SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
264: $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
265: $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
266: $ IWORK, INFO )
267: *
268: * -- LAPACK computational routine --
269: * -- LAPACK is a software package provided by Univ. of Tennessee, --
270: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271: *
272: * .. Scalar Arguments ..
273: INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
274: $ SMLSIZ
275: * ..
276: * .. Array Arguments ..
277: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
278: $ K( * ), PERM( LDGCOL, * )
279: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ),
280: $ DIFL( LDU, * ), DIFR( LDU, * ),
281: $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
282: $ U( LDU, * ), VT( LDU, * ), WORK( * ),
283: $ Z( LDU, * )
284: * ..
285: *
286: * =====================================================================
287: *
288: * .. Parameters ..
289: DOUBLE PRECISION ZERO, ONE
290: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
291: * ..
292: * .. Local Scalars ..
293: INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
294: $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
295: $ NR, NRF, NRP1, SQRE
296: * ..
297: * .. External Subroutines ..
298: EXTERNAL DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
299: * ..
300: * .. Executable Statements ..
301: *
302: * Test the input parameters.
303: *
304: INFO = 0
305: *
306: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
307: INFO = -1
308: ELSE IF( SMLSIZ.LT.3 ) THEN
309: INFO = -2
310: ELSE IF( N.LT.SMLSIZ ) THEN
311: INFO = -3
312: ELSE IF( NRHS.LT.1 ) THEN
313: INFO = -4
314: ELSE IF( LDB.LT.N ) THEN
315: INFO = -6
316: ELSE IF( LDBX.LT.N ) THEN
317: INFO = -8
318: ELSE IF( LDU.LT.N ) THEN
319: INFO = -10
320: ELSE IF( LDGCOL.LT.N ) THEN
321: INFO = -19
322: END IF
323: IF( INFO.NE.0 ) THEN
324: CALL XERBLA( 'DLALSA', -INFO )
325: RETURN
326: END IF
327: *
328: * Book-keeping and setting up the computation tree.
329: *
330: INODE = 1
331: NDIML = INODE + N
332: NDIMR = NDIML + N
333: *
334: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
335: $ IWORK( NDIMR ), SMLSIZ )
336: *
337: * The following code applies back the left singular vector factors.
338: * For applying back the right singular vector factors, go to 50.
339: *
340: IF( ICOMPQ.EQ.1 ) THEN
341: GO TO 50
342: END IF
343: *
344: * The nodes on the bottom level of the tree were solved
345: * by DLASDQ. The corresponding left and right singular vector
346: * matrices are in explicit form. First apply back the left
347: * singular vector matrices.
348: *
349: NDB1 = ( ND+1 ) / 2
350: DO 10 I = NDB1, ND
351: *
352: * IC : center row of each node
353: * NL : number of rows of left subproblem
354: * NR : number of rows of right subproblem
355: * NLF: starting row of the left subproblem
356: * NRF: starting row of the right subproblem
357: *
358: I1 = I - 1
359: IC = IWORK( INODE+I1 )
360: NL = IWORK( NDIML+I1 )
361: NR = IWORK( NDIMR+I1 )
362: NLF = IC - NL
363: NRF = IC + 1
364: CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
365: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
366: CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
367: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
368: 10 CONTINUE
369: *
370: * Next copy the rows of B that correspond to unchanged rows
371: * in the bidiagonal matrix to BX.
372: *
373: DO 20 I = 1, ND
374: IC = IWORK( INODE+I-1 )
375: CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
376: 20 CONTINUE
377: *
378: * Finally go through the left singular vector matrices of all
379: * the other subproblems bottom-up on the tree.
380: *
381: J = 2**NLVL
382: SQRE = 0
383: *
384: DO 40 LVL = NLVL, 1, -1
385: LVL2 = 2*LVL - 1
386: *
387: * find the first node LF and last node LL on
388: * the current level LVL
389: *
390: IF( LVL.EQ.1 ) THEN
391: LF = 1
392: LL = 1
393: ELSE
394: LF = 2**( LVL-1 )
395: LL = 2*LF - 1
396: END IF
397: DO 30 I = LF, LL
398: IM1 = I - 1
399: IC = IWORK( INODE+IM1 )
400: NL = IWORK( NDIML+IM1 )
401: NR = IWORK( NDIMR+IM1 )
402: NLF = IC - NL
403: NRF = IC + 1
404: J = J - 1
405: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
406: $ B( NLF, 1 ), LDB, PERM( NLF, LVL ),
407: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
408: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
409: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
410: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
411: $ INFO )
412: 30 CONTINUE
413: 40 CONTINUE
414: GO TO 90
415: *
416: * ICOMPQ = 1: applying back the right singular vector factors.
417: *
418: 50 CONTINUE
419: *
420: * First now go through the right singular vector matrices of all
421: * the tree nodes top-down.
422: *
423: J = 0
424: DO 70 LVL = 1, NLVL
425: LVL2 = 2*LVL - 1
426: *
427: * Find the first node LF and last node LL on
428: * the current level LVL.
429: *
430: IF( LVL.EQ.1 ) THEN
431: LF = 1
432: LL = 1
433: ELSE
434: LF = 2**( LVL-1 )
435: LL = 2*LF - 1
436: END IF
437: DO 60 I = LL, LF, -1
438: IM1 = I - 1
439: IC = IWORK( INODE+IM1 )
440: NL = IWORK( NDIML+IM1 )
441: NR = IWORK( NDIMR+IM1 )
442: NLF = IC - NL
443: NRF = IC + 1
444: IF( I.EQ.LL ) THEN
445: SQRE = 0
446: ELSE
447: SQRE = 1
448: END IF
449: J = J + 1
450: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
451: $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
452: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
453: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
454: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
455: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
456: $ INFO )
457: 60 CONTINUE
458: 70 CONTINUE
459: *
460: * The nodes on the bottom level of the tree were solved
461: * by DLASDQ. The corresponding right singular vector
462: * matrices are in explicit form. Apply them back.
463: *
464: NDB1 = ( ND+1 ) / 2
465: DO 80 I = NDB1, ND
466: I1 = I - 1
467: IC = IWORK( INODE+I1 )
468: NL = IWORK( NDIML+I1 )
469: NR = IWORK( NDIMR+I1 )
470: NLP1 = NL + 1
471: IF( I.EQ.ND ) THEN
472: NRP1 = NR
473: ELSE
474: NRP1 = NR + 1
475: END IF
476: NLF = IC - NL
477: NRF = IC + 1
478: CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
479: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
480: CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
481: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
482: 80 CONTINUE
483: *
484: 90 CONTINUE
485: *
486: RETURN
487: *
488: * End of DLALSA
489: *
490: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlalsa.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack