*** mlphare.orig.f	Wed Feb 21 11:04:30 2001
--- mlphare.f	Fri Feb 23 14:34:56 2001
***************
*** 4522,4528 ****
        A(3,3) = DBLE(BB(6))
  C
  C          *****************
!       CALL EIGN1M(3,3,A,EIG)
  C          *****************
  C
        WRITE (6,FMT=6000) EIG
--- 4522,4533 ----
        A(3,3) = DBLE(BB(6))
  C
  C          *****************
! C PJX - 23/02/01:
! C The call to EIGN1M has been replaced by the call to EIGN1M_W_LAPACK
! C This routine has the same functionality as EIGNM1 but uses the
! C LAPACK routine DSYEVR to accomplish the task instead
! C      CALL EIGN1M(3,3,A,EIG)
!       CALL EIGN1M_W_LAPACK(3,3,A,EIG)
  C          *****************
  C
        WRITE (6,FMT=6000) EIG
***************
*** 4836,4841 ****
--- 4841,4970 ----
        END
  C
  C
+ C
+ C     =======================================
+       SUBROUTINE EIGN1M_W_LAPACK(N,NDIM,A,EIG)
+ C     =======================================
+ C
+ C     PJX (23/02/01)
+ C
+ C     Find eigenvalues and vectors of a real symmetric
+ C     matrix.
+ C
+ C     Alternative to EIGN1M, using a call to the LAPACK
+ C     routine DSYEVR, to demonstrate use of LAPACK libraries
+ C     in CCP4.
+ C
+ C     Arguments:
+ C     A    (input)  double array (NDIM,N)
+ C                   contains a symmetric NxN matrix
+ C          (output) eigenvectors as columns, ordering same as the
+ C                   eigenvalues in EIG
+ C     NDIM (input)  leading dimension of array A
+ C     N    (input)  dimension of A
+ C     EIG  (output) double array (N)
+ C                   contains the eigenvalues, in ascending order
+ C
+       IMPLICIT NONE
+ C
+ C     ..Parameters
+ C     ..These are used for setting up workspaces needed by DSYEVR
+       INTEGER MAXEIGEN
+       PARAMETER (MAXEIGEN=3)
+       INTEGER MAXEIGEN2
+       PARAMETER (MAXEIGEN2=MAXEIGEN*MAXEIGEN)
+       INTEGER LWORK
+       PARAMETER (LWORK=26*MAXEIGEN)
+       INTEGER LIWORK
+       PARAMETER (LIWORK=10*MAXEIGEN)
+ C
+ C     ..Arguments
+       INTEGER N,NDIM
+       DOUBLE PRECISION A(NDIM,N),EIG(N)
+ C
+ C     ..Local scalars
+       INTEGER   IDO,JDO
+       INTEGER   LDA,LDZ,IL,IU,M,INFO
+       DOUBLE PRECISION VL,VU,ABSTOL
+       CHARACTER JOBZ*1,RANGE*1,UPLO*1
+       LOGICAL   LSYMM
+ C
+ C     ..Local arrays
+       INTEGER   ISUPPZ(2*MAXEIGEN),IWORK(LIWORK)
+       DOUBLE PRECISION WORK(LWORK),VECTORS(MAXEIGEN,MAXEIGEN)
+ C     
+ C     ..External subroutines
+       EXTERNAL DSYEVR,CCPERR
+ C
+ C---- Initial checking
+ C     MAXEIGEN is the maximum number of eigenvalues
+ C     It is used to dimension the work/internal arrays
+       IF (N.GT.MAXEIGEN .OR. NDIM.GT.MAXEIGEN) THEN
+         WRITE(6,FMT=9001)
+  9001   FORMAT(1X,' SUBROUTINE EIGN1M_W_LAPACK',/,
+      +       4X,'MAXEIGEN = ',I4,/,
+      +       4X,'N        = ',I4,/,
+      +       4X,'NDIM     = ',I4,/,
+      +       1X,'MAXEIGEN must be >= N and NDIM')
+         CALL CCPERR(1,' S/R EIGN1M_W_LAPACK: increase MAXEIGEN')
+       END IF
+ C
+ C     Check that A is symmetric
+       LSYMM  = .TRUE.
+       DO IDO = 1,3
+         DO JDO = IDO,3
+           IF (IDO.NE.JDO .AND. (DABS(A(IDO,JDO)-A(JDO,IDO)).GT.1.0E-8))
+      +         LSYMM = .FALSE.
+         END DO
+       END DO
+       IF (.NOT. LSYMM) THEN
+         CALL CCPERR(1,' S/R EIGN1M_W_LAPACK: matrix is not symmetric')
+       END IF
+ C
+ C---- Set up input parameters for DSYEVR
+ C
+ C     JOBZ is V (compute eigenvalues&vectors)
+       JOBZ = 'V'
+ C     RANGE is A (find all eigenvalues/vectors)
+       RANGE = 'A'
+ C     LDA defines the leading dimension of the input array,
+ C     LDZ defines that of the output eigenvector array
+       LDA = NDIM
+       LDZ = MAXEIGEN
+ C     VL,VU and IL,IU not used for RANGE = A (give dummy values)
+       VL = 0.0
+       VU = 0.0
+       IL = 0
+       IU = 0
+ C     ABSTOL sets the absolute tolerance of the eigenvalues
+ C     Set this to zero so DSYEVR calculates its own value
+       ABSTOL = 0.0
+ C     UPLO: A stores upper (U) or lower (L) triangle of
+ C     the original matrix
+ C     In the calling routine A is set as upper-triangular
+       UPLO = 'U'
+ C
+ C     Do the calculations
+       CALL DSYEVR(JOBZ,RANGE,UPLO,N,A,LDA,VL,VU,IL,IU,
+      +     ABSTOL,M,EIG,VECTORS,LDZ,ISUPPZ,WORK,LWORK,IWORK,
+      +     LIWORK,INFO)
+ C
+ C     Check output status for successful completion
+       IF (INFO.NE.0) THEN
+         CALL CCPERR(1,' S/R EIGN1M_W_LAPACK: error from DSYEVR')
+       END IF
+ C
+ C     EIG already contains the eigenvalues in ascending order
+ C     Load the eigenvectors into A before returning
+       DO IDO = 1,NDIM
+         DO JDO = 1,N
+           A(IDO,JDO) = VECTORS(IDO,JDO)
+         END DO
+       END DO
+ C
+ C     Successful completion: return to calling subprogram
+       RETURN
+       END
  C
  C     ===============================
        SUBROUTINE EIGN1M(N,NDIM,A,EIG)