!-------------------------------------------------------------------------------- !M+ ! NAME: ! coordinate_space ! ! PURPOSE: ! Module to decompose a correlation matrix computed from an input ! matrix of vectors into its eigenvectors to allow a user to input ! single vectors and search for the closest match in eigenspace. ! ! CATEGORY: ! Math ! ! LANGUAGE: ! Fortran 90 ! ! CALLING SEQUENCE: ! USE coordinate_space ! ! LANGUAGE: ! Fortran-90 ! ! MODULES: ! type_kinds: Module containing data type kind definitions. ! ! error_handler: Module to define error codes and handle error ! conditions ! ! lapack_utility: Module containing overloaded definitions of the ! various LAPACK routines used. ! ! CONTAINS: ! initialize_coordinate_space: PUBLIC Function to set up and calculate ! the coordinate space of an input matrix. ! ! find_closest_vector: PUBLIC Function to find the closest vector ! in coordinate space of the input ! matrix to the user specified vector. ! ! destroy_coordinate_space: PUBLIC Function to deallocate all the ! matrices used to define the coordinate ! space. ! ! compute_vector_coordinate: PRIVATE function to compute the coordinate ! of the input vector in eignespace. ! ! EXTERNALS: ! None ! ! COMMON BLOCKS: ! None. ! ! SIDE EFFECTS: ! None known. ! ! RESTRICTIONS: ! None. ! ! EXAMPLE: ! Given a MxN matrix of atmospheric temperature profile, temperature, and ! a single temperature profile vector, test_t, of length N, find the closest ! profile in the MxN matrix: ! ! ---------------------------------------------- ! Create coordinate space. Variance cutoff=99.9% ! ---------------------------------------------- ! ! error_status = initialize_coordinate_space( temperature , & ! 99.9, & ! debug = 1 ) ! ! IF ( error_status /= SUCCESS ) THEN ! WRITE( *, '( /5x, "Error calculating coordinate space." )' ) ! STOP ! END IF ! ! ! ------------------------------ ! Match the temperature profiles ! ------------------------------ ! ! error_status = find_closest_vector( test_t, & ! profile_number ) ! ! IF ( error_status /= SUCCESS ) THEN ! WRITE( *, '( /5x, "Error matching temperature vector." )' ) ! STOP ! END IF ! ! WRITE( *, '( /5x, "Matched temperature vector : ", i4 )' ) profile_number ! ! ! ------------------------ ! Cleanup coordinate space ! ------------------------ ! ! error_status = destroy_coordinate_space( ) ! ! IF ( error_status /= SUCCESS ) THEN ! WRITE( *, '( /5x, "Error cleaning up coordinate space." )' ) ! STOP ! END IF ! ! CREATION HISTORY: ! Written by: Paul van Delst, CIMSS/SSEC, 25-Mar-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 coordinate_space ! -------------------- ! Declare modules used ! -------------------- USE type_kinds, ONLY: fp_kind USE error_handler USE lapack_utility ! ----------------------- ! Disable implicit typing ! ----------------------- IMPLICIT NONE ! ------------ ! Visibilities ! ------------ PRIVATE PUBLIC :: initialize_coordinate_space PUBLIC :: find_closest_vector PUBLIC :: destroy_coordinate_space ! ------------------------ ! Module scope shared data ! ------------------------ ! -- Average vector of input MxN matrix, Nx1 REAL( fp_kind ), ALLOCATABLE, DIMENSION( : ) :: average ! -- Vector norm of input MxN matrix, Nx1 REAL( fp_kind ), ALLOCATABLE, DIMENSION( : ) :: vector_norm ! -- Eigenvector analysis arrays for decomposition ! -- of correlation matrix of input matrix ! -- Eigenvalues, Nx1 REAL( fp_kind ), ALLOCATABLE, DIMENSION( : ) :: eigenvalue ! -- Eigenvectors, NxN REAL( fp_kind ), ALLOCATABLE, DIMENSION( :, : ) :: eigenvector ! -- Coordinate space matrix, MxK REAL( fp_kind ), ALLOCATABLE, DIMENSION( :, : ) :: coordinate ! -- The number of eigenvectors retained, K INTEGER :: n_vectors ! -- Logical variable for flagging coordinate ! -- space initialisation LOGICAL :: space_initialized = .FALSE. ! -- Error message for internal writes CHARACTER( 128 ) :: message CONTAINS !-------------------------------------------------------------------------------- !S+ ! NAME: ! initialize_coordinate_space ! ! PURPOSE: ! Function to set up and calculate the coordinate space of an input matrix ! ! CATEGORY: ! Math ! ! LANGUAGE: ! Fortran 90 ! ! CALLING SEQUENCE: ! result = initialize_coordinate_space( matrix, & ! explained_variance = explained_variance, & ! debug = debug & ! message_log = message_log ) ! ! INPUT ARGUMENTS: ! matrix: Input matrix of M instances of vectors ! of length N. ! UNITS: Any ! TYPE: REAL( fp_kind ) ! DIMENSION: M x N ! ATTRIBUTES: INTENT( IN ) ! ! OPTIONAL INPUT ARGUMENTS: ! explained_variance: Tolerance value to select how many ! dimensions or modes (i.e. eigenvectors) ! will be retained in the coordinate space. ! If not specified, the default value is 99.9%. ! UNITS: Percent ! TYPE: REAL( fp_kind ) ! DIMENSION: Scalar ! ATTRIBUTES: OPTIONAL, INTENT( IN ) ! ! debug: Optional scalar argument used to output ! debug information about the computation. ! if debug = 1. If set, debug information ! is written to the standard output (i.e. ! WRITE( *, ... ). The default is no debug output. ! UNITS: None ! TYPE: INTEGER ! DIMENSION: Scalar ! ! 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: ! None ! ! OPTIONAL OUTPUT ARGUMENTS: ! None ! ! FUNCTION RESULT ! Function returns an error flag ! If result = SUCCESS, everything is o.k. ! = FAILURE, if an error occurred ! ! CALLS: ! lapack_syev: LAPACK subroutine to compute all eigenvalues and, ! optionally, eigenvectors of a real symmetric or ! Hermitian matrix. Overloaded to SSYEV and DSYEV ! from the LAPACK F77 distribution. ! SOURCE: lapack module ! ! display_message: Subroutine to output messages ! SOURCE: error_handler module ! ! EXTERNALS: ! None ! ! COMMON BLOCKS: ! None. ! ! SIDE EFFECTS: ! None known. ! ! RESTRICTIONS: ! None. ! ! PROCEDURE: ! Given an input MxN matrix, A, the correlation matrix is determined ! by the following steps: ! ! 1) Calculate the average of the M input vectors of length N, ! ! _ M ! A(n) = sum{ A(i,n) } ! i=1 ! ------------- ! M ! ! 2) Calculate the mean deviation, ! ! dA(m,n) = A(m,n) - A(n) ! ! 3) Calculate the matrix vector norms for every n = 1,2,...,N, ! ________________ ! ^ / M { 2 } ! k(n) = / sum{ dA(i,n) } ! \/ i=1{ } ! ! 4) Normalise the mean deviation vectors, ! ! ~ dA(m,n) ! A(m,n) = --------- ! ^ ! k(n) ! ! 5) Compute the correlation matrix ! ! ~ T ~ ! P = [A] * [A] ! ! 6) Compute the eigenvalues, e, and eigenvectors, E, of the ! correlation matrix, P. ! ! 7) Determine the K eigenvectors required to explain the ! required variance. ! ! 8) Construct the matrix of coordinates in eigenspace, ! ! N ~ ! Y(m,k) = sum( A(m,j) * E(j,k) ) m=1,2,...,M; k=1,2,...,K ! j=1 ! ! where Y(m,k) = MxK matrix of coordinates for the M input ! vectors of length N elements, and ! K = the number of eigenvectors (modes) retained ! for the expansion based upon the explained ! variance tolerance. If all modes are retained ! then K = N. !S- !-------------------------------------------------------------------------------- FUNCTION initialize_coordinate_space ( matrix, & ! Input explained_variance, & ! Optional input debug, & ! Optional input message_log ) & ! Optional input RESULT ( error_status ) !#--------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#--------------------------------------------------------------------------# ! --------- ! Arguments ! --------- REAL( fp_kind ), DIMENSION( :, : ), INTENT( IN ) :: matrix REAL( fp_kind ), OPTIONAL, INTENT( IN ) :: explained_variance INTEGER, OPTIONAL, INTENT( IN ) :: debug CHARACTER( * ), OPTIONAL, INTENT( IN ) :: message_log ! ------ ! Result ! ------ INTEGER :: error_status ! ---------- ! Parameters ! ---------- ! -- Default explained variance tolerance, in %, used to ! -- determine the number of eigenvectors to retain REAL( fp_kind ), PARAMETER :: DEFAULT_VARIANCE_TOLERANCE = 99.9_fp_kind ! -- Profile name CHARACTER( * ), PARAMETER :: ROUTINE_NAME = 'INITIALIZE_COORDINATE_SPACE' ! --------------- ! Local variables ! --------------- ! -- Input matrix dimensions INTEGER :: M, N ! n==nlev, m==nprof ! -- Unit vector matrix, MxN REAL( fp_kind ), DIMENSION( SIZE( matrix, DIM = 1 ), & SIZE( matrix, DIM = 2 ) ) :: unit_vector ! -- Correlation matrix, NxN REAL( fp_kind ), DIMENSION( SIZE( matrix, DIM = 2 ), & SIZE( matrix, DIM = 2 ) ) :: correlation ! -- Eigenvalue sum variables REAL( fp_kind ) :: eigenvalue_sum REAL( fp_kind ) :: eigenvalue_sum_cutoff REAL( fp_kind ) :: variance_tolerance ! -- Error check variables INTEGER :: allocate_status INTEGER :: lapack_status ! -- Loop variables INTEGER :: i, j ! ---------- ! Intrinsics ! ---------- INTRINSIC ABS, & ALLOCATED, & MATMUL, & PRESENT, & REAL, & RESHAPE, & SIZE, & SQRT, & SUM, & TRANSPOSE, & TRIM !#--------------------------------------------------------------------------# !# -- HAS THE COORDINATE SPACE BEEN INITIALIZED? -- # !#--------------------------------------------------------------------------# IF ( space_initialized ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Coordinate space already established.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- GET INPUT MATRIX DIMENSIONS -- # !#--------------------------------------------------------------------------# M = SIZE( matrix, 1 ) N = SIZE( matrix, 2 ) !#--------------------------------------------------------------------------# !# -- ASSIGN THE VARIANCE TOLERANCE -- # !#--------------------------------------------------------------------------# IF ( PRESENT( explained_variance ) ) THEN variance_tolerance = explained_variance ELSE variance_tolerance = DEFAULT_VARIANCE_TOLERANCE END IF !#--------------------------------------------------------------------------# !# -- ALLOCATE ALL THE REQUIRED ARRAYS -- # !#--------------------------------------------------------------------------# ALLOCATE( average( N ), & vector_norm( N ), & eigenvalue( N ), & eigenvector( N, N ), & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error allocating arrays. STAT = ", i5 )' ) & allocate_status error_status = FAILURE CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- CALCULATE AVERAGE VALUE ALONG M DIMENSION -- # !# # !# Calculate average vector: # !# # !# _ M # !# A(n) = sum{ A(i,n) } # !# i=1 # !# ------------- # !# M # !# So: # !# 1) Sum across the #1 (M) dimension # !# 2) Divide by the M dimension # !#--------------------------------------------------------------------------# average(:) = SUM( matrix(:,:), DIM = 1 ) / & ! --------------------------- REAL( M, fp_kind ) !#--------------------------------------------------------------------------# !# -- CALCULATE MEAN DEVIATION MATRIX -- # !# _ # !# Calculate mean deviation, dA(m,n) = A(m,n) - A(n) # !# # !# The RESHAPE intrinsic repeats the average vector M times creating a # !# matrix the same rank and size as the input matrix. # !# # !# Save the result in the unit_vector array. # !#--------------------------------------------------------------------------# unit_vector = matrix - TRANSPOSE( RESHAPE( average, & (/ N, M /), & average ) ) !#--------------------------------------------------------------------------# !# -- CALCULATE THE MATRIX VECTOR NORMS FOR EVERY n = 1,2,...,N -- # !# # !# Calculate vector norms # !# ________________ # !# ^ / M { 2 } # !# k(n) = / sum{ dA(i,n) } # !# \/ i=1{ } # !# # !# Remember that the mean deviations were saved in the unit_vector array. # !# These arrays can get pretty big, particularly if double precision is # !# used, so may as well conserve stack space. # !#--------------------------------------------------------------------------# vector_norm = SQRT( SUM( unit_vector**2, DIM = 1 ) ) !#--------------------------------------------------------------------------# !# -- NORMALISE THE MEAN DEVIATION VECTORS -- # !# # !# Do the normalisation # !# # !# ~ dA(m,n) # !# A(m,n) = --------- # !# ^ # !# k(n) # !# # !# Remember that the mean deviations were saved in the unit_vector array. # !# These arrays can get pretty big, particularly if double precision is # !# used, so may as well conserve stack space. # !#--------------------------------------------------------------------------# unit_vector = unit_vector / TRANSPOSE( RESHAPE( vector_norm, & (/ N, M /), & vector_norm ) ) !#--------------------------------------------------------------------------# !# -- CALCULATE THE CORRELATION MATRIX -- # !# # !# ~ T ~ # !# P = [A] * [A] # !#--------------------------------------------------------------------------# correlation = MATMUL( TRANSPOSE( unit_vector ), & unit_vector ) !#--------------------------------------------------------------------------# !# -- DECOMPOSE THE CORRELATION MATRIX. IT IS A REAL SYMMETRIC -- # !# -- MATRIX SO ALL EIGENVALUES ARE REAL. -- # !# # !# -- **NOTE** -- # !# LAPACK routine S/DSYEV outputs eigenvalues in ASCENDING order # !#--------------------------------------------------------------------------# ! ------------------------------------------------------- ! Call LAPACK routine to determine the eigenvalue vector, ! e, and the eigenvector matrix, E, of the correlation ! matrix, P. ! ------------------------------------------------------- ! -- Assign the correlation matrix to the eigenvector ! -- matrix. LAPACK is lame in that it overwrites input. eigenvector = correlation ! -- Call the LAPACK function error_status = lapack_syev( 'V', & ! INPUT: Compute both eigenvales and vectors 'U', & ! INPUT: Input matrix is stored as upper triangle eigenvector, & ! IN/OUTPUT: Input matrix, output eigenvectors eigenvalue, & ! OUTPUT: Eigenvalues in ascending order lapack_status, & ! OUTPUT: = 0 for normal completion. ! =-i the i-th argument had an illegal value ! =+i the algorithm failed to converge. i ! off-diagonal elements of an intermediate ! tridiagonal form did not converge to zero. direction = 1, & ! Swap to DESCENDING order (default is ASCENDING) message_log = message_log ) ! --------------------------- ! Check for normal completion ! --------------------------- IF ( error_status == FAILURE ) THEN CALL display_message( ROUTINE_NAME, & 'LAPACK eigenvector decomposition error.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- DETERMINE % VARIANCES EXPLAINED -- # !#--------------------------------------------------------------------------# ! ------------------------------------------ ! Determine the % variance explained for set ! tolerance (variance_tolerance) ! ------------------------------------------ eigenvalue_sum_cutoff = SUM( eigenvalue ) * variance_tolerance / 100.0 eigenvalue_loop: DO n_vectors = 1, N ! -- Sum eigenvalues eigenvalue_sum = SUM( eigenvalue( 1 : n_vectors ) ) ! -- Has the required % variance been explained? IF ( eigenvalue_sum > eigenvalue_sum_cutoff ) EXIT eigenvalue_loop END DO eigenvalue_loop ! -- Check if loop went through ALL eigenvalues IF ( n_vectors > N ) n_vectors = N !#--------------------------------------------------------------------------# !# -- CHECK THAT THE CALCULATION WORKED CORRECTLY -- # !# # !# The correlation matrix is given by: # !# # !# ~ T ~ # !# P = [A] * [A] # !# # !# T # !# = [E] * e * [E] # !# # !# Thus, # !# # !# T # !# [E] * P * [E] = e # !# # !# if everything worked. # !# # !# NOTE: If debug code in this routine is removed, then the eigenvector # !# matrix can be used to store the correlations rather than define # !# its own (potentially large) space. The code below is the only # !# place where the correlations and eigenvectors need to be accessed # !# at he same time. Remember that the LAPACK routine used overwrites # !# the input matrix (the correlations) with the eigenvectors. # !#--------------------------------------------------------------------------# decompose_debug: IF ( PRESENT( debug ) ) THEN IF ( debug == 1 ) THEN WRITE( *, '( /5x, "Calculating [E]T * P * [E]..." )' ) ! ------------------ ! Do the matrix math ! ------------------ correlation = MATMUL( TRANSPOSE( eigenvector ), & MATMUL( correlation, eigenvector ) ) ! --------------------------------------------------- ! Output diagonal and off diagonal elements to screen ! --------------------------------------------------- WRITE( *, '( /10x, "Calculated eigenvalues:" )' ) WRITE( *, '( 8( 1x,es13.6 ) )' ) eigenvalue WRITE( *, '( /10x, "Check diagonal (should be same as eigenvalues) : " )' ) WRITE( *, '( 8( 1x,es13.6 ) )' ) ( correlation( i, i ), i = 1, N ) WRITE( *, '( /10x, "Difference (should be zero or less than numerical precision) : " )' ) WRITE( *, '( 8( 1x,es13.6 ) )' ) ( eigenvalue( i ) - correlation( i, i ), i = 1, N ) WRITE( *, '( /10x, "Reverse diagonal (should be zero or less than numerical precision) : " )' ) WRITE( *, '( 8( 1x,es13.6 ) )' ) ( correlation( N-i+1, i ), i = 1, N ) ! -------------------------------------- ! Output number of eigenvectors retained ! -------------------------------------- WRITE( *, '( /10x, "No. of eigenvectors retained: ", i5, " of ", i5 )' ) n_vectors, N END IF END IF decompose_debug !#--------------------------------------------------------------------------# !# -- CALCULATE MATRIX OF CO-ORDINATES -- # !# # !# The matrix of profile coordinates in eigenspace is determined using: # !# # !# N ~ # !# Y(m,k) = sum( A(m,j) * E(j,k) ) m=1,2,...,M; k=1,2,...,K # !# j=1 # !# # !# where Y(m,k) = MxK matrix of coordinates for the M input vectors of # !# length N elements (hence the MxN input matrix), and # !# K = the number of eigenvectors (modes) retained for the # !# expansionbased upon the explained variance tolerance. If # !# all modes are retained then K = N. # !#--------------------------------------------------------------------------# ! -------------------------- ! Allocate coordinate matrix ! -------------------------- ALLOCATE( coordinate( M, n_vectors ), & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Error allocating COORDINATE matrix (MxK).', & error_status, & message_log = message_log ) RETURN END IF ! -------------------------- ! Fill the coordinate matrix ! -------------------------- coordinate = MATMUL( unit_vector( :, : ), & eigenvector( :, 1 : n_vectors ) ) !#--------------------------------------------------------------------------# !# -- NO PROBLEMS HERE -- # !#--------------------------------------------------------------------------# space_initialized = .TRUE. error_status = SUCCESS END FUNCTION initialize_coordinate_space !-------------------------------------------------------------------------------- !S+ ! NAME: ! find_closest_vector ! ! PURPOSE: ! Function to find the closest vector in coordinate space of the input ! matrix to the user specified vector. ! ! CATEGORY: ! Math ! ! LANGUAGE: ! Fortran 90 ! ! CALLING SEQUENCE: ! result = find_closest_vector( vector, & ! vector_number, & ! message_log = message_log ) ! ! INPUT ARGUMENTS: ! vector: Input vector for which a match is required ! UNITS: Same as those of matrix used ! to create coordinate space ! TYPE: REAL( fp_kind ) ! DIMENSION: N ! ATTRIBUTES: INTENT( IN ) ! ! OPTIONAL INPUT ARGUMENTS: ! 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: ! vector_number: The number of the vector in the original ! input matrix that is closest to the input ! vector in eigenspace. ! UNITS: None. ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT( OUT ) ! ! FUNCTION RESULT ! Function returns an error flag ! If result = SUCCESS, everything is o.k. ! = FAILURE, if an error occurred ! ! CALLS: ! compute_vector_coordinate: PRIVATE function to compute the coordinate ! of the input vector in eignespace. ! ! display_message: Subroutine to output messages ! SOURCE: error_handler module ! ! EXTERNALS: ! None ! ! COMMON BLOCKS: ! None. ! ! SIDE EFFECTS: ! None known. ! ! RESTRICTIONS: ! None. ! ! PROCEDURE: ! The closest vector is found by finding minimising the sum of ! the squares of the coordinate differences. The coordinate ! difference matrix is calculated, ! ! coord_difference = coordinate matrix - replicated vector coordinate ! ! 1 ---------> M ! [ c11 c21 c31 ... cM1 ] [ v1 v1 v1 v1 ] ! [ c12 ... ] [ v2 v2 v2 v2 ] ! = [ c13 ... ] - [ v3 v3 v3 .. v3 ] ! [ ... ... ] [ .. .. .. .. ] ! [ c1K cMK ] [ vK vK vK vK ] ! ! [ d11 d21 d31 ... dM1 ] ! [ d12 ... ] ! = [ d13 ... ] ! [ ... ... ] ! [ d1K dMK ] ! ! This array is squared to give the squared difference array ! which is summed along the K dimension to produce a Mx1 vector: ! ! [ K K K ] ! sum_of_squared_differences = [ sum{d1i^2}, sum{d2i^2},... sum{dMi^2} ] ! [ i=1 i=1 i=1 ] ! ! The index of the minimum value of this vector, say the j'th value, ! means that the j'th vector of the original input matrix, (j,1->M), ! is the vector that is the closest in eigenspace to the user input vector. !S- !-------------------------------------------------------------------------------- FUNCTION find_closest_vector ( vector, & ! Input vector_number, & ! Output message_log ) & ! Optional input RESULT ( error_status ) !#--------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#--------------------------------------------------------------------------# ! --------- ! Arguments ! --------- REAL( fp_kind ), DIMENSION( : ), INTENT( IN ) :: vector INTEGER, INTENT( OUT ) :: vector_number CHARACTER( * ), OPTIONAL, INTENT( IN ) :: message_log ! ------ ! Result ! ------ INTEGER :: error_status ! ---------- ! Parameters ! ---------- CHARACTER( * ), PARAMETER :: ROUTINE_NAME = 'FIND_CLOSEST_VECTOR' ! --------------- ! Local variables ! --------------- ! -- Matrix dimensions INTEGER :: M, N ! -- Vector coordinate array, Kx1 REAL( fp_kind ), DIMENSION( n_vectors ) :: vector_coordinate ! -- Vector for sum of squared differences. Can't use ! -- an automatic array as the space may not have been ! -- initialised. REAL( fp_kind ), ALLOCATABLE, DIMENSION( : ) :: sum_vector ! -- Vector required to hold MINLOC result INTEGER, DIMENSION( 1 ) :: number ! -- Error variable INTEGER :: allocate_status ! --------------- ! Intrinsics used ! --------------- INTRINSIC ALLOCATED, & MINLOC, & RESHAPE, & SIZE !#--------------------------------------------------------------------------# !# -- HAS THE COORDINATE SPACE BEEN INITIALISED? -- # !#--------------------------------------------------------------------------# IF ( .NOT. space_initialized ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Coordinate space not established.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- GET INPUT VECTOR DIMENSION AND CHECK STUFF -- # !#--------------------------------------------------------------------------# N = SIZE( vector ) ! -------------------------------- ! Is eigenvector matrix available? ! -------------------------------- IF ( .NOT. ALLOCATED( eigenvector ) ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Eigenvector matrix not available.', & error_status, & message_log = message_log ) RETURN END IF ! ------------------------------------------------ ! Check vector length with that of the eigenvector ! matrix dimension ! ------------------------------------------------ IF ( N /= SIZE( eigenvector, DIM = 1 ) ) THEN error_status = FAILURE WRITE( message, '( "Input vector size is inconsistent. Must have ", i4, " elements." )' ) & SIZE( eigenvector, DIM = 1 ) CALL display_message( ROUTINE_NAME, & message, & error_status, & message_log = message_log ) RETURN END IF ! ------------------------------- ! Is coordinate matrix available? ! ------------------------------- IF ( .NOT. ALLOCATED( coordinate ) ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Coordinate matrix not available.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- COMPUTE THE INPUT VECTOR'S COORDINATES -- # !#--------------------------------------------------------------------------# error_status = compute_vector_coordinate( vector, & vector_coordinate, & message_log = message_log ) IF ( error_status /= SUCCESS ) THEN CALL display_message( ROUTINE_NAME, & 'Input vector coordinate computation failed.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- FIND THE CLOSEST VECTOR IN COORDINATE SPACE -- # !# # !# coord_difference = coordinate matrix - replicated vector # !# coordinate # !# # !# 1 ---------> M # !# [ c11 c21 c31 ... cM1 ] [ v1 v1 v1 v1 ] # !# [ c12 ... ] [ v2 v2 v2 v2 ] # !# = [ c13 ... ] - [ v3 v3 v3 .. v3 ] # !# [ ... ... ] [ .. .. .. .. ] # !# [ c1K cMK ] [ vK vK vK vK ] # !# # !# [ d11 d21 d31 ... dM1 ] # !# [ d12 ... ] # !# = [ d13 ... ] # !# [ ... ... ] # !# [ d1K dMK ] # !# # !# This array is squared to give the squared difference array which is # !# summed along the K dimension to produce a Mx1 vector: # !# # !# [ K K K ] # !# sum_of_squared_differences = [ sum{d1i^2}, sum{d2i^2},... sum{dMi^2} ] # !# [ i=1 i=1 i=1 ] # !# # !# The index of the minimum value of this vector, say the j'th value, # !# means that the j'th vector of the original input matrix, (j,1->N), # !# is the vector that is the closest in the coordinate space to the user # !# input vector. # !# # !#--------------------------------------------------------------------------# ! -------------------------------------------------- ! Get the leading dimension of the coordinate matrix ! -------------------------------------------------- M = SIZE( coordinate, DIM = 1 ) ! ------------------------------------------------------------ ! Allocate an array to hold the sum of the squared differences ! ------------------------------------------------------------ ALLOCATE( sum_vector( M ), STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error allocating SUM_VECTOR array. STAT = ", i5 )' ) & allocate_status error_status = FAILURE CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF ! ---------------------------------------------------------------- ! Calculate the vector of the sum of the squared differences ! ! - The RESHAPE replicates the input vector coordinate ! - The result is differenced from COORDINATE and squared; the **2 ! - The SUM totals up the result along DIM=2, the K dimension ! Note that K==n_vectors. ! ! NOTE: This is done separately as the SGI IRIX f90 compiler has ! a bug in it where MINLOC( SUM(....) ) gives the incorrect ! answer. SGI has reproduced the error and, at least since ! Aug 2001, is working on it. ! ---------------------------------------------------------------- sum_vector = SUM( ( coordinate - TRANSPOSE( RESHAPE( vector_coordinate, & (/ n_vectors, M /), & vector_coordinate ) ) )**2, & DIM = 2 ) ! -------------------------------------------------------- ! Determine the number of the closest vector. MINLOC finds ! the index of the smallest value in the argument. ! -------------------------------------------------------- number = MINLOC( sum_vector ) vector_number = number( 1 ) ! ------------------------------- ! Deallocate the sum vector array ! ------------------------------- DEALLOCATE( sum_vector, STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error deallocating SUM_VECTOR array. STAT = ", i5 )' ) & allocate_status error_status = FAILURE CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- DONE -- # !#--------------------------------------------------------------------------# error_status = success END FUNCTION find_closest_vector !-------------------------------------------------------------------------------- !S+ ! NAME: ! destroy_coordinate_space ! ! PURPOSE: ! Function to deallocate all the matrices used to define ! the coordinate space. ! ! CATEGORY: ! Math ! ! LANGUAGE: ! Fortran 90 ! ! CALLING SEQUENCE: ! result = destroy_coordinate_space( message_log = message_log ) ! ! INPUT ARGUMENTS: ! None. ! ! OPTIONAL INPUT ARGUMENTS: ! 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: ! None. ! ! FUNCTION RESULT ! Function returns an error flag ! If result = SUCCESS, everything is o.k. ! = FAILURE, if an error occurred ! ! CALLS: ! display_message: Subroutine to output messages ! SOURCE: error_handler module ! ! EXTERNALS: ! None ! ! COMMON BLOCKS: ! None. ! ! SIDE EFFECTS: ! None known. ! ! RESTRICTIONS: ! None. !S- !-------------------------------------------------------------------------------- FUNCTION destroy_coordinate_space ( message_log ) & ! Optional input RESULT ( error_status ) !#--------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#--------------------------------------------------------------------------# ! --------- ! Arguments ! --------- CHARACTER( * ), OPTIONAL, INTENT( IN ) :: message_log ! ------ ! Result ! ------ INTEGER :: error_status ! ---------- ! Parameters ! ---------- CHARACTER( * ), PARAMETER :: ROUTINE_NAME = 'DESTROY_COORDINATE_SPACE' ! --------------- ! Local variables ! --------------- INTEGER :: allocate_status ! --------------- ! Intrinsics used ! --------------- INTRINSIC ALLOCATED !#--------------------------------------------------------------------------# !# -- HAS THE COORDINATE SPACE BEEN INITIALISED? -- # !#--------------------------------------------------------------------------# IF ( .NOT. space_initialized ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Coordinate space not established.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- DEALLOCATE ALL THE DYNAMIC, SHARED ARRAYS -- # !#--------------------------------------------------------------------------# DEALLOCATE ( average, & vector_norm, & eigenvalue, & eigenvector, & coordinate, & STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error deallocating arrays. STAT = ", i5 )' ) & allocate_status error_status = FAILURE CALL display_message( ROUTINE_NAME, & TRIM( message ), & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- DONE -- # !#--------------------------------------------------------------------------# space_initialized = .FALSE. error_status = SUCCESS END FUNCTION destroy_coordinate_space !############################################################################## ! -- PRIVATE function to calculate input vector coordinates -- !############################################################################## FUNCTION compute_vector_coordinate ( vector, & ! Input vector_coordinate, & ! Output message_log ) & ! Optional input RESULT ( error_status ) !#--------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#--------------------------------------------------------------------------# ! --------- ! Arguments ! --------- REAL( fp_kind ), DIMENSION( : ), INTENT( IN ) :: vector REAL( fp_kind ), DIMENSION( : ), INTENT( OUT ) :: vector_coordinate CHARACTER( * ), OPTIONAL, INTENT( IN ) :: message_log ! ------ ! Result ! ------ INTEGER :: error_status ! ---------- ! Parameters ! ---------- CHARACTER( * ), PARAMETER :: ROUTINE_NAME = 'COMPUTE_VECTOR_COORDINATE' ! --------------- ! Local variables ! --------------- ! -- Input vector dimension INTEGER :: N ! n==nlev ! -- Normalised vector, Nx1 REAL( fp_kind ), DIMENSION( SIZE( vector ) ) :: normalised_vector ! ---------- ! Intrinsics ! ---------- INTRINSIC MATMUL, & SIZE !#--------------------------------------------------------------------------# !# -- HAS THE COORDINATE SPACE BEEN INITIALISED? -- # !#--------------------------------------------------------------------------# IF ( .NOT. space_initialized ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Coordinate space not established.', & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- GET INPUT VECTOR DIMENSION AND CHECK STUFF -- # !#--------------------------------------------------------------------------# N = SIZE( vector ) ! -------------------------------- ! Is eigenvector matrix available? ! -------------------------------- IF ( .NOT. ALLOCATED( eigenvector ) ) THEN error_status = FAILURE CALL display_message( ROUTINE_NAME, & 'Eigenvector matrix not available.', & error_status, & message_log = message_log ) RETURN END IF ! ------------------------------------------------ ! Check vector length with that of the eigenvector ! matrix dimension ! ------------------------------------------------ IF ( N /= SIZE( eigenvector, DIM = 1 ) ) THEN error_status = FAILURE WRITE( message, '( "Input vector size is inconsistent. Must have ", i4, " elements." )' ) & SIZE( eigenvector, DIM = 1 ) CALL display_message( ROUTINE_NAME, & message, & error_status, & message_log = message_log ) RETURN END IF ! ------------------------------------------- ! Check output vector dimension. It must have ! the same size as the number of eigenvectors ! retained in the coordinate space expansion. ! ------------------------------------------- IF ( SIZE( vector_coordinate ) /= n_vectors ) THEN error_status = FAILURE WRITE( message, '( "Output vector size is inconsistent. Must have ", i4, " elements." )' ) & n_vectors CALL display_message( ROUTINE_NAME, & message, & error_status, & message_log = message_log ) RETURN END IF !#--------------------------------------------------------------------------# !# -- CREATE THE NORMALISED VECTOR, a -- # !# # !# The normalised mean deviation vector is, # !# _ # !# ~ a - A # !# a = ------- # !# ^ # !# k # !# # !#--------------------------------------------------------------------------# normalised_vector(:) = ( vector(:) - average(:) ) / & ! -------------------------- vector_norm(:) !#--------------------------------------------------------------------------# !# -- COMPUTE THE VECTOR COORDINATES -- # !# ~ # !# y = [a] * [E] == 1xN * NxK == 1xK coordinate vector # !#--------------------------------------------------------------------------# ! --------------------------------- ! Initialise the vector coordinates ! --------------------------------- vector_coordinate(:) = 0.0_fp_kind ! ----------------------- ! Compute the corrdinates ! ----------------------- vector_coordinate( 1:n_vectors ) = MATMUL( normalised_vector, & eigenvector( :, 1:n_vectors ) ) !#--------------------------------------------------------------------------# !# -- DONE -- # !#--------------------------------------------------------------------------# error_status = success END FUNCTION compute_vector_coordinate END MODULE coordinate_space !------------------------------------------------------------------------------- ! -- MODIFICATION HISTORY -- !------------------------------------------------------------------------------- ! ! $Id: coordinate_space.f90,v 1.9 2001/09/27 17:48:26 paulv Exp $ ! ! $Date: 2001/09/27 17:48:26 $ ! ! $Revision: 1.9 $ ! ! $Name: Ozone_Match_1-1 $ ! ! $State: Exp $ ! ! $Log: coordinate_space.f90,v $ ! Revision 1.9 2001/09/27 17:48:26 paulv ! - Updated module header documentation. ! - Added subprogram header documentation. ! - No longer use FILE_UTILITY module. ! - Changed name of LAPACK module to LAPACK_UTILITY. ! - Removed all WRITE(*, ... ) statements except for those executed when ! the optional DEBUG argument is set in INITIALIZE_COORDINATE_SPACE. ! - Changed the call to LAPACK_SYEV to reflect the changes made to the ! LAPACK_UTILITY module. The "work" array required by the LAPACK SYEV ! routine is no longer needed in the calling routine. Also, the ! eigenvector/value reversal (i.e. ascending or descending in value) ! is now done within the lapack utility function. ! - In the FIND_CLOSEST changed the method of finding the minimum value from ! ! number = MINLOC( SUM( ( coordinate - TRANSPOSE( RESHAPE( vector_coordinate, & ! (/ n_vectors, M /), & ! vector_coordinate ) ) )**2, & ! DIM = 2 ) ) ! ! vector_number = number( 1 ) ! ! to: ! ! sum_vector = SUM( ( coordinate - TRANSPOSE( RESHAPE( vector_coordinate, & ! (/ n_vectors, M /), & ! vector_coordinate ) ) )**2, & ! DIM = 2 ) ! ! number = MINLOC( sum_vector ) ! vector_number = number( 1 ) ! ! i.e., introduce an intermediate variable SUM_VECTOR. This was done as the ! first method failed on SGI platforms - there is an apparent bug in the ! compiler when one uses MINLOC on an _expression_ (as in the first case) ! as opposed to a vector (as in the second). This change made the code work ! on SGIs. ! - Added "Name" to RCS keyword list. ! ! Revision 1.8 2001/06/20 15:38:53 paulv ! - Added eigenvalue and correlation matrix check difference to debug output ! ! Revision 1.7 2001/06/19 22:15:24 paulv ! - Changed code to use ERROR_HANDLER module. ! - Default type kind used. ! - Cosmetic changes ! ! Revision 1.6 2000/04/06 17:23:58 paulv ! - Removed mean_deviation automatic array from initialize_coordinate_space. ! Mean deviations are stored in the unit_vector matrix to conserve stack space. ! - Removed debug NxN array from initialize_coordinate_space. Debug calcs ! re-use the correlation matrix which is not needed after the eigenvector ! decomposition is done. This also was done to conserve stack space. ! - Updated header documentation. ! ! Revision 1.5 2000/04/06 15:07:05 paulv ! Added GET_LUN() declaration in compute_vector_coordinate for debug ! output. ! ! Revision 1.4 2000/04/05 16:15:47 paulv ! - Renamed compute_coordinate_space to initialize_coordinate_space ! - Renamed cleanup_coordinate_space to destroy_coordinate_space ! - Removed SAVE attributs from common, private arrays ! - Local temporary arrays changed from ALLOCATABLE to automatic ! - Added debug argument to compute_vector_coordinate ! - Added debug argument to find_closest_vector ! ! Revision 1.3 2000/04/03 20:09:31 paulv ! Added documentation header. ! ! Revision 1.2 2000/04/03 18:18:49 paulv ! - Used ONLY clause in the type kinds USE statement to allow easy ! redefinition of the default floating point type used in the calculations. ! - Included lapack module in the compute_coordinate_space routine. This ! overloads the LAPACK call so the FP precision change does not require ! modifying the LAPACK routine call (e.g. from CALL ssyev to CALL dsyev). ! ! Revision 1.1 2000/04/03 15:43:01 paulv ! Initial checked-in version ! !