00001 SUBROUTINE SLAED1( N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, 00002 $ 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 CUTPNT, INFO, LDQ, N 00011 REAL RHO 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER INDXQ( * ), IWORK( * ) 00015 REAL D( * ), Q( LDQ, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * SLAED1 computes the updated eigensystem of a diagonal 00022 * matrix after modification by a rank-one symmetric matrix. This 00023 * routine is used only for the eigenproblem which requires all 00024 * eigenvalues and eigenvectors of a tridiagonal matrix. SLAED7 handles 00025 * the case in which eigenvalues only or eigenvalues and eigenvectors 00026 * of a full symmetric matrix (which was reduced to tridiagonal form) 00027 * are desired. 00028 * 00029 * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) 00030 * 00031 * where Z = Q'u, u is a vector of length N with ones in the 00032 * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. 00033 * 00034 * The eigenvectors of the original matrix are stored in Q, and the 00035 * eigenvalues are in D. The algorithm consists of three stages: 00036 * 00037 * The first stage consists of deflating the size of the problem 00038 * when there are multiple eigenvalues or if there is a zero in 00039 * the Z vector. For each such occurence the dimension of the 00040 * secular equation problem is reduced by one. This stage is 00041 * performed by the routine SLAED2. 00042 * 00043 * The second stage consists of calculating the updated 00044 * eigenvalues. This is done by finding the roots of the secular 00045 * equation via the routine SLAED4 (as called by SLAED3). 00046 * This routine also calculates the eigenvectors of the current 00047 * problem. 00048 * 00049 * The final stage consists of computing the updated eigenvectors 00050 * directly using the updated eigenvalues. The eigenvectors for 00051 * the current problem are multiplied with the eigenvectors from 00052 * the overall problem. 00053 * 00054 * Arguments 00055 * ========= 00056 * 00057 * N (input) INTEGER 00058 * The dimension of the symmetric tridiagonal matrix. N >= 0. 00059 * 00060 * D (input/output) REAL array, dimension (N) 00061 * On entry, the eigenvalues of the rank-1-perturbed matrix. 00062 * On exit, the eigenvalues of the repaired matrix. 00063 * 00064 * Q (input/output) REAL array, dimension (LDQ,N) 00065 * On entry, the eigenvectors of the rank-1-perturbed matrix. 00066 * On exit, the eigenvectors of the repaired tridiagonal matrix. 00067 * 00068 * LDQ (input) INTEGER 00069 * The leading dimension of the array Q. LDQ >= max(1,N). 00070 * 00071 * INDXQ (input/output) INTEGER array, dimension (N) 00072 * On entry, the permutation which separately sorts the two 00073 * subproblems in D into ascending order. 00074 * On exit, the permutation which will reintegrate the 00075 * subproblems back into sorted order, 00076 * i.e. D( INDXQ( I = 1, N ) ) will be in ascending order. 00077 * 00078 * RHO (input) REAL 00079 * The subdiagonal entry used to create the rank-1 modification. 00080 * 00081 * CUTPNT (input) INTEGER 00082 * The location of the last eigenvalue in the leading sub-matrix. 00083 * min(1,N) <= CUTPNT <= N/2. 00084 * 00085 * WORK (workspace) REAL array, dimension (4*N + N**2) 00086 * 00087 * IWORK (workspace) INTEGER array, dimension (4*N) 00088 * 00089 * INFO (output) INTEGER 00090 * = 0: successful exit. 00091 * < 0: if INFO = -i, the i-th argument had an illegal value. 00092 * > 0: if INFO = 1, an eigenvalue did not converge 00093 * 00094 * Further Details 00095 * =============== 00096 * 00097 * Based on contributions by 00098 * Jeff Rutter, Computer Science Division, University of California 00099 * at Berkeley, USA 00100 * Modified by Francoise Tisseur, University of Tennessee. 00101 * 00102 * ===================================================================== 00103 * 00104 * .. Local Scalars .. 00105 INTEGER COLTYP, CPP1, I, IDLMDA, INDX, INDXC, INDXP, 00106 $ IQ2, IS, IW, IZ, K, N1, N2 00107 * .. 00108 * .. External Subroutines .. 00109 EXTERNAL SCOPY, SLAED2, SLAED3, SLAMRG, XERBLA 00110 * .. 00111 * .. Intrinsic Functions .. 00112 INTRINSIC MAX, MIN 00113 * .. 00114 * .. Executable Statements .. 00115 * 00116 * Test the input parameters. 00117 * 00118 INFO = 0 00119 * 00120 IF( N.LT.0 ) THEN 00121 INFO = -1 00122 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 00123 INFO = -4 00124 ELSE IF( MIN( 1, N / 2 ).GT.CUTPNT .OR. ( N / 2 ).LT.CUTPNT ) THEN 00125 INFO = -7 00126 END IF 00127 IF( INFO.NE.0 ) THEN 00128 CALL XERBLA( 'SLAED1', -INFO ) 00129 RETURN 00130 END IF 00131 * 00132 * Quick return if possible 00133 * 00134 IF( N.EQ.0 ) 00135 $ RETURN 00136 * 00137 * The following values are integer pointers which indicate 00138 * the portion of the workspace 00139 * used by a particular array in SLAED2 and SLAED3. 00140 * 00141 IZ = 1 00142 IDLMDA = IZ + N 00143 IW = IDLMDA + N 00144 IQ2 = IW + N 00145 * 00146 INDX = 1 00147 INDXC = INDX + N 00148 COLTYP = INDXC + N 00149 INDXP = COLTYP + N 00150 * 00151 * 00152 * Form the z-vector which consists of the last row of Q_1 and the 00153 * first row of Q_2. 00154 * 00155 CALL SCOPY( CUTPNT, Q( CUTPNT, 1 ), LDQ, WORK( IZ ), 1 ) 00156 CPP1 = CUTPNT + 1 00157 CALL SCOPY( N-CUTPNT, Q( CPP1, CPP1 ), LDQ, WORK( IZ+CUTPNT ), 1 ) 00158 * 00159 * Deflate eigenvalues. 00160 * 00161 CALL SLAED2( K, N, CUTPNT, D, Q, LDQ, INDXQ, RHO, WORK( IZ ), 00162 $ WORK( IDLMDA ), WORK( IW ), WORK( IQ2 ), 00163 $ IWORK( INDX ), IWORK( INDXC ), IWORK( INDXP ), 00164 $ IWORK( COLTYP ), INFO ) 00165 * 00166 IF( INFO.NE.0 ) 00167 $ GO TO 20 00168 * 00169 * Solve Secular Equation. 00170 * 00171 IF( K.NE.0 ) THEN 00172 IS = ( IWORK( COLTYP )+IWORK( COLTYP+1 ) )*CUTPNT + 00173 $ ( IWORK( COLTYP+1 )+IWORK( COLTYP+2 ) )*( N-CUTPNT ) + IQ2 00174 CALL SLAED3( K, N, CUTPNT, D, Q, LDQ, RHO, WORK( IDLMDA ), 00175 $ WORK( IQ2 ), IWORK( INDXC ), IWORK( COLTYP ), 00176 $ WORK( IW ), WORK( IS ), INFO ) 00177 IF( INFO.NE.0 ) 00178 $ GO TO 20 00179 * 00180 * Prepare the INDXQ sorting permutation. 00181 * 00182 N1 = K 00183 N2 = N - K 00184 CALL SLAMRG( N1, N2, D, 1, -1, INDXQ ) 00185 ELSE 00186 DO 10 I = 1, N 00187 INDXQ( I ) = I 00188 10 CONTINUE 00189 END IF 00190 * 00191 20 CONTINUE 00192 RETURN 00193 * 00194 * End of SLAED1 00195 * 00196 END
1.7.2