!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A Module For Computing Matrix functions based on Chebyshev polynomials. MODULE ChebyshevSolversModule USE DataTypesModule, ONLY : NTREAL USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix USE LoggingModule, ONLY : WriteElement, WriteHeader, EnterSubLog, & & ExitSubLog USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, & & DestructMatrixMemoryPool USE PSMatrixAlgebraModule, ONLY : MatrixMultiply, IncrementMatrix, ScaleMatrix USE PSMatrixModule, ONLY : Matrix_ps, FillMatrixIdentity, & & PrintMatrixInformation, ConstructEmptyMatrix, DestructMatrix, & & CopyMatrix USE SolverParametersModule, ONLY : SolverParameters_t, PrintParameters IMPLICIT NONE PRIVATE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> A datatype that represents a Chebyshev polynomial. TYPE, PUBLIC :: ChebyshevPolynomial_t !> Coefficients of the polynomial. REAL(NTREAL), DIMENSION(:), ALLOCATABLE :: coefficients END TYPE ChebyshevPolynomial_t !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Polynomial type PUBLIC :: ConstructPolynomial PUBLIC :: DestructPolynomial PUBLIC :: SetCoefficient !! Solvers PUBLIC :: Compute PUBLIC :: FactorizedCompute !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTERFACE ConstructPolynomial MODULE PROCEDURE ConstructPolynomial_cheby END INTERFACE INTERFACE DestructPolynomial MODULE PROCEDURE DestructPolynomial_cheby END INTERFACE INTERFACE SetCoefficient MODULE PROCEDURE SetCoefficient_cheby END INTERFACE INTERFACE Compute MODULE PROCEDURE Compute_cheby END INTERFACE INTERFACE FactorizedCompute MODULE PROCEDURE FactorizedCompute_cheby END INTERFACE CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Construct a Chebyshev polynomial object. PURE SUBROUTINE ConstructPolynomial_cheby(this, degree) !> The polynomial to construct. TYPE(ChebyshevPolynomial_t), INTENT(INOUT) :: this !> Degree of the polynomial. INTEGER, INTENT(IN) :: degree ALLOCATE(this%coefficients(degree)) this%coefficients = 0 END SUBROUTINE ConstructPolynomial_cheby !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Destruct a polynomial object. PURE SUBROUTINE DestructPolynomial_cheby(this) !> The polynomial to destruct. TYPE(ChebyshevPolynomial_t), INTENT(INOUT) :: this IF (ALLOCATED(this%coefficients)) THEN DEALLOCATE(this%coefficients) END IF END SUBROUTINE DestructPolynomial_cheby !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Set a coefficient of a Chebyshev polynomial. SUBROUTINE SetCoefficient_cheby(this, degree, coefficient) !> The polynomial to set. TYPE(ChebyshevPolynomial_t), INTENT(INOUT) :: this !> Degree for which to set the coefficient. INTEGER, INTENT(IN) :: degree !> Coefficient value. REAL(NTREAL), INTENT(IN) :: coefficient this%coefficients(degree) = coefficient END SUBROUTINE SetCoefficient_cheby !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Compute The Chebyshev Polynomial of the matrix. !> This method uses the standard Chebyshev Polynomial expansion. SUBROUTINE Compute_cheby(InputMat, OutputMat, poly, solver_parameters_in) !> The input matrix TYPE(Matrix_ps), INTENT(IN) :: InputMat !> OutputMat = poly(InputMat) TYPE(Matrix_ps), INTENT(INOUT) :: OutputMat !> The Chebyshev polynomial to compute. TYPE(ChebyshevPolynomial_t), INTENT(IN) :: poly !> Parameters for the solver. TYPE(SolverParameters_t), INTENT(IN), OPTIONAL :: solver_parameters_in !! Handling Solver Parameters TYPE(SolverParameters_t) :: solver_parameters !! Local Matrices TYPE(Matrix_ps) :: Identity TYPE(Matrix_ps) :: BalancedInput TYPE(Matrix_ps) :: Tk TYPE(Matrix_ps) :: Tkminus1 TYPE(Matrix_ps) :: Tkminus2 TYPE(MatrixMemoryPool_p) :: pool !! Local Variables INTEGER :: degree INTEGER :: counter !! Handle The Optional Parameters IF (PRESENT(solver_parameters_in)) THEN solver_parameters = solver_parameters_in ELSE solver_parameters = SolverParameters_t() END IF degree = SIZE(poly%coefficients) IF (solver_parameters%be_verbose) THEN CALL WriteHeader("Chebyshev Solver") CALL EnterSubLog CALL WriteElement(key="Method", value="Standard") CALL WriteElement(key="Degree", value=degree-1) CALL PrintParameters(solver_parameters) END IF !! Initial values for matrices CALL ConstructEmptyMatrix(Identity, InputMat) CALL FillMatrixIdentity(Identity) CALL CopyMatrix(InputMat,BalancedInput) !! Load Balancing Step IF (solver_parameters%do_load_balancing) THEN CALL PermuteMatrix(Identity, Identity, & & solver_parameters%BalancePermutation, memorypool_in=pool) CALL PermuteMatrix(BalancedInput, BalancedInput, & & solver_parameters%BalancePermutation, memorypool_in=pool) END IF !! First Term CALL CopyMatrix(Identity,Tkminus2) IF (degree == 1) THEN CALL CopyMatrix(Tkminus2,OutputMat) CALL ScaleMatrix(OutputMat,poly%coefficients(1)) ELSE CALL CopyMatrix(BalancedInput,Tkminus1) CALL CopyMatrix(Tkminus2,OutputMat) CALL ScaleMatrix(OutputMat,poly%coefficients(1)) CALL IncrementMatrix(Tkminus1,OutputMat, & & alpha_in=poly%coefficients(2)) IF (degree > 2) THEN CALL MatrixMultiply(BalancedInput, Tkminus1, Tk, & & alpha_in=REAL(2.0,NTREAL), & & threshold_in=solver_parameters%threshold, memory_pool_in=pool) CALL IncrementMatrix(Tkminus2,Tk,REAL(-1.0,NTREAL)) CALL IncrementMatrix(Tk, OutputMat, & & alpha_in=poly%coefficients(3)) DO counter = 4, degree CALL CopyMatrix(Tkminus1,Tkminus2) CALL CopyMatrix(Tk,Tkminus1) CALL MatrixMultiply(BalancedInput, Tkminus1, Tk, & & alpha_in=REAL(2.0,NTREAL), & & threshold_in=solver_parameters%threshold, & & memory_pool_in=pool) CALL IncrementMatrix(Tkminus2,Tk, & & REAL(-1.0,NTREAL)) CALL IncrementMatrix(Tk, OutputMat, & & alpha_in=poly%coefficients(counter)) END DO END IF END IF IF (solver_parameters%be_verbose) THEN CALL PrintMatrixInformation(OutputMat) END IF !! Undo Load Balancing Step IF (solver_parameters%do_load_balancing) THEN CALL UndoPermuteMatrix(OutputMat, OutputMat, & & solver_parameters%BalancePermutation, memorypool_in=pool) END IF !! Cleanup IF (solver_parameters%be_verbose) THEN CALL ExitSubLog END IF CALL DestructMatrix(Identity) CALL DestructMatrix(Tk) CALL DestructMatrix(Tkminus1) CALL DestructMatrix(Tkminus2) CALL DestructMatrix(BalancedInput) CALL DestructMatrixMemoryPool(pool) END SUBROUTINE Compute_cheby !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Compute The Chebyshev Polynomial of the matrix. !> This version first factors the Chebyshev Polynomial and computes the !> function using a divide and conquer algorithm. Based on a simplified !> version of the first method in \cite liang2003improved . SUBROUTINE FactorizedCompute_cheby(InputMat, OutputMat, poly, & & solver_parameters_in) !> The input matrix TYPE(Matrix_ps), INTENT(IN) :: InputMat !> OutputMat = poly(InputMat) TYPE(Matrix_ps), INTENT(INOUT) :: OutputMat !> The Chebyshev polynomial to compute. TYPE(ChebyshevPolynomial_t), INTENT(IN) :: poly !> Parameters for the solver. TYPE(SolverParameters_t), INTENT(IN), OPTIONAL :: solver_parameters_in !! Handling Solver Parameters TYPE(SolverParameters_t) :: solver_parameters !! Local Matrices TYPE(Matrix_ps) :: Identity TYPE(Matrix_ps) :: BalancedInput TYPE(Matrix_ps), DIMENSION(:), ALLOCATABLE :: T_Powers TYPE(MatrixMemoryPool_p) :: pool !! Local Variables INTEGER :: degree INTEGER :: log2degree INTEGER :: counter !! Handle The Optional Parameters IF (PRESENT(solver_parameters_in)) THEN solver_parameters = solver_parameters_in ELSE solver_parameters = SolverParameters_t() END IF degree = SIZE(poly%coefficients) IF (solver_parameters%be_verbose) THEN CALL WriteHeader("Chebyshev Solver") CALL EnterSubLog CALL WriteElement(key="Method", value="Recursive") CALL WriteElement(key="Degree", value=degree-1) CALL PrintParameters(solver_parameters) END IF !! Initial values for matrices CALL ConstructEmptyMatrix(Identity, InputMat) CALL FillMatrixIdentity(Identity) CALL CopyMatrix(InputMat,BalancedInput) !! Load Balancing Step IF (solver_parameters%do_load_balancing) THEN CALL PermuteMatrix(Identity, Identity, & & solver_parameters%BalancePermutation, memorypool_in=pool) CALL PermuteMatrix(BalancedInput, BalancedInput, & & solver_parameters%BalancePermutation, memorypool_in=pool) END IF !! Construct The X Powers Array !! First, compute how many powers of two are necessary to compute. log2degree = 1 DO WHILE (2**log2degree .LE. degree) log2degree = log2degree + 1 END DO ALLOCATE(T_Powers(log2degree)) !! Now compute those powers of two CALL CopyMatrix(Identity, T_Powers(1)) IF (degree .EQ. 1) THEN CALL CopyMatrix(T_Powers(1), OutputMat) ELSE CALL CopyMatrix(BalancedInput,T_Powers(2)) DO counter=3,log2degree CALL MatrixMultiply(T_Powers(counter-1), T_Powers(counter-1), & & T_Powers(counter), threshold_in=solver_parameters%threshold, & & alpha_in=REAL(2.0,NTREAL), memory_pool_in=pool) CALL IncrementMatrix(Identity, T_Powers(counter), & & alpha_in=REAL(-1.0,NTREAL)) END DO !! Call Recursive CALL ComputeRecursive(T_Powers, poly, OutputMat, & & pool, 1, solver_parameters) END IF IF (solver_parameters%be_verbose) THEN CALL ExitSubLog CALL PrintMatrixInformation(OutputMat) END IF !! Undo Load Balancing Step IF (solver_parameters%do_load_balancing) THEN CALL UndoPermuteMatrix(OutputMat, OutputMat, & & solver_parameters%BalancePermutation, memorypool_in=pool) END IF !! Cleanup IF (solver_parameters%be_verbose) THEN CALL ExitSubLog END IF DO counter=1,log2degree CALL DestructMatrix(T_Powers(counter)) END DO DEALLOCATE(T_Powers) CALL DestructMatrix(Identity) CALL DestructMatrix(BalancedInput) CALL DestructMatrixMemoryPool(pool) END SUBROUTINE FactorizedCompute_cheby !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> The workhorse routine for the factorized chebyshev computation function. RECURSIVE SUBROUTINE ComputeRecursive(T_Powers, poly, OutputMat, pool, & & depth, solver_parameters) !> The precomputed Chebyshev polynomials. TYPE(Matrix_ps), DIMENSION(:), INTENT(IN) :: T_Powers !> Polynomial coefficients. TYPE(ChebyshevPolynomial_t), INTENT(IN) :: poly !> OutputMat = poly(InputMat) TYPE(Matrix_ps), INTENT(INOUT) :: OutputMat !> The depth of recursion. INTEGER, INTENT(in) :: depth !> Parameters for the solver. TYPE(SolverParameters_t), INTENT(IN) :: solver_parameters !> The memory pool. TYPE(MatrixMemoryPool_p), INTENT(INOUT) :: pool !! Local Data INTEGER :: coefficient_midpoint INTEGER :: left_length, right_length INTEGER :: full_midpoint INTEGER :: counter TYPE(ChebyshevPolynomial_t) :: left_poly TYPE(ChebyshevPolynomial_t) :: right_poly TYPE(Matrix_ps) :: LeftMat TYPE(Matrix_ps) :: RightMat !! First Handle The Base Case IF (SIZE(poly%coefficients) .EQ. 1) THEN CALL CopyMatrix(T_Powers(1), OutputMat) CALL ScaleMatrix(OutputMat, poly%coefficients(1)) ELSE IF (SIZE(poly%coefficients) .EQ. 2) THEN CALL CopyMatrix(T_Powers(1), OutputMat) CALL ScaleMatrix(OutputMat, poly%coefficients(1)) CALL IncrementMatrix(T_Powers(2), OutputMat, & & alpha_in=poly%coefficients(2)) ELSE !! Adjust the coefficients. coefficient_midpoint = SIZE(poly%coefficients)/2 left_length = coefficient_midpoint right_length = SIZE(poly%coefficients) - coefficient_midpoint ALLOCATE(left_poly%coefficients(left_length)) ALLOCATE(right_poly%coefficients(right_length)) left_poly%coefficients(:) = poly%coefficients(:coefficient_midpoint) right_poly%coefficients(:) = poly%coefficients(coefficient_midpoint+1:) DO counter=2,SIZE(left_poly%coefficients) left_poly%coefficients(counter) = left_poly%coefficients(counter) - & & poly%coefficients(SIZE(poly%coefficients) - counter + 2) END DO !! Left recursion CALL ComputeRecursive(T_Powers, left_poly, LeftMat, pool, depth+1, & & solver_parameters) !! Right recursion full_midpoint = SIZE(T_Powers) - depth + 1 CALL ComputeRecursive(T_Powers, right_poly, RightMat, pool, depth+1, & & solver_parameters) !! Sum Together CALL MatrixMultiply(T_Powers(full_midpoint), RightMat, & & OutputMat, threshold_in=solver_parameters%threshold, & & alpha_in=REAL(2.0,NTREAL), memory_pool_in=pool) CALL IncrementMatrix(LeftMat,OutputMat) CALL IncrementMatrix(T_Powers(full_midpoint), & & OutputMat, alpha_in=-1.0*right_poly%coefficients(1)) !! Cleanup DEALLOCATE(left_poly%coefficients) DEALLOCATE(right_poly%coefficients) CALL DestructMatrix(LeftMat) CALL DestructMatrix(RightMat) END IF END SUBROUTINE ComputeRecursive !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE ChebyshevSolversModule