SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER LDA, LDB, M, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DLACPY copies all or part of a two-dimensional matrix A to another * matrix B. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies the part of the matrix A to be copied to B. * = 'U': Upper triangular part * = 'L': Lower triangular part * Otherwise: All of the matrix A * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The m by n matrix A. If UPLO = 'U', only the upper triangle * or trapezoid is accessed; if UPLO = 'L', only the lower * triangle or trapezoid is accessed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (output) DOUBLE PRECISION array, dimension (LDB,N) * On exit, B = A in the locations specified by UPLO. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,M). * * ===================================================================== * * .. Local Scalars .. INTEGER I, J * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * IF( LSAME( UPLO, 'U' ) ) THEN DO 20 J = 1, N DO 10 I = 1, MIN( J, M ) B( I, J ) = A( I, J ) 10 CONTINUE 20 CONTINUE ELSE IF( LSAME( UPLO, 'L' ) ) THEN DO 40 J = 1, N DO 30 I = J, M B( I, J ) = A( I, J ) 30 CONTINUE 40 CONTINUE ELSE DO 60 J = 1, N DO 50 I = 1, M B( I, J ) = A( I, J ) 50 CONTINUE 60 CONTINUE END IF RETURN * * End of DLACPY * END SUBROUTINE DLAE2( A, B, C, RT1, RT2 ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. DOUBLE PRECISION A, B, C, RT1, RT2 * .. * * Purpose * ======= * * DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix * [ A B ] * [ B C ]. * On return, RT1 is the eigenvalue of larger absolute value, and RT2 * is the eigenvalue of smaller absolute value. * * Arguments * ========= * * A (input) DOUBLE PRECISION * The (1,1) element of the 2-by-2 matrix. * * B (input) DOUBLE PRECISION * The (1,2) and (2,1) elements of the 2-by-2 matrix. * * C (input) DOUBLE PRECISION * The (2,2) element of the 2-by-2 matrix. * * RT1 (output) DOUBLE PRECISION * The eigenvalue of larger absolute value. * * RT2 (output) DOUBLE PRECISION * The eigenvalue of smaller absolute value. * * Further Details * =============== * * RT1 is accurate to a few ulps barring over/underflow. * * RT2 may be inaccurate if there is massive cancellation in the * determinant A*C-B*B; higher precision or correctly rounded or * correctly truncated arithmetic would be needed to compute RT2 * accurately in all cases. * * Overflow is possible only if RT1 is within a factor of 5 of overflow. * Underflow is harmless if the input data is 0 or exceeds * underflow_threshold / macheps. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) DOUBLE PRECISION TWO PARAMETER ( TWO = 2.0D0 ) DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) DOUBLE PRECISION HALF PARAMETER ( HALF = 0.5D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION AB, ACMN, ACMX, ADF, DF, RT, SM, TB * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. * * Compute the eigenvalues * SM = A + C DF = A - C ADF = ABS( DF ) TB = B + B AB = ABS( TB ) IF( ABS( A ).GT.ABS( C ) ) THEN ACMX = A ACMN = C ELSE ACMX = C ACMN = A END IF IF( ADF.GT.AB ) THEN RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) ELSE * * Includes case AB=ADF=0 * RT = AB*SQRT( TWO ) END IF IF( SM.LT.ZERO ) THEN RT1 = HALF*( SM-RT ) * * Order of execution important. * To get fully accurate smaller eigenvalue, * next line needs to be executed in higher precision. * RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B ELSE IF( SM.GT.ZERO ) THEN RT1 = HALF*( SM+RT ) * * Order of execution important. * To get fully accurate smaller eigenvalue, * next line needs to be executed in higher precision. * RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B ELSE * * Includes case RT1 = RT2 = 0 * RT1 = HALF*RT RT2 = -HALF*RT END IF RETURN * * End of DLAE2 * END SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, $ WORK, IWORK, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ * .. * .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ), $ WORK( * ) * .. * * Purpose * ======= * * DLAED0 computes all eigenvalues and corresponding eigenvectors of a * symmetric tridiagonal matrix using the divide and conquer method. * * Arguments * ========= * * ICOMPQ (input) INTEGER * = 0: Compute eigenvalues only. * = 1: Compute eigenvectors of original dense symmetric matrix * also. On entry, Q contains the orthogonal matrix used * to reduce the original matrix to tridiagonal form. * = 2: Compute eigenvalues and eigenvectors of tridiagonal * matrix. * * QSIZ (input) INTEGER * The dimension of the orthogonal matrix used to reduce * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the main diagonal of the tridiagonal matrix. * On exit, its eigenvalues. * * E (input) DOUBLE PRECISION array, dimension (N-1) * The off-diagonal elements of the tridiagonal matrix. * On exit, E has been destroyed. * * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) * On entry, Q must contain an N-by-N orthogonal matrix. * If ICOMPQ = 0 Q is not referenced. * If ICOMPQ = 1 On entry, Q is a subset of the columns of the * orthogonal matrix used to reduce the full * matrix to tridiagonal form corresponding to * the subset of the full matrix which is being * decomposed at this time. * If ICOMPQ = 2 On entry, Q will be the identity matrix. * On exit, Q contains the eigenvectors of the * tridiagonal matrix. * * LDQ (input) INTEGER * The leading dimension of the array Q. If eigenvectors are * desired, then LDQ >= max(1,N). In any case, LDQ >= 1. * * QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N) * Referenced only when ICOMPQ = 1. Used to store parts of * the eigenvector matrix when the updating matrix multiplies * take place. * * LDQS (input) INTEGER * The leading dimension of the array QSTORE. If ICOMPQ = 1, * then LDQS >= max(1,N). In any case, LDQS >= 1. * * WORK (workspace) DOUBLE PRECISION array, * If ICOMPQ = 0 or 1, the dimension of WORK must be at least * 1 + 3*N + 2*N*lg N + 2*N**2 * ( lg( N ) = smallest integer k * such that 2^k >= N ) * If ICOMPQ = 2, the dimension of WORK must be at least * 4*N + N**2. * * IWORK (workspace) INTEGER array, * If ICOMPQ = 0 or 1, the dimension of IWORK must be at least * 6 + 6*N + 5*N*lg N. * ( lg( N ) = smallest integer k * such that 2^k >= N ) * If ICOMPQ = 2, the dimension of IWORK must be at least * 3 + 5*N. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: The algorithm failed to compute an eigenvalue while * working on the submatrix lying in rows and columns * INFO/(N+1) through mod(INFO,N+1). * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 ) * .. * .. Local Scalars .. INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM, $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM, $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1, $ SPM2, SUBMAT, SUBPBS, TLVLS DOUBLE PRECISION TEMP * .. * .. External Subroutines .. EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR, $ XERBLA * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, LOG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN INFO = -1 ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -7 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN INFO = -9 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED0', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 ) * * Determine the size and placement of the submatrices, and save in * the leading elements of IWORK. * IWORK( 1 ) = N SUBPBS = 1 TLVLS = 0 10 CONTINUE IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN DO 20 J = SUBPBS, 1, -1 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 IWORK( 2*J-1 ) = IWORK( J ) / 2 20 CONTINUE TLVLS = TLVLS + 1 SUBPBS = 2*SUBPBS GO TO 10 END IF DO 30 J = 2, SUBPBS IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 30 CONTINUE * * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 * using rank-1 modifications (cuts). * SPM1 = SUBPBS - 1 DO 40 I = 1, SPM1 SUBMAT = IWORK( I ) + 1 SMM1 = SUBMAT - 1 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) 40 CONTINUE * INDXQ = 4*N + 3 IF( ICOMPQ.NE.2 ) THEN * * Set up workspaces for eigenvalues only/accumulate new vectors * routine * TEMP = LOG( DBLE( N ) ) / LOG( TWO ) LGN = INT( TEMP ) IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IPRMPT = INDXQ + N + 1 IPERM = IPRMPT + N*LGN IQPTR = IPERM + N*LGN IGIVPT = IQPTR + N + 2 IGIVCL = IGIVPT + N*LGN * IGIVNM = 1 IQ = IGIVNM + 2*N*LGN IWREM = IQ + N**2 + 1 * * Initialize pointers * DO 50 I = 0, SUBPBS IWORK( IPRMPT+I ) = 1 IWORK( IGIVPT+I ) = 1 50 CONTINUE IWORK( IQPTR ) = 1 END IF * * Solve each submatrix eigenproblem at the bottom of the divide and * conquer tree. * CURR = 0 DO 70 I = 0, SPM1 IF( I.EQ.0 ) THEN SUBMAT = 1 MATSIZ = IWORK( 1 ) ELSE SUBMAT = IWORK( I ) + 1 MATSIZ = IWORK( I+1 ) - IWORK( I ) END IF IF( ICOMPQ.EQ.2 ) THEN CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO ) IF( INFO.NE.0 ) $ GO TO 130 ELSE CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK, $ INFO ) IF( INFO.NE.0 ) $ GO TO 130 IF( ICOMPQ.EQ.1 ) THEN CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE, $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+ $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ), $ LDQS ) END IF IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2 CURR = CURR + 1 END IF K = 1 DO 60 J = SUBMAT, IWORK( I+1 ) IWORK( INDXQ+J ) = K K = K + 1 60 CONTINUE 70 CONTINUE * * Successively merge eigensystems of adjacent submatrices * into eigensystem for the corresponding larger matrix. * * while ( SUBPBS > 1 ) * CURLVL = 1 80 CONTINUE IF( SUBPBS.GT.1 ) THEN SPM2 = SUBPBS - 2 DO 90 I = 0, SPM2, 2 IF( I.EQ.0 ) THEN SUBMAT = 1 MATSIZ = IWORK( 2 ) MSD2 = IWORK( 1 ) CURPRB = 0 ELSE SUBMAT = IWORK( I ) + 1 MATSIZ = IWORK( I+2 ) - IWORK( I ) MSD2 = MATSIZ / 2 CURPRB = CURPRB + 1 END IF * * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) * into an eigensystem of size MATSIZ. * DLAED1 is used only for the full eigensystem of a tridiagonal * matrix. * DLAED7 handles the cases in which eigenvalues only or eigenvalues * and eigenvectors of a full symmetric matrix (which was reduced to * tridiagonal form) are desired. * IF( ICOMPQ.EQ.2 ) THEN CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ), $ LDQ, IWORK( INDXQ+SUBMAT ), $ E( SUBMAT+MSD2-1 ), MSD2, WORK, $ IWORK( SUBPBS+1 ), INFO ) ELSE CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB, $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ), $ MSD2, WORK( IQ ), IWORK( IQPTR ), $ IWORK( IPRMPT ), IWORK( IPERM ), $ IWORK( IGIVPT ), IWORK( IGIVCL ), $ WORK( IGIVNM ), WORK( IWREM ), $ IWORK( SUBPBS+1 ), INFO ) END IF IF( INFO.NE.0 ) $ GO TO 130 IWORK( I / 2+1 ) = IWORK( I+2 ) 90 CONTINUE SUBPBS = SUBPBS / 2 CURLVL = CURLVL + 1 GO TO 80 END IF * * end while * * Re-merge the eigenvalues/vectors which were deflated at the final * merge step. * IF( ICOMPQ.EQ.1 ) THEN DO 100 I = 1, N J = IWORK( INDXQ+I ) WORK( I ) = D( J ) CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 ) 100 CONTINUE CALL DCOPY( N, WORK, 1, D, 1 ) ELSE IF( ICOMPQ.EQ.2 ) THEN DO 110 I = 1, N J = IWORK( INDXQ+I ) WORK( I ) = D( J ) CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 ) 110 CONTINUE CALL DCOPY( N, WORK, 1, D, 1 ) CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ ) ELSE DO 120 I = 1, N J = IWORK( INDXQ+I ) WORK( I ) = D( J ) 120 CONTINUE CALL DCOPY( N, WORK, 1, D, 1 ) END IF GO TO 140 * 130 CONTINUE INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 * 140 CONTINUE RETURN * * End of DLAED0 * END SUBROUTINE DLAED1( N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, $ INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER CUTPNT, INFO, LDQ, N DOUBLE PRECISION RHO * .. * .. Array Arguments .. INTEGER INDXQ( * ), IWORK( * ) DOUBLE PRECISION D( * ), Q( LDQ, * ), WORK( * ) * .. * * Purpose * ======= * * DLAED1 computes the updated eigensystem of a diagonal * matrix after modification by a rank-one symmetric matrix. This * routine is used only for the eigenproblem which requires all * eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles * the case in which eigenvalues only or eigenvalues and eigenvectors * of a full symmetric matrix (which was reduced to tridiagonal form) * are desired. * * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) * * where Z = Q'u, u is a vector of length N with ones in the * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. * * The eigenvectors of the original matrix are stored in Q, and the * eigenvalues are in D. The algorithm consists of three stages: * * The first stage consists of deflating the size of the problem * when there are multiple eigenvalues or if there is a zero in * the Z vector. For each such occurence the dimension of the * secular equation problem is reduced by one. This stage is * performed by the routine DLAED2. * * The second stage consists of calculating the updated * eigenvalues. This is done by finding the roots of the secular * equation via the routine DLAED4 (as called by DLAED3). * This routine also calculates the eigenvectors of the current * problem. * * The final stage consists of computing the updated eigenvectors * directly using the updated eigenvalues. The eigenvectors for * the current problem are multiplied with the eigenvectors from * the overall problem. * * Arguments * ========= * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the eigenvalues of the rank-1-perturbed matrix. * On exit, the eigenvalues of the repaired matrix. * * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) * On entry, the eigenvectors of the rank-1-perturbed matrix. * On exit, the eigenvectors of the repaired tridiagonal matrix. * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max(1,N). * * INDXQ (input/output) INTEGER array, dimension (N) * On entry, the permutation which separately sorts the two * subproblems in D into ascending order. * On exit, the permutation which will reintegrate the * subproblems back into sorted order, * i.e. D( INDXQ( I = 1, N ) ) will be in ascending order. * * RHO (input) DOUBLE PRECISION * The subdiagonal entry used to create the rank-1 modification. * * CUTPNT (input) INTEGER * The location of the last eigenvalue in the leading sub-matrix. * min(1,N) <= CUTPNT <= N/2. * * WORK (workspace) DOUBLE PRECISION array, dimension (4*N + N**2) * * IWORK (workspace) INTEGER array, dimension (4*N) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = 1, an eigenvalue did not converge * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * Modified by Francoise Tisseur, University of Tennessee. * * ===================================================================== * * .. Local Scalars .. INTEGER COLTYP, I, IDLMDA, INDX, INDXC, INDXP, IQ2, IS, $ IW, IZ, K, N1, N2, ZPP1 * .. * .. External Subroutines .. EXTERNAL DCOPY, DLAED2, DLAED3, DLAMRG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( N.LT.0 ) THEN INFO = -1 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -4 ELSE IF( MIN( 1, N / 2 ).GT.CUTPNT .OR. ( N / 2 ).LT.CUTPNT ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED1', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * The following values are integer pointers which indicate * the portion of the workspace * used by a particular array in DLAED2 and DLAED3. * IZ = 1 IDLMDA = IZ + N IW = IDLMDA + N IQ2 = IW + N * INDX = 1 INDXC = INDX + N COLTYP = INDXC + N INDXP = COLTYP + N * * * Form the z-vector which consists of the last row of Q_1 and the * first row of Q_2. * CALL DCOPY( CUTPNT, Q( CUTPNT, 1 ), LDQ, WORK( IZ ), 1 ) ZPP1 = CUTPNT + 1 CALL DCOPY( N-CUTPNT, Q( ZPP1, ZPP1 ), LDQ, WORK( IZ+CUTPNT ), 1 ) * * Deflate eigenvalues. * CALL DLAED2( K, N, CUTPNT, D, Q, LDQ, INDXQ, RHO, WORK( IZ ), $ WORK( IDLMDA ), WORK( IW ), WORK( IQ2 ), $ IWORK( INDX ), IWORK( INDXC ), IWORK( INDXP ), $ IWORK( COLTYP ), INFO ) * IF( INFO.NE.0 ) $ GO TO 20 * * Solve Secular Equation. * IF( K.NE.0 ) THEN IS = ( IWORK( COLTYP )+IWORK( COLTYP+1 ) )*CUTPNT + $ ( IWORK( COLTYP+1 )+IWORK( COLTYP+2 ) )*( N-CUTPNT ) + IQ2 CALL DLAED3( K, N, CUTPNT, D, Q, LDQ, RHO, WORK( IDLMDA ), $ WORK( IQ2 ), IWORK( INDXC ), IWORK( COLTYP ), $ WORK( IW ), WORK( IS ), INFO ) IF( INFO.NE.0 ) $ GO TO 20 * * Prepare the INDXQ sorting permutation. * N1 = K N2 = N - K CALL DLAMRG( N1, N2, D, 1, -1, INDXQ ) ELSE DO 10 I = 1, N INDXQ( I ) = I 10 CONTINUE END IF * 20 CONTINUE RETURN * * End of DLAED1 * END SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, $ Q2, INDX, INDXC, INDXP, COLTYP, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, K, LDQ, N, N1 DOUBLE PRECISION RHO * .. * .. Array Arguments .. INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ), $ INDXQ( * ) DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), $ W( * ), Z( * ) * .. * * Purpose * ======= * * DLAED2 merges the two sets of eigenvalues together into a single * sorted set. Then it tries to deflate the size of the problem. * There are two ways in which deflation can occur: when two or more * eigenvalues are close together or if there is a tiny entry in the * Z vector. For each such occurrence the order of the related secular * equation problem is reduced by one. * * Arguments * ========= * * K (output) INTEGER * The number of non-deflated eigenvalues, and the order of the * related secular equation. 0 <= K <=N. * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * N1 (input) INTEGER * The location of the last eigenvalue in the leading sub-matrix. * min(1,N) <= N1 <= N/2. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, D contains the eigenvalues of the two submatrices to * be combined. * On exit, D contains the trailing (N-K) updated eigenvalues * (those which were deflated) sorted into increasing order. * * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) * On entry, Q contains the eigenvectors of two submatrices in * the two square blocks with corners at (1,1), (N1,N1) * and (N1+1, N1+1), (N,N). * On exit, Q contains the trailing (N-K) updated eigenvectors * (those which were deflated) in its last N-K columns. * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max(1,N). * * INDXQ (input/output) INTEGER array, dimension (N) * The permutation which separately sorts the two sub-problems * in D into ascending order. Note that elements in the second * half of this permutation must first have N1 added to their * values. Destroyed on exit. * * RHO (input/output) DOUBLE PRECISION * On entry, the off-diagonal element associated with the rank-1 * cut which originally split the two submatrices which are now * being recombined. * On exit, RHO has been modified to the value required by * DLAED3. * * Z (input) DOUBLE PRECISION array, dimension (N) * On entry, Z contains the updating vector (the last * row of the first sub-eigenvector matrix and the first row of * the second sub-eigenvector matrix). * On exit, the contents of Z have been destroyed by the updating * process. * * DLAMDA (output) DOUBLE PRECISION array, dimension (N) * A copy of the first K eigenvalues which will be used by * DLAED3 to form the secular equation. * * W (output) DOUBLE PRECISION array, dimension (N) * The first k values of the final deflation-altered z-vector * which will be passed to DLAED3. * * Q2 (output) DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2) * A copy of the first K eigenvectors which will be used by * DLAED3 in a matrix multiply (DGEMM) to solve for the new * eigenvectors. * * INDX (workspace) INTEGER array, dimension (N) * The permutation used to sort the contents of DLAMDA into * ascending order. * * INDXC (output) INTEGER array, dimension (N) * The permutation used to arrange the columns of the deflated * Q matrix into three groups: the first group contains non-zero * elements only at and above N1, the second contains * non-zero elements only below N1, and the third is dense. * * INDXP (workspace) INTEGER array, dimension (N) * The permutation used to place deflated values of D at the end * of the array. INDXP(1:K) points to the nondeflated D-values * and INDXP(K+1:N) points to the deflated eigenvalues. * * COLTYP (workspace/output) INTEGER array, dimension (N) * During execution, a label which will indicate which of the * following types a column in the Q2 matrix is: * 1 : non-zero in the upper half only; * 2 : dense; * 3 : non-zero in the lower half only; * 4 : deflated. * On exit, COLTYP(i) is the number of columns of type i, * for i=1 to 4 only. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * Modified by Francoise Tisseur, University of Tennessee. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, $ TWO = 2.0D0, EIGHT = 8.0D0 ) * .. * .. Local Arrays .. INTEGER CTOT( 4 ), PSM( 4 ) * .. * .. Local Scalars .. INTEGER CT, I, IMAX, IQ1, IQ2, J, JMAX, JS, K2, N1P1, $ N2, NJ, PJ DOUBLE PRECISION C, EPS, S, T, TAU, TOL * .. * .. External Functions .. INTEGER IDAMAX DOUBLE PRECISION DLAMCH, DLAPY2 EXTERNAL IDAMAX, DLAMCH, DLAPY2 * .. * .. External Subroutines .. EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -6 ELSE IF( MIN( 1, ( N / 2 ) ).GT.N1 .OR. ( N / 2 ).LT.N1 ) THEN INFO = -3 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED2', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * N2 = N - N1 N1P1 = N1 + 1 * IF( RHO.LT.ZERO ) THEN CALL DSCAL( N2, MONE, Z( N1P1 ), 1 ) END IF * * Normalize z so that norm(z) = 1. Since z is the concatenation of * two normalized vectors, norm2(z) = sqrt(2). * T = ONE / SQRT( TWO ) CALL DSCAL( N, T, Z, 1 ) * * RHO = ABS( norm(z)**2 * RHO ) * RHO = ABS( TWO*RHO ) * * Sort the eigenvalues into increasing order * DO 10 I = N1P1, N INDXQ( I ) = INDXQ( I ) + N1 10 CONTINUE * * re-integrate the deflated parts from the last pass * DO 20 I = 1, N DLAMDA( I ) = D( INDXQ( I ) ) 20 CONTINUE CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDXC ) DO 30 I = 1, N INDX( I ) = INDXQ( INDXC( I ) ) 30 CONTINUE * * Calculate the allowable deflation tolerance * IMAX = IDAMAX( N, Z, 1 ) JMAX = IDAMAX( N, D, 1 ) EPS = DLAMCH( 'Epsilon' ) TOL = EIGHT*EPS*MAX( ABS( D( JMAX ) ), ABS( Z( IMAX ) ) ) * * If the rank-1 modifier is small enough, no more needs to be done * except to reorganize Q so that its columns correspond with the * elements in D. * IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN K = 0 IQ2 = 1 DO 40 J = 1, N I = INDX( J ) CALL DCOPY( N, Q( 1, I ), 1, Q2( IQ2 ), 1 ) DLAMDA( J ) = D( I ) IQ2 = IQ2 + N 40 CONTINUE CALL DLACPY( 'A', N, N, Q2, N, Q, LDQ ) CALL DCOPY( N, DLAMDA, 1, D, 1 ) GO TO 190 END IF * * If there are multiple eigenvalues then the problem deflates. Here * the number of equal eigenvalues are found. As each equal * eigenvalue is found, an elementary reflector is computed to rotate * the corresponding eigensubspace so that the corresponding * components of Z are zero in this new basis. * DO 50 I = 1, N1 COLTYP( I ) = 1 50 CONTINUE DO 60 I = N1P1, N COLTYP( I ) = 3 60 CONTINUE * * K = 0 K2 = N + 1 DO 70 J = 1, N NJ = INDX( J ) IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN * * Deflate due to small z component. * K2 = K2 - 1 COLTYP( NJ ) = 4 INDXP( K2 ) = NJ IF( J.EQ.N ) $ GO TO 100 ELSE PJ = NJ GO TO 80 END IF 70 CONTINUE 80 CONTINUE J = J + 1 NJ = INDX( J ) IF( J.GT.N ) $ GO TO 100 IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN * * Deflate due to small z component. * K2 = K2 - 1 COLTYP( NJ ) = 4 INDXP( K2 ) = NJ ELSE * * Check if eigenvalues are close enough to allow deflation. * S = Z( PJ ) C = Z( NJ ) * * Find sqrt(a**2+b**2) without overflow or * destructive underflow. * TAU = DLAPY2( C, S ) T = D( NJ ) - D( PJ ) C = C / TAU S = -S / TAU IF( ABS( T*C*S ).LE.TOL ) THEN * * Deflation is possible. * Z( NJ ) = TAU Z( PJ ) = ZERO IF( COLTYP( NJ ).NE.COLTYP( PJ ) ) $ COLTYP( NJ ) = 2 COLTYP( PJ ) = 4 CALL DROT( N, Q( 1, PJ ), 1, Q( 1, NJ ), 1, C, S ) T = D( PJ )*C**2 + D( NJ )*S**2 D( NJ ) = D( PJ )*S**2 + D( NJ )*C**2 D( PJ ) = T K2 = K2 - 1 I = 1 90 CONTINUE IF( K2+I.LE.N ) THEN IF( D( PJ ).LT.D( INDXP( K2+I ) ) ) THEN INDXP( K2+I-1 ) = INDXP( K2+I ) INDXP( K2+I ) = PJ I = I + 1 GO TO 90 ELSE INDXP( K2+I-1 ) = PJ END IF ELSE INDXP( K2+I-1 ) = PJ END IF PJ = NJ ELSE K = K + 1 DLAMDA( K ) = D( PJ ) W( K ) = Z( PJ ) INDXP( K ) = PJ PJ = NJ END IF END IF GO TO 80 100 CONTINUE * * Record the last eigenvalue. * K = K + 1 DLAMDA( K ) = D( PJ ) W( K ) = Z( PJ ) INDXP( K ) = PJ * * Count up the total number of the various types of columns, then * form a permutation which positions the four column types into * four uniform groups (although one or more of these groups may be * empty). * DO 110 J = 1, 4 CTOT( J ) = 0 110 CONTINUE DO 120 J = 1, N CT = COLTYP( J ) CTOT( CT ) = CTOT( CT ) + 1 120 CONTINUE * * PSM(*) = Position in SubMatrix (of types 1 through 4) * PSM( 1 ) = 1 PSM( 2 ) = 1 + CTOT( 1 ) PSM( 3 ) = PSM( 2 ) + CTOT( 2 ) PSM( 4 ) = PSM( 3 ) + CTOT( 3 ) K = N - CTOT( 4 ) * * Fill out the INDXC array so that the permutation which it induces * will place all type-1 columns first, all type-2 columns next, * then all type-3's, and finally all type-4's. * DO 130 J = 1, N JS = INDXP( J ) CT = COLTYP( JS ) INDX( PSM( CT ) ) = JS INDXC( PSM( CT ) ) = J PSM( CT ) = PSM( CT ) + 1 130 CONTINUE * * Sort the eigenvalues and corresponding eigenvectors into DLAMDA * and Q2 respectively. The eigenvalues/vectors which were not * deflated go into the first K slots of DLAMDA and Q2 respectively, * while those which were deflated go into the last N - K slots. * I = 1 IQ1 = 1 IQ2 = 1 + ( CTOT( 1 )+CTOT( 2 ) )*N1 DO 140 J = 1, CTOT( 1 ) JS = INDX( I ) CALL DCOPY( N1, Q( 1, JS ), 1, Q2( IQ1 ), 1 ) Z( I ) = D( JS ) I = I + 1 IQ1 = IQ1 + N1 140 CONTINUE * DO 150 J = 1, CTOT( 2 ) JS = INDX( I ) CALL DCOPY( N1, Q( 1, JS ), 1, Q2( IQ1 ), 1 ) CALL DCOPY( N2, Q( N1+1, JS ), 1, Q2( IQ2 ), 1 ) Z( I ) = D( JS ) I = I + 1 IQ1 = IQ1 + N1 IQ2 = IQ2 + N2 150 CONTINUE * DO 160 J = 1, CTOT( 3 ) JS = INDX( I ) CALL DCOPY( N2, Q( N1+1, JS ), 1, Q2( IQ2 ), 1 ) Z( I ) = D( JS ) I = I + 1 IQ2 = IQ2 + N2 160 CONTINUE * IQ1 = IQ2 DO 170 J = 1, CTOT( 4 ) JS = INDX( I ) CALL DCOPY( N, Q( 1, JS ), 1, Q2( IQ2 ), 1 ) IQ2 = IQ2 + N Z( I ) = D( JS ) I = I + 1 170 CONTINUE * * The deflated eigenvalues and their corresponding vectors go back * into the last N - K slots of D and Q respectively. * CALL DLACPY( 'A', N, CTOT( 4 ), Q2( IQ1 ), N, Q( 1, K+1 ), LDQ ) CALL DCOPY( N-K, Z( K+1 ), 1, D( K+1 ), 1 ) * * Copy CTOT into COLTYP for referencing in DLAED3. * DO 180 J = 1, 4 COLTYP( J ) = CTOT( J ) 180 CONTINUE * 190 CONTINUE RETURN * * End of DLAED2 * END SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, $ CTOT, W, S, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, K, LDQ, N, N1 DOUBLE PRECISION RHO * .. * .. Array Arguments .. INTEGER CTOT( * ), INDX( * ) DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), $ S( * ), W( * ) * .. * * Purpose * ======= * * DLAED3 finds the roots of the secular equation, as defined by the * values in D, W, and RHO, between 1 and K. It makes the * appropriate calls to DLAED4 and then updates the eigenvectors by * multiplying the matrix of eigenvectors of the pair of eigensystems * being combined by the matrix of eigenvectors of the K-by-K system * which is solved here. * * This code makes very mild assumptions about floating point * arithmetic. It will work on machines with a guard digit in * add/subtract, or on those binary machines without guard digits * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. * It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. * * Arguments * ========= * * K (input) INTEGER * The number of terms in the rational function to be solved by * DLAED4. K >= 0. * * N (input) INTEGER * The number of rows and columns in the Q matrix. * N >= K (deflation may result in N>K). * * N1 (input) INTEGER * The location of the last eigenvalue in the leading submatrix. * min(1,N) <= N1 <= N/2. * * D (output) DOUBLE PRECISION array, dimension (N) * D(I) contains the updated eigenvalues for * 1 <= I <= K. * * Q (output) DOUBLE PRECISION array, dimension (LDQ,N) * Initially the first K columns are used as workspace. * On output the columns 1 to K contain * the updated eigenvectors. * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max(1,N). * * RHO (input) DOUBLE PRECISION * The value of the parameter in the rank one update equation. * RHO >= 0 required. * * DLAMDA (input/output) DOUBLE PRECISION array, dimension (K) * The first K elements of this array contain the old roots * of the deflated updating problem. These are the poles * of the secular equation. May be changed on output by * having lowest order bit set to zero on Cray X-MP, Cray Y-MP, * Cray-2, or Cray C-90, as described above. * * Q2 (input) DOUBLE PRECISION array, dimension (LDQ2, N) * The first K columns of this matrix contain the non-deflated * eigenvectors for the split problem. * * INDX (input) INTEGER array, dimension (N) * The permutation used to arrange the columns of the deflated * Q matrix into three groups (see DLAED2). * The rows of the eigenvectors found by DLAED4 must be likewise * permuted before the matrix multiply can take place. * * CTOT (input) INTEGER array, dimension (4) * A count of the total number of the various types of columns * in Q, as described in INDX. The fourth column type is any * column which has been deflated. * * W (input/output) DOUBLE PRECISION array, dimension (K) * The first K elements of this array contain the components * of the deflation-adjusted updating vector. Destroyed on * output. * * S (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K * Will contain the eigenvectors of the repaired matrix which * will be multiplied by the previously accumulated eigenvectors * to update the system. * * LDS (input) INTEGER * The leading dimension of S. LDS >= max(1,K). * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = 1, an eigenvalue did not converge * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * Modified by Francoise Tisseur, University of Tennessee. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) * .. * .. Local Scalars .. INTEGER I, II, IQ2, J, N12, N2, N23 DOUBLE PRECISION TEMP * .. * .. External Functions .. DOUBLE PRECISION DLAMC3, DNRM2 EXTERNAL DLAMC3, DNRM2 * .. * .. External Subroutines .. EXTERNAL DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, SIGN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( K.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.K ) THEN INFO = -2 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED3', -INFO ) RETURN END IF * * Quick return if possible * IF( K.EQ.0 ) $ RETURN * * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can * be computed with high relative accuracy (barring over/underflow). * This is a problem on machines without a guard digit in * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), * which on any of these machines zeros out the bottommost * bit of DLAMDA(I) if it is 1; this makes the subsequent * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation * occurs. On binary machines with a guard digit (almost all * machines) it does not change DLAMDA(I) at all. On hexadecimal * and decimal machines with a guard digit, it slightly * changes the bottommost bits of DLAMDA(I). It does not account * for hexadecimal or decimal machines without guard digits * (we know of none). We use a subroutine call to compute * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating * this code. * DO 10 I = 1, K DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) 10 CONTINUE * DO 20 J = 1, K CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) * * If the zero finder fails, the computation is terminated. * IF( INFO.NE.0 ) $ GO TO 120 20 CONTINUE * IF( K.EQ.1 ) $ GO TO 110 IF( K.EQ.2 ) THEN DO 30 J = 1, K W( 1 ) = Q( 1, J ) W( 2 ) = Q( 2, J ) II = INDX( 1 ) Q( 1, J ) = W( II ) II = INDX( 2 ) Q( 2, J ) = W( II ) 30 CONTINUE GO TO 110 END IF * * Compute updated W. * CALL DCOPY( K, W, 1, S, 1 ) * * Initialize W(I) = Q(I,I) * CALL DCOPY( K, Q, LDQ+1, W, 1 ) DO 60 J = 1, K DO 40 I = 1, J - 1 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 40 CONTINUE DO 50 I = J + 1, K W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 50 CONTINUE 60 CONTINUE DO 70 I = 1, K W( I ) = SIGN( SQRT( -W( I ) ), S( I ) ) 70 CONTINUE * * Compute eigenvectors of the modified rank-1 modification. * DO 100 J = 1, K DO 80 I = 1, K S( I ) = W( I ) / Q( I, J ) 80 CONTINUE TEMP = DNRM2( K, S, 1 ) DO 90 I = 1, K II = INDX( I ) Q( I, J ) = S( II ) / TEMP 90 CONTINUE 100 CONTINUE * * Compute the updated eigenvectors. * 110 CONTINUE * N2 = N - N1 N12 = CTOT( 1 ) + CTOT( 2 ) N23 = CTOT( 2 ) + CTOT( 3 ) * CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 ) IQ2 = N1*N12 + 1 IF( N23.NE.0 ) THEN CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23, $ ZERO, Q( N1+1, 1 ), LDQ ) ELSE CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ ) END IF * CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 ) IF( N12.NE.0 ) THEN CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q, $ LDQ ) ELSE CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ ) END IF * * 120 CONTINUE RETURN * * End of DLAED3 * END SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER I, INFO, N DOUBLE PRECISION DLAM, RHO * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), DELTA( * ), Z( * ) * .. * * Purpose * ======= * * This subroutine computes the I-th updated eigenvalue of a symmetric * rank-one modification to a diagonal matrix whose elements are * given in the array d, and that * * D(i) < D(j) for i < j * * and that RHO > 0. This is arranged by the calling routine, and is * no loss in generality. The rank-one modified system is thus * * diag( D ) + RHO * Z * Z_transpose. * * where we assume the Euclidean norm of Z is 1. * * The method consists of approximating the rational functions in the * secular equation by simpler interpolating rational functions. * * Arguments * ========= * * N (input) INTEGER * The length of all arrays. * * I (input) INTEGER * The index of the eigenvalue to be computed. 1 <= I <= N. * * D (input) DOUBLE PRECISION array, dimension (N) * The original eigenvalues. It is assumed that they are in * order, D(I) < D(J) for I < J. * * Z (input) DOUBLE PRECISION array, dimension (N) * The components of the updating vector. * * DELTA (output) DOUBLE PRECISION array, dimension (N) * If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th * component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5 * for detail. The vector DELTA contains the information necessary * to construct the eigenvectors by DLAED3 and DLAED9. * * RHO (input) DOUBLE PRECISION * The scalar in the symmetric updating formula. * * DLAM (output) DOUBLE PRECISION * The computed lambda_I, the I-th updated eigenvalue. * * INFO (output) INTEGER * = 0: successful exit * > 0: if INFO = 1, the updating process failed. * * Internal Parameters * =================== * * Logical variable ORGATI (origin-at-i?) is used for distinguishing * whether D(i) or D(i+1) is treated as the origin. * * ORGATI = .true. origin at i * ORGATI = .false. origin at i+1 * * Logical variable SWTCH3 (switch-for-3-poles?) is for noting * if we are working with THREE poles! * * MAXIT is the maximum number of iterations allowed for each * eigenvalue. * * Further Details * =============== * * Based on contributions by * Ren-Cang Li, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 30 ) DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0, $ TEN = 10.0D0 ) * .. * .. Local Scalars .. LOGICAL ORGATI, SWTCH, SWTCH3 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW, $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI, $ RHOINV, TAU, TEMP, TEMP1, W * .. * .. Local Arrays .. DOUBLE PRECISION ZZ( 3 ) * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. External Subroutines .. EXTERNAL DLAED5, DLAED6 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * * Since this routine is called in an inner loop, we do no argument * checking. * * Quick return for N=1 and 2. * INFO = 0 IF( N.EQ.1 ) THEN * * Presumably, I=1 upon entry * DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 ) DELTA( 1 ) = ONE RETURN END IF IF( N.EQ.2 ) THEN CALL DLAED5( I, D, Z, DELTA, RHO, DLAM ) RETURN END IF * * Compute machine epsilon * EPS = DLAMCH( 'Epsilon' ) RHOINV = ONE / RHO * * The case I = N * IF( I.EQ.N ) THEN * * Initialize some basic variables * II = N - 1 NITER = 1 * * Calculate initial guess * MIDPT = RHO / TWO * * If ||Z||_2 is not one, then TEMP should be set to * RHO * ||Z||_2^2 / TWO * DO 10 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 10 CONTINUE * PSI = ZERO DO 20 J = 1, N - 2 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 20 CONTINUE * C = RHOINV + PSI W = C + Z( II )*Z( II ) / DELTA( II ) + $ Z( N )*Z( N ) / DELTA( N ) * IF( W.LE.ZERO ) THEN TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) + $ Z( N )*Z( N ) / RHO IF( C.LE.TEMP ) THEN TAU = RHO ELSE DEL = D( N ) - D( N-1 ) A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) B = Z( N )*Z( N )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) ELSE TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) END IF END IF * * It can be proved that * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO * DLTLB = MIDPT DLTUB = RHO ELSE DEL = D( N ) - D( N-1 ) A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) B = Z( N )*Z( N )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) ELSE TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) END IF * * It can be proved that * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 * DLTLB = ZERO DLTUB = MIDPT END IF * DO 30 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - TAU 30 CONTINUE * * Evaluate PSI and the derivative DPSI * DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 40 J = 1, II TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 40 CONTINUE ERRETM = ABS( ERRETM ) * * Evaluate PHI and the derivative DPHI * TEMP = Z( N ) / DELTA( N ) PHI = Z( N )*TEMP DPHI = TEMP*TEMP ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + $ ABS( TAU )*( DPSI+DPHI ) * W = RHOINV + PHI + PSI * * Test for convergence * IF( ABS( W ).LE.EPS*ERRETM ) THEN DLAM = D( I ) + TAU GO TO 250 END IF * IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF * * Calculate the new step * NITER = NITER + 1 C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI A = ( DELTA( N-1 )+DELTA( N ) )*W - $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) B = DELTA( N-1 )*DELTA( N )*W IF( C.LT.ZERO ) $ C = ABS( C ) IF( C.EQ.ZERO ) THEN * ETA = B/A * ETA = RHO - TAU ETA = DLTUB - TAU ELSE IF( A.GE.ZERO ) THEN ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF * * Note, eta should be positive if w is negative, and * eta should be negative otherwise. However, * if for some reason caused by roundoff, eta*w > 0, * we simply use one Newton step instead. This way * will guarantee eta*w < 0. * IF( W*ETA.GT.ZERO ) $ ETA = -W / ( DPSI+DPHI ) TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF DO 50 J = 1, N DELTA( J ) = DELTA( J ) - ETA 50 CONTINUE * TAU = TAU + ETA * * Evaluate PSI and the derivative DPSI * DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 60 J = 1, II TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 60 CONTINUE ERRETM = ABS( ERRETM ) * * Evaluate PHI and the derivative DPHI * TEMP = Z( N ) / DELTA( N ) PHI = Z( N )*TEMP DPHI = TEMP*TEMP ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + $ ABS( TAU )*( DPSI+DPHI ) * W = RHOINV + PHI + PSI * * Main loop to update the values of the array DELTA * ITER = NITER + 1 * DO 90 NITER = ITER, MAXIT * * Test for convergence * IF( ABS( W ).LE.EPS*ERRETM ) THEN DLAM = D( I ) + TAU GO TO 250 END IF * IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF * * Calculate the new step * C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI A = ( DELTA( N-1 )+DELTA( N ) )*W - $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) B = DELTA( N-1 )*DELTA( N )*W IF( A.GE.ZERO ) THEN ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF * * Note, eta should be positive if w is negative, and * eta should be negative otherwise. However, * if for some reason caused by roundoff, eta*w > 0, * we simply use one Newton step instead. This way * will guarantee eta*w < 0. * IF( W*ETA.GT.ZERO ) $ ETA = -W / ( DPSI+DPHI ) TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF DO 70 J = 1, N DELTA( J ) = DELTA( J ) - ETA 70 CONTINUE * TAU = TAU + ETA * * Evaluate PSI and the derivative DPSI * DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 80 J = 1, II TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 80 CONTINUE ERRETM = ABS( ERRETM ) * * Evaluate PHI and the derivative DPHI * TEMP = Z( N ) / DELTA( N ) PHI = Z( N )*TEMP DPHI = TEMP*TEMP ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + $ ABS( TAU )*( DPSI+DPHI ) * W = RHOINV + PHI + PSI 90 CONTINUE * * Return with INFO = 1, NITER = MAXIT and not converged * INFO = 1 DLAM = D( I ) + TAU GO TO 250 * * End for the case I = N * ELSE * * The case for I < N * NITER = 1 IP1 = I + 1 * * Calculate initial guess * DEL = D( IP1 ) - D( I ) MIDPT = DEL / TWO DO 100 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 100 CONTINUE * PSI = ZERO DO 110 J = 1, I - 1 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 110 CONTINUE * PHI = ZERO DO 120 J = N, I + 2, -1 PHI = PHI + Z( J )*Z( J ) / DELTA( J ) 120 CONTINUE C = RHOINV + PSI + PHI W = C + Z( I )*Z( I ) / DELTA( I ) + $ Z( IP1 )*Z( IP1 ) / DELTA( IP1 ) * IF( W.GT.ZERO ) THEN * * d(i)< the ith eigenvalue < (d(i)+d(i+1))/2 * * We choose d(i) as origin. * ORGATI = .TRUE. A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 ) B = Z( I )*Z( I )*DEL IF( A.GT.ZERO ) THEN TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) ELSE TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) END IF DLTLB = ZERO DLTUB = MIDPT ELSE * * (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1) * * We choose d(i+1) as origin. * ORGATI = .FALSE. A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 ) B = Z( IP1 )*Z( IP1 )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) ) ELSE TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C ) END IF DLTLB = -MIDPT DLTUB = ZERO END IF * IF( ORGATI ) THEN DO 130 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - TAU 130 CONTINUE ELSE DO 140 J = 1, N DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU 140 CONTINUE END IF IF( ORGATI ) THEN II = I ELSE II = I + 1 END IF IIM1 = II - 1 IIP1 = II + 1 * * Evaluate PSI and the derivative DPSI * DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 150 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 150 CONTINUE ERRETM = ABS( ERRETM ) * * Evaluate PHI and the derivative DPHI * DPHI = ZERO PHI = ZERO DO 160 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 160 CONTINUE * W = RHOINV + PHI + PSI * * W is the value of the secular function with * its ii-th element removed. * SWTCH3 = .FALSE. IF( ORGATI ) THEN IF( W.LT.ZERO ) $ SWTCH3 = .TRUE. ELSE IF( W.GT.ZERO ) $ SWTCH3 = .TRUE. END IF IF( II.EQ.1 .OR. II.EQ.N ) $ SWTCH3 = .FALSE. * TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = W + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + $ THREE*ABS( TEMP ) + ABS( TAU )*DW * * Test for convergence * IF( ABS( W ).LE.EPS*ERRETM ) THEN IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF GO TO 250 END IF * IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF * * Calculate the new step * NITER = NITER + 1 IF( .NOT.SWTCH3 ) THEN IF( ORGATI ) THEN C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )* $ ( Z( I ) / DELTA( I ) )**2 ELSE C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* $ ( Z( IP1 ) / DELTA( IP1 ) )**2 END IF A = ( DELTA( I )+DELTA( IP1 ) )*W - $ DELTA( I )*DELTA( IP1 )*DW B = DELTA( I )*DELTA( IP1 )*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN IF( ORGATI ) THEN A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )* $ ( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )* $ ( DPSI+DPHI ) END IF END IF ETA = B / A ELSE IF( A.LE.ZERO ) THEN ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF ELSE * * Interpolation using THREE most relevant poles * TEMP = RHOINV + PSI + PHI IF( ORGATI ) THEN TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* $ ( ( DPSI-TEMP1 )+DPHI ) ELSE TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* $ ( DPSI+( DPHI-TEMP1 ) ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF ZZ( 2 ) = Z( II )*Z( II ) CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, $ INFO ) IF( INFO.NE.0 ) $ GO TO 250 END IF * * Note, eta should be positive if w is negative, and * eta should be negative otherwise. However, * if for some reason caused by roundoff, eta*w > 0, * we simply use one Newton step instead. This way * will guarantee eta*w < 0. * IF( W*ETA.GE.ZERO ) $ ETA = -W / DW TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF * PREW = W * DO 180 J = 1, N DELTA( J ) = DELTA( J ) - ETA 180 CONTINUE * * Evaluate PSI and the derivative DPSI * DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 190 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 190 CONTINUE ERRETM = ABS( ERRETM ) * * Evaluate PHI and the derivative DPHI * DPHI = ZERO PHI = ZERO DO 200 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 200 CONTINUE * TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = RHOINV + PHI + PSI + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + $ THREE*ABS( TEMP ) + ABS( TAU+ETA )*DW * SWTCH = .FALSE. IF( ORGATI ) THEN IF( -W.GT.ABS( PREW ) / TEN ) $ SWTCH = .TRUE. ELSE IF( W.GT.ABS( PREW ) / TEN ) $ SWTCH = .TRUE. END IF * TAU = TAU + ETA * * Main loop to update the values of the array DELTA * ITER = NITER + 1 * DO 240 NITER = ITER, MAXIT * * Test for convergence * IF( ABS( W ).LE.EPS*ERRETM ) THEN IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF GO TO 250 END IF * IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF * * Calculate the new step * IF( .NOT.SWTCH3 ) THEN IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN C = W - DELTA( IP1 )*DW - $ ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2 ELSE C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* $ ( Z( IP1 ) / DELTA( IP1 ) )**2 END IF ELSE TEMP = Z( II ) / DELTA( II ) IF( ORGATI ) THEN DPSI = DPSI + TEMP*TEMP ELSE DPHI = DPHI + TEMP*TEMP END IF C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI END IF A = ( DELTA( I )+DELTA( IP1 ) )*W - $ DELTA( I )*DELTA( IP1 )*DW B = DELTA( I )*DELTA( IP1 )*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN A = Z( I )*Z( I ) + DELTA( IP1 )* $ DELTA( IP1 )*( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + $ DELTA( I )*DELTA( I )*( DPSI+DPHI ) END IF ELSE A = DELTA( I )*DELTA( I )*DPSI + $ DELTA( IP1 )*DELTA( IP1 )*DPHI END IF END IF ETA = B / A ELSE IF( A.LE.ZERO ) THEN ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF ELSE * * Interpolation using THREE most relevant poles * TEMP = RHOINV + PSI + PHI IF( SWTCH ) THEN C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI ELSE IF( ORGATI ) THEN TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* $ ( ( DPSI-TEMP1 )+DPHI ) ELSE TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* $ ( DPSI+( DPHI-TEMP1 ) ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF END IF CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, $ INFO ) IF( INFO.NE.0 ) $ GO TO 250 END IF * * Note, eta should be positive if w is negative, and * eta should be negative otherwise. However, * if for some reason caused by roundoff, eta*w > 0, * we simply use one Newton step instead. This way * will guarantee eta*w < 0. * IF( W*ETA.GE.ZERO ) $ ETA = -W / DW TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF * DO 210 J = 1, N DELTA( J ) = DELTA( J ) - ETA 210 CONTINUE * TAU = TAU + ETA PREW = W * * Evaluate PSI and the derivative DPSI * DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 220 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 220 CONTINUE ERRETM = ABS( ERRETM ) * * Evaluate PHI and the derivative DPHI * DPHI = ZERO PHI = ZERO DO 230 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 230 CONTINUE * TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = RHOINV + PHI + PSI + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + $ THREE*ABS( TEMP ) + ABS( TAU )*DW IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN ) $ SWTCH = .NOT.SWTCH * 240 CONTINUE * * Return with INFO = 1, NITER = MAXIT and not converged * INFO = 1 IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF * END IF * 250 CONTINUE * RETURN * * End of DLAED4 * END SUBROUTINE DLAED5( I, D, Z, DELTA, RHO, DLAM ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER I DOUBLE PRECISION DLAM, RHO * .. * .. Array Arguments .. DOUBLE PRECISION D( 2 ), DELTA( 2 ), Z( 2 ) * .. * * Purpose * ======= * * This subroutine computes the I-th eigenvalue of a symmetric rank-one * modification of a 2-by-2 diagonal matrix * * diag( D ) + RHO * Z * transpose(Z) . * * The diagonal elements in the array D are assumed to satisfy * * D(i) < D(j) for i < j . * * We also assume RHO > 0 and that the Euclidean norm of the vector * Z is one. * * Arguments * ========= * * I (input) INTEGER * The index of the eigenvalue to be computed. I = 1 or I = 2. * * D (input) DOUBLE PRECISION array, dimension (2) * The original eigenvalues. We assume D(1) < D(2). * * Z (input) DOUBLE PRECISION array, dimension (2) * The components of the updating vector. * * DELTA (output) DOUBLE PRECISION array, dimension (2) * The vector DELTA contains the information necessary * to construct the eigenvectors. * * RHO (input) DOUBLE PRECISION * The scalar in the symmetric updating formula. * * DLAM (output) DOUBLE PRECISION * The computed lambda_I, the I-th updated eigenvalue. * * Further Details * =============== * * Based on contributions by * Ren-Cang Li, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, FOUR PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ FOUR = 4.0D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION B, C, DEL, TAU, TEMP, W * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. * DEL = D( 2 ) - D( 1 ) IF( I.EQ.1 ) THEN W = ONE + TWO*RHO*( Z( 2 )*Z( 2 )-Z( 1 )*Z( 1 ) ) / DEL IF( W.GT.ZERO ) THEN B = DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) C = RHO*Z( 1 )*Z( 1 )*DEL * * B > ZERO, always * TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) ) DLAM = D( 1 ) + TAU DELTA( 1 ) = -Z( 1 ) / TAU DELTA( 2 ) = Z( 2 ) / ( DEL-TAU ) ELSE B = -DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) C = RHO*Z( 2 )*Z( 2 )*DEL IF( B.GT.ZERO ) THEN TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) ) ELSE TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO END IF DLAM = D( 2 ) + TAU DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) DELTA( 2 ) = -Z( 2 ) / TAU END IF TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) DELTA( 1 ) = DELTA( 1 ) / TEMP DELTA( 2 ) = DELTA( 2 ) / TEMP ELSE * * Now I=2 * B = -DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) ) C = RHO*Z( 2 )*Z( 2 )*DEL IF( B.GT.ZERO ) THEN TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO ELSE TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) ) END IF DLAM = D( 2 ) + TAU DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) DELTA( 2 ) = -Z( 2 ) / TAU TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) DELTA( 1 ) = DELTA( 1 ) / TEMP DELTA( 2 ) = DELTA( 2 ) / TEMP END IF RETURN * * End OF DLAED5 * END SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) * * -- LAPACK routine (version 3.1.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * February 2007 * * .. Scalar Arguments .. LOGICAL ORGATI INTEGER INFO, KNITER DOUBLE PRECISION FINIT, RHO, TAU * .. * .. Array Arguments .. DOUBLE PRECISION D( 3 ), Z( 3 ) * .. * * Purpose * ======= * * DLAED6 computes the positive or negative root (closest to the origin) * of * z(1) z(2) z(3) * f(x) = rho + --------- + ---------- + --------- * d(1)-x d(2)-x d(3)-x * * It is assumed that * * if ORGATI = .true. the root is between d(2) and d(3); * otherwise it is between d(1) and d(2) * * This routine will be called by DLAED4 when necessary. In most cases, * the root sought is the smallest in magnitude, though it might not be * in some extremely rare situations. * * Arguments * ========= * * KNITER (input) INTEGER * Refer to DLAED4 for its significance. * * ORGATI (input) LOGICAL * If ORGATI is true, the needed root is between d(2) and * d(3); otherwise it is between d(1) and d(2). See * DLAED4 for further details. * * RHO (input) DOUBLE PRECISION * Refer to the equation f(x) above. * * D (input) DOUBLE PRECISION array, dimension (3) * D satisfies d(1) < d(2) < d(3). * * Z (input) DOUBLE PRECISION array, dimension (3) * Each of the elements in z must be positive. * * FINIT (input) DOUBLE PRECISION * The value of f at 0. It is more accurate than the one * evaluated inside this routine (if someone wants to do * so). * * TAU (output) DOUBLE PRECISION * The root of the equation f(x). * * INFO (output) INTEGER * = 0: successful exit * > 0: if INFO = 1, failure to converge * * Further Details * =============== * * 30/06/99: Based on contributions by * Ren-Cang Li, Computer Science Division, University of California * at Berkeley, USA * * 10/02/03: This version has a few statements commented out for thread * safety (machine parameters are computed on each entry). SJH. * * 05/10/06: Modified from a new version of Ren-Cang Li, use * Gragg-Thornton-Warner cubic convergent scheme for better stability. * * ===================================================================== * * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 40 ) DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 ) * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Local Arrays .. DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 ) * .. * .. Local Scalars .. LOGICAL SCALE INTEGER I, ITER, NITER DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F, $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1, $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, $ LBD, UBD * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT * .. * .. Executable Statements .. * INFO = 0 * IF( ORGATI ) THEN LBD = D(2) UBD = D(3) ELSE LBD = D(1) UBD = D(2) END IF IF( FINIT .LT. ZERO )THEN LBD = ZERO ELSE UBD = ZERO END IF * NITER = 1 TAU = ZERO IF( KNITER.EQ.2 ) THEN IF( ORGATI ) THEN TEMP = ( D( 3 )-D( 2 ) ) / TWO C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP ) A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 ) B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 ) ELSE TEMP = ( D( 1 )-D( 2 ) ) / TWO C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP ) A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 ) B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 ) END IF TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) A = A / TEMP B = B / TEMP C = C / TEMP IF( C.EQ.ZERO ) THEN TAU = B / A ELSE IF( A.LE.ZERO ) THEN TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF IF( TAU .LT. LBD .OR. TAU .GT. UBD ) $ TAU = ( LBD+UBD )/TWO IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN TAU = ZERO ELSE TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) + $ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) + $ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) ) IF( TEMP .LE. ZERO )THEN LBD = TAU ELSE UBD = TAU END IF IF( ABS( FINIT ).LE.ABS( TEMP ) ) $ TAU = ZERO END IF END IF * * get machine parameters for possible scaling to avoid overflow * * modified by Sven: parameters SMALL1, SMINV1, SMALL2, * SMINV2, EPS are not SAVEd anymore between one call to the * others but recomputed at each call * EPS = DLAMCH( 'Epsilon' ) BASE = DLAMCH( 'Base' ) SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) / $ THREE ) ) SMINV1 = ONE / SMALL1 SMALL2 = SMALL1*SMALL1 SMINV2 = SMINV1*SMINV1 * * Determine if scaling of inputs necessary to avoid overflow * when computing 1/TEMP**3 * IF( ORGATI ) THEN TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) ) ELSE TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) ) END IF SCALE = .FALSE. IF( TEMP.LE.SMALL1 ) THEN SCALE = .TRUE. IF( TEMP.LE.SMALL2 ) THEN * * Scale up by power of radix nearest 1/SAFMIN**(2/3) * SCLFAC = SMINV2 SCLINV = SMALL2 ELSE * * Scale up by power of radix nearest 1/SAFMIN**(1/3) * SCLFAC = SMINV1 SCLINV = SMALL1 END IF * * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1) * DO 10 I = 1, 3 DSCALE( I ) = D( I )*SCLFAC ZSCALE( I ) = Z( I )*SCLFAC 10 CONTINUE TAU = TAU*SCLFAC LBD = LBD*SCLFAC UBD = UBD*SCLFAC ELSE * * Copy D and Z to DSCALE and ZSCALE * DO 20 I = 1, 3 DSCALE( I ) = D( I ) ZSCALE( I ) = Z( I ) 20 CONTINUE END IF * FC = ZERO DF = ZERO DDF = ZERO DO 30 I = 1, 3 TEMP = ONE / ( DSCALE( I )-TAU ) TEMP1 = ZSCALE( I )*TEMP TEMP2 = TEMP1*TEMP TEMP3 = TEMP2*TEMP FC = FC + TEMP1 / DSCALE( I ) DF = DF + TEMP2 DDF = DDF + TEMP3 30 CONTINUE F = FINIT + TAU*FC * IF( ABS( F ).LE.ZERO ) $ GO TO 60 IF( F .LE. ZERO )THEN LBD = TAU ELSE UBD = TAU END IF * * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent * scheme * * It is not hard to see that * * 1) Iterations will go up monotonically * if FINIT < 0; * * 2) Iterations will go down monotonically * if FINIT > 0. * ITER = NITER + 1 * DO 50 NITER = ITER, MAXIT * IF( ORGATI ) THEN TEMP1 = DSCALE( 2 ) - TAU TEMP2 = DSCALE( 3 ) - TAU ELSE TEMP1 = DSCALE( 1 ) - TAU TEMP2 = DSCALE( 2 ) - TAU END IF A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF B = TEMP1*TEMP2*F C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) A = A / TEMP B = B / TEMP C = C / TEMP IF( C.EQ.ZERO ) THEN ETA = B / A ELSE IF( A.LE.ZERO ) THEN ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF IF( F*ETA.GE.ZERO ) THEN ETA = -F / DF END IF * TAU = TAU + ETA IF( TAU .LT. LBD .OR. TAU .GT. UBD ) $ TAU = ( LBD + UBD )/TWO * FC = ZERO ERRETM = ZERO DF = ZERO DDF = ZERO DO 40 I = 1, 3 TEMP = ONE / ( DSCALE( I )-TAU ) TEMP1 = ZSCALE( I )*TEMP TEMP2 = TEMP1*TEMP TEMP3 = TEMP2*TEMP TEMP4 = TEMP1 / DSCALE( I ) FC = FC + TEMP4 ERRETM = ERRETM + ABS( TEMP4 ) DF = DF + TEMP2 DDF = DDF + TEMP3 40 CONTINUE F = FINIT + TAU*FC ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + $ ABS( TAU )*DF IF( ABS( F ).LE.EPS*ERRETM ) $ GO TO 60 IF( F .LE. ZERO )THEN LBD = TAU ELSE UBD = TAU END IF 50 CONTINUE INFO = 1 60 CONTINUE * * Undo scaling * IF( SCALE ) $ TAU = TAU*SCLINV RETURN * * End of DLAED6 * END SUBROUTINE DLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, $ LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, $ PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, $ INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N, $ QSIZ, TLVLS DOUBLE PRECISION RHO * .. * .. Array Arguments .. INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ), $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * ) DOUBLE PRECISION D( * ), GIVNUM( 2, * ), Q( LDQ, * ), $ QSTORE( * ), WORK( * ) * .. * * Purpose * ======= * * DLAED7 computes the updated eigensystem of a diagonal * matrix after modification by a rank-one symmetric matrix. This * routine is used only for the eigenproblem which requires all * eigenvalues and optionally eigenvectors of a dense symmetric matrix * that has been reduced to tridiagonal form. DLAED1 handles * the case in which all eigenvalues and eigenvectors of a symmetric * tridiagonal matrix are desired. * * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) * * where Z = Q'u, u is a vector of length N with ones in the * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. * * The eigenvectors of the original matrix are stored in Q, and the * eigenvalues are in D. The algorithm consists of three stages: * * The first stage consists of deflating the size of the problem * when there are multiple eigenvalues or if there is a zero in * the Z vector. For each such occurence the dimension of the * secular equation problem is reduced by one. This stage is * performed by the routine DLAED8. * * The second stage consists of calculating the updated * eigenvalues. This is done by finding the roots of the secular * equation via the routine DLAED4 (as called by DLAED9). * This routine also calculates the eigenvectors of the current * problem. * * The final stage consists of computing the updated eigenvectors * directly using the updated eigenvalues. The eigenvectors for * the current problem are multiplied with the eigenvectors from * the overall problem. * * Arguments * ========= * * ICOMPQ (input) INTEGER * = 0: Compute eigenvalues only. * = 1: Compute eigenvectors of original dense symmetric matrix * also. On entry, Q contains the orthogonal matrix used * to reduce the original matrix to tridiagonal form. * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * QSIZ (input) INTEGER * The dimension of the orthogonal matrix used to reduce * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. * * TLVLS (input) INTEGER * The total number of merging levels in the overall divide and * conquer tree. * * CURLVL (input) INTEGER * The current level in the overall merge routine, * 0 <= CURLVL <= TLVLS. * * CURPBM (input) INTEGER * The current problem in the current level in the overall * merge routine (counting from upper left to lower right). * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the eigenvalues of the rank-1-perturbed matrix. * On exit, the eigenvalues of the repaired matrix. * * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) * On entry, the eigenvectors of the rank-1-perturbed matrix. * On exit, the eigenvectors of the repaired tridiagonal matrix. * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max(1,N). * * INDXQ (output) INTEGER array, dimension (N) * The permutation which will reintegrate the subproblem just * solved back into sorted order, i.e., D( INDXQ( I = 1, N ) ) * will be in ascending order. * * RHO (input) DOUBLE PRECISION * The subdiagonal element used to create the rank-1 * modification. * * CUTPNT (input) INTEGER * Contains the location of the last eigenvalue in the leading * sub-matrix. min(1,N) <= CUTPNT <= N. * * QSTORE (input/output) DOUBLE PRECISION array, dimension (N**2+1) * Stores eigenvectors of submatrices encountered during * divide and conquer, packed together. QPTR points to * beginning of the submatrices. * * QPTR (input/output) INTEGER array, dimension (N+2) * List of indices pointing to beginning of submatrices stored * in QSTORE. The submatrices are numbered starting at the * bottom left of the divide and conquer tree, from left to * right and bottom to top. * * PRMPTR (input) INTEGER array, dimension (N lg N) * Contains a list of pointers which indicate where in PERM a * level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) * indicates the size of the permutation and also the size of * the full, non-deflated problem. * * PERM (input) INTEGER array, dimension (N lg N) * Contains the permutations (from deflation and sorting) to be * applied to each eigenblock. * * GIVPTR (input) INTEGER array, dimension (N lg N) * Contains a list of pointers which indicate where in GIVCOL a * level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) * indicates the number of Givens rotations. * * GIVCOL (input) INTEGER array, dimension (2, N lg N) * Each pair of numbers indicates a pair of columns to take place * in a Givens rotation. * * GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N) * Each number indicates the S value to be used in the * corresponding Givens rotation. * * WORK (workspace) DOUBLE PRECISION array, dimension (3*N+QSIZ*N) * * IWORK (workspace) INTEGER array, dimension (4*N) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = 1, an eigenvalue did not converge * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) * .. * .. Local Scalars .. INTEGER COLTYP, CURR, I, IDLMDA, INDX, INDXC, INDXP, $ IQ2, IS, IW, IZ, K, LDQ2, N1, N2, PTR * .. * .. External Subroutines .. EXTERNAL DGEMM, DLAED8, DLAED9, DLAEDA, DLAMRG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN INFO = -4 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -9 ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN INFO = -12 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED7', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * The following values are for bookkeeping purposes only. They are * integer pointers which indicate the portion of the workspace * used by a particular array in DLAED8 and DLAED9. * IF( ICOMPQ.EQ.1 ) THEN LDQ2 = QSIZ ELSE LDQ2 = N END IF * IZ = 1 IDLMDA = IZ + N IW = IDLMDA + N IQ2 = IW + N IS = IQ2 + N*LDQ2 * INDX = 1 INDXC = INDX + N COLTYP = INDXC + N INDXP = COLTYP + N * * Form the z-vector which consists of the last row of Q_1 and the * first row of Q_2. * PTR = 1 + 2**TLVLS DO 10 I = 1, CURLVL - 1 PTR = PTR + 2**( TLVLS-I ) 10 CONTINUE CURR = PTR + CURPBM CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, $ GIVCOL, GIVNUM, QSTORE, QPTR, WORK( IZ ), $ WORK( IZ+N ), INFO ) * * When solving the final problem, we no longer need the stored data, * so we will overwrite the data from this level onto the previously * used storage space. * IF( CURLVL.EQ.TLVLS ) THEN QPTR( CURR ) = 1 PRMPTR( CURR ) = 1 GIVPTR( CURR ) = 1 END IF * * Sort and Deflate eigenvalues. * CALL DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, $ WORK( IZ ), WORK( IDLMDA ), WORK( IQ2 ), LDQ2, $ WORK( IW ), PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ), $ GIVCOL( 1, GIVPTR( CURR ) ), $ GIVNUM( 1, GIVPTR( CURR ) ), IWORK( INDXP ), $ IWORK( INDX ), INFO ) PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR ) * * Solve Secular Equation. * IF( K.NE.0 ) THEN CALL DLAED9( K, 1, K, N, D, WORK( IS ), K, RHO, WORK( IDLMDA ), $ WORK( IW ), QSTORE( QPTR( CURR ) ), K, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( ICOMPQ.EQ.1 ) THEN CALL DGEMM( 'N', 'N', QSIZ, K, K, ONE, WORK( IQ2 ), LDQ2, $ QSTORE( QPTR( CURR ) ), K, ZERO, Q, LDQ ) END IF QPTR( CURR+1 ) = QPTR( CURR ) + K**2 * * Prepare the INDXQ sorting permutation. * N1 = K N2 = N - K CALL DLAMRG( N1, N2, D, 1, -1, INDXQ ) ELSE QPTR( CURR+1 ) = QPTR( CURR ) DO 20 I = 1, N INDXQ( I ) = I 20 CONTINUE END IF * 30 CONTINUE RETURN * * End of DLAED7 * END SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, $ CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, $ GIVCOL, GIVNUM, INDXP, INDX, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER CUTPNT, GIVPTR, ICOMPQ, INFO, K, LDQ, LDQ2, N, $ QSIZ DOUBLE PRECISION RHO * .. * .. Array Arguments .. INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), $ INDXQ( * ), PERM( * ) DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) * .. * * Purpose * ======= * * DLAED8 merges the two sets of eigenvalues together into a single * sorted set. Then it tries to deflate the size of the problem. * There are two ways in which deflation can occur: when two or more * eigenvalues are close together or if there is a tiny element in the * Z vector. For each such occurrence the order of the related secular * equation problem is reduced by one. * * Arguments * ========= * * ICOMPQ (input) INTEGER * = 0: Compute eigenvalues only. * = 1: Compute eigenvectors of original dense symmetric matrix * also. On entry, Q contains the orthogonal matrix used * to reduce the original matrix to tridiagonal form. * * K (output) INTEGER * The number of non-deflated eigenvalues, and the order of the * related secular equation. * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * QSIZ (input) INTEGER * The dimension of the orthogonal matrix used to reduce * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the eigenvalues of the two submatrices to be * combined. On exit, the trailing (N-K) updated eigenvalues * (those which were deflated) sorted into increasing order. * * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) * If ICOMPQ = 0, Q is not referenced. Otherwise, * on entry, Q contains the eigenvectors of the partially solved * system which has been previously updated in matrix * multiplies with other partially solved eigensystems. * On exit, Q contains the trailing (N-K) updated eigenvectors * (those which were deflated) in its last N-K columns. * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max(1,N). * * INDXQ (input) INTEGER array, dimension (N) * The permutation which separately sorts the two sub-problems * in D into ascending order. Note that elements in the second * half of this permutation must first have CUTPNT added to * their values in order to be accurate. * * RHO (input/output) DOUBLE PRECISION * On entry, the off-diagonal element associated with the rank-1 * cut which originally split the two submatrices which are now * being recombined. * On exit, RHO has been modified to the value required by * DLAED3. * * CUTPNT (input) INTEGER * The location of the last eigenvalue in the leading * sub-matrix. min(1,N) <= CUTPNT <= N. * * Z (input) DOUBLE PRECISION array, dimension (N) * On entry, Z contains the updating vector (the last row of * the first sub-eigenvector matrix and the first row of the * second sub-eigenvector matrix). * On exit, the contents of Z are destroyed by the updating * process. * * DLAMDA (output) DOUBLE PRECISION array, dimension (N) * A copy of the first K eigenvalues which will be used by * DLAED3 to form the secular equation. * * Q2 (output) DOUBLE PRECISION array, dimension (LDQ2,N) * If ICOMPQ = 0, Q2 is not referenced. Otherwise, * a copy of the first K eigenvectors which will be used by * DLAED7 in a matrix multiply (DGEMM) to update the new * eigenvectors. * * LDQ2 (input) INTEGER * The leading dimension of the array Q2. LDQ2 >= max(1,N). * * W (output) DOUBLE PRECISION array, dimension (N) * The first k values of the final deflation-altered z-vector and * will be passed to DLAED3. * * PERM (output) INTEGER array, dimension (N) * The permutations (from deflation and sorting) to be applied * to each eigenblock. * * GIVPTR (output) INTEGER * The number of Givens rotations which took place in this * subproblem. * * GIVCOL (output) INTEGER array, dimension (2, N) * Each pair of numbers indicates a pair of columns to take place * in a Givens rotation. * * GIVNUM (output) DOUBLE PRECISION array, dimension (2, N) * Each number indicates the S value to be used in the * corresponding Givens rotation. * * INDXP (workspace) INTEGER array, dimension (N) * The permutation used to place deflated values of D at the end * of the array. INDXP(1:K) points to the nondeflated D-values * and INDXP(K+1:N) points to the deflated eigenvalues. * * INDX (workspace) INTEGER array, dimension (N) * The permutation used to sort the contents of D into ascending * order. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, $ TWO = 2.0D0, EIGHT = 8.0D0 ) * .. * .. Local Scalars .. * INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2 DOUBLE PRECISION C, EPS, S, T, TAU, TOL * .. * .. External Functions .. INTEGER IDAMAX DOUBLE PRECISION DLAMCH, DLAPY2 EXTERNAL IDAMAX, DLAMCH, DLAPY2 * .. * .. External Subroutines .. EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN INFO = -4 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN INFO = -7 ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN INFO = -10 ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN INFO = -14 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED8', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * N1 = CUTPNT N2 = N - N1 N1P1 = N1 + 1 * IF( RHO.LT.ZERO ) THEN CALL DSCAL( N2, MONE, Z( N1P1 ), 1 ) END IF * * Normalize z so that norm(z) = 1 * T = ONE / SQRT( TWO ) DO 10 J = 1, N INDX( J ) = J 10 CONTINUE CALL DSCAL( N, T, Z, 1 ) RHO = ABS( TWO*RHO ) * * Sort the eigenvalues into increasing order * DO 20 I = CUTPNT + 1, N INDXQ( I ) = INDXQ( I ) + CUTPNT 20 CONTINUE DO 30 I = 1, N DLAMDA( I ) = D( INDXQ( I ) ) W( I ) = Z( INDXQ( I ) ) 30 CONTINUE I = 1 J = CUTPNT + 1 CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) DO 40 I = 1, N D( I ) = DLAMDA( INDX( I ) ) Z( I ) = W( INDX( I ) ) 40 CONTINUE * * Calculate the allowable deflation tolerence * IMAX = IDAMAX( N, Z, 1 ) JMAX = IDAMAX( N, D, 1 ) EPS = DLAMCH( 'Epsilon' ) TOL = EIGHT*EPS*ABS( D( JMAX ) ) * * If the rank-1 modifier is small enough, no more needs to be done * except to reorganize Q so that its columns correspond with the * elements in D. * IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN K = 0 IF( ICOMPQ.EQ.0 ) THEN DO 50 J = 1, N PERM( J ) = INDXQ( INDX( J ) ) 50 CONTINUE ELSE DO 60 J = 1, N PERM( J ) = INDXQ( INDX( J ) ) CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 60 CONTINUE CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), $ LDQ ) END IF RETURN END IF * * If there are multiple eigenvalues then the problem deflates. Here * the number of equal eigenvalues are found. As each equal * eigenvalue is found, an elementary reflector is computed to rotate * the corresponding eigensubspace so that the corresponding * components of Z are zero in this new basis. * K = 0 GIVPTR = 0 K2 = N + 1 DO 70 J = 1, N IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN * * Deflate due to small z component. * K2 = K2 - 1 INDXP( K2 ) = J IF( J.EQ.N ) $ GO TO 110 ELSE JLAM = J GO TO 80 END IF 70 CONTINUE 80 CONTINUE J = J + 1 IF( J.GT.N ) $ GO TO 100 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN * * Deflate due to small z component. * K2 = K2 - 1 INDXP( K2 ) = J ELSE * * Check if eigenvalues are close enough to allow deflation. * S = Z( JLAM ) C = Z( J ) * * Find sqrt(a**2+b**2) without overflow or * destructive underflow. * TAU = DLAPY2( C, S ) T = D( J ) - D( JLAM ) C = C / TAU S = -S / TAU IF( ABS( T*C*S ).LE.TOL ) THEN * * Deflation is possible. * Z( J ) = TAU Z( JLAM ) = ZERO * * Record the appropriate Givens rotation * GIVPTR = GIVPTR + 1 GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) ) GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) ) GIVNUM( 1, GIVPTR ) = C GIVNUM( 2, GIVPTR ) = S IF( ICOMPQ.EQ.1 ) THEN CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1, $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S ) END IF T = D( JLAM )*C*C + D( J )*S*S D( J ) = D( JLAM )*S*S + D( J )*C*C D( JLAM ) = T K2 = K2 - 1 I = 1 90 CONTINUE IF( K2+I.LE.N ) THEN IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN INDXP( K2+I-1 ) = INDXP( K2+I ) INDXP( K2+I ) = JLAM I = I + 1 GO TO 90 ELSE INDXP( K2+I-1 ) = JLAM END IF ELSE INDXP( K2+I-1 ) = JLAM END IF JLAM = J ELSE K = K + 1 W( K ) = Z( JLAM ) DLAMDA( K ) = D( JLAM ) INDXP( K ) = JLAM JLAM = J END IF END IF GO TO 80 100 CONTINUE * * Record the last eigenvalue. * K = K + 1 W( K ) = Z( JLAM ) DLAMDA( K ) = D( JLAM ) INDXP( K ) = JLAM * 110 CONTINUE * * Sort the eigenvalues and corresponding eigenvectors into DLAMDA * and Q2 respectively. The eigenvalues/vectors which were not * deflated go into the first K slots of DLAMDA and Q2 respectively, * while those which were deflated go into the last N - K slots. * IF( ICOMPQ.EQ.0 ) THEN DO 120 J = 1, N JP = INDXP( J ) DLAMDA( J ) = D( JP ) PERM( J ) = INDXQ( INDX( JP ) ) 120 CONTINUE ELSE DO 130 J = 1, N JP = INDXP( J ) DLAMDA( J ) = D( JP ) PERM( J ) = INDXQ( INDX( JP ) ) CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 130 CONTINUE END IF * * The deflated eigenvalues and their corresponding vectors go back * into the last N - K slots of D and Q respectively. * IF( K.LT.N ) THEN IF( ICOMPQ.EQ.0 ) THEN CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) ELSE CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, $ Q( 1, K+1 ), LDQ ) END IF END IF * RETURN * * End of DLAED8 * END SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, $ S, LDS, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N DOUBLE PRECISION RHO * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), $ W( * ) * .. * * Purpose * ======= * * DLAED9 finds the roots of the secular equation, as defined by the * values in D, Z, and RHO, between KSTART and KSTOP. It makes the * appropriate calls to DLAED4 and then stores the new matrix of * eigenvectors for use in calculating the next level of Z vectors. * * Arguments * ========= * * K (input) INTEGER * The number of terms in the rational function to be solved by * DLAED4. K >= 0. * * KSTART (input) INTEGER * KSTOP (input) INTEGER * The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP * are to be computed. 1 <= KSTART <= KSTOP <= K. * * N (input) INTEGER * The number of rows and columns in the Q matrix. * N >= K (delation may result in N > K). * * D (output) DOUBLE PRECISION array, dimension (N) * D(I) contains the updated eigenvalues * for KSTART <= I <= KSTOP. * * Q (workspace) DOUBLE PRECISION array, dimension (LDQ,N) * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max( 1, N ). * * RHO (input) DOUBLE PRECISION * The value of the parameter in the rank one update equation. * RHO >= 0 required. * * DLAMDA (input) DOUBLE PRECISION array, dimension (K) * The first K elements of this array contain the old roots * of the deflated updating problem. These are the poles * of the secular equation. * * W (input) DOUBLE PRECISION array, dimension (K) * The first K elements of this array contain the components * of the deflation-adjusted updating vector. * * S (output) DOUBLE PRECISION array, dimension (LDS, K) * Will contain the eigenvectors of the repaired matrix which * will be stored for subsequent Z vector calculation and * multiplied by the previously accumulated eigenvectors * to update the system. * * LDS (input) INTEGER * The leading dimension of S. LDS >= max( 1, K ). * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = 1, an eigenvalue did not converge * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Local Scalars .. INTEGER I, J DOUBLE PRECISION TEMP * .. * .. External Functions .. DOUBLE PRECISION DLAMC3, DNRM2 EXTERNAL DLAMC3, DNRM2 * .. * .. External Subroutines .. EXTERNAL DCOPY, DLAED4, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, SIGN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( K.LT.0 ) THEN INFO = -1 ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN INFO = -2 ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) ) $ THEN INFO = -3 ELSE IF( N.LT.K ) THEN INFO = -4 ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN INFO = -7 ELSE IF( LDS.LT.MAX( 1, K ) ) THEN INFO = -12 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAED9', -INFO ) RETURN END IF * * Quick return if possible * IF( K.EQ.0 ) $ RETURN * * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can * be computed with high relative accuracy (barring over/underflow). * This is a problem on machines without a guard digit in * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), * which on any of these machines zeros out the bottommost * bit of DLAMDA(I) if it is 1; this makes the subsequent * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation * occurs. On binary machines with a guard digit (almost all * machines) it does not change DLAMDA(I) at all. On hexadecimal * and decimal machines with a guard digit, it slightly * changes the bottommost bits of DLAMDA(I). It does not account * for hexadecimal or decimal machines without guard digits * (we know of none). We use a subroutine call to compute * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating * this code. * DO 10 I = 1, N DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) 10 CONTINUE * DO 20 J = KSTART, KSTOP CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) * * If the zero finder fails, the computation is terminated. * IF( INFO.NE.0 ) $ GO TO 120 20 CONTINUE * IF( K.EQ.1 .OR. K.EQ.2 ) THEN DO 40 I = 1, K DO 30 J = 1, K S( J, I ) = Q( J, I ) 30 CONTINUE 40 CONTINUE GO TO 120 END IF * * Compute updated W. * CALL DCOPY( K, W, 1, S, 1 ) * * Initialize W(I) = Q(I,I) * CALL DCOPY( K, Q, LDQ+1, W, 1 ) DO 70 J = 1, K DO 50 I = 1, J - 1 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 50 CONTINUE DO 60 I = J + 1, K W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 60 CONTINUE 70 CONTINUE DO 80 I = 1, K W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) ) 80 CONTINUE * * Compute eigenvectors of the modified rank-1 modification. * DO 110 J = 1, K DO 90 I = 1, K Q( I, J ) = W( I ) / Q( I, J ) 90 CONTINUE TEMP = DNRM2( K, Q( 1, J ), 1 ) DO 100 I = 1, K S( I, J ) = Q( I, J ) / TEMP 100 CONTINUE 110 CONTINUE * 120 CONTINUE RETURN * * End of DLAED9 * END SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER CURLVL, CURPBM, INFO, N, TLVLS * .. * .. Array Arguments .. INTEGER GIVCOL( 2, * ), GIVPTR( * ), PERM( * ), $ PRMPTR( * ), QPTR( * ) DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) * .. * * Purpose * ======= * * DLAEDA computes the Z vector corresponding to the merge step in the * CURLVLth step of the merge process with TLVLS steps for the CURPBMth * problem. * * Arguments * ========= * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * TLVLS (input) INTEGER * The total number of merging levels in the overall divide and * conquer tree. * * CURLVL (input) INTEGER * The current level in the overall merge routine, * 0 <= curlvl <= tlvls. * * CURPBM (input) INTEGER * The current problem in the current level in the overall * merge routine (counting from upper left to lower right). * * PRMPTR (input) INTEGER array, dimension (N lg N) * Contains a list of pointers which indicate where in PERM a * level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) * indicates the size of the permutation and incidentally the * size of the full, non-deflated problem. * * PERM (input) INTEGER array, dimension (N lg N) * Contains the permutations (from deflation and sorting) to be * applied to each eigenblock. * * GIVPTR (input) INTEGER array, dimension (N lg N) * Contains a list of pointers which indicate where in GIVCOL a * level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) * indicates the number of Givens rotations. * * GIVCOL (input) INTEGER array, dimension (2, N lg N) * Each pair of numbers indicates a pair of columns to take place * in a Givens rotation. * * GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N) * Each number indicates the S value to be used in the * corresponding Givens rotation. * * Q (input) DOUBLE PRECISION array, dimension (N**2) * Contains the square eigenblocks from previous levels, the * starting positions for blocks are given by QPTR. * * QPTR (input) INTEGER array, dimension (N+2) * Contains a list of pointers which indicate where in Q an * eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates * the size of the block. * * Z (output) DOUBLE PRECISION array, dimension (N) * On output this vector contains the updating vector (the last * row of the first sub-eigenvector matrix and the first row of * the second sub-eigenvector matrix). * * ZTEMP (workspace) DOUBLE PRECISION array, dimension (N) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, HALF, ONE PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2, $ PTR, ZPTR1 * .. * .. External Subroutines .. EXTERNAL DCOPY, DGEMV, DROT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, INT, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( N.LT.0 ) THEN INFO = -1 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLAEDA', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Determine location of first number in second half. * MID = N / 2 + 1 * * Gather last/first rows of appropriate eigenblocks into center of Z * PTR = 1 * * Determine location of lowest level subproblem in the full storage * scheme * CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1 * * Determine size of these matrices. We add HALF to the value of * the SQRT in case the machine underestimates one of these square * roots. * BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) ) DO 10 K = 1, MID - BSIZ1 - 1 Z( K ) = ZERO 10 CONTINUE CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1, $ Z( MID-BSIZ1 ), 1 ) CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 ) DO 20 K = MID + BSIZ2, N Z( K ) = ZERO 20 CONTINUE * * Loop thru remaining levels 1 -> CURLVL applying the Givens * rotations and permutation and then multiplying the center matrices * against the current Z. * PTR = 2**TLVLS + 1 DO 70 K = 1, CURLVL - 1 CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1 PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR ) PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 ) ZPTR1 = MID - PSIZ1 * * Apply Givens at CURR and CURR+1 * DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1 CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1, $ Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ), $ GIVNUM( 2, I ) ) 30 CONTINUE DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1 CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1, $ Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ), $ GIVNUM( 2, I ) ) 40 CONTINUE PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR ) PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 ) DO 50 I = 0, PSIZ1 - 1 ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 ) 50 CONTINUE DO 60 I = 0, PSIZ2 - 1 ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 ) 60 CONTINUE * * Multiply Blocks at CURR and CURR+1 * * Determine size of these matrices. We add HALF to the value of * the SQRT in case the machine underestimates one of these * square roots. * BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+ $ 1 ) ) ) ) IF( BSIZ1.GT.0 ) THEN CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ), $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 ) END IF CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ), $ 1 ) IF( BSIZ2.GT.0 ) THEN CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ), $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 ) END IF CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1, $ Z( MID+BSIZ2 ), 1 ) * PTR = PTR + 2**( TLVLS-K ) 70 CONTINUE * RETURN * * End of DLAEDA * END SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 * .. * * Purpose * ======= * * DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix * [ A B ] * [ B C ]. * On return, RT1 is the eigenvalue of larger absolute value, RT2 is the * eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right * eigenvector for RT1, giving the decomposition * * [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] * [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. * * Arguments * ========= * * A (input) DOUBLE PRECISION * The (1,1) element of the 2-by-2 matrix. * * B (input) DOUBLE PRECISION * The (1,2) element and the conjugate of the (2,1) element of * the 2-by-2 matrix. * * C (input) DOUBLE PRECISION * The (2,2) element of the 2-by-2 matrix. * * RT1 (output) DOUBLE PRECISION * The eigenvalue of larger absolute value. * * RT2 (output) DOUBLE PRECISION * The eigenvalue of smaller absolute value. * * CS1 (output) DOUBLE PRECISION * SN1 (output) DOUBLE PRECISION * The vector (CS1, SN1) is a unit right eigenvector for RT1. * * Further Details * =============== * * RT1 is accurate to a few ulps barring over/underflow. * * RT2 may be inaccurate if there is massive cancellation in the * determinant A*C-B*B; higher precision or correctly rounded or * correctly truncated arithmetic would be needed to compute RT2 * accurately in all cases. * * CS1 and SN1 are accurate to a few ulps barring over/underflow. * * Overflow is possible only if RT1 is within a factor of 5 of overflow. * Underflow is harmless if the input data is 0 or exceeds * underflow_threshold / macheps. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) DOUBLE PRECISION TWO PARAMETER ( TWO = 2.0D0 ) DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) DOUBLE PRECISION HALF PARAMETER ( HALF = 0.5D0 ) * .. * .. Local Scalars .. INTEGER SGN1, SGN2 DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, $ TB, TN * .. * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. * * Compute the eigenvalues * SM = A + C DF = A - C ADF = ABS( DF ) TB = B + B AB = ABS( TB ) IF( ABS( A ).GT.ABS( C ) ) THEN ACMX = A ACMN = C ELSE ACMX = C ACMN = A END IF IF( ADF.GT.AB ) THEN RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) ELSE * * Includes case AB=ADF=0 * RT = AB*SQRT( TWO ) END IF IF( SM.LT.ZERO ) THEN RT1 = HALF*( SM-RT ) SGN1 = -1 * * Order of execution important. * To get fully accurate smaller eigenvalue, * next line needs to be executed in higher precision. * RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B ELSE IF( SM.GT.ZERO ) THEN RT1 = HALF*( SM+RT ) SGN1 = 1 * * Order of execution important. * To get fully accurate smaller eigenvalue, * next line needs to be executed in higher precision. * RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B ELSE * * Includes case RT1 = RT2 = 0 * RT1 = HALF*RT RT2 = -HALF*RT SGN1 = 1 END IF * * Compute the eigenvector * IF( DF.GE.ZERO ) THEN CS = DF + RT SGN2 = 1 ELSE CS = DF - RT SGN2 = -1 END IF ACS = ABS( CS ) IF( ACS.GT.AB ) THEN CT = -TB / CS SN1 = ONE / SQRT( ONE+CT*CT ) CS1 = CT*SN1 ELSE IF( AB.EQ.ZERO ) THEN CS1 = ONE SN1 = ZERO ELSE TN = -CS / TB CS1 = ONE / SQRT( ONE+TN*TN ) SN1 = TN*CS1 END IF END IF IF( SGN1.EQ.SGN2 ) THEN TN = CS1 CS1 = -SN1 SN1 = TN END IF RETURN * * End of DLAEV2 * END SUBROUTINE DLAMRG( N1, N2, A, DTRD1, DTRD2, INDEX ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER DTRD1, DTRD2, N1, N2 * .. * .. Array Arguments .. INTEGER INDEX( * ) DOUBLE PRECISION A( * ) * .. * * Purpose * ======= * * DLAMRG will create a permutation list which will merge the elements * of A (which is composed of two independently sorted sets) into a * single set which is sorted in ascending order. * * Arguments * ========= * * N1 (input) INTEGER * N2 (input) INTEGER * These arguements contain the respective lengths of the two * sorted lists to be merged. * * A (input) DOUBLE PRECISION array, dimension (N1+N2) * The first N1 elements of A contain a list of numbers which * are sorted in either ascending or descending order. Likewise * for the final N2 elements. * * DTRD1 (input) INTEGER * DTRD2 (input) INTEGER * These are the strides to be taken through the array A. * Allowable strides are 1 and -1. They indicate whether a * subset of A is sorted in ascending (DTRDx = 1) or descending * (DTRDx = -1) order. * * INDEX (output) INTEGER array, dimension (N1+N2) * On exit this array will contain a permutation such that * if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be * sorted in ascending order. * * ===================================================================== * * .. Local Scalars .. INTEGER I, IND1, IND2, N1SV, N2SV * .. * .. Executable Statements .. * N1SV = N1 N2SV = N2 IF( DTRD1.GT.0 ) THEN IND1 = 1 ELSE IND1 = N1 END IF IF( DTRD2.GT.0 ) THEN IND2 = 1 + N1 ELSE IND2 = N1 + N2 END IF I = 1 * while ( (N1SV > 0) & (N2SV > 0) ) 10 CONTINUE IF( N1SV.GT.0 .AND. N2SV.GT.0 ) THEN IF( A( IND1 ).LE.A( IND2 ) ) THEN INDEX( I ) = IND1 I = I + 1 IND1 = IND1 + DTRD1 N1SV = N1SV - 1 ELSE INDEX( I ) = IND2 I = I + 1 IND2 = IND2 + DTRD2 N2SV = N2SV - 1 END IF GO TO 10 END IF * end while IF( N1SV.EQ.0 ) THEN DO 20 N1SV = 1, N2SV INDEX( I ) = IND2 I = I + 1 IND2 = IND2 + DTRD2 20 CONTINUE ELSE * N2SV .EQ. 0 DO 30 N2SV = 1, N1SV INDEX( I ) = IND1 I = I + 1 IND1 = IND1 + DTRD1 30 CONTINUE END IF * RETURN * * End of DLAMRG * END DOUBLE PRECISION FUNCTION DLANST( NORM, N, D, E ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER NORM INTEGER N * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), E( * ) * .. * * Purpose * ======= * * DLANST returns the value of the one norm, or the Frobenius norm, or * the infinity norm, or the element of largest absolute value of a * real symmetric tridiagonal matrix A. * * Description * =========== * * DLANST returns the value * * DLANST = ( max(abs(A(i,j))), NORM = 'M' or 'm' * ( * ( norm1(A), NORM = '1', 'O' or 'o' * ( * ( normI(A), NORM = 'I' or 'i' * ( * ( normF(A), NORM = 'F', 'f', 'E' or 'e' * * where norm1 denotes the one norm of a matrix (maximum column sum), * normI denotes the infinity norm of a matrix (maximum row sum) and * normF denotes the Frobenius norm of a matrix (square root of sum of * squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. * * Arguments * ========= * * NORM (input) CHARACTER*1 * Specifies the value to be returned in DLANST as described * above. * * N (input) INTEGER * The order of the matrix A. N >= 0. When N = 0, DLANST is * set to zero. * * D (input) DOUBLE PRECISION array, dimension (N) * The diagonal elements of A. * * E (input) DOUBLE PRECISION array, dimension (N-1) * The (n-1) sub-diagonal or super-diagonal elements of A. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER I DOUBLE PRECISION ANORM, SCALE, SUM * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLASSQ * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. * .. Executable Statements .. * IF( N.LE.0 ) THEN ANORM = ZERO ELSE IF( LSAME( NORM, 'M' ) ) THEN * * Find max(abs(A(i,j))). * ANORM = ABS( D( N ) ) DO 10 I = 1, N - 1 ANORM = MAX( ANORM, ABS( D( I ) ) ) ANORM = MAX( ANORM, ABS( E( I ) ) ) 10 CONTINUE ELSE IF( LSAME( NORM, 'O' ) .OR. NORM.EQ.'1' .OR. $ LSAME( NORM, 'I' ) ) THEN * * Find norm1(A). * IF( N.EQ.1 ) THEN ANORM = ABS( D( 1 ) ) ELSE ANORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), $ ABS( E( N-1 ) )+ABS( D( N ) ) ) DO 20 I = 2, N - 1 ANORM = MAX( ANORM, ABS( D( I ) )+ABS( E( I ) )+ $ ABS( E( I-1 ) ) ) 20 CONTINUE END IF ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN * * Find normF(A). * SCALE = ZERO SUM = ONE IF( N.GT.1 ) THEN CALL DLASSQ( N-1, E, 1, SCALE, SUM ) SUM = 2*SUM END IF CALL DLASSQ( N, D, 1, SCALE, SUM ) ANORM = SCALE*SQRT( SUM ) END IF * DLANST = ANORM RETURN * * End of DLANST * END DOUBLE PRECISION FUNCTION DLAPY2( X, Y ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. DOUBLE PRECISION X, Y * .. * * Purpose * ======= * * DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary * overflow. * * Arguments * ========= * * X (input) DOUBLE PRECISION * Y (input) DOUBLE PRECISION * X and Y specify the values x and y. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION W, XABS, YABS, Z * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * XABS = ABS( X ) YABS = ABS( Y ) W = MAX( XABS, YABS ) Z = MIN( XABS, YABS ) IF( Z.EQ.ZERO ) THEN DLAPY2 = W ELSE DLAPY2 = W*SQRT( ONE+( Z / W )**2 ) END IF RETURN * * End of DLAPY2 * END SUBROUTINE DLARTG( F, G, CS, SN, R ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. DOUBLE PRECISION CS, F, G, R, SN * .. * * Purpose * ======= * * DLARTG generate a plane rotation so that * * [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. * [ -SN CS ] [ G ] [ 0 ] * * This is a slower, more accurate version of the BLAS1 routine DROTG, * with the following other differences: * F and G are unchanged on return. * If G=0, then CS=1 and SN=0. * If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any * floating point operations (saves work in DBDSQR when * there are zeros on the diagonal). * * If F exceeds G in magnitude, CS will be positive. * * Arguments * ========= * * F (input) DOUBLE PRECISION * The first component of vector to be rotated. * * G (input) DOUBLE PRECISION * The second component of vector to be rotated. * * CS (output) DOUBLE PRECISION * The cosine of the rotation. * * SN (output) DOUBLE PRECISION * The sine of the rotation. * * R (output) DOUBLE PRECISION * The nonzero component of the rotated vector. * * This version has a few statements commented out for thread safety * (machine parameters are computed on each entry). 10 feb 03, SJH. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) DOUBLE PRECISION TWO PARAMETER ( TWO = 2.0D0 ) * .. * .. Local Scalars .. * LOGICAL FIRST INTEGER COUNT, I DOUBLE PRECISION EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, SQRT * .. * .. Save statement .. * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 * .. * .. Data statements .. * DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * * IF( FIRST ) THEN SAFMIN = DLAMCH( 'S' ) EPS = DLAMCH( 'E' ) SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / $ LOG( DLAMCH( 'B' ) ) / TWO ) SAFMX2 = ONE / SAFMN2 * FIRST = .FALSE. * END IF IF( G.EQ.ZERO ) THEN CS = ONE SN = ZERO R = F ELSE IF( F.EQ.ZERO ) THEN CS = ZERO SN = ONE R = G ELSE F1 = F G1 = G SCALE = MAX( ABS( F1 ), ABS( G1 ) ) IF( SCALE.GE.SAFMX2 ) THEN COUNT = 0 10 CONTINUE COUNT = COUNT + 1 F1 = F1*SAFMN2 G1 = G1*SAFMN2 SCALE = MAX( ABS( F1 ), ABS( G1 ) ) IF( SCALE.GE.SAFMX2 ) $ GO TO 10 R = SQRT( F1**2+G1**2 ) CS = F1 / R SN = G1 / R DO 20 I = 1, COUNT R = R*SAFMX2 20 CONTINUE ELSE IF( SCALE.LE.SAFMN2 ) THEN COUNT = 0 30 CONTINUE COUNT = COUNT + 1 F1 = F1*SAFMX2 G1 = G1*SAFMX2 SCALE = MAX( ABS( F1 ), ABS( G1 ) ) IF( SCALE.LE.SAFMN2 ) $ GO TO 30 R = SQRT( F1**2+G1**2 ) CS = F1 / R SN = G1 / R DO 40 I = 1, COUNT R = R*SAFMN2 40 CONTINUE ELSE R = SQRT( F1**2+G1**2 ) CS = F1 / R SN = G1 / R END IF IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN CS = -CS SN = -SN R = -R END IF END IF RETURN * * End of DLARTG * END SUBROUTINE DLASCL( TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER TYPE INTEGER INFO, KL, KU, LDA, M, N DOUBLE PRECISION CFROM, CTO * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DLASCL multiplies the M by N real matrix A by the real scalar * CTO/CFROM. This is done without over/underflow as long as the final * result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that * A may be full, upper triangular, lower triangular, upper Hessenberg, * or banded. * * Arguments * ========= * * TYPE (input) CHARACTER*1 * TYPE indices the storage type of the input matrix. * = 'G': A is a full matrix. * = 'L': A is a lower triangular matrix. * = 'U': A is an upper triangular matrix. * = 'H': A is an upper Hessenberg matrix. * = 'B': A is a symmetric band matrix with lower bandwidth KL * and upper bandwidth KU and with the only the lower * half stored. * = 'Q': A is a symmetric band matrix with lower bandwidth KL * and upper bandwidth KU and with the only the upper * half stored. * = 'Z': A is a band matrix with lower bandwidth KL and upper * bandwidth KU. * * KL (input) INTEGER * The lower bandwidth of A. Referenced only if TYPE = 'B', * 'Q' or 'Z'. * * KU (input) INTEGER * The upper bandwidth of A. Referenced only if TYPE = 'B', * 'Q' or 'Z'. * * CFROM (input) DOUBLE PRECISION * CTO (input) DOUBLE PRECISION * The matrix A is multiplied by CTO/CFROM. A(I,J) is computed * without over/underflow if the final result CTO*A(I,J)/CFROM * can be represented without over/underflow. CFROM must be * nonzero. * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * The matrix to be multiplied by CTO/CFROM. See TYPE for the * storage type. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * INFO (output) INTEGER * 0 - successful exit * <0 - if INFO = -i, the i-th argument had an illegal value. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. LOGICAL DONE INTEGER I, ITYPE, J, K1, K2, K3, K4 DOUBLE PRECISION BIGNUM, CFROM1, CFROMC, CTO1, CTOC, MUL, SMLNUM * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH EXTERNAL LSAME, DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 * IF( LSAME( TYPE, 'G' ) ) THEN ITYPE = 0 ELSE IF( LSAME( TYPE, 'L' ) ) THEN ITYPE = 1 ELSE IF( LSAME( TYPE, 'U' ) ) THEN ITYPE = 2 ELSE IF( LSAME( TYPE, 'H' ) ) THEN ITYPE = 3 ELSE IF( LSAME( TYPE, 'B' ) ) THEN ITYPE = 4 ELSE IF( LSAME( TYPE, 'Q' ) ) THEN ITYPE = 5 ELSE IF( LSAME( TYPE, 'Z' ) ) THEN ITYPE = 6 ELSE ITYPE = -1 END IF * IF( ITYPE.EQ.-1 ) THEN INFO = -1 ELSE IF( CFROM.EQ.ZERO ) THEN INFO = -4 ELSE IF( M.LT.0 ) THEN INFO = -6 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.4 .AND. N.NE.M ) .OR. $ ( ITYPE.EQ.5 .AND. N.NE.M ) ) THEN INFO = -7 ELSE IF( ITYPE.LE.3 .AND. LDA.LT.MAX( 1, M ) ) THEN INFO = -9 ELSE IF( ITYPE.GE.4 ) THEN IF( KL.LT.0 .OR. KL.GT.MAX( M-1, 0 ) ) THEN INFO = -2 ELSE IF( KU.LT.0 .OR. KU.GT.MAX( N-1, 0 ) .OR. $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. KL.NE.KU ) ) $ THEN INFO = -3 ELSE IF( ( ITYPE.EQ.4 .AND. LDA.LT.KL+1 ) .OR. $ ( ITYPE.EQ.5 .AND. LDA.LT.KU+1 ) .OR. $ ( ITYPE.EQ.6 .AND. LDA.LT.2*KL+KU+1 ) ) THEN INFO = -9 END IF END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLASCL', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 .OR. M.EQ.0 ) $ RETURN * * Get machine parameters * SMLNUM = DLAMCH( 'S' ) BIGNUM = ONE / SMLNUM * CFROMC = CFROM CTOC = CTO * 10 CONTINUE CFROM1 = CFROMC*SMLNUM CTO1 = CTOC / BIGNUM IF( ABS( CFROM1 ).GT.ABS( CTOC ) .AND. CTOC.NE.ZERO ) THEN MUL = SMLNUM DONE = .FALSE. CFROMC = CFROM1 ELSE IF( ABS( CTO1 ).GT.ABS( CFROMC ) ) THEN MUL = BIGNUM DONE = .FALSE. CTOC = CTO1 ELSE MUL = CTOC / CFROMC DONE = .TRUE. END IF * IF( ITYPE.EQ.0 ) THEN * * Full matrix * DO 30 J = 1, N DO 20 I = 1, M A( I, J ) = A( I, J )*MUL 20 CONTINUE 30 CONTINUE * ELSE IF( ITYPE.EQ.1 ) THEN * * Lower triangular matrix * DO 50 J = 1, N DO 40 I = J, M A( I, J ) = A( I, J )*MUL 40 CONTINUE 50 CONTINUE * ELSE IF( ITYPE.EQ.2 ) THEN * * Upper triangular matrix * DO 70 J = 1, N DO 60 I = 1, MIN( J, M ) A( I, J ) = A( I, J )*MUL 60 CONTINUE 70 CONTINUE * ELSE IF( ITYPE.EQ.3 ) THEN * * Upper Hessenberg matrix * DO 90 J = 1, N DO 80 I = 1, MIN( J+1, M ) A( I, J ) = A( I, J )*MUL 80 CONTINUE 90 CONTINUE * ELSE IF( ITYPE.EQ.4 ) THEN * * Lower half of a symmetric band matrix * K3 = KL + 1 K4 = N + 1 DO 110 J = 1, N DO 100 I = 1, MIN( K3, K4-J ) A( I, J ) = A( I, J )*MUL 100 CONTINUE 110 CONTINUE * ELSE IF( ITYPE.EQ.5 ) THEN * * Upper half of a symmetric band matrix * K1 = KU + 2 K3 = KU + 1 DO 130 J = 1, N DO 120 I = MAX( K1-J, 1 ), K3 A( I, J ) = A( I, J )*MUL 120 CONTINUE 130 CONTINUE * ELSE IF( ITYPE.EQ.6 ) THEN * * Band matrix * K1 = KL + KU + 2 K2 = KL + 1 K3 = 2*KL + KU + 1 K4 = KL + KU + 1 + M DO 150 J = 1, N DO 140 I = MAX( K1-J, K2 ), MIN( K3, K4-J ) A( I, J ) = A( I, J )*MUL 140 CONTINUE 150 CONTINUE * END IF * IF( .NOT.DONE ) $ GO TO 10 * RETURN * * End of DLASCL * END SUBROUTINE DLASET( UPLO, M, N, ALPHA, BETA, A, LDA ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER LDA, M, N DOUBLE PRECISION ALPHA, BETA * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DLASET initializes an m-by-n matrix A to BETA on the diagonal and * ALPHA on the offdiagonals. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies the part of the matrix A to be set. * = 'U': Upper triangular part is set; the strictly lower * triangular part of A is not changed. * = 'L': Lower triangular part is set; the strictly upper * triangular part of A is not changed. * Otherwise: All of the matrix A is set. * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * ALPHA (input) DOUBLE PRECISION * The constant to which the offdiagonal elements are to be set. * * BETA (input) DOUBLE PRECISION * The constant to which the diagonal elements are to be set. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On exit, the leading m-by-n submatrix of A is set as follows: * * if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n, * if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n, * otherwise, A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j, * * and, for all UPLO, A(i,i) = BETA, 1<=i<=min(m,n). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * ===================================================================== * * .. Local Scalars .. INTEGER I, J * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * IF( LSAME( UPLO, 'U' ) ) THEN * * Set the strictly upper triangular or trapezoidal part of the * array to ALPHA. * DO 20 J = 2, N DO 10 I = 1, MIN( J-1, M ) A( I, J ) = ALPHA 10 CONTINUE 20 CONTINUE * ELSE IF( LSAME( UPLO, 'L' ) ) THEN * * Set the strictly lower triangular or trapezoidal part of the * array to ALPHA. * DO 40 J = 1, MIN( M, N ) DO 30 I = J + 1, M A( I, J ) = ALPHA 30 CONTINUE 40 CONTINUE * ELSE * * Set the leading m-by-n submatrix to ALPHA. * DO 60 J = 1, N DO 50 I = 1, M A( I, J ) = ALPHA 50 CONTINUE 60 CONTINUE END IF * * Set the first min(M,N) diagonal elements to BETA. * DO 70 I = 1, MIN( M, N ) A( I, I ) = BETA 70 CONTINUE * RETURN * * End of DLASET * END SUBROUTINE DLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER DIRECT, PIVOT, SIDE INTEGER LDA, M, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), C( * ), S( * ) * .. * * Purpose * ======= * * DLASR applies a sequence of plane rotations to a real matrix A, * from either the left or the right. * * When SIDE = 'L', the transformation takes the form * * A := P*A * * and when SIDE = 'R', the transformation takes the form * * A := A*P**T * * where P is an orthogonal matrix consisting of a sequence of z plane * rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', * and P**T is the transpose of P. * * When DIRECT = 'F' (Forward sequence), then * * P = P(z-1) * ... * P(2) * P(1) * * and when DIRECT = 'B' (Backward sequence), then * * P = P(1) * P(2) * ... * P(z-1) * * where P(k) is a plane rotation matrix defined by the 2-by-2 rotation * * R(k) = ( c(k) s(k) ) * = ( -s(k) c(k) ). * * When PIVOT = 'V' (Variable pivot), the rotation is performed * for the plane (k,k+1), i.e., P(k) has the form * * P(k) = ( 1 ) * ( ... ) * ( 1 ) * ( c(k) s(k) ) * ( -s(k) c(k) ) * ( 1 ) * ( ... ) * ( 1 ) * * where R(k) appears as a rank-2 modification to the identity matrix in * rows and columns k and k+1. * * When PIVOT = 'T' (Top pivot), the rotation is performed for the * plane (1,k+1), so P(k) has the form * * P(k) = ( c(k) s(k) ) * ( 1 ) * ( ... ) * ( 1 ) * ( -s(k) c(k) ) * ( 1 ) * ( ... ) * ( 1 ) * * where R(k) appears in rows and columns 1 and k+1. * * Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is * performed for the plane (k,z), giving P(k) the form * * P(k) = ( 1 ) * ( ... ) * ( 1 ) * ( c(k) s(k) ) * ( 1 ) * ( ... ) * ( 1 ) * ( -s(k) c(k) ) * * where R(k) appears in rows and columns k and z. The rotations are * performed without ever forming P(k) explicitly. * * Arguments * ========= * * SIDE (input) CHARACTER*1 * Specifies whether the plane rotation matrix P is applied to * A on the left or the right. * = 'L': Left, compute A := P*A * = 'R': Right, compute A:= A*P**T * * PIVOT (input) CHARACTER*1 * Specifies the plane for which P(k) is a plane rotation * matrix. * = 'V': Variable pivot, the plane (k,k+1) * = 'T': Top pivot, the plane (1,k+1) * = 'B': Bottom pivot, the plane (k,z) * * DIRECT (input) CHARACTER*1 * Specifies whether P is a forward or backward sequence of * plane rotations. * = 'F': Forward, P = P(z-1)*...*P(2)*P(1) * = 'B': Backward, P = P(1)*P(2)*...*P(z-1) * * M (input) INTEGER * The number of rows of the matrix A. If m <= 1, an immediate * return is effected. * * N (input) INTEGER * The number of columns of the matrix A. If n <= 1, an * immediate return is effected. * * C (input) DOUBLE PRECISION array, dimension * (M-1) if SIDE = 'L' * (N-1) if SIDE = 'R' * The cosines c(k) of the plane rotations. * * S (input) DOUBLE PRECISION array, dimension * (M-1) if SIDE = 'L' * (N-1) if SIDE = 'R' * The sines s(k) of the plane rotations. The 2-by-2 plane * rotation part of the matrix P(k), R(k), has the form * R(k) = ( c(k) s(k) ) * ( -s(k) c(k) ). * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * The M-by-N matrix A. On exit, A is overwritten by P*A if * SIDE = 'R' or by A*P**T if SIDE = 'L'. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER I, INFO, J DOUBLE PRECISION CTEMP, STEMP, TEMP * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters * INFO = 0 IF( .NOT.( LSAME( SIDE, 'L' ) .OR. LSAME( SIDE, 'R' ) ) ) THEN INFO = 1 ELSE IF( .NOT.( LSAME( PIVOT, 'V' ) .OR. LSAME( PIVOT, $ 'T' ) .OR. LSAME( PIVOT, 'B' ) ) ) THEN INFO = 2 ELSE IF( .NOT.( LSAME( DIRECT, 'F' ) .OR. LSAME( DIRECT, 'B' ) ) ) $ THEN INFO = 3 ELSE IF( M.LT.0 ) THEN INFO = 4 ELSE IF( N.LT.0 ) THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = 9 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLASR ', INFO ) RETURN END IF * * Quick return if possible * IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) $ RETURN IF( LSAME( SIDE, 'L' ) ) THEN * * Form P * A * IF( LSAME( PIVOT, 'V' ) ) THEN IF( LSAME( DIRECT, 'F' ) ) THEN DO 20 J = 1, M - 1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 10 I = 1, N TEMP = A( J+1, I ) A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I ) A( J, I ) = STEMP*TEMP + CTEMP*A( J, I ) 10 CONTINUE END IF 20 CONTINUE ELSE IF( LSAME( DIRECT, 'B' ) ) THEN DO 40 J = M - 1, 1, -1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 30 I = 1, N TEMP = A( J+1, I ) A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I ) A( J, I ) = STEMP*TEMP + CTEMP*A( J, I ) 30 CONTINUE END IF 40 CONTINUE END IF ELSE IF( LSAME( PIVOT, 'T' ) ) THEN IF( LSAME( DIRECT, 'F' ) ) THEN DO 60 J = 2, M CTEMP = C( J-1 ) STEMP = S( J-1 ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 50 I = 1, N TEMP = A( J, I ) A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I ) A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I ) 50 CONTINUE END IF 60 CONTINUE ELSE IF( LSAME( DIRECT, 'B' ) ) THEN DO 80 J = M, 2, -1 CTEMP = C( J-1 ) STEMP = S( J-1 ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 70 I = 1, N TEMP = A( J, I ) A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I ) A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I ) 70 CONTINUE END IF 80 CONTINUE END IF ELSE IF( LSAME( PIVOT, 'B' ) ) THEN IF( LSAME( DIRECT, 'F' ) ) THEN DO 100 J = 1, M - 1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 90 I = 1, N TEMP = A( J, I ) A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP 90 CONTINUE END IF 100 CONTINUE ELSE IF( LSAME( DIRECT, 'B' ) ) THEN DO 120 J = M - 1, 1, -1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 110 I = 1, N TEMP = A( J, I ) A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP 110 CONTINUE END IF 120 CONTINUE END IF END IF ELSE IF( LSAME( SIDE, 'R' ) ) THEN * * Form A * P' * IF( LSAME( PIVOT, 'V' ) ) THEN IF( LSAME( DIRECT, 'F' ) ) THEN DO 140 J = 1, N - 1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 130 I = 1, M TEMP = A( I, J+1 ) A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J ) A( I, J ) = STEMP*TEMP + CTEMP*A( I, J ) 130 CONTINUE END IF 140 CONTINUE ELSE IF( LSAME( DIRECT, 'B' ) ) THEN DO 160 J = N - 1, 1, -1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 150 I = 1, M TEMP = A( I, J+1 ) A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J ) A( I, J ) = STEMP*TEMP + CTEMP*A( I, J ) 150 CONTINUE END IF 160 CONTINUE END IF ELSE IF( LSAME( PIVOT, 'T' ) ) THEN IF( LSAME( DIRECT, 'F' ) ) THEN DO 180 J = 2, N CTEMP = C( J-1 ) STEMP = S( J-1 ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 170 I = 1, M TEMP = A( I, J ) A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 ) A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 ) 170 CONTINUE END IF 180 CONTINUE ELSE IF( LSAME( DIRECT, 'B' ) ) THEN DO 200 J = N, 2, -1 CTEMP = C( J-1 ) STEMP = S( J-1 ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 190 I = 1, M TEMP = A( I, J ) A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 ) A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 ) 190 CONTINUE END IF 200 CONTINUE END IF ELSE IF( LSAME( PIVOT, 'B' ) ) THEN IF( LSAME( DIRECT, 'F' ) ) THEN DO 220 J = 1, N - 1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 210 I = 1, M TEMP = A( I, J ) A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP 210 CONTINUE END IF 220 CONTINUE ELSE IF( LSAME( DIRECT, 'B' ) ) THEN DO 240 J = N - 1, 1, -1 CTEMP = C( J ) STEMP = S( J ) IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN DO 230 I = 1, M TEMP = A( I, J ) A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP 230 CONTINUE END IF 240 CONTINUE END IF END IF END IF * RETURN * * End of DLASR * END SUBROUTINE DLASRT( ID, N, D, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER ID INTEGER INFO, N * .. * .. Array Arguments .. DOUBLE PRECISION D( * ) * .. * * Purpose * ======= * * Sort the numbers in D in increasing order (if ID = 'I') or * in decreasing order (if ID = 'D' ). * * Use Quick Sort, reverting to Insertion sort on arrays of * size <= 20. Dimension of STACK limits N to about 2**32. * * Arguments * ========= * * ID (input) CHARACTER*1 * = 'I': sort D in increasing order; * = 'D': sort D in decreasing order. * * N (input) INTEGER * The length of the array D. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the array to be sorted. * On exit, D has been sorted into increasing order * (D(1) <= ... <= D(N) ) or into decreasing order * (D(1) >= ... >= D(N) ), depending on ID. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== * * .. Parameters .. INTEGER SELECT PARAMETER ( SELECT = 20 ) * .. * .. Local Scalars .. INTEGER DIR, ENDD, I, J, START, STKPNT DOUBLE PRECISION D1, D2, D3, DMNMX, TMP * .. * .. Local Arrays .. INTEGER STACK( 2, 32 ) * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * * Test the input paramters. * INFO = 0 DIR = -1 IF( LSAME( ID, 'D' ) ) THEN DIR = 0 ELSE IF( LSAME( ID, 'I' ) ) THEN DIR = 1 END IF IF( DIR.EQ.-1 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DLASRT', -INFO ) RETURN END IF * * Quick return if possible * IF( N.LE.1 ) $ RETURN * STKPNT = 1 STACK( 1, 1 ) = 1 STACK( 2, 1 ) = N 10 CONTINUE START = STACK( 1, STKPNT ) ENDD = STACK( 2, STKPNT ) STKPNT = STKPNT - 1 IF( ENDD-START.LE.SELECT .AND. ENDD-START.GT.0 ) THEN * * Do Insertion sort on D( START:ENDD ) * IF( DIR.EQ.0 ) THEN * * Sort into decreasing order * DO 30 I = START + 1, ENDD DO 20 J = I, START + 1, -1 IF( D( J ).GT.D( J-1 ) ) THEN DMNMX = D( J ) D( J ) = D( J-1 ) D( J-1 ) = DMNMX ELSE GO TO 30 END IF 20 CONTINUE 30 CONTINUE * ELSE * * Sort into increasing order * DO 50 I = START + 1, ENDD DO 40 J = I, START + 1, -1 IF( D( J ).LT.D( J-1 ) ) THEN DMNMX = D( J ) D( J ) = D( J-1 ) D( J-1 ) = DMNMX ELSE GO TO 50 END IF 40 CONTINUE 50 CONTINUE * END IF * ELSE IF( ENDD-START.GT.SELECT ) THEN * * Partition D( START:ENDD ) and stack parts, largest one first * * Choose partition entry as median of 3 * D1 = D( START ) D2 = D( ENDD ) I = ( START+ENDD ) / 2 D3 = D( I ) IF( D1.LT.D2 ) THEN IF( D3.LT.D1 ) THEN DMNMX = D1 ELSE IF( D3.LT.D2 ) THEN DMNMX = D3 ELSE DMNMX = D2 END IF ELSE IF( D3.LT.D2 ) THEN DMNMX = D2 ELSE IF( D3.LT.D1 ) THEN DMNMX = D3 ELSE DMNMX = D1 END IF END IF * IF( DIR.EQ.0 ) THEN * * Sort into decreasing order * I = START - 1 J = ENDD + 1 60 CONTINUE 70 CONTINUE J = J - 1 IF( D( J ).LT.DMNMX ) $ GO TO 70 80 CONTINUE I = I + 1 IF( D( I ).GT.DMNMX ) $ GO TO 80 IF( I.LT.J ) THEN TMP = D( I ) D( I ) = D( J ) D( J ) = TMP GO TO 60 END IF IF( J-START.GT.ENDD-J-1 ) THEN STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = START STACK( 2, STKPNT ) = J STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = J + 1 STACK( 2, STKPNT ) = ENDD ELSE STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = J + 1 STACK( 2, STKPNT ) = ENDD STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = START STACK( 2, STKPNT ) = J END IF ELSE * * Sort into increasing order * I = START - 1 J = ENDD + 1 90 CONTINUE 100 CONTINUE J = J - 1 IF( D( J ).GT.DMNMX ) $ GO TO 100 110 CONTINUE I = I + 1 IF( D( I ).LT.DMNMX ) $ GO TO 110 IF( I.LT.J ) THEN TMP = D( I ) D( I ) = D( J ) D( J ) = TMP GO TO 90 END IF IF( J-START.GT.ENDD-J-1 ) THEN STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = START STACK( 2, STKPNT ) = J STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = J + 1 STACK( 2, STKPNT ) = ENDD ELSE STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = J + 1 STACK( 2, STKPNT ) = ENDD STKPNT = STKPNT + 1 STACK( 1, STKPNT ) = START STACK( 2, STKPNT ) = J END IF END IF END IF IF( STKPNT.GT.0 ) $ GO TO 10 RETURN * * End of DLASRT * END SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INCX, N DOUBLE PRECISION SCALE, SUMSQ * .. * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * Purpose * ======= * * DLASSQ returns the values scl and smsq such that * * ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, * * where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is * assumed to be non-negative and scl returns the value * * scl = max( scale, abs( x( i ) ) ). * * scale and sumsq must be supplied in SCALE and SUMSQ and * scl and smsq are overwritten on SCALE and SUMSQ respectively. * * The routine makes only one pass through the vector x. * * Arguments * ========= * * N (input) INTEGER * The number of elements to be used from the vector X. * * X (input) DOUBLE PRECISION array, dimension (N) * The vector for which a scaled sum of squares is computed. * x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. * * INCX (input) INTEGER * The increment between successive values of the vector X. * INCX > 0. * * SCALE (input/output) DOUBLE PRECISION * On entry, the value scale in the equation above. * On exit, SCALE is overwritten with scl , the scaling factor * for the sum of squares. * * SUMSQ (input/output) DOUBLE PRECISION * On entry, the value sumsq in the equation above. * On exit, SUMSQ is overwritten with smsq , the basic sum of * squares from which scl has been factored out. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER IX DOUBLE PRECISION ABSXI * .. * .. Intrinsic Functions .. INTRINSIC ABS * .. * .. Executable Statements .. * IF( N.GT.0 ) THEN DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX IF( X( IX ).NE.ZERO ) THEN ABSXI = ABS( X( IX ) ) IF( SCALE.LT.ABSXI ) THEN SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 SCALE = ABSXI ELSE SUMSQ = SUMSQ + ( ABSXI / SCALE )**2 END IF END IF 10 CONTINUE END IF RETURN * * End of DLASSQ * END SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, $ LIWORK, INFO ) * * -- LAPACK driver routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER COMPZ INTEGER INFO, LDZ, LIWORK, LWORK, N * .. * .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) * .. * * Purpose * ======= * * DSTEDC computes all eigenvalues and, optionally, eigenvectors of a * symmetric tridiagonal matrix using the divide and conquer method. * The eigenvectors of a full or band real symmetric matrix can also be * found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this * matrix to tridiagonal form. * * This code makes very mild assumptions about floating point * arithmetic. It will work on machines with a guard digit in * add/subtract, or on those binary machines without guard digits * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. * It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. See DLAED3 for details. * * Arguments * ========= * * COMPZ (input) CHARACTER*1 * = 'N': Compute eigenvalues only. * = 'I': Compute eigenvectors of tridiagonal matrix also. * = 'V': Compute eigenvectors of original dense symmetric * matrix also. On entry, Z contains the orthogonal * matrix used to reduce the original matrix to * tridiagonal form. * * N (input) INTEGER * The dimension of the symmetric tridiagonal matrix. N >= 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the diagonal elements of the tridiagonal matrix. * On exit, if INFO = 0, the eigenvalues in ascending order. * * E (input/output) DOUBLE PRECISION array, dimension (N-1) * On entry, the subdiagonal elements of the tridiagonal matrix. * On exit, E has been destroyed. * * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) * On entry, if COMPZ = 'V', then Z contains the orthogonal * matrix used in the reduction to tridiagonal form. * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the * orthonormal eigenvectors of the original symmetric matrix, * and if COMPZ = 'I', Z contains the orthonormal eigenvectors * of the symmetric tridiagonal matrix. * If COMPZ = 'N', then Z is not referenced. * * LDZ (input) INTEGER * The leading dimension of the array Z. LDZ >= 1. * If eigenvectors are desired, then LDZ >= max(1,N). * * WORK (workspace/output) DOUBLE PRECISION array, * dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. * If COMPZ = 'N' or N <= 1 then LWORK must be at least 1. * If COMPZ = 'V' and N > 1 then LWORK must be at least * ( 1 + 3*N + 2*N*lg N + 3*N**2 ), * where lg( N ) = smallest integer k such * that 2**k >= N. * If COMPZ = 'I' and N > 1 then LWORK must be at least * ( 1 + 4*N + N**2 ). * Note that for COMPZ = 'I' or 'V', then if N is less than or * equal to the minimum divide size, usually 25, then LWORK need * only be max(1,2*(N-1)). * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. * * LIWORK (input) INTEGER * The dimension of the array IWORK. * If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1. * If COMPZ = 'V' and N > 1 then LIWORK must be at least * ( 6 + 6*N + 5*N*lg N ). * If COMPZ = 'I' and N > 1 then LIWORK must be at least * ( 3 + 5*N ). * Note that for COMPZ = 'I' or 'V', then if N is less than or * equal to the minimum divide size, usually 25, then LIWORK * need only be 1. * * If LIWORK = -1, then a workspace query is assumed; the * routine only calculates the optimal size of the IWORK array, * returns this value as the first entry of the IWORK array, and * no error message related to LIWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: The algorithm failed to compute an eigenvalue while * working on the submatrix lying in rows and columns * INFO/(N+1) through mod(INFO,N+1). * * Further Details * =============== * * Based on contributions by * Jeff Rutter, Computer Science Division, University of California * at Berkeley, USA * Modified by Francoise Tisseur, University of Tennessee. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW DOUBLE PRECISION EPS, ORGNRM, P, TINY * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLANST EXTERNAL LSAME, ILAENV, DLAMCH, DLANST * .. * .. External Subroutines .. EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT, $ DSTEQR, DSTERF, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) * IF( LSAME( COMPZ, 'N' ) ) THEN ICOMPZ = 0 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN ICOMPZ = 1 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN ICOMPZ = 2 ELSE ICOMPZ = -1 END IF IF( ICOMPZ.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( ( LDZ.LT.1 ) .OR. $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN INFO = -6 END IF * IF( INFO.EQ.0 ) THEN * * Compute the workspace requirements * SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 ) IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN LIWMIN = 1 LWMIN = 1 ELSE IF( N.LE.SMLSIZ ) THEN LIWMIN = 1 LWMIN = 2*( N - 1 ) ELSE LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) ) IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IF( 2**LGN.LT.N ) $ LGN = LGN + 1 IF( ICOMPZ.EQ.1 ) THEN LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2 LIWMIN = 6 + 6*N + 5*N*LGN ELSE IF( ICOMPZ.EQ.2 ) THEN LWMIN = 1 + 4*N + N**2 LIWMIN = 3 + 5*N END IF END IF WORK( 1 ) = LWMIN IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN INFO = -8 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN INFO = -10 END IF END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSTEDC', -INFO ) RETURN ELSE IF (LQUERY) THEN RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN IF( N.EQ.1 ) THEN IF( ICOMPZ.NE.0 ) $ Z( 1, 1 ) = ONE RETURN END IF * * If the following conditional clause is removed, then the routine * will use the Divide and Conquer routine to compute only the * eigenvalues, which requires (3N + 3N**2) real workspace and * (2 + 5N + 2N lg(N)) integer workspace. * Since on many architectures DSTERF is much faster than any other * algorithm for finding eigenvalues only, it is used here * as the default. If the conditional clause is removed, then * information on the size of workspace needs to be changed. * * If COMPZ = 'N', use DSTERF to compute the eigenvalues. * IF( ICOMPZ.EQ.0 ) THEN CALL DSTERF( N, D, E, INFO ) GO TO 50 END IF * * If N is smaller than the minimum divide size (SMLSIZ+1), then * solve the problem with another solver. * IF( N.LE.SMLSIZ ) THEN * CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) * ELSE * * If COMPZ = 'V', the Z matrix must be stored elsewhere for later * use. * IF( ICOMPZ.EQ.1 ) THEN STOREZ = 1 + N*N ELSE STOREZ = 1 END IF * IF( ICOMPZ.EQ.2 ) THEN CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) END IF * * Scale. * ORGNRM = DLANST( 'M', N, D, E ) IF( ORGNRM.EQ.ZERO ) $ GO TO 50 * EPS = DLAMCH( 'Epsilon' ) * START = 1 * * while ( START <= N ) * 10 CONTINUE IF( START.LE.N ) THEN * * Let FINISH be the position of the next subdiagonal entry * such that E( FINISH ) <= TINY or FINISH = N if no such * subdiagonal exists. The matrix identified by the elements * between START and FINISH constitutes an independent * sub-problem. * FINISH = START 20 CONTINUE IF( FINISH.LT.N ) THEN TINY = EPS*SQRT( ABS( D( FINISH ) ) )* $ SQRT( ABS( D( FINISH+1 ) ) ) IF( ABS( E( FINISH ) ).GT.TINY ) THEN FINISH = FINISH + 1 GO TO 20 END IF END IF * * (Sub) Problem determined. Compute its size and solve it. * M = FINISH - START + 1 IF( M.EQ.1 ) THEN START = FINISH + 1 GO TO 10 END IF IF( M.GT.SMLSIZ ) THEN * * Scale. * ORGNRM = DLANST( 'M', M, D( START ), E( START ) ) CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M, $ INFO ) CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ), $ M-1, INFO ) * IF( ICOMPZ.EQ.1 ) THEN STRTRW = 1 ELSE STRTRW = START END IF CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ), $ Z( STRTRW, START ), LDZ, WORK( 1 ), N, $ WORK( STOREZ ), IWORK, INFO ) IF( INFO.NE.0 ) THEN INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + $ MOD( INFO, ( M+1 ) ) + START - 1 GO TO 50 END IF * * Scale back. * CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M, $ INFO ) * ELSE IF( ICOMPZ.EQ.1 ) THEN * * Since QR won't update a Z matrix which is larger than * the length of D, we must solve the sub-problem in a * workspace and then multiply back into Z. * CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M, $ WORK( M*M+1 ), INFO ) CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ, $ WORK( STOREZ ), N ) CALL DGEMM( 'N', 'N', N, M, M, ONE, $ WORK( STOREZ ), N, WORK, M, ZERO, $ Z( 1, START ), LDZ ) ELSE IF( ICOMPZ.EQ.2 ) THEN CALL DSTEQR( 'I', M, D( START ), E( START ), $ Z( START, START ), LDZ, WORK, INFO ) ELSE CALL DSTERF( M, D( START ), E( START ), INFO ) END IF IF( INFO.NE.0 ) THEN INFO = START*( N+1 ) + FINISH GO TO 50 END IF END IF * START = FINISH + 1 GO TO 10 END IF * * endwhile * * If the problem split any number of times, then the eigenvalues * will not be properly ordered. Here we permute the eigenvalues * (and the associated eigenvectors) into ascending order. * IF( M.NE.N ) THEN IF( ICOMPZ.EQ.0 ) THEN * * Use Quick Sort * CALL DLASRT( 'I', N, D, INFO ) * ELSE * * Use Selection Sort to minimize swaps of eigenvectors * DO 40 II = 2, N I = II - 1 K = I P = D( I ) DO 30 J = II, N IF( D( J ).LT.P ) THEN K = J P = D( J ) END IF 30 CONTINUE IF( K.NE.I ) THEN D( K ) = D( I ) D( I ) = P CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) END IF 40 CONTINUE END IF END IF END IF * 50 CONTINUE WORK( 1 ) = LWMIN IWORK( 1 ) = LIWMIN * RETURN * * End of DSTEDC * END SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER COMPZ INTEGER INFO, LDZ, N * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) * .. * * Purpose * ======= * * DSTEQR computes all eigenvalues and, optionally, eigenvectors of a * symmetric tridiagonal matrix using the implicit QL or QR method. * The eigenvectors of a full or band symmetric matrix can also be found * if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to * tridiagonal form. * * Arguments * ========= * * COMPZ (input) CHARACTER*1 * = 'N': Compute eigenvalues only. * = 'V': Compute eigenvalues and eigenvectors of the original * symmetric matrix. On entry, Z must contain the * orthogonal matrix used to reduce the original matrix * to tridiagonal form. * = 'I': Compute eigenvalues and eigenvectors of the * tridiagonal matrix. Z is initialized to the identity * matrix. * * N (input) INTEGER * The order of the matrix. N >= 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the diagonal elements of the tridiagonal matrix. * On exit, if INFO = 0, the eigenvalues in ascending order. * * E (input/output) DOUBLE PRECISION array, dimension (N-1) * On entry, the (n-1) subdiagonal elements of the tridiagonal * matrix. * On exit, E has been destroyed. * * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) * On entry, if COMPZ = 'V', then Z contains the orthogonal * matrix used in the reduction to tridiagonal form. * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the * orthonormal eigenvectors of the original symmetric matrix, * and if COMPZ = 'I', Z contains the orthonormal eigenvectors * of the symmetric tridiagonal matrix. * If COMPZ = 'N', then Z is not referenced. * * LDZ (input) INTEGER * The leading dimension of the array Z. LDZ >= 1, and if * eigenvectors are desired, then LDZ >= max(1,N). * * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) * If COMPZ = 'N', then WORK is not referenced. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: the algorithm has failed to find all the eigenvalues in * a total of 30*N iterations; if INFO = i, then i * elements of E have not converged to zero; on exit, D * and E contain the elements of a symmetric tridiagonal * matrix which is orthogonally similar to the original * matrix. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, THREE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ THREE = 3.0D0 ) INTEGER MAXIT PARAMETER ( MAXIT = 30 ) * .. * .. Local Scalars .. INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, $ NM1, NMAXIT DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2 * .. * .. External Subroutines .. EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR, $ DLASRT, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SIGN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( LSAME( COMPZ, 'N' ) ) THEN ICOMPZ = 0 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN ICOMPZ = 1 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN ICOMPZ = 2 ELSE ICOMPZ = -1 END IF IF( ICOMPZ.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, $ N ) ) ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSTEQR', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * IF( N.EQ.1 ) THEN IF( ICOMPZ.EQ.2 ) $ Z( 1, 1 ) = ONE RETURN END IF * * Determine the unit roundoff and over/underflow thresholds. * EPS = DLAMCH( 'E' ) EPS2 = EPS**2 SAFMIN = DLAMCH( 'S' ) SAFMAX = ONE / SAFMIN SSFMAX = SQRT( SAFMAX ) / THREE SSFMIN = SQRT( SAFMIN ) / EPS2 * * Compute the eigenvalues and eigenvectors of the tridiagonal * matrix. * IF( ICOMPZ.EQ.2 ) $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) * NMAXIT = N*MAXIT JTOT = 0 * * Determine where the matrix splits and choose QL or QR iteration * for each block, according to whether top or bottom diagonal * element is smaller. * L1 = 1 NM1 = N - 1 * 10 CONTINUE IF( L1.GT.N ) $ GO TO 160 IF( L1.GT.1 ) $ E( L1-1 ) = ZERO IF( L1.LE.NM1 ) THEN DO 20 M = L1, NM1 TST = ABS( E( M ) ) IF( TST.EQ.ZERO ) $ GO TO 30 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ $ 1 ) ) ) )*EPS ) THEN E( M ) = ZERO GO TO 30 END IF 20 CONTINUE END IF M = N * 30 CONTINUE L = L1 LSV = L LEND = M LENDSV = LEND L1 = M + 1 IF( LEND.EQ.L ) $ GO TO 10 * * Scale submatrix in rows and columns L to LEND * ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) ISCALE = 0 IF( ANORM.EQ.ZERO ) $ GO TO 10 IF( ANORM.GT.SSFMAX ) THEN ISCALE = 1 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, $ INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, $ INFO ) ELSE IF( ANORM.LT.SSFMIN ) THEN ISCALE = 2 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, $ INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, $ INFO ) END IF * * Choose between QL and QR iteration * IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN LEND = LSV L = LENDSV END IF * IF( LEND.GT.L ) THEN * * QL Iteration * * Look for small subdiagonal element. * 40 CONTINUE IF( L.NE.LEND ) THEN LENDM1 = LEND - 1 DO 50 M = L, LENDM1 TST = ABS( E( M ) )**2 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ $ SAFMIN )GO TO 60 50 CONTINUE END IF * M = LEND * 60 CONTINUE IF( M.LT.LEND ) $ E( M ) = ZERO P = D( L ) IF( M.EQ.L ) $ GO TO 80 * * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 * to compute its eigensystem. * IF( M.EQ.L+1 ) THEN IF( ICOMPZ.GT.0 ) THEN CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) WORK( L ) = C WORK( N-1+L ) = S CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ), $ WORK( N-1+L ), Z( 1, L ), LDZ ) ELSE CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) END IF D( L ) = RT1 D( L+1 ) = RT2 E( L ) = ZERO L = L + 2 IF( L.LE.LEND ) $ GO TO 40 GO TO 140 END IF * IF( JTOT.EQ.NMAXIT ) $ GO TO 140 JTOT = JTOT + 1 * * Form shift. * G = ( D( L+1 )-P ) / ( TWO*E( L ) ) R = DLAPY2( G, ONE ) G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) * S = ONE C = ONE P = ZERO * * Inner loop * MM1 = M - 1 DO 70 I = MM1, L, -1 F = S*E( I ) B = C*E( I ) CALL DLARTG( G, F, C, S, R ) IF( I.NE.M-1 ) $ E( I+1 ) = R G = D( I+1 ) - P R = ( D( I )-G )*S + TWO*C*B P = S*R D( I+1 ) = G + P G = C*R - B * * If eigenvectors are desired, then save rotations. * IF( ICOMPZ.GT.0 ) THEN WORK( I ) = C WORK( N-1+I ) = -S END IF * 70 CONTINUE * * If eigenvectors are desired, then apply saved rotations. * IF( ICOMPZ.GT.0 ) THEN MM = M - L + 1 CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), $ Z( 1, L ), LDZ ) END IF * D( L ) = D( L ) - P E( L ) = G GO TO 40 * * Eigenvalue found. * 80 CONTINUE D( L ) = P * L = L + 1 IF( L.LE.LEND ) $ GO TO 40 GO TO 140 * ELSE * * QR Iteration * * Look for small superdiagonal element. * 90 CONTINUE IF( L.NE.LEND ) THEN LENDP1 = LEND + 1 DO 100 M = L, LENDP1, -1 TST = ABS( E( M-1 ) )**2 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ $ SAFMIN )GO TO 110 100 CONTINUE END IF * M = LEND * 110 CONTINUE IF( M.GT.LEND ) $ E( M-1 ) = ZERO P = D( L ) IF( M.EQ.L ) $ GO TO 130 * * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 * to compute its eigensystem. * IF( M.EQ.L-1 ) THEN IF( ICOMPZ.GT.0 ) THEN CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) WORK( M ) = C WORK( N-1+M ) = S CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ), $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) ELSE CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) END IF D( L-1 ) = RT1 D( L ) = RT2 E( L-1 ) = ZERO L = L - 2 IF( L.GE.LEND ) $ GO TO 90 GO TO 140 END IF * IF( JTOT.EQ.NMAXIT ) $ GO TO 140 JTOT = JTOT + 1 * * Form shift. * G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) R = DLAPY2( G, ONE ) G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) * S = ONE C = ONE P = ZERO * * Inner loop * LM1 = L - 1 DO 120 I = M, LM1 F = S*E( I ) B = C*E( I ) CALL DLARTG( G, F, C, S, R ) IF( I.NE.M ) $ E( I-1 ) = R G = D( I ) - P R = ( D( I+1 )-G )*S + TWO*C*B P = S*R D( I ) = G + P G = C*R - B * * If eigenvectors are desired, then save rotations. * IF( ICOMPZ.GT.0 ) THEN WORK( I ) = C WORK( N-1+I ) = S END IF * 120 CONTINUE * * If eigenvectors are desired, then apply saved rotations. * IF( ICOMPZ.GT.0 ) THEN MM = L - M + 1 CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), $ Z( 1, M ), LDZ ) END IF * D( L ) = D( L ) - P E( LM1 ) = G GO TO 90 * * Eigenvalue found. * 130 CONTINUE D( L ) = P * L = L - 1 IF( L.GE.LEND ) $ GO TO 90 GO TO 140 * END IF * * Undo scaling if necessary * 140 CONTINUE IF( ISCALE.EQ.1 ) THEN CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, $ D( LSV ), N, INFO ) CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), $ N, INFO ) ELSE IF( ISCALE.EQ.2 ) THEN CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, $ D( LSV ), N, INFO ) CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), $ N, INFO ) END IF * * Check for no convergence to an eigenvalue after a total * of N*MAXIT iterations. * IF( JTOT.LT.NMAXIT ) $ GO TO 10 DO 150 I = 1, N - 1 IF( E( I ).NE.ZERO ) $ INFO = INFO + 1 150 CONTINUE GO TO 190 * * Order eigenvalues and eigenvectors. * 160 CONTINUE IF( ICOMPZ.EQ.0 ) THEN * * Use Quick Sort * CALL DLASRT( 'I', N, D, INFO ) * ELSE * * Use Selection Sort to minimize swaps of eigenvectors * DO 180 II = 2, N I = II - 1 K = I P = D( I ) DO 170 J = II, N IF( D( J ).LT.P ) THEN K = J P = D( J ) END IF 170 CONTINUE IF( K.NE.I ) THEN D( K ) = D( I ) D( I ) = P CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) END IF 180 CONTINUE END IF * 190 CONTINUE RETURN * * End of DSTEQR * END SUBROUTINE DSTERF( N, D, E, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, N * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), E( * ) * .. * * Purpose * ======= * * DSTERF computes all eigenvalues of a symmetric tridiagonal matrix * using the Pal-Walker-Kahan variant of the QL or QR algorithm. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix. N >= 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the n diagonal elements of the tridiagonal matrix. * On exit, if INFO = 0, the eigenvalues in ascending order. * * E (input/output) DOUBLE PRECISION array, dimension (N-1) * On entry, the (n-1) subdiagonal elements of the tridiagonal * matrix. * On exit, E has been destroyed. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: the algorithm failed to find all of the eigenvalues in * a total of 30*N iterations; if INFO = i, then i * elements of E have not converged to zero. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, THREE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ THREE = 3.0D0 ) INTEGER MAXIT PARAMETER ( MAXIT = 30 ) * .. * .. Local Scalars .. INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M, $ NMAXIT DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, $ SIGMA, SSFMAX, SSFMIN * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 EXTERNAL DLAMCH, DLANST, DLAPY2 * .. * .. External Subroutines .. EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, SIGN, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * * Quick return if possible * IF( N.LT.0 ) THEN INFO = -1 CALL XERBLA( 'DSTERF', -INFO ) RETURN END IF IF( N.LE.1 ) $ RETURN * * Determine the unit roundoff for this environment. * EPS = DLAMCH( 'E' ) EPS2 = EPS**2 SAFMIN = DLAMCH( 'S' ) SAFMAX = ONE / SAFMIN SSFMAX = SQRT( SAFMAX ) / THREE SSFMIN = SQRT( SAFMIN ) / EPS2 * * Compute the eigenvalues of the tridiagonal matrix. * NMAXIT = N*MAXIT SIGMA = ZERO JTOT = 0 * * Determine where the matrix splits and choose QL or QR iteration * for each block, according to whether top or bottom diagonal * element is smaller. * L1 = 1 * 10 CONTINUE IF( L1.GT.N ) $ GO TO 170 IF( L1.GT.1 ) $ E( L1-1 ) = ZERO DO 20 M = L1, N - 1 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ $ 1 ) ) ) )*EPS ) THEN E( M ) = ZERO GO TO 30 END IF 20 CONTINUE M = N * 30 CONTINUE L = L1 LSV = L LEND = M LENDSV = LEND L1 = M + 1 IF( LEND.EQ.L ) $ GO TO 10 * * Scale submatrix in rows and columns L to LEND * ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) ISCALE = 0 IF( ANORM.GT.SSFMAX ) THEN ISCALE = 1 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, $ INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, $ INFO ) ELSE IF( ANORM.LT.SSFMIN ) THEN ISCALE = 2 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, $ INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, $ INFO ) END IF * DO 40 I = L, LEND - 1 E( I ) = E( I )**2 40 CONTINUE * * Choose between QL and QR iteration * IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN LEND = LSV L = LENDSV END IF * IF( LEND.GE.L ) THEN * * QL Iteration * * Look for small subdiagonal element. * 50 CONTINUE IF( L.NE.LEND ) THEN DO 60 M = L, LEND - 1 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) $ GO TO 70 60 CONTINUE END IF M = LEND * 70 CONTINUE IF( M.LT.LEND ) $ E( M ) = ZERO P = D( L ) IF( M.EQ.L ) $ GO TO 90 * * If remaining matrix is 2 by 2, use DLAE2 to compute its * eigenvalues. * IF( M.EQ.L+1 ) THEN RTE = SQRT( E( L ) ) CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) D( L ) = RT1 D( L+1 ) = RT2 E( L ) = ZERO L = L + 2 IF( L.LE.LEND ) $ GO TO 50 GO TO 150 END IF * IF( JTOT.EQ.NMAXIT ) $ GO TO 150 JTOT = JTOT + 1 * * Form shift. * RTE = SQRT( E( L ) ) SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) R = DLAPY2( SIGMA, ONE ) SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) * C = ONE S = ZERO GAMMA = D( M ) - SIGMA P = GAMMA*GAMMA * * Inner loop * DO 80 I = M - 1, L, -1 BB = E( I ) R = P + BB IF( I.NE.M-1 ) $ E( I+1 ) = S*R OLDC = C C = P / R S = BB / R OLDGAM = GAMMA ALPHA = D( I ) GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) IF( C.NE.ZERO ) THEN P = ( GAMMA*GAMMA ) / C ELSE P = OLDC*BB END IF 80 CONTINUE * E( L ) = S*P D( L ) = SIGMA + GAMMA GO TO 50 * * Eigenvalue found. * 90 CONTINUE D( L ) = P * L = L + 1 IF( L.LE.LEND ) $ GO TO 50 GO TO 150 * ELSE * * QR Iteration * * Look for small superdiagonal element. * 100 CONTINUE DO 110 M = L, LEND + 1, -1 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) $ GO TO 120 110 CONTINUE M = LEND * 120 CONTINUE IF( M.GT.LEND ) $ E( M-1 ) = ZERO P = D( L ) IF( M.EQ.L ) $ GO TO 140 * * If remaining matrix is 2 by 2, use DLAE2 to compute its * eigenvalues. * IF( M.EQ.L-1 ) THEN RTE = SQRT( E( L-1 ) ) CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) D( L ) = RT1 D( L-1 ) = RT2 E( L-1 ) = ZERO L = L - 2 IF( L.GE.LEND ) $ GO TO 100 GO TO 150 END IF * IF( JTOT.EQ.NMAXIT ) $ GO TO 150 JTOT = JTOT + 1 * * Form shift. * RTE = SQRT( E( L-1 ) ) SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) R = DLAPY2( SIGMA, ONE ) SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) * C = ONE S = ZERO GAMMA = D( M ) - SIGMA P = GAMMA*GAMMA * * Inner loop * DO 130 I = M, L - 1 BB = E( I ) R = P + BB IF( I.NE.M ) $ E( I-1 ) = S*R OLDC = C C = P / R S = BB / R OLDGAM = GAMMA ALPHA = D( I+1 ) GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM D( I ) = OLDGAM + ( ALPHA-GAMMA ) IF( C.NE.ZERO ) THEN P = ( GAMMA*GAMMA ) / C ELSE P = OLDC*BB END IF 130 CONTINUE * E( L-1 ) = S*P D( L ) = SIGMA + GAMMA GO TO 100 * * Eigenvalue found. * 140 CONTINUE D( L ) = P * L = L - 1 IF( L.GE.LEND ) $ GO TO 100 GO TO 150 * END IF * * Undo scaling if necessary * 150 CONTINUE IF( ISCALE.EQ.1 ) $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, $ D( LSV ), N, INFO ) IF( ISCALE.EQ.2 ) $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, $ D( LSV ), N, INFO ) * * Check for no convergence to an eigenvalue after a total * of N*MAXIT iterations. * IF( JTOT.LT.NMAXIT ) $ GO TO 10 DO 160 I = 1, N - 1 IF( E( I ).NE.ZERO ) $ INFO = INFO + 1 160 CONTINUE GO TO 180 * * Sort eigenvalues in increasing order. * 170 CONTINUE CALL DLASRT( 'I', N, D, INFO ) * 180 CONTINUE RETURN * * End of DSTERF * END SUBROUTINE DSTEVD( JOBZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, $ LIWORK, INFO ) * * -- LAPACK driver routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER JOBZ INTEGER INFO, LDZ, LIWORK, LWORK, N * .. * .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) * .. * * Purpose * ======= * * DSTEVD computes all eigenvalues and, optionally, eigenvectors of a * real symmetric tridiagonal matrix. If eigenvectors are desired, it * uses a divide and conquer algorithm. * * The divide and conquer algorithm makes very mild assumptions about * floating point arithmetic. It will work on machines with a guard * digit in add/subtract, or on those binary machines without guard * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or * Cray-2. It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. * * Arguments * ========= * * JOBZ (input) CHARACTER*1 * = 'N': Compute eigenvalues only; * = 'V': Compute eigenvalues and eigenvectors. * * N (input) INTEGER * The order of the matrix. N >= 0. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the n diagonal elements of the tridiagonal matrix * A. * On exit, if INFO = 0, the eigenvalues in ascending order. * * E (input/output) DOUBLE PRECISION array, dimension (N-1) * On entry, the (n-1) subdiagonal elements of the tridiagonal * matrix A, stored in elements 1 to N-1 of E. * On exit, the contents of E are destroyed. * * Z (output) DOUBLE PRECISION array, dimension (LDZ, N) * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal * eigenvectors of the matrix A, with the i-th column of Z * holding the eigenvector associated with D(i). * If JOBZ = 'N', then Z is not referenced. * * LDZ (input) INTEGER * The leading dimension of the array Z. LDZ >= 1, and if * JOBZ = 'V', LDZ >= max(1,N). * * WORK (workspace/output) DOUBLE PRECISION array, * dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. * If JOBZ = 'N' or N <= 1 then LWORK must be at least 1. * If JOBZ = 'V' and N > 1 then LWORK must be at least * ( 1 + 4*N + N**2 ). * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal sizes of the WORK and IWORK * arrays, returns these values as the first entries of the WORK * and IWORK arrays, and no error message related to LWORK or * LIWORK is issued by XERBLA. * * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. * * LIWORK (input) INTEGER * The dimension of the array IWORK. * If JOBZ = 'N' or N <= 1 then LIWORK must be at least 1. * If JOBZ = 'V' and N > 1 then LIWORK must be at least 3+5*N. * * If LIWORK = -1, then a workspace query is assumed; the * routine only calculates the optimal sizes of the WORK and * IWORK arrays, returns these values as the first entries of * the WORK and IWORK arrays, and no error message related to * LWORK or LIWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, the algorithm failed to converge; i * off-diagonal elements of E did not converge to zero. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. LOGICAL LQUERY, WANTZ INTEGER ISCALE, LIWMIN, LWMIN DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, $ TNRM * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH, DLANST EXTERNAL LSAME, DLAMCH, DLANST * .. * .. External Subroutines .. EXTERNAL DSCAL, DSTEDC, DSTERF, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC SQRT * .. * .. Executable Statements .. * * Test the input parameters. * WANTZ = LSAME( JOBZ, 'V' ) LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) * INFO = 0 LIWMIN = 1 LWMIN = 1 IF( N.GT.1 .AND. WANTZ ) THEN LWMIN = 1 + 4*N + N**2 LIWMIN = 3 + 5*N END IF * IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN INFO = -6 END IF * IF( INFO.EQ.0 ) THEN WORK( 1 ) = LWMIN IWORK( 1 ) = LIWMIN * IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN INFO = -8 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN INFO = -10 END IF END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSTEVD', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * IF( N.EQ.1 ) THEN IF( WANTZ ) $ Z( 1, 1 ) = ONE RETURN END IF * * Get machine constants. * SAFMIN = DLAMCH( 'Safe minimum' ) EPS = DLAMCH( 'Precision' ) SMLNUM = SAFMIN / EPS BIGNUM = ONE / SMLNUM RMIN = SQRT( SMLNUM ) RMAX = SQRT( BIGNUM ) * * Scale matrix to allowable range, if necessary. * ISCALE = 0 TNRM = DLANST( 'M', N, D, E ) IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN ISCALE = 1 SIGMA = RMIN / TNRM ELSE IF( TNRM.GT.RMAX ) THEN ISCALE = 1 SIGMA = RMAX / TNRM END IF IF( ISCALE.EQ.1 ) THEN CALL DSCAL( N, SIGMA, D, 1 ) CALL DSCAL( N-1, SIGMA, E( 1 ), 1 ) END IF * * For eigenvalues only, call DSTERF. For eigenvalues and * eigenvectors, call DSTEDC. * IF( .NOT.WANTZ ) THEN CALL DSTERF( N, D, E, INFO ) ELSE CALL DSTEDC( 'I', N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, $ INFO ) END IF * * If matrix was scaled, then rescale eigenvalues appropriately. * IF( ISCALE.EQ.1 ) $ CALL DSCAL( N, ONE / SIGMA, D, 1 ) * WORK( 1 ) = LWMIN IWORK( 1 ) = LIWMIN * RETURN * * End of DSTEVD * END