Annotation of rpl/lapack/lapack/zlaed0.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
! 2: $ IWORK, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER INFO, LDQ, LDQS, N, QSIZ
! 11: * ..
! 12: * .. Array Arguments ..
! 13: INTEGER IWORK( * )
! 14: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
! 15: COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * Using the divide and conquer method, ZLAED0 computes all eigenvalues
! 22: * of a symmetric tridiagonal matrix which is one diagonal block of
! 23: * those from reducing a dense or band Hermitian matrix and
! 24: * corresponding eigenvectors of the dense or band matrix.
! 25: *
! 26: * Arguments
! 27: * =========
! 28: *
! 29: * QSIZ (input) INTEGER
! 30: * The dimension of the unitary matrix used to reduce
! 31: * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
! 32: *
! 33: * N (input) INTEGER
! 34: * The dimension of the symmetric tridiagonal matrix. N >= 0.
! 35: *
! 36: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 37: * On entry, the diagonal elements of the tridiagonal matrix.
! 38: * On exit, the eigenvalues in ascending order.
! 39: *
! 40: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
! 41: * On entry, the off-diagonal elements of the tridiagonal matrix.
! 42: * On exit, E has been destroyed.
! 43: *
! 44: * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
! 45: * On entry, Q must contain an QSIZ x N matrix whose columns
! 46: * unitarily orthonormal. It is a part of the unitary matrix
! 47: * that reduces the full dense Hermitian matrix to a
! 48: * (reducible) symmetric tridiagonal matrix.
! 49: *
! 50: * LDQ (input) INTEGER
! 51: * The leading dimension of the array Q. LDQ >= max(1,N).
! 52: *
! 53: * IWORK (workspace) INTEGER array,
! 54: * the dimension of IWORK must be at least
! 55: * 6 + 6*N + 5*N*lg N
! 56: * ( lg( N ) = smallest integer k
! 57: * such that 2^k >= N )
! 58: *
! 59: * RWORK (workspace) DOUBLE PRECISION array,
! 60: * dimension (1 + 3*N + 2*N*lg N + 3*N**2)
! 61: * ( lg( N ) = smallest integer k
! 62: * such that 2^k >= N )
! 63: *
! 64: * QSTORE (workspace) COMPLEX*16 array, dimension (LDQS, N)
! 65: * Used to store parts of
! 66: * the eigenvector matrix when the updating matrix multiplies
! 67: * take place.
! 68: *
! 69: * LDQS (input) INTEGER
! 70: * The leading dimension of the array QSTORE.
! 71: * LDQS >= max(1,N).
! 72: *
! 73: * INFO (output) INTEGER
! 74: * = 0: successful exit.
! 75: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 76: * > 0: The algorithm failed to compute an eigenvalue while
! 77: * working on the submatrix lying in rows and columns
! 78: * INFO/(N+1) through mod(INFO,N+1).
! 79: *
! 80: * =====================================================================
! 81: *
! 82: * Warning: N could be as big as QSIZ!
! 83: *
! 84: * .. Parameters ..
! 85: DOUBLE PRECISION TWO
! 86: PARAMETER ( TWO = 2.D+0 )
! 87: * ..
! 88: * .. Local Scalars ..
! 89: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
! 90: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
! 91: $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
! 92: $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
! 93: DOUBLE PRECISION TEMP
! 94: * ..
! 95: * .. External Subroutines ..
! 96: EXTERNAL DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
! 97: * ..
! 98: * .. External Functions ..
! 99: INTEGER ILAENV
! 100: EXTERNAL ILAENV
! 101: * ..
! 102: * .. Intrinsic Functions ..
! 103: INTRINSIC ABS, DBLE, INT, LOG, MAX
! 104: * ..
! 105: * .. Executable Statements ..
! 106: *
! 107: * Test the input parameters.
! 108: *
! 109: INFO = 0
! 110: *
! 111: * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
! 112: * INFO = -1
! 113: * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
! 114: * $ THEN
! 115: IF( QSIZ.LT.MAX( 0, N ) ) THEN
! 116: INFO = -1
! 117: ELSE IF( N.LT.0 ) THEN
! 118: INFO = -2
! 119: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
! 120: INFO = -6
! 121: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
! 122: INFO = -8
! 123: END IF
! 124: IF( INFO.NE.0 ) THEN
! 125: CALL XERBLA( 'ZLAED0', -INFO )
! 126: RETURN
! 127: END IF
! 128: *
! 129: * Quick return if possible
! 130: *
! 131: IF( N.EQ.0 )
! 132: $ RETURN
! 133: *
! 134: SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
! 135: *
! 136: * Determine the size and placement of the submatrices, and save in
! 137: * the leading elements of IWORK.
! 138: *
! 139: IWORK( 1 ) = N
! 140: SUBPBS = 1
! 141: TLVLS = 0
! 142: 10 CONTINUE
! 143: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
! 144: DO 20 J = SUBPBS, 1, -1
! 145: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
! 146: IWORK( 2*J-1 ) = IWORK( J ) / 2
! 147: 20 CONTINUE
! 148: TLVLS = TLVLS + 1
! 149: SUBPBS = 2*SUBPBS
! 150: GO TO 10
! 151: END IF
! 152: DO 30 J = 2, SUBPBS
! 153: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
! 154: 30 CONTINUE
! 155: *
! 156: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
! 157: * using rank-1 modifications (cuts).
! 158: *
! 159: SPM1 = SUBPBS - 1
! 160: DO 40 I = 1, SPM1
! 161: SUBMAT = IWORK( I ) + 1
! 162: SMM1 = SUBMAT - 1
! 163: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
! 164: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
! 165: 40 CONTINUE
! 166: *
! 167: INDXQ = 4*N + 3
! 168: *
! 169: * Set up workspaces for eigenvalues only/accumulate new vectors
! 170: * routine
! 171: *
! 172: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
! 173: LGN = INT( TEMP )
! 174: IF( 2**LGN.LT.N )
! 175: $ LGN = LGN + 1
! 176: IF( 2**LGN.LT.N )
! 177: $ LGN = LGN + 1
! 178: IPRMPT = INDXQ + N + 1
! 179: IPERM = IPRMPT + N*LGN
! 180: IQPTR = IPERM + N*LGN
! 181: IGIVPT = IQPTR + N + 2
! 182: IGIVCL = IGIVPT + N*LGN
! 183: *
! 184: IGIVNM = 1
! 185: IQ = IGIVNM + 2*N*LGN
! 186: IWREM = IQ + N**2 + 1
! 187: * Initialize pointers
! 188: DO 50 I = 0, SUBPBS
! 189: IWORK( IPRMPT+I ) = 1
! 190: IWORK( IGIVPT+I ) = 1
! 191: 50 CONTINUE
! 192: IWORK( IQPTR ) = 1
! 193: *
! 194: * Solve each submatrix eigenproblem at the bottom of the divide and
! 195: * conquer tree.
! 196: *
! 197: CURR = 0
! 198: DO 70 I = 0, SPM1
! 199: IF( I.EQ.0 ) THEN
! 200: SUBMAT = 1
! 201: MATSIZ = IWORK( 1 )
! 202: ELSE
! 203: SUBMAT = IWORK( I ) + 1
! 204: MATSIZ = IWORK( I+1 ) - IWORK( I )
! 205: END IF
! 206: LL = IQ - 1 + IWORK( IQPTR+CURR )
! 207: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
! 208: $ RWORK( LL ), MATSIZ, RWORK, INFO )
! 209: CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
! 210: $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
! 211: $ RWORK( IWREM ) )
! 212: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
! 213: CURR = CURR + 1
! 214: IF( INFO.GT.0 ) THEN
! 215: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
! 216: RETURN
! 217: END IF
! 218: K = 1
! 219: DO 60 J = SUBMAT, IWORK( I+1 )
! 220: IWORK( INDXQ+J ) = K
! 221: K = K + 1
! 222: 60 CONTINUE
! 223: 70 CONTINUE
! 224: *
! 225: * Successively merge eigensystems of adjacent submatrices
! 226: * into eigensystem for the corresponding larger matrix.
! 227: *
! 228: * while ( SUBPBS > 1 )
! 229: *
! 230: CURLVL = 1
! 231: 80 CONTINUE
! 232: IF( SUBPBS.GT.1 ) THEN
! 233: SPM2 = SUBPBS - 2
! 234: DO 90 I = 0, SPM2, 2
! 235: IF( I.EQ.0 ) THEN
! 236: SUBMAT = 1
! 237: MATSIZ = IWORK( 2 )
! 238: MSD2 = IWORK( 1 )
! 239: CURPRB = 0
! 240: ELSE
! 241: SUBMAT = IWORK( I ) + 1
! 242: MATSIZ = IWORK( I+2 ) - IWORK( I )
! 243: MSD2 = MATSIZ / 2
! 244: CURPRB = CURPRB + 1
! 245: END IF
! 246: *
! 247: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
! 248: * into an eigensystem of size MATSIZ. ZLAED7 handles the case
! 249: * when the eigenvectors of a full or band Hermitian matrix (which
! 250: * was reduced to tridiagonal form) are desired.
! 251: *
! 252: * I am free to use Q as a valuable working space until Loop 150.
! 253: *
! 254: CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
! 255: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
! 256: $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
! 257: $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
! 258: $ IWORK( IPERM ), IWORK( IGIVPT ),
! 259: $ IWORK( IGIVCL ), RWORK( IGIVNM ),
! 260: $ Q( 1, SUBMAT ), RWORK( IWREM ),
! 261: $ IWORK( SUBPBS+1 ), INFO )
! 262: IF( INFO.GT.0 ) THEN
! 263: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
! 264: RETURN
! 265: END IF
! 266: IWORK( I / 2+1 ) = IWORK( I+2 )
! 267: 90 CONTINUE
! 268: SUBPBS = SUBPBS / 2
! 269: CURLVL = CURLVL + 1
! 270: GO TO 80
! 271: END IF
! 272: *
! 273: * end while
! 274: *
! 275: * Re-merge the eigenvalues/vectors which were deflated at the final
! 276: * merge step.
! 277: *
! 278: DO 100 I = 1, N
! 279: J = IWORK( INDXQ+I )
! 280: RWORK( I ) = D( J )
! 281: CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
! 282: 100 CONTINUE
! 283: CALL DCOPY( N, RWORK, 1, D, 1 )
! 284: *
! 285: RETURN
! 286: *
! 287: * End of ZLAED0
! 288: *
! 289: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaed0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack