!------------------------------------------------------------------------------ !M+ ! NAME: ! lapack_utility ! ! PURPOSE: ! Module to provide Fortran 90 interfaces to some LAPACK F77 routines. ! (added to as additional LAPACK functionality is required.) ! ! CATEGORY: ! Math ! ! LANGUAGE: ! Fortran-90 ! ! CALLING SEQUENCE: ! USE lapack_utility ! ! MODULES: ! type_kinds: Module containing data type kind definitions. ! ! error_handler: Module to define error codes and handle error ! conditions ! ! CONTAINS: ! lapack_syev: Function to compute all eigenvalues and, optionally, ! eigenvectors of a real symmetric matrix. Is a wrapper ! for the LAPACK F77 SSYEV and DSYEV subroutines. ! ! EXTERNALS: ! LAPACK library routines. ! ! COMMON BLOCKS: ! None. ! ! SIDE EFFECTS: ! None known. ! ! RESTRICTIONS: ! None. ! ! INCLUDES: ! None. ! ! CREATION HISTORY: ! Written by: Paul van Delst, CIMSS/SSEC, 03-Apr-2000 ! paul.vandelst@ssec.wisc.edu ! ! Copyright (C) 2000 Paul van Delst ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation; either version 2 ! of the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program; if not, write to the Free Software ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. !M- !------------------------------------------------------------------------------ MODULE lapack_utility ! -------------------- ! Declare module usage ! -------------------- USE type_kinds USE error_handler ! ----------------------- ! Disable implicit typing ! ----------------------- IMPLICIT NONE ! ------------------ ! Explicit visibility ! ------------------- PRIVATE PUBLIC :: lapack_syev ! --------------------- ! Procedure overloading ! --------------------- INTERFACE lapack_syev MODULE PROCEDURE lapack_ssyev MODULE PROCEDURE lapack_dsyev END INTERFACE ! lapack_syev ! ---------------- ! Module variables ! ---------------- CHARACTER( 128 ) :: message CONTAINS !-------------------------------------------------------------------------------- !S+ ! NAME: ! lapack_syev ! ! PURPOSE: ! Function to compute all eigenvalues and, optionally, eigenvectors of a ! real symmetric matrix. Is a wrapper for the LAPACK F77 SSYEV and DSYEV ! subroutines. ! ! CATEGORY: ! Math ! ! LANGUAGE: ! Fortran 90 ! ! CALLING SEQUENCE: ! result = lapack_syev( jobz, & ! Input ! uplo, & ! Input ! a, & ! In/Output ! w, & ! Output ! info, & ! Output ! direction = direction, & ! Optional input ! message_log = message_log ) & ! Optional Input ! ! INPUT ARGUMENTS: ! jobz: Single character describing the action required. ! If = 'N': Compute eigenvalues only; ! = 'V': Compute eigenvalues and eigenvectors. ! UNITS: None ! TYPE: CHARACTER( 1 ) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT( IN ) ! ! uplo: Single character describing which triangle of the ! input matrix is stored. ! If = 'U': Upper triangle of input matrix A is stored; ! = 'L': Lower triangle of input matrix A is stored. ! UNITS: None ! TYPE: CHARACTER( 1 ) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT( IN ) ! ! a: Floating point symmetric array for which the ! eigenvalues/vectors are required. ! If UPLO = 'U', the leading N-by-N upper triangular ! part of A contains the upper triangular part of the ! matrix A. If UPLO = 'L', the leading N-by-N lower ! triangular part of A contains the lower triangular ! part of the matrix A. On exit, if JOBZ = 'V', then ! if INFO = 0, A contains the orthonormal eigenvectors ! of the matrix A. If JOBZ = 'N', then on exit the lower ! triangle (if UPLO='L') or the upper triangle ! (if UPLO='U') of A, including the diagonal, is destroyed. ! NOTE: Regardless of values of JOBZ and UPLO, the contents ! of A are overwritten upon output. ! UNITS: Depends ! TYPE: REAL( fp_kind ) ! DIMENSION: M x N, where M > or = MAX( 1,N ) ! ATTRIBUTES: INTENT( IN OUT ) ! ! OPTIONAL INPUT ARGUMENTS: ! direction: Optional scalar argument used to determine the ordering ! of the output eigenvalues/vectors. ! If = 0, eigenvalues (and corresponding eigenvectors) are ! in ASCENDING order. ! = 1, eigenvalues (and corresponding eigenvectors) are ! in DESCENDING order. ! The default action, i.e. if the argument is not specified, ! is to return the eigenvalues in ASCENDING order. ! UNITS: None ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT( IN OUT ) ! ! message_log: Character string specifying a filename in which any ! messages will be logged. If not specified, or if an ! error occurs opening the log file, the default action ! is to output messages to the screen. ! UNITS: None ! TYPE: Character ! DIMENSION: Scalar ! ATTRIBUTES: OPTIONAL, INTENT( IN ) ! ! OUTPUT ARGUMENTS: ! w: Floating point vector continaing the eigenvalues. ! UNITS: Depends ! TYPE: REAL( fp_kind ) ! DIMENSION: N ! ATTRIBUTES: INTENT( OUT ) ! ! info: Integer indicating if decomposition within the LAPACK ! routines was successful. ! If = 0: successful exit ! < 0: if INFO = -i, the i-th argument (to the LAPACK ! routine) had an illegal value ! > 0: if INFO = i, the algorithm failed to converge; i ! off-diagonal elements of an intermediate tridiagonal ! form did not converge to zero. ! UNITS: None ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT( OUT ) ! ! OPTIONAL OUTPUT ARGUMENTS: ! None ! ! FUNCTION RESULT ! Function returns an error flag ! If result = SUCCESS, everything is o.k. ! = FAILURE, if an error occurred ! ! CALLS: ! SSYEV or DSYEV: Single and Double LAPACK subroutine to compute ! all eigenvalues and, optionally, eigenvectors of ! a real symmetric matrix. The functions are ! overloaded in this subprogram. ! SOURCE: LAPACK library. ! ! display_message: Subroutine to output messages ! SOURCE: error_handler module ! ! EXTERNALS: ! LAPACK library routines. ! ! COMMON BLOCKS: ! None. ! ! SIDE EFFECTS: ! Input matrix, A, is overwritten on output. ! ! RESTRICTIONS: ! See LAPACK routine SSYEV/DSYEV.. ! ! PROCEDURE: ! The following procedure details what occurs within the LAPACK ! SSYEV/DSYEV routines, not this one. ! ! The input symmetric/Hermitian matrix, A, is first reduced to a real ! symmetric tridiagonal form, T, by an orthogonal/unitary similarity ! transformation. This is done via the LAPACK routines SSYTRD/DSYTRD. ! ! If only eigenvalues are required, they are calculated using the ! Pal-Walker-Kahan variant of the QL or QR algorithm. This is done ! via the LAPACK routines SSTERF/DSTERF. ! ! If both eigenvalues and vectors are required, first the orthogonal ! matrix Q is determined, which is defined as the product of n-1 ! elementary reflectors of order N, as returned by SSYTRD/DSYTRD. This ! is done via the LAPACK routines SORGTR/DORGTR. Second, the eigenvalues ! and vectors are calculated for the symmetric tridiagonal matrix using ! the generated orthogonal matrix by the implicit QL or QR method. This ! is done via the LAPACK routines SSTEQR/DSTEQR. !S- !-------------------------------------------------------------------------------- FUNCTION lapack_ssyev( jobz, & ! Input uplo, & ! Input a, & ! In/Output w, & ! Output info, & ! Output direction, & ! Optional input message_log ) & ! Optional Input RESULT( error_status ) !#--------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#--------------------------------------------------------------------------# ! --------- ! Arguments ! --------- CHARACTER( 1 ), INTENT( IN ) :: jobz CHARACTER( 1 ), INTENT( IN ) :: uplo REAL( Single ), DIMENSION( :, : ), INTENT( IN OUT ) :: a REAL( Single ), DIMENSION( : ), INTENT( OUT ) :: w INTEGER, INTENT( OUT ) :: info INTEGER, OPTIONAL, INTENT( IN ) :: direction CHARACTER( * ), OPTIONAL, INTENT( IN ) :: message_log ! --------------- ! Function result ! --------------- INTEGER :: error_status ! ---------------- ! Local parameters ! ---------------- CHARACTER( * ), PARAMETER :: ROUTINE_NAME = 'f90_ssyev' CHARACTER( * ), PARAMETER, DIMENSION( 8 ) :: SSYEV_ARGS = & (/ 'JOBZ ', 'UPLO ', 'N ', 'A ', 'LDA ', 'W ', 'WORK ', 'LWORK' /) ! --------------- ! Local variables ! --------------- INTEGER :: allocate_status INTEGER :: i INTEGER :: n ! Order of matrix A INTEGER :: lda ! The leading dimension of A, MAX(1,N) INTEGER :: lwork ! The length of the work array REAL( Single ), DIMENSION( : ), ALLOCATABLE :: work REAL( Single ), DIMENSION( 2 ) :: dummy_work ! -------------------------- ! Externals (lapack library) ! -------------------------- EXTERNAL ssyev ! ---------- ! Intrinsics ! ---------- INTRINSIC SIZE, & TRIM write(*,*) 'Entered lapack_ssyev' !#--------------------------------------------------------------------------# !# -- GET DIMENSIONS AND CALL SSYEV TO DETERMINE OPTIMUM WORK ARRAY SIZE -- # !#--------------------------------------------------------------------------# ! -- Dimensions of A n = SIZE( a, DIM = 2 ) lda = SIZE( a, DIM = 1 ) ! -- Call SSYEV to determine optimum LWORK lwork = -1 CALL ssyev( jobz, uplo, n, a, lda, w, dummy_work, lwork, info ) write(*,*) 'Finished first ssyev in lapack ssyev' IF ( info /= 0 ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Illegal value for SSYEV argument '//& TRIM( SSYEV_ARGS( ABS( info ) ) ), & error_status, & message_log = message_log ) RETURN END IF ! -- Assign optimum LWORK lwork = dummy_work( 1 ) !#--------------------------------------------------------------------------# !# -- ALLOCATE WORK ARRAY -- # !#--------------------------------------------------------------------------# ALLOCATE( work( lwork ), & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN error_status = FAILURE WRITE( message, '( "Error allocating WORK array. STAT = ", i5 )' ) & allocate_status CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- CALL THE LAPACK ROUTINE SSYEV -- # !#--------------------------------------------------------------------------# CALL ssyev( jobz, uplo, n, a, lda, w, work, lwork, info ) !#--------------------------------------------------------------------------# !# -- CHECK FOR ANY ERRORS -- # !#--------------------------------------------------------------------------# ! -- Set successful return flag error_status = SUCCESS ! -- Check info output IF ( info /= 0 ) THEN error_status = WARNING WRITE( message, '( i10 )' ) info message = TRIM( message )//' off-diagonal elements of an intermediate '//& 'tridiagonal form did not converge to zero.' CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) END IF !#--------------------------------------------------------------------------# !# -- REVERSE ORDER OF ARRAYS IF REQUIRED -- # !#--------------------------------------------------------------------------# IF ( PRESENT( direction ) ) THEN IF ( direction == 1 ) THEN ! ---------------------------------------------------------- ! Reverse order of eigenVECTORS from ascending -> descending ! ---------------------------------------------------------- ! -- Is work array large enough? IF ( lwork < lda ) THEN DEALLOCATE( work ) ALLOCATE( work( lda ), & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN error_status = WARNING WRITE( message, '( "Error allocating eigenvector swap array - No swap done. STAT = ", i5 )' ) & allocate_status CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF END IF ! -- Swap the eigenvectors DO i = 1, n/2 work = a( :, i ) a( :, i ) = a( :, n-i+1 ) a( :, n-i+1 ) = work END DO ! ---------------------------------------------------------- ! Reverse order of eigenVALUES from ascending -> descending ! ---------------------------------------------------------- work( 1:n ) = w DO i = 1, n w( i ) = work( n - i + 1 ) END DO END IF END IF !#--------------------------------------------------------------------------# !# -- DONE -- # !#--------------------------------------------------------------------------# DEALLOCATE( work ) END FUNCTION lapack_ssyev FUNCTION lapack_dsyev( jobz, & ! Input uplo, & ! Input a, & ! In/Output w, & ! Output info, & ! Output direction, & ! Optional input message_log ) & ! Optional Input RESULT( error_status ) !#--------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#--------------------------------------------------------------------------# ! --------- ! Arguments ! --------- CHARACTER( 1 ), INTENT( IN ) :: jobz CHARACTER( 1 ), INTENT( IN ) :: uplo REAL( Double ), DIMENSION( :, : ), INTENT( IN OUT ) :: a REAL( Double ), DIMENSION( : ), INTENT( OUT ) :: w INTEGER, INTENT( OUT ) :: info INTEGER, OPTIONAL, INTENT( IN ) :: direction CHARACTER( * ), OPTIONAL, INTENT( IN ) :: message_log ! --------------- ! Function result ! --------------- INTEGER :: error_status ! ---------------- ! Local parameters ! ---------------- CHARACTER( * ), PARAMETER :: ROUTINE_NAME = 'F90_DSYEV' CHARACTER( * ), PARAMETER, DIMENSION( 8 ) :: DSYEV_ARGS = & (/ 'JOBZ ', 'UPLO ', 'N ', 'A ', 'LDA ', 'W ', 'WORK ', 'LWORK' /) ! --------------- ! Local variables ! --------------- INTEGER :: allocate_status INTEGER :: i INTEGER :: n ! Order of matrix A INTEGER :: lda ! The leading dimension of A, MAX(1,N) INTEGER :: lwork ! The length of the work array REAL( Double ), DIMENSION( : ), ALLOCATABLE :: work REAL( Double ), DIMENSION( 2 ) :: dummy_work ! -------------------------- ! Externals (lapack library) ! -------------------------- EXTERNAL dsyev ! ---------- ! Intrinsics ! ---------- INTRINSIC SIZE, & TRIM !#--------------------------------------------------------------------------# !# -- GET DIMENSIONS AND CALL DSYEV TO DETERMINE OPTIMUM WORK ARRAY SIZE -- # !#--------------------------------------------------------------------------# ! -- Dimensions of A n = SIZE( a, DIM = 2 ) lda = SIZE( a, DIM = 1 ) ! -- Call DSYEV to determine optimum LWORK lwork = -1 CALL dsyev( jobz, uplo, n, a, lda, w, dummy_work, lwork, info ) IF ( info /= 0 ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Illegal value for DSYEV argument '//& TRIM( DSYEV_ARGS( ABS( info ) ) ), & error_status, & message_log = message_log ) RETURN END IF ! -- Assign optimum LWORK lwork = dummy_work( 1 ) !#--------------------------------------------------------------------------# !# -- ALLOCATE WORK ARRAY -- # !#--------------------------------------------------------------------------# ALLOCATE( work( lwork ), & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN error_status = FAILURE WRITE( message, '( "Error allocating WORK array. STAT = ", i5 )' ) & allocate_status CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- CALL THE LAPACK ROUTINE DSYEV -- # !#--------------------------------------------------------------------------# CALL dsyev( jobz, uplo, n, a, lda, w, work, lwork, info ) !#--------------------------------------------------------------------------# !# -- CHECK FOR ANY ERRORS -- # !#--------------------------------------------------------------------------# ! -- Set successful return flag error_status = SUCCESS ! -- Check info output IF ( info /= 0 ) THEN error_status = WARNING WRITE( message, '( i10 )' ) info message = TRIM( message )//' off-diagonal elements of an intermediate '//& 'tridiagonal form did not converge to zero.' CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) END IF !#--------------------------------------------------------------------------# !# -- REVERSE ORDER OF ARRAYS IF REQUIRED -- # !#--------------------------------------------------------------------------# IF ( PRESENT( direction ) ) THEN IF ( direction == 1 ) THEN ! ---------------------------------------------------------- ! Reverse order of eigenVECTORS from ascending -> descending ! ---------------------------------------------------------- ! -- Is work array large enough? IF ( lwork < lda ) THEN DEALLOCATE( work ) ALLOCATE( work( lda ), & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN error_status = WARNING WRITE( message, '( "Error allocating eigenvector swap array - No swap done. STAT = ", i5 )' ) & allocate_status CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF END IF ! -- Swap the eigenvectors DO i = 1, n/2 work( 1:lda ) = a( :, i ) a( :, i ) = a( :, n-i+1 ) a( :, n-i+1 ) = work( 1:lda ) END DO ! ---------------------------------------------------------- ! Reverse order of eigenVALUES from ascending -> descending ! ---------------------------------------------------------- work( 1:n ) = w DO i = 1, n w( i ) = work( n - i + 1 ) END DO END IF END IF !#--------------------------------------------------------------------------# !# -- DONE -- # !#--------------------------------------------------------------------------# DEALLOCATE( work ) END FUNCTION lapack_dsyev END MODULE lapack_utility !------------------------------------------------------------------------------- ! -- MODIFICATION HISTORY -- !------------------------------------------------------------------------------- ! ! $Id: lapack_utility.f90,v 1.2 2001/09/28 16:57:14 paulv Exp $ ! ! $Date: 2001/09/28 16:57:14 $ ! ! $Revision: 1.2 $ ! ! $Name: Ozone_Match_1-1 $ ! ! $State: Exp $ ! ! $Log: lapack_utility.f90,v $ ! Revision 1.2 2001/09/28 16:57:14 paulv ! - Added suprogram header documentation. ! ! Revision 1.1 2001/09/26 17:07:13 paulv ! - Renamed module from LAPACK. ! - Retooled SYEV functions to eliminate the need for the work array. ! ! !