c
c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
C
C
C
      PROGRAM EXAMPLE_LAPACK_2
C

C
C
C ***** TO DEMONSTRATE THE USE OF THE LAPACK EIGEN-SOLVER.
C
C ***** THIS PROGRAM ILLUSTRATES THE FOLLOWING:
C
C       (1) USE OF THE SIMPLE DRIVER ROUTINE TO OBTAIN THE EIGENVALUES OF
C           THE TEST PROBLEM STORED AS A SYMMETRIC BANDED MATRIX.
C       (2) USE OF THE SINGULAR VALUE DECOMPOSITION ROUTINE TO OBTAIN THE 
C           SINGULAR VALUES OF THE TEST PROBLEM STORED AS A FULL MATRIX.
C
C ***** THE TEST PROBLEM USED IS THE TEST PROBLEM USED IN THE CXML 
C       SPARSE SOLVERS. THIS IS A SYMMETRIC POSITIVE DEFINITE MATRIX.
C
C ***** SAMPLE OUTPUT FOLLOWS:
C
C   ILLUSTRATING THE USE OF LAPACK EIGEN-SOLVER AND SINGULAR VALUE DECOMPOSITION
C   ROUTINES.
C      TEST MATRIX IS SYMMETRIC POSITIVE DEFINITE
C      ORDER OF THE MATRIX:   16
C          DETERMINING THE EIGENVALUES, WITH THE MATRIX STORED IN A SYMMETRIC
C           BANDED FORM
C          EIGENVALUES IN ASCENDING ORDER:
C              0.76393202E+00
C              0.17639320E+01
C              0.17639320E+01
C              0.27639320E+01
C              0.30000000E+01
C              0.30000000E+01
C              0.40000000E+01
C              0.40000000E+01
C              0.40000000E+01
C              0.40000000E+01
C              0.50000000E+01
C              0.50000000E+01
C              0.52360680E+01
C              0.62360680E+01
C              0.62360680E+01
C              0.72360680E+01
C          DETERMINING THE SINGULAR VALUES, WITH THE MATRIX STORED IN A FULL FORM
C          SINGULAR VALUES IN DESCENDING ORDER:
C              0.72360680E+01
C              0.62360680E+01
C              0.62360680E+01
C              0.52360680E+01
C              0.50000000E+01
C              0.50000000E+01
C              0.40000000E+01
C              0.40000000E+01
C              0.40000000E+01
C              0.40000000E+01
C              0.30000000E+01
C              0.30000000E+01
C              0.27639320E+01
C              0.17639320E+01
C              0.17639320E+01
C              0.76393202E+00
C      MAXIMUM ERROR BETWEEN THE SINGULAR VALUES AND THE EIGENVALUES:
C       0.62172489E-14
C
C
C
C      IMPLICIT REAL*8 (A-H, O-Z)
C
      PARAMETER (NMAX = 100)
      PARAMETER (NMAX_BAND = 12)
C
      REAL*8 ABND(NMAX_BAND, NMAX), A(NMAX, NMAX), 
     $       EIGEN(NMAX), SINGV(NMAX), 
     $       DUM(1,1), RWORK(5*NMAX)
C
      INTEGER LWORK, I, NX, NY, NXNY, IDUM
C
C ***** SETUP OUTPUT FILE
C
      IOUNIT = 7
      OPEN (UNIT=IOUNIT, FILE='OUTPUT.DATA', STATUS='UNKNOWN')
      REWIND IOUNIT
C
      WRITE (IOUNIT, 101)
C
C ***** SET UP THE PROBLEM SIZE
C
      NX = 4
      NY = 4
      NXNY = NX*NY
      WRITE (IOUNIT, 102) NXNY
C 
C ***** GENERATE THE MATRIX AND STORE THE UPPER TRIANGULAR PART IN 
C       THE LAPACK BANDED FORM.
C
      IBW = NX+1
      DO J = 1, NXNY
          DO I = 1, IBW
              ABND(I,J) = 0.0D0
          END DO
      END DO
C
      DO J = 1, NXNY
          ABND(IBW,J) = 4.0D0
      END DO
C
      DO J = 1, NY
          DO I = 2, NX
              K = (J-1) * NX + I
              ABND(IBW-1,K) = -1.0D0
          END DO
      END DO
C
      DO J = IBW, NXNY
          ABND(1,J) = -1.0D0
      END DO
C
C ***** CALL THE ROUTINE FOR EIGENVALUE DETERMINATION. 
C
      LDUM = 1
      CALL DSBEV ( 'N', 'U', NXNY, NX, ABND, NMAX_BAND,  
     $              EIGEN, DUM, LDUM, RWORK, INFO)
C
C ***** CHECK THAT THE ROUTINE COMPLETED WITHOUT ERROR
C
      IF (INFO.NE.0) THEN
          WRITE (IOUNIT, 103) INFO
          GO TO 999
      END IF
C
C ***** PRINT OUT THE EIGENVALUES
C
      WRITE (IOUNIT, 105)
      WRITE (IOUNIT, 106) (EIGEN(I), I=1, NXNY)
C 
C ***** GENERATE THE MATRIX AND STORE IT IN THE FULL FORM FOR THE 
C       SINGULAR VALUE DECOMPOSITION ROUTINE.
C
      DO J = 1, NXNY
          DO I = 1, NXNY
              A(I,J) = 0.0D0
          END DO
          A(J,J) = 4.0D0
      END DO
C
      DO J = 1, NY
          DO I = 2, NX
              K = (J-1) * NX + I
              A(K,K-1) = -1.0D0
              A(K-1,K) = -1.0D0
          END DO
      END DO
C
      DO J = NX+1, NXNY
          A(J-NX,J) = -1.0D0
          A(J,J-NX) = -1.0D0
      END DO
C
C ***** CALL THE ROUTINE TO OBTAIN THE SINGULAR VALUES.
C
      LWORK = 5 * NMAX
      CALL DGESVD ( 'N', 'N', NXNY, NXNY, A, NMAX,
     $              SINGV, DUM, LDUM, DUM, LDUM, RWORK, LWORK, INFO)
C
C ***** CHECK THAT THE ROUTINE COMPLETED WITHOUT ERROR
C
      IF (INFO.NE.0) THEN
          WRITE (IOUNIT, 104) INFO
          GO TO 999
      END IF
C
C ***** PRINT OUT THE SINGULAR VALUES
C
      WRITE (IOUNIT, 107)
      WRITE (IOUNIT, 106) (SINGV(I), I=1, NXNY)
C
C ***** GENERATE THE MAXIMUM ERROR BETWEEN THE SINGULAR VALUES AND THE 
C       EIGENVALUES. NOTE THAT THE MATRIX IS SYMMETRIC POSITIVE DEFINITE.
C
      ERROR = 0.0D0
      DO I = 1, NXNY
          ERROR = MAX ( ERROR, ABS(EIGEN(I)-SINGV(NXNY-I+1)) )
      ENDDO
C
      WRITE (IOUNIT, 108) ERROR
C
  999 CONTINUE
C
  101 FORMAT (/,2X, 'ILLUSTRATING THE USE OF LAPACK EIGEN-SOLVER ', 
     $              'AND SINGULAR VALUE DECOMPOSITION',
     $        /,2X, 'ROUTINES.',
     $        //,5X, 'TEST MATRIX IS SYMMETRIC POSITIVE DEFINITE',/)
  102 FORMAT (/,5X, 'ORDER OF THE MATRIX:',I5)
  103 FORMAT (/,9X, 'ERROR IN THE EIGENVALUE ROUTINE. INFO:',I10)
  104 FORMAT (/,9X, 'ERROR IN THE SVD ROUTINE. INFO:',I10)
  105 FORMAT (///,9X, 'DETERMINING THE EIGENVALUES, WITH THE MATRIX ', 
     $              'STORED IN A SYMMETRIC'
     $        /,9X, ' BANDED FORM',
     $        //,9X, 'EIGENVALUES IN ASCENDING ORDER:',/)
  106 FORMAT (/,12X, E15.8)
  107 FORMAT (///,9X, 'DETERMINING THE SINGULAR VALUES, WITH THE ', 
     $              'MATRIX STORED IN A FULL FORM',
     $        //,9X, 'SINGULAR VALUES IN DESCENDING ORDER:',/)
  108 FORMAT (///,5X, 'MAXIMUM ERROR BETWEEN THE SINGULAR VALUES AND ',
     $              'THE EIGENVALUES:',
     $        /,5X,E15.8)
C
      STOP
      END
