1: SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
2: $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
3: $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
4: $ IWORK, INFO )
5: *
6: * -- LAPACK routine (version 3.2) --
7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9: * November 2006
10: *
11: * .. Scalar Arguments ..
12: INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
13: $ SMLSIZ
14: * ..
15: * .. Array Arguments ..
16: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
17: $ K( * ), PERM( LDGCOL, * )
18: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ),
19: $ DIFL( LDU, * ), DIFR( LDU, * ),
20: $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
21: $ U( LDU, * ), VT( LDU, * ), WORK( * ),
22: $ Z( LDU, * )
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * DLALSA is an itermediate step in solving the least squares problem
29: * by computing the SVD of the coefficient matrix in compact form (The
30: * singular vectors are computed as products of simple orthorgonal
31: * matrices.).
32: *
33: * If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
34: * matrix of an upper bidiagonal matrix to the right hand side; and if
35: * ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
36: * right hand side. The singular vector matrices were generated in
37: * compact form by DLALSA.
38: *
39: * Arguments
40: * =========
41: *
42: *
43: * ICOMPQ (input) INTEGER
44: * Specifies whether the left or the right singular vector
45: * matrix is involved.
46: * = 0: Left singular vector matrix
47: * = 1: Right singular vector matrix
48: *
49: * SMLSIZ (input) INTEGER
50: * The maximum size of the subproblems at the bottom of the
51: * computation tree.
52: *
53: * N (input) INTEGER
54: * The row and column dimensions of the upper bidiagonal matrix.
55: *
56: * NRHS (input) INTEGER
57: * The number of columns of B and BX. NRHS must be at least 1.
58: *
59: * B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
60: * On input, B contains the right hand sides of the least
61: * squares problem in rows 1 through M.
62: * On output, B contains the solution X in rows 1 through N.
63: *
64: * LDB (input) INTEGER
65: * The leading dimension of B in the calling subprogram.
66: * LDB must be at least max(1,MAX( M, N ) ).
67: *
68: * BX (output) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
69: * On exit, the result of applying the left or right singular
70: * vector matrix to B.
71: *
72: * LDBX (input) INTEGER
73: * The leading dimension of BX.
74: *
75: * U (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
76: * On entry, U contains the left singular vector matrices of all
77: * subproblems at the bottom level.
78: *
79: * LDU (input) INTEGER, LDU = > N.
80: * The leading dimension of arrays U, VT, DIFL, DIFR,
81: * POLES, GIVNUM, and Z.
82: *
83: * VT (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
84: * On entry, VT' contains the right singular vector matrices of
85: * all subproblems at the bottom level.
86: *
87: * K (input) INTEGER array, dimension ( N ).
88: *
89: * DIFL (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
90: * where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
91: *
92: * DIFR (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
93: * On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
94: * distances between singular values on the I-th level and
95: * singular values on the (I -1)-th level, and DIFR(*, 2 * I)
96: * record the normalizing factors of the right singular vectors
97: * matrices of subproblems on I-th level.
98: *
99: * Z (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
100: * On entry, Z(1, I) contains the components of the deflation-
101: * adjusted updating row vector for subproblems on the I-th
102: * level.
103: *
104: * POLES (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
105: * On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
106: * singular values involved in the secular equations on the I-th
107: * level.
108: *
109: * GIVPTR (input) INTEGER array, dimension ( N ).
110: * On entry, GIVPTR( I ) records the number of Givens
111: * rotations performed on the I-th problem on the computation
112: * tree.
113: *
114: * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
115: * On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
116: * locations of Givens rotations performed on the I-th level on
117: * the computation tree.
118: *
119: * LDGCOL (input) INTEGER, LDGCOL = > N.
120: * The leading dimension of arrays GIVCOL and PERM.
121: *
122: * PERM (input) INTEGER array, dimension ( LDGCOL, NLVL ).
123: * On entry, PERM(*, I) records permutations done on the I-th
124: * level of the computation tree.
125: *
126: * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
127: * On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
128: * values of Givens rotations performed on the I-th level on the
129: * computation tree.
130: *
131: * C (input) DOUBLE PRECISION array, dimension ( N ).
132: * On entry, if the I-th subproblem is not square,
133: * C( I ) contains the C-value of a Givens rotation related to
134: * the right null space of the I-th subproblem.
135: *
136: * S (input) DOUBLE PRECISION array, dimension ( N ).
137: * On entry, if the I-th subproblem is not square,
138: * S( I ) contains the S-value of a Givens rotation related to
139: * the right null space of the I-th subproblem.
140: *
141: * WORK (workspace) DOUBLE PRECISION array.
142: * The dimension must be at least N.
143: *
144: * IWORK (workspace) INTEGER array.
145: * The dimension must be at least 3 * N
146: *
147: * INFO (output) INTEGER
148: * = 0: successful exit.
149: * < 0: if INFO = -i, the i-th argument had an illegal value.
150: *
151: * Further Details
152: * ===============
153: *
154: * Based on contributions by
155: * Ming Gu and Ren-Cang Li, Computer Science Division, University of
156: * California at Berkeley, USA
157: * Osni Marques, LBNL/NERSC, USA
158: *
159: * =====================================================================
160: *
161: * .. Parameters ..
162: DOUBLE PRECISION ZERO, ONE
163: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
164: * ..
165: * .. Local Scalars ..
166: INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
167: $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
168: $ NR, NRF, NRP1, SQRE
169: * ..
170: * .. External Subroutines ..
171: EXTERNAL DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
172: * ..
173: * .. Executable Statements ..
174: *
175: * Test the input parameters.
176: *
177: INFO = 0
178: *
179: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
180: INFO = -1
181: ELSE IF( SMLSIZ.LT.3 ) THEN
182: INFO = -2
183: ELSE IF( N.LT.SMLSIZ ) THEN
184: INFO = -3
185: ELSE IF( NRHS.LT.1 ) THEN
186: INFO = -4
187: ELSE IF( LDB.LT.N ) THEN
188: INFO = -6
189: ELSE IF( LDBX.LT.N ) THEN
190: INFO = -8
191: ELSE IF( LDU.LT.N ) THEN
192: INFO = -10
193: ELSE IF( LDGCOL.LT.N ) THEN
194: INFO = -19
195: END IF
196: IF( INFO.NE.0 ) THEN
197: CALL XERBLA( 'DLALSA', -INFO )
198: RETURN
199: END IF
200: *
201: * Book-keeping and setting up the computation tree.
202: *
203: INODE = 1
204: NDIML = INODE + N
205: NDIMR = NDIML + N
206: *
207: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
208: $ IWORK( NDIMR ), SMLSIZ )
209: *
210: * The following code applies back the left singular vector factors.
211: * For applying back the right singular vector factors, go to 50.
212: *
213: IF( ICOMPQ.EQ.1 ) THEN
214: GO TO 50
215: END IF
216: *
217: * The nodes on the bottom level of the tree were solved
218: * by DLASDQ. The corresponding left and right singular vector
219: * matrices are in explicit form. First apply back the left
220: * singular vector matrices.
221: *
222: NDB1 = ( ND+1 ) / 2
223: DO 10 I = NDB1, ND
224: *
225: * IC : center row of each node
226: * NL : number of rows of left subproblem
227: * NR : number of rows of right subproblem
228: * NLF: starting row of the left subproblem
229: * NRF: starting row of the right subproblem
230: *
231: I1 = I - 1
232: IC = IWORK( INODE+I1 )
233: NL = IWORK( NDIML+I1 )
234: NR = IWORK( NDIMR+I1 )
235: NLF = IC - NL
236: NRF = IC + 1
237: CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
238: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
239: CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
240: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
241: 10 CONTINUE
242: *
243: * Next copy the rows of B that correspond to unchanged rows
244: * in the bidiagonal matrix to BX.
245: *
246: DO 20 I = 1, ND
247: IC = IWORK( INODE+I-1 )
248: CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
249: 20 CONTINUE
250: *
251: * Finally go through the left singular vector matrices of all
252: * the other subproblems bottom-up on the tree.
253: *
254: J = 2**NLVL
255: SQRE = 0
256: *
257: DO 40 LVL = NLVL, 1, -1
258: LVL2 = 2*LVL - 1
259: *
260: * find the first node LF and last node LL on
261: * the current level LVL
262: *
263: IF( LVL.EQ.1 ) THEN
264: LF = 1
265: LL = 1
266: ELSE
267: LF = 2**( LVL-1 )
268: LL = 2*LF - 1
269: END IF
270: DO 30 I = LF, LL
271: IM1 = I - 1
272: IC = IWORK( INODE+IM1 )
273: NL = IWORK( NDIML+IM1 )
274: NR = IWORK( NDIMR+IM1 )
275: NLF = IC - NL
276: NRF = IC + 1
277: J = J - 1
278: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
279: $ B( NLF, 1 ), LDB, PERM( NLF, LVL ),
280: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
281: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
282: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
283: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
284: $ INFO )
285: 30 CONTINUE
286: 40 CONTINUE
287: GO TO 90
288: *
289: * ICOMPQ = 1: applying back the right singular vector factors.
290: *
291: 50 CONTINUE
292: *
293: * First now go through the right singular vector matrices of all
294: * the tree nodes top-down.
295: *
296: J = 0
297: DO 70 LVL = 1, NLVL
298: LVL2 = 2*LVL - 1
299: *
300: * Find the first node LF and last node LL on
301: * the current level LVL.
302: *
303: IF( LVL.EQ.1 ) THEN
304: LF = 1
305: LL = 1
306: ELSE
307: LF = 2**( LVL-1 )
308: LL = 2*LF - 1
309: END IF
310: DO 60 I = LL, LF, -1
311: IM1 = I - 1
312: IC = IWORK( INODE+IM1 )
313: NL = IWORK( NDIML+IM1 )
314: NR = IWORK( NDIMR+IM1 )
315: NLF = IC - NL
316: NRF = IC + 1
317: IF( I.EQ.LL ) THEN
318: SQRE = 0
319: ELSE
320: SQRE = 1
321: END IF
322: J = J + 1
323: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
324: $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
325: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
326: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
327: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
328: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
329: $ INFO )
330: 60 CONTINUE
331: 70 CONTINUE
332: *
333: * The nodes on the bottom level of the tree were solved
334: * by DLASDQ. The corresponding right singular vector
335: * matrices are in explicit form. Apply them back.
336: *
337: NDB1 = ( ND+1 ) / 2
338: DO 80 I = NDB1, ND
339: I1 = I - 1
340: IC = IWORK( INODE+I1 )
341: NL = IWORK( NDIML+I1 )
342: NR = IWORK( NDIMR+I1 )
343: NLP1 = NL + 1
344: IF( I.EQ.ND ) THEN
345: NRP1 = NR
346: ELSE
347: NRP1 = NR + 1
348: END IF
349: NLF = IC - NL
350: NRF = IC + 1
351: CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
352: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
353: CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
354: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
355: 80 CONTINUE
356: *
357: 90 CONTINUE
358: *
359: RETURN
360: *
361: * End of DLALSA
362: *
363: 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