LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slaed3.f
Go to the documentation of this file.
00001 *> \brief \b SLAED3
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLAED3 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed3.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed3.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed3.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
00022 *                          CTOT, W, S, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, K, LDQ, N, N1
00026 *       REAL               RHO
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            CTOT( * ), INDX( * )
00030 *       REAL               D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
00031 *      $                   S( * ), W( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> SLAED3 finds the roots of the secular equation, as defined by the
00041 *> values in D, W, and RHO, between 1 and K.  It makes the
00042 *> appropriate calls to SLAED4 and then updates the eigenvectors by
00043 *> multiplying the matrix of eigenvectors of the pair of eigensystems
00044 *> being combined by the matrix of eigenvectors of the K-by-K system
00045 *> which is solved here.
00046 *>
00047 *> This code makes very mild assumptions about floating point
00048 *> arithmetic. It will work on machines with a guard digit in
00049 *> add/subtract, or on those binary machines without guard digits
00050 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
00051 *> It could conceivably fail on hexadecimal or decimal machines
00052 *> without guard digits, but we know of none.
00053 *> \endverbatim
00054 *
00055 *  Arguments:
00056 *  ==========
00057 *
00058 *> \param[in] K
00059 *> \verbatim
00060 *>          K is INTEGER
00061 *>          The number of terms in the rational function to be solved by
00062 *>          SLAED4.  K >= 0.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] N
00066 *> \verbatim
00067 *>          N is INTEGER
00068 *>          The number of rows and columns in the Q matrix.
00069 *>          N >= K (deflation may result in N>K).
00070 *> \endverbatim
00071 *>
00072 *> \param[in] N1
00073 *> \verbatim
00074 *>          N1 is INTEGER
00075 *>          The location of the last eigenvalue in the leading submatrix.
00076 *>          min(1,N) <= N1 <= N/2.
00077 *> \endverbatim
00078 *>
00079 *> \param[out] D
00080 *> \verbatim
00081 *>          D is REAL array, dimension (N)
00082 *>          D(I) contains the updated eigenvalues for
00083 *>          1 <= I <= K.
00084 *> \endverbatim
00085 *>
00086 *> \param[out] Q
00087 *> \verbatim
00088 *>          Q is REAL array, dimension (LDQ,N)
00089 *>          Initially the first K columns are used as workspace.
00090 *>          On output the columns 1 to K contain
00091 *>          the updated eigenvectors.
00092 *> \endverbatim
00093 *>
00094 *> \param[in] LDQ
00095 *> \verbatim
00096 *>          LDQ is INTEGER
00097 *>          The leading dimension of the array Q.  LDQ >= max(1,N).
00098 *> \endverbatim
00099 *>
00100 *> \param[in] RHO
00101 *> \verbatim
00102 *>          RHO is REAL
00103 *>          The value of the parameter in the rank one update equation.
00104 *>          RHO >= 0 required.
00105 *> \endverbatim
00106 *>
00107 *> \param[in,out] DLAMDA
00108 *> \verbatim
00109 *>          DLAMDA is REAL array, dimension (K)
00110 *>          The first K elements of this array contain the old roots
00111 *>          of the deflated updating problem.  These are the poles
00112 *>          of the secular equation. May be changed on output by
00113 *>          having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
00114 *>          Cray-2, or Cray C-90, as described above.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] Q2
00118 *> \verbatim
00119 *>          Q2 is REAL array, dimension (LDQ2, N)
00120 *>          The first K columns of this matrix contain the non-deflated
00121 *>          eigenvectors for the split problem.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] INDX
00125 *> \verbatim
00126 *>          INDX is INTEGER array, dimension (N)
00127 *>          The permutation used to arrange the columns of the deflated
00128 *>          Q matrix into three groups (see SLAED2).
00129 *>          The rows of the eigenvectors found by SLAED4 must be likewise
00130 *>          permuted before the matrix multiply can take place.
00131 *> \endverbatim
00132 *>
00133 *> \param[in] CTOT
00134 *> \verbatim
00135 *>          CTOT is INTEGER array, dimension (4)
00136 *>          A count of the total number of the various types of columns
00137 *>          in Q, as described in INDX.  The fourth column type is any
00138 *>          column which has been deflated.
00139 *> \endverbatim
00140 *>
00141 *> \param[in,out] W
00142 *> \verbatim
00143 *>          W is REAL array, dimension (K)
00144 *>          The first K elements of this array contain the components
00145 *>          of the deflation-adjusted updating vector. Destroyed on
00146 *>          output.
00147 *> \endverbatim
00148 *>
00149 *> \param[out] S
00150 *> \verbatim
00151 *>          S is REAL array, dimension (N1 + 1)*K
00152 *>          Will contain the eigenvectors of the repaired matrix which
00153 *>          will be multiplied by the previously accumulated eigenvectors
00154 *>          to update the system.
00155 *> \endverbatim
00156 *>
00157 *> \param[out] INFO
00158 *> \verbatim
00159 *>          INFO is INTEGER
00160 *>          = 0:  successful exit.
00161 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00162 *>          > 0:  if INFO = 1, an eigenvalue did not converge
00163 *> \endverbatim
00164 *
00165 *  Authors:
00166 *  ========
00167 *
00168 *> \author Univ. of Tennessee 
00169 *> \author Univ. of California Berkeley 
00170 *> \author Univ. of Colorado Denver 
00171 *> \author NAG Ltd. 
00172 *
00173 *> \date November 2011
00174 *
00175 *> \ingroup auxOTHERcomputational
00176 *
00177 *> \par Contributors:
00178 *  ==================
00179 *>
00180 *> Jeff Rutter, Computer Science Division, University of California
00181 *> at Berkeley, USA \n
00182 *>  Modified by Francoise Tisseur, University of Tennessee
00183 *>
00184 *  =====================================================================
00185       SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
00186      $                   CTOT, W, S, INFO )
00187 *
00188 *  -- LAPACK computational routine (version 3.4.0) --
00189 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00190 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00191 *     November 2011
00192 *
00193 *     .. Scalar Arguments ..
00194       INTEGER            INFO, K, LDQ, N, N1
00195       REAL               RHO
00196 *     ..
00197 *     .. Array Arguments ..
00198       INTEGER            CTOT( * ), INDX( * )
00199       REAL               D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
00200      $                   S( * ), W( * )
00201 *     ..
00202 *
00203 *  =====================================================================
00204 *
00205 *     .. Parameters ..
00206       REAL               ONE, ZERO
00207       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
00208 *     ..
00209 *     .. Local Scalars ..
00210       INTEGER            I, II, IQ2, J, N12, N2, N23
00211       REAL               TEMP
00212 *     ..
00213 *     .. External Functions ..
00214       REAL               SLAMC3, SNRM2
00215       EXTERNAL           SLAMC3, SNRM2
00216 *     ..
00217 *     .. External Subroutines ..
00218       EXTERNAL           SCOPY, SGEMM, SLACPY, SLAED4, SLASET, XERBLA
00219 *     ..
00220 *     .. Intrinsic Functions ..
00221       INTRINSIC          MAX, SIGN, SQRT
00222 *     ..
00223 *     .. Executable Statements ..
00224 *
00225 *     Test the input parameters.
00226 *
00227       INFO = 0
00228 *
00229       IF( K.LT.0 ) THEN
00230          INFO = -1
00231       ELSE IF( N.LT.K ) THEN
00232          INFO = -2
00233       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00234          INFO = -6
00235       END IF
00236       IF( INFO.NE.0 ) THEN
00237          CALL XERBLA( 'SLAED3', -INFO )
00238          RETURN
00239       END IF
00240 *
00241 *     Quick return if possible
00242 *
00243       IF( K.EQ.0 )
00244      $   RETURN
00245 *
00246 *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
00247 *     be computed with high relative accuracy (barring over/underflow).
00248 *     This is a problem on machines without a guard digit in
00249 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
00250 *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
00251 *     which on any of these machines zeros out the bottommost
00252 *     bit of DLAMDA(I) if it is 1; this makes the subsequent
00253 *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
00254 *     occurs. On binary machines with a guard digit (almost all
00255 *     machines) it does not change DLAMDA(I) at all. On hexadecimal
00256 *     and decimal machines with a guard digit, it slightly
00257 *     changes the bottommost bits of DLAMDA(I). It does not account
00258 *     for hexadecimal or decimal machines without guard digits
00259 *     (we know of none). We use a subroutine call to compute
00260 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
00261 *     this code.
00262 *
00263       DO 10 I = 1, K
00264          DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
00265    10 CONTINUE
00266 *
00267       DO 20 J = 1, K
00268          CALL SLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
00269 *
00270 *        If the zero finder fails, the computation is terminated.
00271 *
00272          IF( INFO.NE.0 )
00273      $      GO TO 120
00274    20 CONTINUE
00275 *
00276       IF( K.EQ.1 )
00277      $   GO TO 110
00278       IF( K.EQ.2 ) THEN
00279          DO 30 J = 1, K
00280             W( 1 ) = Q( 1, J )
00281             W( 2 ) = Q( 2, J )
00282             II = INDX( 1 )
00283             Q( 1, J ) = W( II )
00284             II = INDX( 2 )
00285             Q( 2, J ) = W( II )
00286    30    CONTINUE
00287          GO TO 110
00288       END IF
00289 *
00290 *     Compute updated W.
00291 *
00292       CALL SCOPY( K, W, 1, S, 1 )
00293 *
00294 *     Initialize W(I) = Q(I,I)
00295 *
00296       CALL SCOPY( K, Q, LDQ+1, W, 1 )
00297       DO 60 J = 1, K
00298          DO 40 I = 1, J - 1
00299             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00300    40    CONTINUE
00301          DO 50 I = J + 1, K
00302             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00303    50    CONTINUE
00304    60 CONTINUE
00305       DO 70 I = 1, K
00306          W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
00307    70 CONTINUE
00308 *
00309 *     Compute eigenvectors of the modified rank-1 modification.
00310 *
00311       DO 100 J = 1, K
00312          DO 80 I = 1, K
00313             S( I ) = W( I ) / Q( I, J )
00314    80    CONTINUE
00315          TEMP = SNRM2( K, S, 1 )
00316          DO 90 I = 1, K
00317             II = INDX( I )
00318             Q( I, J ) = S( II ) / TEMP
00319    90    CONTINUE
00320   100 CONTINUE
00321 *
00322 *     Compute the updated eigenvectors.
00323 *
00324   110 CONTINUE
00325 *
00326       N2 = N - N1
00327       N12 = CTOT( 1 ) + CTOT( 2 )
00328       N23 = CTOT( 2 ) + CTOT( 3 )
00329 *
00330       CALL SLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
00331       IQ2 = N1*N12 + 1
00332       IF( N23.NE.0 ) THEN
00333          CALL SGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
00334      $               ZERO, Q( N1+1, 1 ), LDQ )
00335       ELSE
00336          CALL SLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
00337       END IF
00338 *
00339       CALL SLACPY( 'A', N12, K, Q, LDQ, S, N12 )
00340       IF( N12.NE.0 ) THEN
00341          CALL SGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
00342      $               LDQ )
00343       ELSE
00344          CALL SLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
00345       END IF
00346 *
00347 *
00348   120 CONTINUE
00349       RETURN
00350 *
00351 *     End of SLAED3
00352 *
00353       END
 All Files Functions