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: *> \date June 2017
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 DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
266: $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
267: $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
268: $ IWORK, INFO )
269: *
270: * -- LAPACK computational routine (version 3.7.1) --
271: * -- LAPACK is a software package provided by Univ. of Tennessee, --
272: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
273: * June 2017
274: *
275: * .. Scalar Arguments ..
276: INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
277: $ SMLSIZ
278: * ..
279: * .. Array Arguments ..
280: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
281: $ K( * ), PERM( LDGCOL, * )
282: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ),
283: $ DIFL( LDU, * ), DIFR( LDU, * ),
284: $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
285: $ U( LDU, * ), VT( LDU, * ), WORK( * ),
286: $ Z( LDU, * )
287: * ..
288: *
289: * =====================================================================
290: *
291: * .. Parameters ..
292: DOUBLE PRECISION ZERO, ONE
293: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
294: * ..
295: * .. Local Scalars ..
296: INTEGER I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
297: $ ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
298: $ NR, NRF, NRP1, SQRE
299: * ..
300: * .. External Subroutines ..
301: EXTERNAL DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
302: * ..
303: * .. Executable Statements ..
304: *
305: * Test the input parameters.
306: *
307: INFO = 0
308: *
309: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
310: INFO = -1
311: ELSE IF( SMLSIZ.LT.3 ) THEN
312: INFO = -2
313: ELSE IF( N.LT.SMLSIZ ) THEN
314: INFO = -3
315: ELSE IF( NRHS.LT.1 ) THEN
316: INFO = -4
317: ELSE IF( LDB.LT.N ) THEN
318: INFO = -6
319: ELSE IF( LDBX.LT.N ) THEN
320: INFO = -8
321: ELSE IF( LDU.LT.N ) THEN
322: INFO = -10
323: ELSE IF( LDGCOL.LT.N ) THEN
324: INFO = -19
325: END IF
326: IF( INFO.NE.0 ) THEN
327: CALL XERBLA( 'DLALSA', -INFO )
328: RETURN
329: END IF
330: *
331: * Book-keeping and setting up the computation tree.
332: *
333: INODE = 1
334: NDIML = INODE + N
335: NDIMR = NDIML + N
336: *
337: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
338: $ IWORK( NDIMR ), SMLSIZ )
339: *
340: * The following code applies back the left singular vector factors.
341: * For applying back the right singular vector factors, go to 50.
342: *
343: IF( ICOMPQ.EQ.1 ) THEN
344: GO TO 50
345: END IF
346: *
347: * The nodes on the bottom level of the tree were solved
348: * by DLASDQ. The corresponding left and right singular vector
349: * matrices are in explicit form. First apply back the left
350: * singular vector matrices.
351: *
352: NDB1 = ( ND+1 ) / 2
353: DO 10 I = NDB1, ND
354: *
355: * IC : center row of each node
356: * NL : number of rows of left subproblem
357: * NR : number of rows of right subproblem
358: * NLF: starting row of the left subproblem
359: * NRF: starting row of the right subproblem
360: *
361: I1 = I - 1
362: IC = IWORK( INODE+I1 )
363: NL = IWORK( NDIML+I1 )
364: NR = IWORK( NDIMR+I1 )
365: NLF = IC - NL
366: NRF = IC + 1
367: CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
368: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
369: CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
370: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
371: 10 CONTINUE
372: *
373: * Next copy the rows of B that correspond to unchanged rows
374: * in the bidiagonal matrix to BX.
375: *
376: DO 20 I = 1, ND
377: IC = IWORK( INODE+I-1 )
378: CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
379: 20 CONTINUE
380: *
381: * Finally go through the left singular vector matrices of all
382: * the other subproblems bottom-up on the tree.
383: *
384: J = 2**NLVL
385: SQRE = 0
386: *
387: DO 40 LVL = NLVL, 1, -1
388: LVL2 = 2*LVL - 1
389: *
390: * find the first node LF and last node LL on
391: * the current level LVL
392: *
393: IF( LVL.EQ.1 ) THEN
394: LF = 1
395: LL = 1
396: ELSE
397: LF = 2**( LVL-1 )
398: LL = 2*LF - 1
399: END IF
400: DO 30 I = LF, LL
401: IM1 = I - 1
402: IC = IWORK( INODE+IM1 )
403: NL = IWORK( NDIML+IM1 )
404: NR = IWORK( NDIMR+IM1 )
405: NLF = IC - NL
406: NRF = IC + 1
407: J = J - 1
408: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
409: $ B( NLF, 1 ), LDB, PERM( NLF, LVL ),
410: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
411: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
412: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
413: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
414: $ INFO )
415: 30 CONTINUE
416: 40 CONTINUE
417: GO TO 90
418: *
419: * ICOMPQ = 1: applying back the right singular vector factors.
420: *
421: 50 CONTINUE
422: *
423: * First now go through the right singular vector matrices of all
424: * the tree nodes top-down.
425: *
426: J = 0
427: DO 70 LVL = 1, NLVL
428: LVL2 = 2*LVL - 1
429: *
430: * Find the first node LF and last node LL on
431: * the current level LVL.
432: *
433: IF( LVL.EQ.1 ) THEN
434: LF = 1
435: LL = 1
436: ELSE
437: LF = 2**( LVL-1 )
438: LL = 2*LF - 1
439: END IF
440: DO 60 I = LL, LF, -1
441: IM1 = I - 1
442: IC = IWORK( INODE+IM1 )
443: NL = IWORK( NDIML+IM1 )
444: NR = IWORK( NDIMR+IM1 )
445: NLF = IC - NL
446: NRF = IC + 1
447: IF( I.EQ.LL ) THEN
448: SQRE = 0
449: ELSE
450: SQRE = 1
451: END IF
452: J = J + 1
453: CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
454: $ BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
455: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
456: $ GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
457: $ DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
458: $ Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
459: $ INFO )
460: 60 CONTINUE
461: 70 CONTINUE
462: *
463: * The nodes on the bottom level of the tree were solved
464: * by DLASDQ. The corresponding right singular vector
465: * matrices are in explicit form. Apply them back.
466: *
467: NDB1 = ( ND+1 ) / 2
468: DO 80 I = NDB1, ND
469: I1 = I - 1
470: IC = IWORK( INODE+I1 )
471: NL = IWORK( NDIML+I1 )
472: NR = IWORK( NDIMR+I1 )
473: NLP1 = NL + 1
474: IF( I.EQ.ND ) THEN
475: NRP1 = NR
476: ELSE
477: NRP1 = NR + 1
478: END IF
479: NLF = IC - NL
480: NRF = IC + 1
481: CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
482: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
483: CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
484: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
485: 80 CONTINUE
486: *
487: 90 CONTINUE
488: *
489: RETURN
490: *
491: * End of DLALSA
492: *
493: 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