LAPACK 3.3.0

dlaed3.f

Go to the documentation of this file.
00001       SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
00002      $                   CTOT, W, S, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, K, LDQ, N, N1
00011       DOUBLE PRECISION   RHO
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            CTOT( * ), INDX( * )
00015       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
00016      $                   S( * ), W( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DLAED3 finds the roots of the secular equation, as defined by the
00023 *  values in D, W, and RHO, between 1 and K.  It makes the
00024 *  appropriate calls to DLAED4 and then updates the eigenvectors by
00025 *  multiplying the matrix of eigenvectors of the pair of eigensystems
00026 *  being combined by the matrix of eigenvectors of the K-by-K system
00027 *  which is solved here.
00028 *
00029 *  This code makes very mild assumptions about floating point
00030 *  arithmetic. It will work on machines with a guard digit in
00031 *  add/subtract, or on those binary machines without guard digits
00032 *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
00033 *  It could conceivably fail on hexadecimal or decimal machines
00034 *  without guard digits, but we know of none.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  K       (input) INTEGER
00040 *          The number of terms in the rational function to be solved by
00041 *          DLAED4.  K >= 0.
00042 *
00043 *  N       (input) INTEGER
00044 *          The number of rows and columns in the Q matrix.
00045 *          N >= K (deflation may result in N>K).
00046 *
00047 *  N1      (input) INTEGER
00048 *          The location of the last eigenvalue in the leading submatrix.
00049 *          min(1,N) <= N1 <= N/2.
00050 *
00051 *  D       (output) DOUBLE PRECISION array, dimension (N)
00052 *          D(I) contains the updated eigenvalues for
00053 *          1 <= I <= K.
00054 *
00055 *  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)
00056 *          Initially the first K columns are used as workspace.
00057 *          On output the columns 1 to K contain
00058 *          the updated eigenvectors.
00059 *
00060 *  LDQ     (input) INTEGER
00061 *          The leading dimension of the array Q.  LDQ >= max(1,N).
00062 *
00063 *  RHO     (input) DOUBLE PRECISION
00064 *          The value of the parameter in the rank one update equation.
00065 *          RHO >= 0 required.
00066 *
00067 *  DLAMDA  (input/output) DOUBLE PRECISION array, dimension (K)
00068 *          The first K elements of this array contain the old roots
00069 *          of the deflated updating problem.  These are the poles
00070 *          of the secular equation. May be changed on output by
00071 *          having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
00072 *          Cray-2, or Cray C-90, as described above.
00073 *
00074 *  Q2      (input) DOUBLE PRECISION array, dimension (LDQ2, N)
00075 *          The first K columns of this matrix contain the non-deflated
00076 *          eigenvectors for the split problem.
00077 *
00078 *  INDX    (input) INTEGER array, dimension (N)
00079 *          The permutation used to arrange the columns of the deflated
00080 *          Q matrix into three groups (see DLAED2).
00081 *          The rows of the eigenvectors found by DLAED4 must be likewise
00082 *          permuted before the matrix multiply can take place.
00083 *
00084 *  CTOT    (input) INTEGER array, dimension (4)
00085 *          A count of the total number of the various types of columns
00086 *          in Q, as described in INDX.  The fourth column type is any
00087 *          column which has been deflated.
00088 *
00089 *  W       (input/output) DOUBLE PRECISION array, dimension (K)
00090 *          The first K elements of this array contain the components
00091 *          of the deflation-adjusted updating vector. Destroyed on
00092 *          output.
00093 *
00094 *  S       (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K
00095 *          Will contain the eigenvectors of the repaired matrix which
00096 *          will be multiplied by the previously accumulated eigenvectors
00097 *          to update the system.
00098 *
00099 *  LDS     (input) INTEGER
00100 *          The leading dimension of S.  LDS >= max(1,K).
00101 *
00102 *  INFO    (output) INTEGER
00103 *          = 0:  successful exit.
00104 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00105 *          > 0:  if INFO = 1, an eigenvalue did not converge
00106 *
00107 *  Further Details
00108 *  ===============
00109 *
00110 *  Based on contributions by
00111 *     Jeff Rutter, Computer Science Division, University of California
00112 *     at Berkeley, USA
00113 *  Modified by Francoise Tisseur, University of Tennessee.
00114 *
00115 *  =====================================================================
00116 *
00117 *     .. Parameters ..
00118       DOUBLE PRECISION   ONE, ZERO
00119       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
00120 *     ..
00121 *     .. Local Scalars ..
00122       INTEGER            I, II, IQ2, J, N12, N2, N23
00123       DOUBLE PRECISION   TEMP
00124 *     ..
00125 *     .. External Functions ..
00126       DOUBLE PRECISION   DLAMC3, DNRM2
00127       EXTERNAL           DLAMC3, DNRM2
00128 *     ..
00129 *     .. External Subroutines ..
00130       EXTERNAL           DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA
00131 *     ..
00132 *     .. Intrinsic Functions ..
00133       INTRINSIC          MAX, SIGN, SQRT
00134 *     ..
00135 *     .. Executable Statements ..
00136 *
00137 *     Test the input parameters.
00138 *
00139       INFO = 0
00140 *
00141       IF( K.LT.0 ) THEN
00142          INFO = -1
00143       ELSE IF( N.LT.K ) THEN
00144          INFO = -2
00145       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00146          INFO = -6
00147       END IF
00148       IF( INFO.NE.0 ) THEN
00149          CALL XERBLA( 'DLAED3', -INFO )
00150          RETURN
00151       END IF
00152 *
00153 *     Quick return if possible
00154 *
00155       IF( K.EQ.0 )
00156      $   RETURN
00157 *
00158 *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
00159 *     be computed with high relative accuracy (barring over/underflow).
00160 *     This is a problem on machines without a guard digit in
00161 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
00162 *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
00163 *     which on any of these machines zeros out the bottommost
00164 *     bit of DLAMDA(I) if it is 1; this makes the subsequent
00165 *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
00166 *     occurs. On binary machines with a guard digit (almost all
00167 *     machines) it does not change DLAMDA(I) at all. On hexadecimal
00168 *     and decimal machines with a guard digit, it slightly
00169 *     changes the bottommost bits of DLAMDA(I). It does not account
00170 *     for hexadecimal or decimal machines without guard digits
00171 *     (we know of none). We use a subroutine call to compute
00172 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
00173 *     this code.
00174 *
00175       DO 10 I = 1, K
00176          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
00177    10 CONTINUE
00178 *
00179       DO 20 J = 1, K
00180          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
00181 *
00182 *        If the zero finder fails, the computation is terminated.
00183 *
00184          IF( INFO.NE.0 )
00185      $      GO TO 120
00186    20 CONTINUE
00187 *
00188       IF( K.EQ.1 )
00189      $   GO TO 110
00190       IF( K.EQ.2 ) THEN
00191          DO 30 J = 1, K
00192             W( 1 ) = Q( 1, J )
00193             W( 2 ) = Q( 2, J )
00194             II = INDX( 1 )
00195             Q( 1, J ) = W( II )
00196             II = INDX( 2 )
00197             Q( 2, J ) = W( II )
00198    30    CONTINUE
00199          GO TO 110
00200       END IF
00201 *
00202 *     Compute updated W.
00203 *
00204       CALL DCOPY( K, W, 1, S, 1 )
00205 *
00206 *     Initialize W(I) = Q(I,I)
00207 *
00208       CALL DCOPY( K, Q, LDQ+1, W, 1 )
00209       DO 60 J = 1, K
00210          DO 40 I = 1, J - 1
00211             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00212    40    CONTINUE
00213          DO 50 I = J + 1, K
00214             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00215    50    CONTINUE
00216    60 CONTINUE
00217       DO 70 I = 1, K
00218          W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
00219    70 CONTINUE
00220 *
00221 *     Compute eigenvectors of the modified rank-1 modification.
00222 *
00223       DO 100 J = 1, K
00224          DO 80 I = 1, K
00225             S( I ) = W( I ) / Q( I, J )
00226    80    CONTINUE
00227          TEMP = DNRM2( K, S, 1 )
00228          DO 90 I = 1, K
00229             II = INDX( I )
00230             Q( I, J ) = S( II ) / TEMP
00231    90    CONTINUE
00232   100 CONTINUE
00233 *
00234 *     Compute the updated eigenvectors.
00235 *
00236   110 CONTINUE
00237 *
00238       N2 = N - N1
00239       N12 = CTOT( 1 ) + CTOT( 2 )
00240       N23 = CTOT( 2 ) + CTOT( 3 )
00241 *
00242       CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
00243       IQ2 = N1*N12 + 1
00244       IF( N23.NE.0 ) THEN
00245          CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
00246      $               ZERO, Q( N1+1, 1 ), LDQ )
00247       ELSE
00248          CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
00249       END IF
00250 *
00251       CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 )
00252       IF( N12.NE.0 ) THEN
00253          CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
00254      $               LDQ )
00255       ELSE
00256          CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
00257       END IF
00258 *
00259 *
00260   120 CONTINUE
00261       RETURN
00262 *
00263 *     End of DLAED3
00264 *
00265       END
 All Files Functions