Annotation of rpl/lapack/lapack/dlaed0.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
! 2: $ WORK, 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 ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
! 11: * ..
! 12: * .. Array Arguments ..
! 13: INTEGER IWORK( * )
! 14: DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
! 15: $ WORK( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DLAED0 computes all eigenvalues and corresponding eigenvectors of a
! 22: * symmetric tridiagonal matrix using the divide and conquer method.
! 23: *
! 24: * Arguments
! 25: * =========
! 26: *
! 27: * ICOMPQ (input) INTEGER
! 28: * = 0: Compute eigenvalues only.
! 29: * = 1: Compute eigenvectors of original dense symmetric matrix
! 30: * also. On entry, Q contains the orthogonal matrix used
! 31: * to reduce the original matrix to tridiagonal form.
! 32: * = 2: Compute eigenvalues and eigenvectors of tridiagonal
! 33: * matrix.
! 34: *
! 35: * QSIZ (input) INTEGER
! 36: * The dimension of the orthogonal matrix used to reduce
! 37: * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
! 38: *
! 39: * N (input) INTEGER
! 40: * The dimension of the symmetric tridiagonal matrix. N >= 0.
! 41: *
! 42: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 43: * On entry, the main diagonal of the tridiagonal matrix.
! 44: * On exit, its eigenvalues.
! 45: *
! 46: * E (input) DOUBLE PRECISION array, dimension (N-1)
! 47: * The off-diagonal elements of the tridiagonal matrix.
! 48: * On exit, E has been destroyed.
! 49: *
! 50: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
! 51: * On entry, Q must contain an N-by-N orthogonal matrix.
! 52: * If ICOMPQ = 0 Q is not referenced.
! 53: * If ICOMPQ = 1 On entry, Q is a subset of the columns of the
! 54: * orthogonal matrix used to reduce the full
! 55: * matrix to tridiagonal form corresponding to
! 56: * the subset of the full matrix which is being
! 57: * decomposed at this time.
! 58: * If ICOMPQ = 2 On entry, Q will be the identity matrix.
! 59: * On exit, Q contains the eigenvectors of the
! 60: * tridiagonal matrix.
! 61: *
! 62: * LDQ (input) INTEGER
! 63: * The leading dimension of the array Q. If eigenvectors are
! 64: * desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
! 65: *
! 66: * QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
! 67: * Referenced only when ICOMPQ = 1. Used to store parts of
! 68: * the eigenvector matrix when the updating matrix multiplies
! 69: * take place.
! 70: *
! 71: * LDQS (input) INTEGER
! 72: * The leading dimension of the array QSTORE. If ICOMPQ = 1,
! 73: * then LDQS >= max(1,N). In any case, LDQS >= 1.
! 74: *
! 75: * WORK (workspace) DOUBLE PRECISION array,
! 76: * If ICOMPQ = 0 or 1, the dimension of WORK must be at least
! 77: * 1 + 3*N + 2*N*lg N + 2*N**2
! 78: * ( lg( N ) = smallest integer k
! 79: * such that 2^k >= N )
! 80: * If ICOMPQ = 2, the dimension of WORK must be at least
! 81: * 4*N + N**2.
! 82: *
! 83: * IWORK (workspace) INTEGER array,
! 84: * If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
! 85: * 6 + 6*N + 5*N*lg N.
! 86: * ( lg( N ) = smallest integer k
! 87: * such that 2^k >= N )
! 88: * If ICOMPQ = 2, the dimension of IWORK must be at least
! 89: * 3 + 5*N.
! 90: *
! 91: * INFO (output) INTEGER
! 92: * = 0: successful exit.
! 93: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 94: * > 0: The algorithm failed to compute an eigenvalue while
! 95: * working on the submatrix lying in rows and columns
! 96: * INFO/(N+1) through mod(INFO,N+1).
! 97: *
! 98: * Further Details
! 99: * ===============
! 100: *
! 101: * Based on contributions by
! 102: * Jeff Rutter, Computer Science Division, University of California
! 103: * at Berkeley, USA
! 104: *
! 105: * =====================================================================
! 106: *
! 107: * .. Parameters ..
! 108: DOUBLE PRECISION ZERO, ONE, TWO
! 109: PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
! 110: * ..
! 111: * .. Local Scalars ..
! 112: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
! 113: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
! 114: $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
! 115: $ SPM2, SUBMAT, SUBPBS, TLVLS
! 116: DOUBLE PRECISION TEMP
! 117: * ..
! 118: * .. External Subroutines ..
! 119: EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
! 120: $ XERBLA
! 121: * ..
! 122: * .. External Functions ..
! 123: INTEGER ILAENV
! 124: EXTERNAL ILAENV
! 125: * ..
! 126: * .. Intrinsic Functions ..
! 127: INTRINSIC ABS, DBLE, INT, LOG, MAX
! 128: * ..
! 129: * .. Executable Statements ..
! 130: *
! 131: * Test the input parameters.
! 132: *
! 133: INFO = 0
! 134: *
! 135: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
! 136: INFO = -1
! 137: ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
! 138: INFO = -2
! 139: ELSE IF( N.LT.0 ) THEN
! 140: INFO = -3
! 141: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
! 142: INFO = -7
! 143: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
! 144: INFO = -9
! 145: END IF
! 146: IF( INFO.NE.0 ) THEN
! 147: CALL XERBLA( 'DLAED0', -INFO )
! 148: RETURN
! 149: END IF
! 150: *
! 151: * Quick return if possible
! 152: *
! 153: IF( N.EQ.0 )
! 154: $ RETURN
! 155: *
! 156: SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
! 157: *
! 158: * Determine the size and placement of the submatrices, and save in
! 159: * the leading elements of IWORK.
! 160: *
! 161: IWORK( 1 ) = N
! 162: SUBPBS = 1
! 163: TLVLS = 0
! 164: 10 CONTINUE
! 165: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
! 166: DO 20 J = SUBPBS, 1, -1
! 167: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
! 168: IWORK( 2*J-1 ) = IWORK( J ) / 2
! 169: 20 CONTINUE
! 170: TLVLS = TLVLS + 1
! 171: SUBPBS = 2*SUBPBS
! 172: GO TO 10
! 173: END IF
! 174: DO 30 J = 2, SUBPBS
! 175: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
! 176: 30 CONTINUE
! 177: *
! 178: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
! 179: * using rank-1 modifications (cuts).
! 180: *
! 181: SPM1 = SUBPBS - 1
! 182: DO 40 I = 1, SPM1
! 183: SUBMAT = IWORK( I ) + 1
! 184: SMM1 = SUBMAT - 1
! 185: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
! 186: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
! 187: 40 CONTINUE
! 188: *
! 189: INDXQ = 4*N + 3
! 190: IF( ICOMPQ.NE.2 ) THEN
! 191: *
! 192: * Set up workspaces for eigenvalues only/accumulate new vectors
! 193: * routine
! 194: *
! 195: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
! 196: LGN = INT( TEMP )
! 197: IF( 2**LGN.LT.N )
! 198: $ LGN = LGN + 1
! 199: IF( 2**LGN.LT.N )
! 200: $ LGN = LGN + 1
! 201: IPRMPT = INDXQ + N + 1
! 202: IPERM = IPRMPT + N*LGN
! 203: IQPTR = IPERM + N*LGN
! 204: IGIVPT = IQPTR + N + 2
! 205: IGIVCL = IGIVPT + N*LGN
! 206: *
! 207: IGIVNM = 1
! 208: IQ = IGIVNM + 2*N*LGN
! 209: IWREM = IQ + N**2 + 1
! 210: *
! 211: * Initialize pointers
! 212: *
! 213: DO 50 I = 0, SUBPBS
! 214: IWORK( IPRMPT+I ) = 1
! 215: IWORK( IGIVPT+I ) = 1
! 216: 50 CONTINUE
! 217: IWORK( IQPTR ) = 1
! 218: END IF
! 219: *
! 220: * Solve each submatrix eigenproblem at the bottom of the divide and
! 221: * conquer tree.
! 222: *
! 223: CURR = 0
! 224: DO 70 I = 0, SPM1
! 225: IF( I.EQ.0 ) THEN
! 226: SUBMAT = 1
! 227: MATSIZ = IWORK( 1 )
! 228: ELSE
! 229: SUBMAT = IWORK( I ) + 1
! 230: MATSIZ = IWORK( I+1 ) - IWORK( I )
! 231: END IF
! 232: IF( ICOMPQ.EQ.2 ) THEN
! 233: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
! 234: $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
! 235: IF( INFO.NE.0 )
! 236: $ GO TO 130
! 237: ELSE
! 238: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
! 239: $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
! 240: $ INFO )
! 241: IF( INFO.NE.0 )
! 242: $ GO TO 130
! 243: IF( ICOMPQ.EQ.1 ) THEN
! 244: CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
! 245: $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
! 246: $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
! 247: $ LDQS )
! 248: END IF
! 249: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
! 250: CURR = CURR + 1
! 251: END IF
! 252: K = 1
! 253: DO 60 J = SUBMAT, IWORK( I+1 )
! 254: IWORK( INDXQ+J ) = K
! 255: K = K + 1
! 256: 60 CONTINUE
! 257: 70 CONTINUE
! 258: *
! 259: * Successively merge eigensystems of adjacent submatrices
! 260: * into eigensystem for the corresponding larger matrix.
! 261: *
! 262: * while ( SUBPBS > 1 )
! 263: *
! 264: CURLVL = 1
! 265: 80 CONTINUE
! 266: IF( SUBPBS.GT.1 ) THEN
! 267: SPM2 = SUBPBS - 2
! 268: DO 90 I = 0, SPM2, 2
! 269: IF( I.EQ.0 ) THEN
! 270: SUBMAT = 1
! 271: MATSIZ = IWORK( 2 )
! 272: MSD2 = IWORK( 1 )
! 273: CURPRB = 0
! 274: ELSE
! 275: SUBMAT = IWORK( I ) + 1
! 276: MATSIZ = IWORK( I+2 ) - IWORK( I )
! 277: MSD2 = MATSIZ / 2
! 278: CURPRB = CURPRB + 1
! 279: END IF
! 280: *
! 281: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
! 282: * into an eigensystem of size MATSIZ.
! 283: * DLAED1 is used only for the full eigensystem of a tridiagonal
! 284: * matrix.
! 285: * DLAED7 handles the cases in which eigenvalues only or eigenvalues
! 286: * and eigenvectors of a full symmetric matrix (which was reduced to
! 287: * tridiagonal form) are desired.
! 288: *
! 289: IF( ICOMPQ.EQ.2 ) THEN
! 290: CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
! 291: $ LDQ, IWORK( INDXQ+SUBMAT ),
! 292: $ E( SUBMAT+MSD2-1 ), MSD2, WORK,
! 293: $ IWORK( SUBPBS+1 ), INFO )
! 294: ELSE
! 295: CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
! 296: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
! 297: $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
! 298: $ MSD2, WORK( IQ ), IWORK( IQPTR ),
! 299: $ IWORK( IPRMPT ), IWORK( IPERM ),
! 300: $ IWORK( IGIVPT ), IWORK( IGIVCL ),
! 301: $ WORK( IGIVNM ), WORK( IWREM ),
! 302: $ IWORK( SUBPBS+1 ), INFO )
! 303: END IF
! 304: IF( INFO.NE.0 )
! 305: $ GO TO 130
! 306: IWORK( I / 2+1 ) = IWORK( I+2 )
! 307: 90 CONTINUE
! 308: SUBPBS = SUBPBS / 2
! 309: CURLVL = CURLVL + 1
! 310: GO TO 80
! 311: END IF
! 312: *
! 313: * end while
! 314: *
! 315: * Re-merge the eigenvalues/vectors which were deflated at the final
! 316: * merge step.
! 317: *
! 318: IF( ICOMPQ.EQ.1 ) THEN
! 319: DO 100 I = 1, N
! 320: J = IWORK( INDXQ+I )
! 321: WORK( I ) = D( J )
! 322: CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
! 323: 100 CONTINUE
! 324: CALL DCOPY( N, WORK, 1, D, 1 )
! 325: ELSE IF( ICOMPQ.EQ.2 ) THEN
! 326: DO 110 I = 1, N
! 327: J = IWORK( INDXQ+I )
! 328: WORK( I ) = D( J )
! 329: CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
! 330: 110 CONTINUE
! 331: CALL DCOPY( N, WORK, 1, D, 1 )
! 332: CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
! 333: ELSE
! 334: DO 120 I = 1, N
! 335: J = IWORK( INDXQ+I )
! 336: WORK( I ) = D( J )
! 337: 120 CONTINUE
! 338: CALL DCOPY( N, WORK, 1, D, 1 )
! 339: END IF
! 340: GO TO 140
! 341: *
! 342: 130 CONTINUE
! 343: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
! 344: *
! 345: 140 CONTINUE
! 346: RETURN
! 347: *
! 348: * End of DLAED0
! 349: *
! 350: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack