00001 SUBROUTINE SLASQ1( N, D, E, WORK, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * 00005 * -- Contributed by Osni Marques of the Lawrence Berkeley National -- 00006 * -- Laboratory and Beresford Parlett of the Univ. of California at -- 00007 * -- Berkeley -- 00008 * -- November 2008 -- 00009 * 00010 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00011 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00012 * 00013 * .. Scalar Arguments .. 00014 INTEGER INFO, N 00015 * .. 00016 * .. Array Arguments .. 00017 REAL D( * ), E( * ), WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * SLASQ1 computes the singular values of a real N-by-N bidiagonal 00024 * matrix with diagonal D and off-diagonal E. The singular values 00025 * are computed to high relative accuracy, in the absence of 00026 * denormalization, underflow and overflow. The algorithm was first 00027 * presented in 00028 * 00029 * "Accurate singular values and differential qd algorithms" by K. V. 00030 * Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, 00031 * 1994, 00032 * 00033 * and the present implementation is described in "An implementation of 00034 * the dqds Algorithm (Positive Case)", LAPACK Working Note. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * N (input) INTEGER 00040 * The number of rows and columns in the matrix. N >= 0. 00041 * 00042 * D (input/output) REAL array, dimension (N) 00043 * On entry, D contains the diagonal elements of the 00044 * bidiagonal matrix whose SVD is desired. On normal exit, 00045 * D contains the singular values in decreasing order. 00046 * 00047 * E (input/output) REAL array, dimension (N) 00048 * On entry, elements E(1:N-1) contain the off-diagonal elements 00049 * of the bidiagonal matrix whose SVD is desired. 00050 * On exit, E is overwritten. 00051 * 00052 * WORK (workspace) REAL array, dimension (4*N) 00053 * 00054 * INFO (output) INTEGER 00055 * = 0: successful exit 00056 * < 0: if INFO = -i, the i-th argument had an illegal value 00057 * > 0: the algorithm failed 00058 * = 1, a split was marked by a positive value in E 00059 * = 2, current block of Z not diagonalized after 30*N 00060 * iterations (in inner while loop) 00061 * = 3, termination criterion of outer while loop not met 00062 * (program created more than N unreduced blocks) 00063 * 00064 * ===================================================================== 00065 * 00066 * .. Parameters .. 00067 REAL ZERO 00068 PARAMETER ( ZERO = 0.0E0 ) 00069 * .. 00070 * .. Local Scalars .. 00071 INTEGER I, IINFO 00072 REAL EPS, SCALE, SAFMIN, SIGMN, SIGMX 00073 * .. 00074 * .. External Subroutines .. 00075 EXTERNAL SCOPY, SLAS2, SLASCL, SLASQ2, SLASRT, XERBLA 00076 * .. 00077 * .. External Functions .. 00078 REAL SLAMCH 00079 EXTERNAL SLAMCH 00080 * .. 00081 * .. Intrinsic Functions .. 00082 INTRINSIC ABS, MAX, SQRT 00083 * .. 00084 * .. Executable Statements .. 00085 * 00086 INFO = 0 00087 IF( N.LT.0 ) THEN 00088 INFO = -2 00089 CALL XERBLA( 'SLASQ1', -INFO ) 00090 RETURN 00091 ELSE IF( N.EQ.0 ) THEN 00092 RETURN 00093 ELSE IF( N.EQ.1 ) THEN 00094 D( 1 ) = ABS( D( 1 ) ) 00095 RETURN 00096 ELSE IF( N.EQ.2 ) THEN 00097 CALL SLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX ) 00098 D( 1 ) = SIGMX 00099 D( 2 ) = SIGMN 00100 RETURN 00101 END IF 00102 * 00103 * Estimate the largest singular value. 00104 * 00105 SIGMX = ZERO 00106 DO 10 I = 1, N - 1 00107 D( I ) = ABS( D( I ) ) 00108 SIGMX = MAX( SIGMX, ABS( E( I ) ) ) 00109 10 CONTINUE 00110 D( N ) = ABS( D( N ) ) 00111 * 00112 * Early return if SIGMX is zero (matrix is already diagonal). 00113 * 00114 IF( SIGMX.EQ.ZERO ) THEN 00115 CALL SLASRT( 'D', N, D, IINFO ) 00116 RETURN 00117 END IF 00118 * 00119 DO 20 I = 1, N 00120 SIGMX = MAX( SIGMX, D( I ) ) 00121 20 CONTINUE 00122 * 00123 * Copy D and E into WORK (in the Z format) and scale (squaring the 00124 * input data makes scaling by a power of the radix pointless). 00125 * 00126 EPS = SLAMCH( 'Precision' ) 00127 SAFMIN = SLAMCH( 'Safe minimum' ) 00128 SCALE = SQRT( EPS / SAFMIN ) 00129 CALL SCOPY( N, D, 1, WORK( 1 ), 2 ) 00130 CALL SCOPY( N-1, E, 1, WORK( 2 ), 2 ) 00131 CALL SLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1, 00132 $ IINFO ) 00133 * 00134 * Compute the q's and e's. 00135 * 00136 DO 30 I = 1, 2*N - 1 00137 WORK( I ) = WORK( I )**2 00138 30 CONTINUE 00139 WORK( 2*N ) = ZERO 00140 * 00141 CALL SLASQ2( N, WORK, INFO ) 00142 * 00143 IF( INFO.EQ.0 ) THEN 00144 DO 40 I = 1, N 00145 D( I ) = SQRT( WORK( I ) ) 00146 40 CONTINUE 00147 CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO ) 00148 END IF 00149 * 00150 RETURN 00151 * 00152 * End of SLASQ1 00153 * 00154 END
1.7.2