Annotation of rpl/lapack/lapack/dlals0.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLALS0
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLALS0 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlals0.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlals0.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlals0.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
! 22: * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
! 23: * POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
! 27: * $ LDGNUM, NL, NR, NRHS, SQRE
! 28: * DOUBLE PRECISION C, S
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * INTEGER GIVCOL( LDGCOL, * ), PERM( * )
! 32: * DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ),
! 33: * $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
! 34: * $ POLES( LDGNUM, * ), WORK( * ), Z( * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: *> DLALS0 applies back the multiplying factors of either the left or the
! 44: *> right singular vector matrix of a diagonal matrix appended by a row
! 45: *> to the right hand side matrix B in solving the least squares problem
! 46: *> using the divide-and-conquer SVD approach.
! 47: *>
! 48: *> For the left singular vector matrix, three types of orthogonal
! 49: *> matrices are involved:
! 50: *>
! 51: *> (1L) Givens rotations: the number of such rotations is GIVPTR; the
! 52: *> pairs of columns/rows they were applied to are stored in GIVCOL;
! 53: *> and the C- and S-values of these rotations are stored in GIVNUM.
! 54: *>
! 55: *> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
! 56: *> row, and for J=2:N, PERM(J)-th row of B is to be moved to the
! 57: *> J-th row.
! 58: *>
! 59: *> (3L) The left singular vector matrix of the remaining matrix.
! 60: *>
! 61: *> For the right singular vector matrix, four types of orthogonal
! 62: *> matrices are involved:
! 63: *>
! 64: *> (1R) The right singular vector matrix of the remaining matrix.
! 65: *>
! 66: *> (2R) If SQRE = 1, one extra Givens rotation to generate the right
! 67: *> null space.
! 68: *>
! 69: *> (3R) The inverse transformation of (2L).
! 70: *>
! 71: *> (4R) The inverse transformation of (1L).
! 72: *> \endverbatim
! 73: *
! 74: * Arguments:
! 75: * ==========
! 76: *
! 77: *> \param[in] ICOMPQ
! 78: *> \verbatim
! 79: *> ICOMPQ is INTEGER
! 80: *> Specifies whether singular vectors are to be computed in
! 81: *> factored form:
! 82: *> = 0: Left singular vector matrix.
! 83: *> = 1: Right singular vector matrix.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[in] NL
! 87: *> \verbatim
! 88: *> NL is INTEGER
! 89: *> The row dimension of the upper block. NL >= 1.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] NR
! 93: *> \verbatim
! 94: *> NR is INTEGER
! 95: *> The row dimension of the lower block. NR >= 1.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in] SQRE
! 99: *> \verbatim
! 100: *> SQRE is INTEGER
! 101: *> = 0: the lower block is an NR-by-NR square matrix.
! 102: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
! 103: *>
! 104: *> The bidiagonal matrix has row dimension N = NL + NR + 1,
! 105: *> and column dimension M = N + SQRE.
! 106: *> \endverbatim
! 107: *>
! 108: *> \param[in] NRHS
! 109: *> \verbatim
! 110: *> NRHS is INTEGER
! 111: *> The number of columns of B and BX. NRHS must be at least 1.
! 112: *> \endverbatim
! 113: *>
! 114: *> \param[in,out] B
! 115: *> \verbatim
! 116: *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS )
! 117: *> On input, B contains the right hand sides of the least
! 118: *> squares problem in rows 1 through M. On output, B contains
! 119: *> the solution X in rows 1 through N.
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[in] LDB
! 123: *> \verbatim
! 124: *> LDB is INTEGER
! 125: *> The leading dimension of B. LDB must be at least
! 126: *> max(1,MAX( M, N ) ).
! 127: *> \endverbatim
! 128: *>
! 129: *> \param[out] BX
! 130: *> \verbatim
! 131: *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS )
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[in] LDBX
! 135: *> \verbatim
! 136: *> LDBX is INTEGER
! 137: *> The leading dimension of BX.
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[in] PERM
! 141: *> \verbatim
! 142: *> PERM is INTEGER array, dimension ( N )
! 143: *> The permutations (from deflation and sorting) applied
! 144: *> to the two blocks.
! 145: *> \endverbatim
! 146: *>
! 147: *> \param[in] GIVPTR
! 148: *> \verbatim
! 149: *> GIVPTR is INTEGER
! 150: *> The number of Givens rotations which took place in this
! 151: *> subproblem.
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[in] GIVCOL
! 155: *> \verbatim
! 156: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
! 157: *> Each pair of numbers indicates a pair of rows/columns
! 158: *> involved in a Givens rotation.
! 159: *> \endverbatim
! 160: *>
! 161: *> \param[in] LDGCOL
! 162: *> \verbatim
! 163: *> LDGCOL is INTEGER
! 164: *> The leading dimension of GIVCOL, must be at least N.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[in] GIVNUM
! 168: *> \verbatim
! 169: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 170: *> Each number indicates the C or S value used in the
! 171: *> corresponding Givens rotation.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[in] LDGNUM
! 175: *> \verbatim
! 176: *> LDGNUM is INTEGER
! 177: *> The leading dimension of arrays DIFR, POLES and
! 178: *> GIVNUM, must be at least K.
! 179: *> \endverbatim
! 180: *>
! 181: *> \param[in] POLES
! 182: *> \verbatim
! 183: *> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 184: *> On entry, POLES(1:K, 1) contains the new singular
! 185: *> values obtained from solving the secular equation, and
! 186: *> POLES(1:K, 2) is an array containing the poles in the secular
! 187: *> equation.
! 188: *> \endverbatim
! 189: *>
! 190: *> \param[in] DIFL
! 191: *> \verbatim
! 192: *> DIFL is DOUBLE PRECISION array, dimension ( K ).
! 193: *> On entry, DIFL(I) is the distance between I-th updated
! 194: *> (undeflated) singular value and the I-th (undeflated) old
! 195: *> singular value.
! 196: *> \endverbatim
! 197: *>
! 198: *> \param[in] DIFR
! 199: *> \verbatim
! 200: *> DIFR is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
! 201: *> On entry, DIFR(I, 1) contains the distances between I-th
! 202: *> updated (undeflated) singular value and the I+1-th
! 203: *> (undeflated) old singular value. And DIFR(I, 2) is the
! 204: *> normalizing factor for the I-th right singular vector.
! 205: *> \endverbatim
! 206: *>
! 207: *> \param[in] Z
! 208: *> \verbatim
! 209: *> Z is DOUBLE PRECISION array, dimension ( K )
! 210: *> Contain the components of the deflation-adjusted updating row
! 211: *> vector.
! 212: *> \endverbatim
! 213: *>
! 214: *> \param[in] K
! 215: *> \verbatim
! 216: *> K is INTEGER
! 217: *> Contains the dimension of the non-deflated matrix,
! 218: *> This is the order of the related secular equation. 1 <= K <=N.
! 219: *> \endverbatim
! 220: *>
! 221: *> \param[in] C
! 222: *> \verbatim
! 223: *> C is DOUBLE PRECISION
! 224: *> C contains garbage if SQRE =0 and the C-value of a Givens
! 225: *> rotation related to the right null space if SQRE = 1.
! 226: *> \endverbatim
! 227: *>
! 228: *> \param[in] S
! 229: *> \verbatim
! 230: *> S is DOUBLE PRECISION
! 231: *> S contains garbage if SQRE =0 and the S-value of a Givens
! 232: *> rotation related to the right null space if SQRE = 1.
! 233: *> \endverbatim
! 234: *>
! 235: *> \param[out] WORK
! 236: *> \verbatim
! 237: *> WORK is DOUBLE PRECISION array, dimension ( K )
! 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 DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
268: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
269: $ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
270: *
1.8 ! bertrand 271: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 272: * -- LAPACK is a software package provided by Univ. of Tennessee, --
273: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 274: * November 2011
1.1 bertrand 275: *
276: * .. Scalar Arguments ..
277: INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
278: $ LDGNUM, NL, NR, NRHS, SQRE
279: DOUBLE PRECISION C, S
280: * ..
281: * .. Array Arguments ..
282: INTEGER GIVCOL( LDGCOL, * ), PERM( * )
283: DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ),
284: $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
285: $ POLES( LDGNUM, * ), WORK( * ), Z( * )
286: * ..
287: *
288: * =====================================================================
289: *
290: * .. Parameters ..
291: DOUBLE PRECISION ONE, ZERO, NEGONE
292: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
293: * ..
294: * .. Local Scalars ..
295: INTEGER I, J, M, N, NLP1
296: DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
297: * ..
298: * .. External Subroutines ..
299: EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DROT, DSCAL,
300: $ XERBLA
301: * ..
302: * .. External Functions ..
303: DOUBLE PRECISION DLAMC3, DNRM2
304: EXTERNAL DLAMC3, DNRM2
305: * ..
306: * .. Intrinsic Functions ..
307: INTRINSIC MAX
308: * ..
309: * .. Executable Statements ..
310: *
311: * Test the input parameters.
312: *
313: INFO = 0
314: *
315: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
316: INFO = -1
317: ELSE IF( NL.LT.1 ) THEN
318: INFO = -2
319: ELSE IF( NR.LT.1 ) THEN
320: INFO = -3
321: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
322: INFO = -4
323: END IF
324: *
325: N = NL + NR + 1
326: *
327: IF( NRHS.LT.1 ) THEN
328: INFO = -5
329: ELSE IF( LDB.LT.N ) THEN
330: INFO = -7
331: ELSE IF( LDBX.LT.N ) THEN
332: INFO = -9
333: ELSE IF( GIVPTR.LT.0 ) THEN
334: INFO = -11
335: ELSE IF( LDGCOL.LT.N ) THEN
336: INFO = -13
337: ELSE IF( LDGNUM.LT.N ) THEN
338: INFO = -15
339: ELSE IF( K.LT.1 ) THEN
340: INFO = -20
341: END IF
342: IF( INFO.NE.0 ) THEN
343: CALL XERBLA( 'DLALS0', -INFO )
344: RETURN
345: END IF
346: *
347: M = N + SQRE
348: NLP1 = NL + 1
349: *
350: IF( ICOMPQ.EQ.0 ) THEN
351: *
352: * Apply back orthogonal transformations from the left.
353: *
354: * Step (1L): apply back the Givens rotations performed.
355: *
356: DO 10 I = 1, GIVPTR
357: CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
358: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
359: $ GIVNUM( I, 1 ) )
360: 10 CONTINUE
361: *
362: * Step (2L): permute rows of B.
363: *
364: CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
365: DO 20 I = 2, N
366: CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
367: 20 CONTINUE
368: *
369: * Step (3L): apply the inverse of the left singular vector
370: * matrix to BX.
371: *
372: IF( K.EQ.1 ) THEN
373: CALL DCOPY( NRHS, BX, LDBX, B, LDB )
374: IF( Z( 1 ).LT.ZERO ) THEN
375: CALL DSCAL( NRHS, NEGONE, B, LDB )
376: END IF
377: ELSE
378: DO 50 J = 1, K
379: DIFLJ = DIFL( J )
380: DJ = POLES( J, 1 )
381: DSIGJ = -POLES( J, 2 )
382: IF( J.LT.K ) THEN
383: DIFRJ = -DIFR( J, 1 )
384: DSIGJP = -POLES( J+1, 2 )
385: END IF
386: IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
387: $ THEN
388: WORK( J ) = ZERO
389: ELSE
390: WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
391: $ ( POLES( J, 2 )+DJ )
392: END IF
393: DO 30 I = 1, J - 1
394: IF( ( Z( I ).EQ.ZERO ) .OR.
395: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
396: WORK( I ) = ZERO
397: ELSE
398: WORK( I ) = POLES( I, 2 )*Z( I ) /
399: $ ( DLAMC3( POLES( I, 2 ), DSIGJ )-
400: $ DIFLJ ) / ( POLES( I, 2 )+DJ )
401: END IF
402: 30 CONTINUE
403: DO 40 I = J + 1, K
404: IF( ( Z( I ).EQ.ZERO ) .OR.
405: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
406: WORK( I ) = ZERO
407: ELSE
408: WORK( I ) = POLES( I, 2 )*Z( I ) /
409: $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+
410: $ DIFRJ ) / ( POLES( I, 2 )+DJ )
411: END IF
412: 40 CONTINUE
413: WORK( 1 ) = NEGONE
414: TEMP = DNRM2( K, WORK, 1 )
415: CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
416: $ B( J, 1 ), LDB )
417: CALL DLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
418: $ LDB, INFO )
419: 50 CONTINUE
420: END IF
421: *
422: * Move the deflated rows of BX to B also.
423: *
424: IF( K.LT.MAX( M, N ) )
425: $ CALL DLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
426: $ B( K+1, 1 ), LDB )
427: ELSE
428: *
429: * Apply back the right orthogonal transformations.
430: *
431: * Step (1R): apply back the new right singular vector matrix
432: * to B.
433: *
434: IF( K.EQ.1 ) THEN
435: CALL DCOPY( NRHS, B, LDB, BX, LDBX )
436: ELSE
437: DO 80 J = 1, K
438: DSIGJ = POLES( J, 2 )
439: IF( Z( J ).EQ.ZERO ) THEN
440: WORK( J ) = ZERO
441: ELSE
442: WORK( J ) = -Z( J ) / DIFL( J ) /
443: $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
444: END IF
445: DO 60 I = 1, J - 1
446: IF( Z( J ).EQ.ZERO ) THEN
447: WORK( I ) = ZERO
448: ELSE
449: WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
450: $ 2 ) )-DIFR( I, 1 ) ) /
451: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
452: END IF
453: 60 CONTINUE
454: DO 70 I = J + 1, K
455: IF( Z( J ).EQ.ZERO ) THEN
456: WORK( I ) = ZERO
457: ELSE
458: WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
459: $ 2 ) )-DIFL( I ) ) /
460: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
461: END IF
462: 70 CONTINUE
463: CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
464: $ BX( J, 1 ), LDBX )
465: 80 CONTINUE
466: END IF
467: *
468: * Step (2R): if SQRE = 1, apply back the rotation that is
469: * related to the right null space of the subproblem.
470: *
471: IF( SQRE.EQ.1 ) THEN
472: CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
473: CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
474: END IF
475: IF( K.LT.MAX( M, N ) )
476: $ CALL DLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
477: $ LDBX )
478: *
479: * Step (3R): permute rows of B.
480: *
481: CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
482: IF( SQRE.EQ.1 ) THEN
483: CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
484: END IF
485: DO 90 I = 2, N
486: CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
487: 90 CONTINUE
488: *
489: * Step (4R): apply back the Givens rotations performed.
490: *
491: DO 100 I = GIVPTR, 1, -1
492: CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
493: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
494: $ -GIVNUM( I, 1 ) )
495: 100 CONTINUE
496: END IF
497: *
498: RETURN
499: *
500: * End of DLALS0
501: *
502: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlals0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack