Annotation of rpl/lapack/lapack/dlalsa.f, revision 1.1
1.1 ! bertrand 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