Annotation of rpl/lapack/lapack/dlalsa.f, revision 1.9
1.9 ! bertrand 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: * =====================================================================
1.1 bertrand 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: *
1.9 ! bertrand 272: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 273: * -- LAPACK is a software package provided by Univ. of Tennessee, --
274: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 275: * November 2011
1.1 bertrand 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>
![[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