Loading src/CMakeLists.txt +1 −0 Original line number Diff line number Diff line Loading @@ -32,6 +32,7 @@ TRIBITS_ADD_LIBRARY(Utils Futility_DBC.f90 UnitTest.f90 IntrType.f90 Hash.f90 Search.f90 Constants_Conversion.f90 ExtendedMath.f90 Loading src/Geom_Graph.f90 +227 −119 Original line number Diff line number Diff line Loading @@ -27,6 +27,7 @@ PRIVATE PUBLIC :: GraphType PUBLIC :: DAGraphType PUBLIC :: OPERATOR(==) PUBLIC :: ASSIGNMENT(=) !PUBLIC :: OPERATOR(+) INTEGER(SIK),PARAMETER :: GRAPH_NULL_EDGE=0 Loading @@ -35,6 +36,8 @@ INTEGER(SIK),PARAMETER :: GRAPH_QUADRATIC_EDGE=-1 !> @brief a Directed Acyclic Graph Type TYPE :: DAGraphType !> Logical indicating if the graph is sorted by ascending node value LOGICAL(SBK) :: sorted=.FALSE. !> The number of nodes on the graph INTEGER(SIK) :: n=0 !> Nodes (Size n) Loading @@ -50,9 +53,13 @@ TYPE :: DAGraphType !> @copybrief Geom_Graph::clear_DAGraphType !> @copydetails Geom_Graph::clear_DAGraphType PROCEDURE,PASS :: clear => clear_DAGraphType !> @copybrief Geom_Graph::defineEdge_DAGraphType !> @copydetails Geom_Graph::defineEdge_DAGraphType PROCEDURE,PASS :: insertNode=> insertNode_DAGraphType !> @copybrief Geom_Graph::insertNodeSorted_DAGraphType !> @copydetails Geom_Graph::insertNodeSorted_DAGraphType PROCEDURE,PASS,PRIVATE :: insertNodeSorted => insertNodeSorted_DAGraphType !> @copybrief Geom_Graph::insertNodeAtIndex_DAGraphType !> @copydetails Geom_Graph::insertNodeAtIndex_DAGraphType PROCEDURE,PASS,PRIVATE :: insertNodeAtIndex => insertNodeAtIndex_DAGraphType GENERIC :: insertNode => insertNodeSorted,insertNodeAtIndex !> @copybrief Geom_Graph::removeNode_DAGraphType !> @copydetails Geom_Graph::removeNode_DAGraphType PROCEDURE,PASS :: removeNode=> removeNode_DAGraphType Loading Loading @@ -225,6 +232,12 @@ INTERFACE OPERATOR(==) MODULE PROCEDURE isequal_graphType ENDINTERFACE INTERFACE ASSIGNMENT(=) !> @copybrief Geom_Graph::assign_DAGraphType !> @copydetails Geom_Graph::assign_DAGraphType MODULE PROCEDURE assign_DAGraphType ENDINTERFACE !INTERFACE OPERATOR(+) ! !> @copybrief Geom_Graph::add_graphType ! !> @copydetails Geom_Graph::add_graphType Loading Loading @@ -2209,19 +2222,21 @@ ENDFUNCTION isequal_GraphType SUBROUTINE init_DAGraphType(this,n,nodes) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: n INTEGER(SIK),ALLOCATABLE,INTENT(IN) :: nodes(:) INTEGER(SIK),INTENT(IN) :: nodes(:) REQUIRE(.NOT.ALLOCATED(this%nodes)) REQUIRE(SIZE(nodes) == n .OR. n == 0) IF(.NOT. ALLOCATED(this%nodes) .AND. ALLOCATED(nodes)) THEN IF(n > 0) THEN IF(SIZE(nodes) == n) THEN this%n=n ALLOCATE(this%nodes(n)) ALLOCATE(this%edgeMatrix(n,n)) IF(n > 0) THEN this%nodes=nodes this%edgeMatrix=0 ENDIF ENDIF CALL sort(this%nodes) this%sorted=.TRUE. ENDIF this%edgeMatrix=0 ENDSUBROUTINE init_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2234,35 +2249,60 @@ SUBROUTINE clear_DAGraphType(this) this%n=0 IF(ALLOCATED(this%nodes)) DEALLOCATE(this%nodes) IF(ALLOCATED(this%edgeMatrix)) DEALLOCATE(this%edgeMatrix) this%sorted=.FALSE. ENDSUBROUTINE clear_DAGraphType ! !------------------------------------------------------------------------------- !> @brief !> @param !> SUBROUTINE insertNode_DAGraphType(this,ID,ind) SUBROUTINE insertNodeSorted_DAGraphType(this,ID) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ID INTEGER(SIK),INTENT(IN),OPTIONAL :: ind INTEGER(SIK) :: index INTEGER(SIK) :: index,a,b INTEGER(SIK),ALLOCATABLE :: tmpNodes(:),tmpEdges(:,:) index=1 IF(PRESENT(ind)) index=ind IF(ALLOCATED(this%edgeMatrix)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) IF(ALL(this%nodes /= ID)) THEN !Store values temporarily ALLOCATE(tmpNodes(this%n)) ALLOCATE(tmpEdges(this%n,this%n)) tmpNodes=this%nodes tmpEdges=this%edgeMatrix CALL MOVE_ALLOC(this%nodes,tmpNodes) CALL MOVE_ALLOC(this%edgeMatrix,tmpEdges) !Resize arrays on graph type DEALLOCATE(this%nodes) DEALLOCATE(this%edgeMatrix) this%n=this%n+1 ALLOCATE(this%nodes(this%n)) ALLOCATE(this%edgeMatrix(this%n,this%n)) this%edgeMatrix=0 !Find the location for this index. If the graph is sorted, then we'll insert !the value at the location that maintains the sorting. Otherwise, we'll just !stick it at the beginning !Also handle some special cases up front to simplify the binary search for sorted insertion IF(this%n == 0) THEN index = 1 ELSEIF(this%sorted) THEN a = 1 b = this%n IF(ID < this%nodes(a)) THEN index = 1 ELSEIF(ID > this%nodes(b)) THEN index = this%n ELSE index = (a + b)/2 DO WHILE(.TRUE.) IF(ID < this%nodes(index-1)) THEN b = index ELSEIF(ID > this%nodes(index)) THEN a = index ELSE EXIT ENDIF index = (a + b)/2 ENDDO ENDIF ELSE index = 1 ENDIF !reassign old data and assign new node this%nodes(1:index-1)=tmpNodes(1:index-1) this%nodes(index)=ID Loading @@ -2271,17 +2311,55 @@ SUBROUTINE insertNode_DAGraphType(this,ID,ind) this%edgeMatrix(index+1:,1:index-1)=tmpEdges(index:,1:index-1) this%edgeMatrix(1:index-1,index+1:)=tmpEdges(1:index-1,index:) this%edgeMatrix(index+1:,index+1:)=tmpEdges(index:,index:) DEALLOCATE(tmpNodes) DEALLOCATE(tmpEdges) ENDIF ELSE ALLOCATE(this%edgeMatrix(1,1)) ALLOCATE(this%nodes(1)) ENDSUBROUTINE insertNodeSorted_DAGraphType ! !------------------------------------------------------------------------------- !> @brief !> @param !> SUBROUTINE insertNodeAtIndex_DAGraphType(this,ID,index) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ID INTEGER(SIK),INTENT(IN) :: index INTEGER(SIK),ALLOCATABLE :: tmpNodes(:),tmpEdges(:,:) REQUIRE(ALLOCATED(this%edgeMatrix)) IF(ALL(this%nodes /= ID)) THEN !Store values temporarily CALL MOVE_ALLOC(this%nodes,tmpNodes) CALL MOVE_ALLOC(this%edgeMatrix,tmpEdges) !Resize arrays on graph type this%n=this%n+1 ALLOCATE(this%nodes(this%n)) ALLOCATE(this%edgeMatrix(this%n,this%n)) this%edgeMatrix=0 this%nodes=ID this%n=1 !reassign old data and assign new node this%nodes(1:index-1)=tmpNodes(1:index-1) this%nodes(index)=ID this%nodes(index+1:this%n)=tmpNodes(index:) this%edgeMatrix(1:index-1,1:index-1)=tmpEdges(1:index-1,1:index-1) this%edgeMatrix(index+1:,1:index-1)=tmpEdges(index:,1:index-1) this%edgeMatrix(1:index-1,index+1:)=tmpEdges(1:index-1,index:) this%edgeMatrix(index+1:,index+1:)=tmpEdges(index:,index:) !Update the sorted status IF(this%sorted) THEN IF(index > 1) THEN IF(this%nodes(index) <= this%nodes(index-1)) THEN this%sorted=.FALSE. ENDIF ENDIF IF(index < this%n) THEN IF(this%nodes(index) >= this%nodes(index+1)) THEN this%sorted=.FALSE. ENDIF ENDIF ENDIF ENDSUBROUTINE insertNode_DAGraphType ENDIF ENDSUBROUTINE insertNodeAtIndex_DAGraphType ! !------------------------------------------------------------------------------- !> @brief Loading @@ -2291,17 +2369,11 @@ SUBROUTINE removeNode_DAGraphType(this,ID,ind) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN),OPTIONAL :: ID INTEGER(SIK),INTENT(IN),OPTIONAL :: ind INTEGER(SIK) :: i,index INTEGER(SIK) :: index index=0 IF(PRESENT(ID)) THEN !Find the index to remove DO i=1,this%n IF(ID == this%nodes(i)) THEN index=i EXIT ENDIF ENDDO index = this%getIndex(ID) !If the ID was found, remove the index. IF(index > 0) CALL removeNodeByIndex_DAGraphType(this,index) ELSEIF(PRESENT(ind)) THEN Loading @@ -2316,24 +2388,20 @@ ENDSUBROUTINE removeNode_DAGraphType SUBROUTINE removeNodeByIndex_DAGraphType(this,ind) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ind INTEGER(SIK) :: i INTEGER(SIK),ALLOCATABLE :: tmpN(:),tmpEdge(:,:) IF(ALLOCATED(this%edgeMatrix)) THEN IF((1 <= ind) .AND. (ind <= this%n)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) REQUIRE(1 <= ind) REQUIRE(ind <= this%n) !Store values temporarily i=this%n ALLOCATE(tmpN(i)) ALLOCATE(tmpEdge(i,i)) tmpN=this%nodes tmpEdge=this%edgeMatrix CALL MOVE_ALLOC(this%nodes,tmpN) CALL MOVE_ALLOC(this%edgeMatrix,tmpEdge) !Resize arrays on graph type DEALLOCATE(this%nodes) DEALLOCATE(this%edgeMatrix) this%n=this%n-1 IF(this%n > 0) THEN ALLOCATE(this%nodes(this%n)) ALLOCATE(this%edgeMatrix(this%n,this%n)) IF(this%n > 0) THEN this%edgeMatrix=0 !reassign old data and assign new node this%nodes(1:ind-1)=tmpN(1:ind-1) Loading @@ -2345,10 +2413,7 @@ SUBROUTINE removeNodeByIndex_DAGraphType(this,ind) this%edgeMatrix(ind:,ind:)=tmpEdge(ind+1:,ind+1:) ENDIF ENDIF !CALL demalloc(tmpint) DEALLOCATE(tmpEdge) ENDIF ENDIF ENDSUBROUTINE removeNodeByIndex_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2365,9 +2430,9 @@ FUNCTION isStartNode_DAGraphType(this,ID,ind) RESULT(bool) bool=.FALSE. IF(PRESENT(ID)) THEN index=this%getIndex(ID) IF(ALL(this%edgeMatrix(:,index) == 0)) bool=.TRUE. bool = (ALL(this%edgeMatrix(:,index) == 0)) ELSEIF(PRESENT(ind)) THEN IF(ALL(this%edgeMatrix(:,ind) == 0)) bool=.TRUE. bool = (ALL(this%edgeMatrix(:,ind) == 0)) ENDIF ENDFUNCTION isStartNode_DAGraphType ! Loading @@ -2383,11 +2448,13 @@ SUBROUTINE getNextStartNode_DAGraphType(this,old,ID) ID=0 DO i=1,this%n IF(this%isStartNode(IND=i) .AND. ALL(this%nodes(i) /= old)) THEN IF(this%isStartNode(IND=i)) THEN IF(ALL(this%nodes(i) /= old)) THEN ID=this%nodes(i) old(i)=ID EXIT ENDIF ENDIF ENDDO ENDSUBROUTINE getNextStartNode_DAGraphType ! Loading @@ -2401,11 +2468,12 @@ SUBROUTINE defineEdge_DAGraphType(this,fromID,toID) INTEGER(SIK),INTENT(IN) :: toID INTEGER(SIK) :: fromInd,toInd IF(ALLOCATED(this%edgeMatrix)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) fromInd=this%getIndex(fromID) toInd=this%getIndex(toID) IF((fromInd > 0) .AND. (toInd > 0)) this%edgeMatrix(fromInd,toInd)=1 ENDIF ENDSUBROUTINE defineEdge_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2418,31 +2486,67 @@ SUBROUTINE removeEdge_DAGraphType(this,fromID,toID) INTEGER(SIK),INTENT(IN) :: toID INTEGER(SIK) :: fromInd,toInd IF(ALLOCATED(this%edgeMatrix)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) fromInd=this%getIndex(fromID) toInd=this%getIndex(toID) IF((fromInd > 0) .AND. (toInd > 0)) this%edgeMatrix(fromInd,toInd)=0 ENDIF ENDSUBROUTINE removeEdge_DAGraphType ! !------------------------------------------------------------------------------- !> @brief !> @param !> @brief retrieves the index of a node in a @c DAGraphType !> @param this the graph object !> @param ID the node ID to search for !> @returns ind the index of node @c ID !> !> Returns 0 if the ID could not be found. If the graph is sorted, a binary !> search is used. Otherwise, a linear search is used. !> FUNCTION getIndex_DAGraphType(this,ID) RESULT(ind) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ID INTEGER(SIK) :: i,ind INTEGER(SIK) :: ind ! INTEGER(SIK) :: a,b REQUIRE(ALLOCATED(this%nodes)) ind=0 IF(ALLOCATED(this%edgeMatrix)) THEN DO i=1,this%n IF(this%nodes(i) == ID) THEN ind=i IF(this%sorted) THEN a = 1 b = this%n IF(ID == this%nodes(a)) THEN ind = 1 ELSEIF(ID == this%nodes(b)) THEN ind = this%n ELSE ind = (a + b)/2 DO WHILE(.TRUE.) IF(ID < this%nodes(ind)) THEN b = ind ELSEIF(ID > this%nodes(ind)) THEN a = ind ELSE EXIT ENDIF IF(a == b) EXIT ind = (a + b)/2 ENDDO IF(this%nodes(ind) /= ID) THEN ind = 0 ENDIF ENDIF ELSE DO a=1,this%n IF(this%nodes(a) == ID) THEN ind=a EXIT ENDIF ENDDO IF(a > this%n) ind=0 ENDIF ENDFUNCTION getIndex_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2459,18 +2563,22 @@ SUBROUTINE KATS_DAGraphType(this) ALLOCATE(oldIDs(this%n)) oldIDs=0 CALL this%getNextStartNode(oldIDs,ID) sortedGraph%n=0 ALLOCATE(sortedGraph%nodes(sortedGraph%n)) ALLOCATE(sortedGraph%edgeMatrix(sortedGraph%n,sortedGraph%n)) DO WHILE(ID /= 0) !Add the node to the sorted graph CALL sortedGraph%insertNode(ID,1) !Remove all edges from that node DO i=1,this%n CALL this%removeEdge(ID,this%nodes(i)) ENDDO !Get the next "start node CALL this%getNextStartNode(oldIDs,ID) ENDDO IF(ALL(this%edgeMatrix == 0)) THEN SELECTTYPE(this); TYPE IS(DAGraphtype) this=sortedGraph ENDSELECT ENDIF DEALLOCATE(oldIDs) ENDSUBROUTINE KATS_DAGraphType Loading @@ -2482,13 +2590,13 @@ ENDSUBROUTINE KATS_DAGraphType SUBROUTINE assign_DAGraphType(g0,g1) CLASS(DAGraphType),INTENT(INOUT) :: g0 CLASS(DAGraphType),INTENT(IN) :: g1 CALL clear_DAGraphType(g0) IF(g1%n > 0) THEN CALL g0%clear() g0%n=g1%n ALLOCATE(g0%nodes(g1%n)) ALLOCATE(g0%edgeMatrix(g1%n,g1%n)) g0%nodes=g1%nodes g0%edgeMatrix=g1%edgeMatrix IF(g1%n > 0) THEN g0%nodes(:)=g1%nodes g0%edgeMatrix(:,:)=g1%edgeMatrix ENDIF ENDSUBROUTINE assign_DAGraphType ! Loading src/Hash.f90 0 → 100644 +53 −0 Original line number Diff line number Diff line !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! Futility Development Group ! ! All rights reserved. ! ! ! ! Futility is a jointly-maintained, open-source project between the University ! ! of Michigan and Oak Ridge National Laboratory. The copyright and license ! ! can be found in LICENSE.txt in the head directory of this repository. ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! !> @brief Provides hashing algorithms !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! MODULE HashModule USE IntrType IMPLICIT NONE PRIVATE PUBLIC :: stringHash CONTAINS ! !------------------------------------------------------------------------------- !> @brief Hashes a string to an integer value !> @param string the string to hash !> @param p the rolling polynomial base !> @param m the modulus !> !> This algorithm is based on https://cp-algorithms.com/string/string-hashing.html !> !> It allows for ASCII values 32 and up. Recommended value for p is 97 (a prime !> number close to 95, which is the number of "normal" ASCII values expected in the !> 32-126 range). Recommended value for m is a large prime number such as 10**9 + 9. !> However, other arguments may be valid for certain applications. !> PURE FUNCTION stringHash(string,p,m) RESULT(hash) CHARACTER(LEN=*),INTENT(IN) :: string INTEGER(SLK),INTENT(IN) :: p INTEGER(SLK),INTENT(IN) :: m INTEGER(SLK) :: hash ! INTEGER(SIK) :: i INTEGER(SLK) :: p_pow hash = 0 p_pow = 1 DO i=1,LEN(string) !1-31 are special characters that we disallow hash = MOD(hash + (IACHAR(string(i:i))-31_SIK)*p_pow,m) p_pow = MOD(p_pow * p,m) ENDDO !i ENDFUNCTION stringHash ! ENDMODULE HashModule No newline at end of file src/Interpolators.f90 +124 −48 Original line number Diff line number Diff line Loading @@ -47,25 +47,27 @@ CONTAINS !> @param table contains data to be interpolated between defined at points !> given in labels !> @param point point at which interpolation needs to be done to return a value !> @param lextrap logical indicating if extrapolation should be performed if called for !> @returns Interpolant interpolated data found through linear interpolation to the point point FUNCTION Interp_1D(labels,table,point) RESULT(Interpolant) FUNCTION Interp_1D(labels,table,point,lextrap) RESULT(Interpolant) REAL(SRK),INTENT(IN) :: labels(:) REAL(SRK),INTENT(IN) :: table(:) REAL(SRK),INTENT(IN) :: point REAL(SRK) :: Interpolant,f(2) INTEGER(SIK) :: i,N_i,i_p LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK) :: Interpolant,f(2),t(2) LOGICAL(SBK) :: canExtrapolate REQUIRE(SIZE(labels) > 1) REQUIRE(SIZE(labels) == SIZE(table)) REQUIRE(isStrictlyIncDec(labels)) CALL Get_points_and_weights(labels,point,f,i_p,N_i) canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Get interpolated value Interpolant=0.0_SRK DO i=MAX(i_p-1,1),MIN(i_p,N_i) Interpolant=Interpolant+f(i-i_p+2)*table(i) ENDDO CALL Get_points_and_weights(labels,table,point,f,t,canExtrapolate) Interpolant=f(1)*t(1)+f(2)*t(2) ENDFUNCTION Interp_1D ! Loading @@ -77,14 +79,17 @@ ENDFUNCTION Interp_1D !> @param table contains data to be interpolated between defined at points !> given in labels1 and labels2 !> @param point point at which interpolation needs to be done to return a value !> @param lextrap logical indicating if extrapolation should be performed if called for !> @returns Interpolant interpolated data found through linear interpolation to the point point FUNCTION Interp_2D(labels1,labels2,table,point) RESULT(Interpolant) FUNCTION Interp_2D(labels1,labels2,table,point,lextrap) RESULT(Interpolant) REAL(SRK),INTENT(IN) :: labels1(:) REAL(SRK),INTENT(IN) :: labels2(:) REAL(SRK),INTENT(IN) :: table(:,:) REAL(SRK),INTENT(IN) :: point(:) REAL(SRK) :: Interpolant,f1(2),f2(2) INTEGER(SIK) :: i,j,N_i,N_j,i_p,j_p LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK) :: Interpolant,f(2),t(2),Interp1D(2) INTEGER(SIK) :: j,i_p(2),j_p(2) LOGICAL(SBK) :: canExtrapolate REQUIRE(SIZE(labels1) > 1) REQUIRE(SIZE(labels2) > 1) Loading @@ -93,16 +98,21 @@ FUNCTION Interp_2D(labels1,labels2,table,point) RESULT(Interpolant) REQUIRE(isStrictlyIncDec(labels1)) REQUIRE(isStrictlyIncDec(labels2)) CALL Get_points_and_weights(labels1,point(1),f1,i_p,N_i) CALL Get_points_and_weights(labels2,point(2),f2,j_p,N_j) !Get interpolated value Interpolant=0.0_SRK DO j=MAX(j_p-1,1),MIN(j_p,N_j) DO i=MAX(i_p-1,1),MIN(i_p,N_i) Interpolant=Interpolant+f1(i-i_p+2)*f2(j-j_p+2)*table(i,j) ENDDO canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Get interval in labels around the point CALL Get_interval(labels1,point(1),i_p) CALL Get_interval(labels2,point(2),j_p) !Interpolate w.r.t labels 1 DO j=1,2 CALL Get_points_and_weights(labels1,table(:,j_p(j)),point(1),f,t,canExtrapolate) Interp1D(j)=f(1)*t(1)+f(2)*t(2) ENDDO !Interpolate w.r.t labels 2 CALL Get_points_and_weights(labels2(j_p(1):j_p(2)),Interp1D,point(2),f,t,canExtrapolate) Interpolant=f(1)*t(1)+f(2)*t(2) ENDFUNCTION Interp_2D ! Loading @@ -115,16 +125,19 @@ ENDFUNCTION Interp_2D !> @param table contains data to be interpolated between defined at points !> given in labels1, labels2, and labels 3 !> @param point point at which interpolation needs to be done to return a value !> @param lextrap logical indicating if extrapolation should be performed if called for !> @returns Interpolant interpolated data found through linear interpolation to the point point !> FUNCTION Interp_3D(labels1,labels2,labels3,table,point) RESULT(Interpolant) FUNCTION Interp_3D(labels1,labels2,labels3,table,point,lextrap) RESULT(Interpolant) REAL(SRK),INTENT(IN) :: labels1(:) REAL(SRK),INTENT(IN) :: labels2(:) REAL(SRK),INTENT(IN) :: labels3(:) REAL(SRK),INTENT(IN) :: table(:,:,:) REAL(SRK),INTENT(IN) :: point(:) REAL(SRK) :: Interpolant,f1(2),f2(2),f3(2) INTEGER(SIK) :: i,j,k,N_i,N_j,N_k,i_p,j_p,k_p LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK) :: Interpolant,f(2),t(2),Interp2D(2,2),Interp1D(2) INTEGER(SIK) :: j,k,i_p(2),j_p(2),k_p(2) LOGICAL(SBK) :: canExtrapolate REQUIRE(SIZE(labels1) > 1) REQUIRE(SIZE(labels2) > 1) Loading @@ -136,41 +149,58 @@ FUNCTION Interp_3D(labels1,labels2,labels3,table,point) RESULT(Interpolant) REQUIRE(isStrictlyIncDec(labels2)) REQUIRE(isStrictlyIncDec(labels3)) CALL Get_points_and_weights(labels1,point(1),f1,i_p,N_i) CALL Get_points_and_weights(labels2,point(2),f2,j_p,N_j) CALL Get_points_and_weights(labels3,point(3),f3,k_p,N_k) !Get interpolated value Interpolant=0.0_SRK DO k=MAX(k_p-1,1),MIN(k_p,N_k) DO j=MAX(j_p-1,1),MIN(j_p,N_j) DO i=MAX(i_p-1,1),MIN(i_p,N_i) Interpolant=Interpolant+f1(i-i_p+2)*f2(j-j_p+2)*f3(k-k_p+2)*table(i,j,k) canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Get interval in labels around the point CALL Get_interval(labels1,point(1),i_p) CALL Get_interval(labels2,point(2),j_p) CALL Get_interval(labels3,point(3),k_p) !Interpolate w.r.t labels 1 DO k=1,2 DO j=1,2 CALL Get_points_and_weights(labels1,table(:,j_p(j),k_p(k)),point(1),f,t,canExtrapolate) Interp2D(j,k)=f(1)*t(1)+f(2)*t(2) ENDDO ENDDO !Interpolate w.r.t labels 2 DO k=1,2 CALL Get_points_and_weights(labels2(j_p(1):j_p(2)),Interp2D(:,k),point(2),f,t,canExtrapolate) Interp1D(k)=f(1)*t(1)+f(2)*t(2) ENDDO !Interpolate w.r.t labels 3 CALL Get_points_and_weights(labels3(k_p(1):k_p(2)),Interp1D,point(3),f,t,canExtrapolate) Interpolant=f(1)*t(1)+f(2)*t(2) ENDFUNCTION Interp_3D ! !------------------------------------------------------------------------------- !> @brief This routine takes a 1D array of table labels (the points along a single !> dimension at which table entries are defined) and the interpolation point !> coordinate along that dimension and returns the index of the label entry !> directly ahead of that point along with the weights of that point and the !> preceeding point. !> coordinate along that dimension and returns the weights and table !> values for linear interpolation !> @param labels axis labels at which data points in table are defined !> @param table the 1-D table that contains values associated with labels !> @param point point at which interpolation needs to be done to return a value !> @param f weights associated with the points on either side of point !> @param i_p index of data point directly ahead of the interpolation point in labels !> @param N_i number of entries of labels !> @param t table associated with the table values on either side of point !> @param lextrap logical to perform extrapolation if neccessary !> SUBROUTINE Get_points_and_weights(labels,point,f,i_p,N_i) REAL(SRK),INTENT(IN) :: labels(:) SUBROUTINE Get_points_and_weights(labels,table,point,f,t,lextrap) REAL(SRK),INTENT(IN) :: labels(:),table(:) REAL(SRK),INTENT(IN) :: point REAL(SRK),INTENT(OUT) :: f(2) INTEGER(SIK),INTENT(OUT) :: i_p,N_i LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK),INTENT(OUT) :: f(2),t(2) INTEGER(SIK) :: i_p,N_i LOGICAL(SBK) :: canExtrapolate N_i=SIZE(labels) canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Determine labels index(s) and weights f(:)=1.0_SRK IF(labels(1) < labels(2)) THEN Loading @@ -178,16 +208,62 @@ SUBROUTINE Get_points_and_weights(labels,point,f,i_p,N_i) i_p=getFirstGreater(labels,point) IF(i_p > 1 .AND. i_p < N_i+1) THEN f(1)=(labels(i_p)-point)/(labels(i_p)-labels(i_p-1)) f(2)=1.0_SRK-f(1) t(1)=table(i_p-1) t(2)=table(i_p) ELSEIF(i_p == 1) THEN t(1)=table(1) t(2)=table(2) IF(canExtrapolate) f(1)=(labels(2)-point)/(labels(2)-labels(1)) ELSEIF(i_p == N_i+1) THEN t(1)=table(N_i) t(2)=table(N_i-1) IF(canExtrapolate) f(1)=(point-labels(N_i-1))/(labels(N_i)-labels(N_i-1)) ENDIF f(2)=1.0_SRK-f(1) ELSE !Descending order i_p=getFirstGreaterEqual(labels,point) IF(i_p > 1 .AND. i_p < N_i+1) THEN f(1)=(point-labels(i_p))/(labels(i_p-1)-labels(i_p)) f(2)=1.0_SRK-f(1) t(1)=table(i_p-1) t(2)=table(i_p) f(2)=(labels(i_p-1)-point)/(labels(i_p-1)-labels(i_p)) ELSEIF(i_p == 1) THEN t(1)=table(2) t(2)=table(1) IF(canExtrapolate) f(2)=(labels(2)-point)/(labels(2)-labels(1)) ELSEIF(i_p == N_i+1) THEN t(1)=table(N_i-1) t(2)=table(N_i) IF(canExtrapolate) f(2)=(point-labels(N_i-1))/(labels(N_i)-labels(N_i-1)) ENDIF f(1)=1.0_SRK-f(2) ENDIF ENDSUBROUTINE Get_points_and_weights ! !------------------------------------------------------------------------------- !> @brief This routine gets the interval in labels that contains point. !> If the point lays outside of labels then the closest interval !> in labels is returned. !> @param labels axis labels at which data points in table are defined !> @param point point at which interpolation needs to be done to return a value !> @param i_p index of the interval around point in labels !> SUBROUTINE Get_interval(labels,point,i_p) REAL(SRK),INTENT(IN) :: labels(:) REAL(SRK),INTENT(IN) :: point INTEGER(SIK),INTENT(OUT) :: i_p(2) IF(labels(1) < labels(2)) THEN !Ascending order i_p(2)=getFirstGreater(labels,point) ELSE !Descending order i_p(2)=getFirstGreaterEqual(labels,point) ENDIF IF(i_p(2) == 1) i_p(2)=i_p(2)+1 IF(i_p(2) > SIZE(labels)) i_p(2)=i_p(2)-1 i_p(1)=i_p(2)-1 ENDSUBROUTINE Get_interval ! ENDMODULE InterpolatorsModule src/ParameterLists.f90 +969 −516 File changed.Preview size limit exceeded, changes collapsed. Show changes Loading
src/CMakeLists.txt +1 −0 Original line number Diff line number Diff line Loading @@ -32,6 +32,7 @@ TRIBITS_ADD_LIBRARY(Utils Futility_DBC.f90 UnitTest.f90 IntrType.f90 Hash.f90 Search.f90 Constants_Conversion.f90 ExtendedMath.f90 Loading
src/Geom_Graph.f90 +227 −119 Original line number Diff line number Diff line Loading @@ -27,6 +27,7 @@ PRIVATE PUBLIC :: GraphType PUBLIC :: DAGraphType PUBLIC :: OPERATOR(==) PUBLIC :: ASSIGNMENT(=) !PUBLIC :: OPERATOR(+) INTEGER(SIK),PARAMETER :: GRAPH_NULL_EDGE=0 Loading @@ -35,6 +36,8 @@ INTEGER(SIK),PARAMETER :: GRAPH_QUADRATIC_EDGE=-1 !> @brief a Directed Acyclic Graph Type TYPE :: DAGraphType !> Logical indicating if the graph is sorted by ascending node value LOGICAL(SBK) :: sorted=.FALSE. !> The number of nodes on the graph INTEGER(SIK) :: n=0 !> Nodes (Size n) Loading @@ -50,9 +53,13 @@ TYPE :: DAGraphType !> @copybrief Geom_Graph::clear_DAGraphType !> @copydetails Geom_Graph::clear_DAGraphType PROCEDURE,PASS :: clear => clear_DAGraphType !> @copybrief Geom_Graph::defineEdge_DAGraphType !> @copydetails Geom_Graph::defineEdge_DAGraphType PROCEDURE,PASS :: insertNode=> insertNode_DAGraphType !> @copybrief Geom_Graph::insertNodeSorted_DAGraphType !> @copydetails Geom_Graph::insertNodeSorted_DAGraphType PROCEDURE,PASS,PRIVATE :: insertNodeSorted => insertNodeSorted_DAGraphType !> @copybrief Geom_Graph::insertNodeAtIndex_DAGraphType !> @copydetails Geom_Graph::insertNodeAtIndex_DAGraphType PROCEDURE,PASS,PRIVATE :: insertNodeAtIndex => insertNodeAtIndex_DAGraphType GENERIC :: insertNode => insertNodeSorted,insertNodeAtIndex !> @copybrief Geom_Graph::removeNode_DAGraphType !> @copydetails Geom_Graph::removeNode_DAGraphType PROCEDURE,PASS :: removeNode=> removeNode_DAGraphType Loading Loading @@ -225,6 +232,12 @@ INTERFACE OPERATOR(==) MODULE PROCEDURE isequal_graphType ENDINTERFACE INTERFACE ASSIGNMENT(=) !> @copybrief Geom_Graph::assign_DAGraphType !> @copydetails Geom_Graph::assign_DAGraphType MODULE PROCEDURE assign_DAGraphType ENDINTERFACE !INTERFACE OPERATOR(+) ! !> @copybrief Geom_Graph::add_graphType ! !> @copydetails Geom_Graph::add_graphType Loading Loading @@ -2209,19 +2222,21 @@ ENDFUNCTION isequal_GraphType SUBROUTINE init_DAGraphType(this,n,nodes) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: n INTEGER(SIK),ALLOCATABLE,INTENT(IN) :: nodes(:) INTEGER(SIK),INTENT(IN) :: nodes(:) REQUIRE(.NOT.ALLOCATED(this%nodes)) REQUIRE(SIZE(nodes) == n .OR. n == 0) IF(.NOT. ALLOCATED(this%nodes) .AND. ALLOCATED(nodes)) THEN IF(n > 0) THEN IF(SIZE(nodes) == n) THEN this%n=n ALLOCATE(this%nodes(n)) ALLOCATE(this%edgeMatrix(n,n)) IF(n > 0) THEN this%nodes=nodes this%edgeMatrix=0 ENDIF ENDIF CALL sort(this%nodes) this%sorted=.TRUE. ENDIF this%edgeMatrix=0 ENDSUBROUTINE init_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2234,35 +2249,60 @@ SUBROUTINE clear_DAGraphType(this) this%n=0 IF(ALLOCATED(this%nodes)) DEALLOCATE(this%nodes) IF(ALLOCATED(this%edgeMatrix)) DEALLOCATE(this%edgeMatrix) this%sorted=.FALSE. ENDSUBROUTINE clear_DAGraphType ! !------------------------------------------------------------------------------- !> @brief !> @param !> SUBROUTINE insertNode_DAGraphType(this,ID,ind) SUBROUTINE insertNodeSorted_DAGraphType(this,ID) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ID INTEGER(SIK),INTENT(IN),OPTIONAL :: ind INTEGER(SIK) :: index INTEGER(SIK) :: index,a,b INTEGER(SIK),ALLOCATABLE :: tmpNodes(:),tmpEdges(:,:) index=1 IF(PRESENT(ind)) index=ind IF(ALLOCATED(this%edgeMatrix)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) IF(ALL(this%nodes /= ID)) THEN !Store values temporarily ALLOCATE(tmpNodes(this%n)) ALLOCATE(tmpEdges(this%n,this%n)) tmpNodes=this%nodes tmpEdges=this%edgeMatrix CALL MOVE_ALLOC(this%nodes,tmpNodes) CALL MOVE_ALLOC(this%edgeMatrix,tmpEdges) !Resize arrays on graph type DEALLOCATE(this%nodes) DEALLOCATE(this%edgeMatrix) this%n=this%n+1 ALLOCATE(this%nodes(this%n)) ALLOCATE(this%edgeMatrix(this%n,this%n)) this%edgeMatrix=0 !Find the location for this index. If the graph is sorted, then we'll insert !the value at the location that maintains the sorting. Otherwise, we'll just !stick it at the beginning !Also handle some special cases up front to simplify the binary search for sorted insertion IF(this%n == 0) THEN index = 1 ELSEIF(this%sorted) THEN a = 1 b = this%n IF(ID < this%nodes(a)) THEN index = 1 ELSEIF(ID > this%nodes(b)) THEN index = this%n ELSE index = (a + b)/2 DO WHILE(.TRUE.) IF(ID < this%nodes(index-1)) THEN b = index ELSEIF(ID > this%nodes(index)) THEN a = index ELSE EXIT ENDIF index = (a + b)/2 ENDDO ENDIF ELSE index = 1 ENDIF !reassign old data and assign new node this%nodes(1:index-1)=tmpNodes(1:index-1) this%nodes(index)=ID Loading @@ -2271,17 +2311,55 @@ SUBROUTINE insertNode_DAGraphType(this,ID,ind) this%edgeMatrix(index+1:,1:index-1)=tmpEdges(index:,1:index-1) this%edgeMatrix(1:index-1,index+1:)=tmpEdges(1:index-1,index:) this%edgeMatrix(index+1:,index+1:)=tmpEdges(index:,index:) DEALLOCATE(tmpNodes) DEALLOCATE(tmpEdges) ENDIF ELSE ALLOCATE(this%edgeMatrix(1,1)) ALLOCATE(this%nodes(1)) ENDSUBROUTINE insertNodeSorted_DAGraphType ! !------------------------------------------------------------------------------- !> @brief !> @param !> SUBROUTINE insertNodeAtIndex_DAGraphType(this,ID,index) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ID INTEGER(SIK),INTENT(IN) :: index INTEGER(SIK),ALLOCATABLE :: tmpNodes(:),tmpEdges(:,:) REQUIRE(ALLOCATED(this%edgeMatrix)) IF(ALL(this%nodes /= ID)) THEN !Store values temporarily CALL MOVE_ALLOC(this%nodes,tmpNodes) CALL MOVE_ALLOC(this%edgeMatrix,tmpEdges) !Resize arrays on graph type this%n=this%n+1 ALLOCATE(this%nodes(this%n)) ALLOCATE(this%edgeMatrix(this%n,this%n)) this%edgeMatrix=0 this%nodes=ID this%n=1 !reassign old data and assign new node this%nodes(1:index-1)=tmpNodes(1:index-1) this%nodes(index)=ID this%nodes(index+1:this%n)=tmpNodes(index:) this%edgeMatrix(1:index-1,1:index-1)=tmpEdges(1:index-1,1:index-1) this%edgeMatrix(index+1:,1:index-1)=tmpEdges(index:,1:index-1) this%edgeMatrix(1:index-1,index+1:)=tmpEdges(1:index-1,index:) this%edgeMatrix(index+1:,index+1:)=tmpEdges(index:,index:) !Update the sorted status IF(this%sorted) THEN IF(index > 1) THEN IF(this%nodes(index) <= this%nodes(index-1)) THEN this%sorted=.FALSE. ENDIF ENDIF IF(index < this%n) THEN IF(this%nodes(index) >= this%nodes(index+1)) THEN this%sorted=.FALSE. ENDIF ENDIF ENDIF ENDSUBROUTINE insertNode_DAGraphType ENDIF ENDSUBROUTINE insertNodeAtIndex_DAGraphType ! !------------------------------------------------------------------------------- !> @brief Loading @@ -2291,17 +2369,11 @@ SUBROUTINE removeNode_DAGraphType(this,ID,ind) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN),OPTIONAL :: ID INTEGER(SIK),INTENT(IN),OPTIONAL :: ind INTEGER(SIK) :: i,index INTEGER(SIK) :: index index=0 IF(PRESENT(ID)) THEN !Find the index to remove DO i=1,this%n IF(ID == this%nodes(i)) THEN index=i EXIT ENDIF ENDDO index = this%getIndex(ID) !If the ID was found, remove the index. IF(index > 0) CALL removeNodeByIndex_DAGraphType(this,index) ELSEIF(PRESENT(ind)) THEN Loading @@ -2316,24 +2388,20 @@ ENDSUBROUTINE removeNode_DAGraphType SUBROUTINE removeNodeByIndex_DAGraphType(this,ind) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ind INTEGER(SIK) :: i INTEGER(SIK),ALLOCATABLE :: tmpN(:),tmpEdge(:,:) IF(ALLOCATED(this%edgeMatrix)) THEN IF((1 <= ind) .AND. (ind <= this%n)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) REQUIRE(1 <= ind) REQUIRE(ind <= this%n) !Store values temporarily i=this%n ALLOCATE(tmpN(i)) ALLOCATE(tmpEdge(i,i)) tmpN=this%nodes tmpEdge=this%edgeMatrix CALL MOVE_ALLOC(this%nodes,tmpN) CALL MOVE_ALLOC(this%edgeMatrix,tmpEdge) !Resize arrays on graph type DEALLOCATE(this%nodes) DEALLOCATE(this%edgeMatrix) this%n=this%n-1 IF(this%n > 0) THEN ALLOCATE(this%nodes(this%n)) ALLOCATE(this%edgeMatrix(this%n,this%n)) IF(this%n > 0) THEN this%edgeMatrix=0 !reassign old data and assign new node this%nodes(1:ind-1)=tmpN(1:ind-1) Loading @@ -2345,10 +2413,7 @@ SUBROUTINE removeNodeByIndex_DAGraphType(this,ind) this%edgeMatrix(ind:,ind:)=tmpEdge(ind+1:,ind+1:) ENDIF ENDIF !CALL demalloc(tmpint) DEALLOCATE(tmpEdge) ENDIF ENDIF ENDSUBROUTINE removeNodeByIndex_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2365,9 +2430,9 @@ FUNCTION isStartNode_DAGraphType(this,ID,ind) RESULT(bool) bool=.FALSE. IF(PRESENT(ID)) THEN index=this%getIndex(ID) IF(ALL(this%edgeMatrix(:,index) == 0)) bool=.TRUE. bool = (ALL(this%edgeMatrix(:,index) == 0)) ELSEIF(PRESENT(ind)) THEN IF(ALL(this%edgeMatrix(:,ind) == 0)) bool=.TRUE. bool = (ALL(this%edgeMatrix(:,ind) == 0)) ENDIF ENDFUNCTION isStartNode_DAGraphType ! Loading @@ -2383,11 +2448,13 @@ SUBROUTINE getNextStartNode_DAGraphType(this,old,ID) ID=0 DO i=1,this%n IF(this%isStartNode(IND=i) .AND. ALL(this%nodes(i) /= old)) THEN IF(this%isStartNode(IND=i)) THEN IF(ALL(this%nodes(i) /= old)) THEN ID=this%nodes(i) old(i)=ID EXIT ENDIF ENDIF ENDDO ENDSUBROUTINE getNextStartNode_DAGraphType ! Loading @@ -2401,11 +2468,12 @@ SUBROUTINE defineEdge_DAGraphType(this,fromID,toID) INTEGER(SIK),INTENT(IN) :: toID INTEGER(SIK) :: fromInd,toInd IF(ALLOCATED(this%edgeMatrix)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) fromInd=this%getIndex(fromID) toInd=this%getIndex(toID) IF((fromInd > 0) .AND. (toInd > 0)) this%edgeMatrix(fromInd,toInd)=1 ENDIF ENDSUBROUTINE defineEdge_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2418,31 +2486,67 @@ SUBROUTINE removeEdge_DAGraphType(this,fromID,toID) INTEGER(SIK),INTENT(IN) :: toID INTEGER(SIK) :: fromInd,toInd IF(ALLOCATED(this%edgeMatrix)) THEN REQUIRE(ALLOCATED(this%edgeMatrix)) fromInd=this%getIndex(fromID) toInd=this%getIndex(toID) IF((fromInd > 0) .AND. (toInd > 0)) this%edgeMatrix(fromInd,toInd)=0 ENDIF ENDSUBROUTINE removeEdge_DAGraphType ! !------------------------------------------------------------------------------- !> @brief !> @param !> @brief retrieves the index of a node in a @c DAGraphType !> @param this the graph object !> @param ID the node ID to search for !> @returns ind the index of node @c ID !> !> Returns 0 if the ID could not be found. If the graph is sorted, a binary !> search is used. Otherwise, a linear search is used. !> FUNCTION getIndex_DAGraphType(this,ID) RESULT(ind) CLASS(DAGraphType),INTENT(INOUT) :: this INTEGER(SIK),INTENT(IN) :: ID INTEGER(SIK) :: i,ind INTEGER(SIK) :: ind ! INTEGER(SIK) :: a,b REQUIRE(ALLOCATED(this%nodes)) ind=0 IF(ALLOCATED(this%edgeMatrix)) THEN DO i=1,this%n IF(this%nodes(i) == ID) THEN ind=i IF(this%sorted) THEN a = 1 b = this%n IF(ID == this%nodes(a)) THEN ind = 1 ELSEIF(ID == this%nodes(b)) THEN ind = this%n ELSE ind = (a + b)/2 DO WHILE(.TRUE.) IF(ID < this%nodes(ind)) THEN b = ind ELSEIF(ID > this%nodes(ind)) THEN a = ind ELSE EXIT ENDIF IF(a == b) EXIT ind = (a + b)/2 ENDDO IF(this%nodes(ind) /= ID) THEN ind = 0 ENDIF ENDIF ELSE DO a=1,this%n IF(this%nodes(a) == ID) THEN ind=a EXIT ENDIF ENDDO IF(a > this%n) ind=0 ENDIF ENDFUNCTION getIndex_DAGraphType ! !------------------------------------------------------------------------------- Loading @@ -2459,18 +2563,22 @@ SUBROUTINE KATS_DAGraphType(this) ALLOCATE(oldIDs(this%n)) oldIDs=0 CALL this%getNextStartNode(oldIDs,ID) sortedGraph%n=0 ALLOCATE(sortedGraph%nodes(sortedGraph%n)) ALLOCATE(sortedGraph%edgeMatrix(sortedGraph%n,sortedGraph%n)) DO WHILE(ID /= 0) !Add the node to the sorted graph CALL sortedGraph%insertNode(ID,1) !Remove all edges from that node DO i=1,this%n CALL this%removeEdge(ID,this%nodes(i)) ENDDO !Get the next "start node CALL this%getNextStartNode(oldIDs,ID) ENDDO IF(ALL(this%edgeMatrix == 0)) THEN SELECTTYPE(this); TYPE IS(DAGraphtype) this=sortedGraph ENDSELECT ENDIF DEALLOCATE(oldIDs) ENDSUBROUTINE KATS_DAGraphType Loading @@ -2482,13 +2590,13 @@ ENDSUBROUTINE KATS_DAGraphType SUBROUTINE assign_DAGraphType(g0,g1) CLASS(DAGraphType),INTENT(INOUT) :: g0 CLASS(DAGraphType),INTENT(IN) :: g1 CALL clear_DAGraphType(g0) IF(g1%n > 0) THEN CALL g0%clear() g0%n=g1%n ALLOCATE(g0%nodes(g1%n)) ALLOCATE(g0%edgeMatrix(g1%n,g1%n)) g0%nodes=g1%nodes g0%edgeMatrix=g1%edgeMatrix IF(g1%n > 0) THEN g0%nodes(:)=g1%nodes g0%edgeMatrix(:,:)=g1%edgeMatrix ENDIF ENDSUBROUTINE assign_DAGraphType ! Loading
src/Hash.f90 0 → 100644 +53 −0 Original line number Diff line number Diff line !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! Futility Development Group ! ! All rights reserved. ! ! ! ! Futility is a jointly-maintained, open-source project between the University ! ! of Michigan and Oak Ridge National Laboratory. The copyright and license ! ! can be found in LICENSE.txt in the head directory of this repository. ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! !> @brief Provides hashing algorithms !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! MODULE HashModule USE IntrType IMPLICIT NONE PRIVATE PUBLIC :: stringHash CONTAINS ! !------------------------------------------------------------------------------- !> @brief Hashes a string to an integer value !> @param string the string to hash !> @param p the rolling polynomial base !> @param m the modulus !> !> This algorithm is based on https://cp-algorithms.com/string/string-hashing.html !> !> It allows for ASCII values 32 and up. Recommended value for p is 97 (a prime !> number close to 95, which is the number of "normal" ASCII values expected in the !> 32-126 range). Recommended value for m is a large prime number such as 10**9 + 9. !> However, other arguments may be valid for certain applications. !> PURE FUNCTION stringHash(string,p,m) RESULT(hash) CHARACTER(LEN=*),INTENT(IN) :: string INTEGER(SLK),INTENT(IN) :: p INTEGER(SLK),INTENT(IN) :: m INTEGER(SLK) :: hash ! INTEGER(SIK) :: i INTEGER(SLK) :: p_pow hash = 0 p_pow = 1 DO i=1,LEN(string) !1-31 are special characters that we disallow hash = MOD(hash + (IACHAR(string(i:i))-31_SIK)*p_pow,m) p_pow = MOD(p_pow * p,m) ENDDO !i ENDFUNCTION stringHash ! ENDMODULE HashModule No newline at end of file
src/Interpolators.f90 +124 −48 Original line number Diff line number Diff line Loading @@ -47,25 +47,27 @@ CONTAINS !> @param table contains data to be interpolated between defined at points !> given in labels !> @param point point at which interpolation needs to be done to return a value !> @param lextrap logical indicating if extrapolation should be performed if called for !> @returns Interpolant interpolated data found through linear interpolation to the point point FUNCTION Interp_1D(labels,table,point) RESULT(Interpolant) FUNCTION Interp_1D(labels,table,point,lextrap) RESULT(Interpolant) REAL(SRK),INTENT(IN) :: labels(:) REAL(SRK),INTENT(IN) :: table(:) REAL(SRK),INTENT(IN) :: point REAL(SRK) :: Interpolant,f(2) INTEGER(SIK) :: i,N_i,i_p LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK) :: Interpolant,f(2),t(2) LOGICAL(SBK) :: canExtrapolate REQUIRE(SIZE(labels) > 1) REQUIRE(SIZE(labels) == SIZE(table)) REQUIRE(isStrictlyIncDec(labels)) CALL Get_points_and_weights(labels,point,f,i_p,N_i) canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Get interpolated value Interpolant=0.0_SRK DO i=MAX(i_p-1,1),MIN(i_p,N_i) Interpolant=Interpolant+f(i-i_p+2)*table(i) ENDDO CALL Get_points_and_weights(labels,table,point,f,t,canExtrapolate) Interpolant=f(1)*t(1)+f(2)*t(2) ENDFUNCTION Interp_1D ! Loading @@ -77,14 +79,17 @@ ENDFUNCTION Interp_1D !> @param table contains data to be interpolated between defined at points !> given in labels1 and labels2 !> @param point point at which interpolation needs to be done to return a value !> @param lextrap logical indicating if extrapolation should be performed if called for !> @returns Interpolant interpolated data found through linear interpolation to the point point FUNCTION Interp_2D(labels1,labels2,table,point) RESULT(Interpolant) FUNCTION Interp_2D(labels1,labels2,table,point,lextrap) RESULT(Interpolant) REAL(SRK),INTENT(IN) :: labels1(:) REAL(SRK),INTENT(IN) :: labels2(:) REAL(SRK),INTENT(IN) :: table(:,:) REAL(SRK),INTENT(IN) :: point(:) REAL(SRK) :: Interpolant,f1(2),f2(2) INTEGER(SIK) :: i,j,N_i,N_j,i_p,j_p LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK) :: Interpolant,f(2),t(2),Interp1D(2) INTEGER(SIK) :: j,i_p(2),j_p(2) LOGICAL(SBK) :: canExtrapolate REQUIRE(SIZE(labels1) > 1) REQUIRE(SIZE(labels2) > 1) Loading @@ -93,16 +98,21 @@ FUNCTION Interp_2D(labels1,labels2,table,point) RESULT(Interpolant) REQUIRE(isStrictlyIncDec(labels1)) REQUIRE(isStrictlyIncDec(labels2)) CALL Get_points_and_weights(labels1,point(1),f1,i_p,N_i) CALL Get_points_and_weights(labels2,point(2),f2,j_p,N_j) !Get interpolated value Interpolant=0.0_SRK DO j=MAX(j_p-1,1),MIN(j_p,N_j) DO i=MAX(i_p-1,1),MIN(i_p,N_i) Interpolant=Interpolant+f1(i-i_p+2)*f2(j-j_p+2)*table(i,j) ENDDO canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Get interval in labels around the point CALL Get_interval(labels1,point(1),i_p) CALL Get_interval(labels2,point(2),j_p) !Interpolate w.r.t labels 1 DO j=1,2 CALL Get_points_and_weights(labels1,table(:,j_p(j)),point(1),f,t,canExtrapolate) Interp1D(j)=f(1)*t(1)+f(2)*t(2) ENDDO !Interpolate w.r.t labels 2 CALL Get_points_and_weights(labels2(j_p(1):j_p(2)),Interp1D,point(2),f,t,canExtrapolate) Interpolant=f(1)*t(1)+f(2)*t(2) ENDFUNCTION Interp_2D ! Loading @@ -115,16 +125,19 @@ ENDFUNCTION Interp_2D !> @param table contains data to be interpolated between defined at points !> given in labels1, labels2, and labels 3 !> @param point point at which interpolation needs to be done to return a value !> @param lextrap logical indicating if extrapolation should be performed if called for !> @returns Interpolant interpolated data found through linear interpolation to the point point !> FUNCTION Interp_3D(labels1,labels2,labels3,table,point) RESULT(Interpolant) FUNCTION Interp_3D(labels1,labels2,labels3,table,point,lextrap) RESULT(Interpolant) REAL(SRK),INTENT(IN) :: labels1(:) REAL(SRK),INTENT(IN) :: labels2(:) REAL(SRK),INTENT(IN) :: labels3(:) REAL(SRK),INTENT(IN) :: table(:,:,:) REAL(SRK),INTENT(IN) :: point(:) REAL(SRK) :: Interpolant,f1(2),f2(2),f3(2) INTEGER(SIK) :: i,j,k,N_i,N_j,N_k,i_p,j_p,k_p LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK) :: Interpolant,f(2),t(2),Interp2D(2,2),Interp1D(2) INTEGER(SIK) :: j,k,i_p(2),j_p(2),k_p(2) LOGICAL(SBK) :: canExtrapolate REQUIRE(SIZE(labels1) > 1) REQUIRE(SIZE(labels2) > 1) Loading @@ -136,41 +149,58 @@ FUNCTION Interp_3D(labels1,labels2,labels3,table,point) RESULT(Interpolant) REQUIRE(isStrictlyIncDec(labels2)) REQUIRE(isStrictlyIncDec(labels3)) CALL Get_points_and_weights(labels1,point(1),f1,i_p,N_i) CALL Get_points_and_weights(labels2,point(2),f2,j_p,N_j) CALL Get_points_and_weights(labels3,point(3),f3,k_p,N_k) !Get interpolated value Interpolant=0.0_SRK DO k=MAX(k_p-1,1),MIN(k_p,N_k) DO j=MAX(j_p-1,1),MIN(j_p,N_j) DO i=MAX(i_p-1,1),MIN(i_p,N_i) Interpolant=Interpolant+f1(i-i_p+2)*f2(j-j_p+2)*f3(k-k_p+2)*table(i,j,k) canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Get interval in labels around the point CALL Get_interval(labels1,point(1),i_p) CALL Get_interval(labels2,point(2),j_p) CALL Get_interval(labels3,point(3),k_p) !Interpolate w.r.t labels 1 DO k=1,2 DO j=1,2 CALL Get_points_and_weights(labels1,table(:,j_p(j),k_p(k)),point(1),f,t,canExtrapolate) Interp2D(j,k)=f(1)*t(1)+f(2)*t(2) ENDDO ENDDO !Interpolate w.r.t labels 2 DO k=1,2 CALL Get_points_and_weights(labels2(j_p(1):j_p(2)),Interp2D(:,k),point(2),f,t,canExtrapolate) Interp1D(k)=f(1)*t(1)+f(2)*t(2) ENDDO !Interpolate w.r.t labels 3 CALL Get_points_and_weights(labels3(k_p(1):k_p(2)),Interp1D,point(3),f,t,canExtrapolate) Interpolant=f(1)*t(1)+f(2)*t(2) ENDFUNCTION Interp_3D ! !------------------------------------------------------------------------------- !> @brief This routine takes a 1D array of table labels (the points along a single !> dimension at which table entries are defined) and the interpolation point !> coordinate along that dimension and returns the index of the label entry !> directly ahead of that point along with the weights of that point and the !> preceeding point. !> coordinate along that dimension and returns the weights and table !> values for linear interpolation !> @param labels axis labels at which data points in table are defined !> @param table the 1-D table that contains values associated with labels !> @param point point at which interpolation needs to be done to return a value !> @param f weights associated with the points on either side of point !> @param i_p index of data point directly ahead of the interpolation point in labels !> @param N_i number of entries of labels !> @param t table associated with the table values on either side of point !> @param lextrap logical to perform extrapolation if neccessary !> SUBROUTINE Get_points_and_weights(labels,point,f,i_p,N_i) REAL(SRK),INTENT(IN) :: labels(:) SUBROUTINE Get_points_and_weights(labels,table,point,f,t,lextrap) REAL(SRK),INTENT(IN) :: labels(:),table(:) REAL(SRK),INTENT(IN) :: point REAL(SRK),INTENT(OUT) :: f(2) INTEGER(SIK),INTENT(OUT) :: i_p,N_i LOGICAL(SBK),INTENT(IN),OPTIONAL :: lextrap REAL(SRK),INTENT(OUT) :: f(2),t(2) INTEGER(SIK) :: i_p,N_i LOGICAL(SBK) :: canExtrapolate N_i=SIZE(labels) canExtrapolate=.FALSE. IF(PRESENT(lextrap)) THEN canExtrapolate=lextrap ENDIF !Determine labels index(s) and weights f(:)=1.0_SRK IF(labels(1) < labels(2)) THEN Loading @@ -178,16 +208,62 @@ SUBROUTINE Get_points_and_weights(labels,point,f,i_p,N_i) i_p=getFirstGreater(labels,point) IF(i_p > 1 .AND. i_p < N_i+1) THEN f(1)=(labels(i_p)-point)/(labels(i_p)-labels(i_p-1)) f(2)=1.0_SRK-f(1) t(1)=table(i_p-1) t(2)=table(i_p) ELSEIF(i_p == 1) THEN t(1)=table(1) t(2)=table(2) IF(canExtrapolate) f(1)=(labels(2)-point)/(labels(2)-labels(1)) ELSEIF(i_p == N_i+1) THEN t(1)=table(N_i) t(2)=table(N_i-1) IF(canExtrapolate) f(1)=(point-labels(N_i-1))/(labels(N_i)-labels(N_i-1)) ENDIF f(2)=1.0_SRK-f(1) ELSE !Descending order i_p=getFirstGreaterEqual(labels,point) IF(i_p > 1 .AND. i_p < N_i+1) THEN f(1)=(point-labels(i_p))/(labels(i_p-1)-labels(i_p)) f(2)=1.0_SRK-f(1) t(1)=table(i_p-1) t(2)=table(i_p) f(2)=(labels(i_p-1)-point)/(labels(i_p-1)-labels(i_p)) ELSEIF(i_p == 1) THEN t(1)=table(2) t(2)=table(1) IF(canExtrapolate) f(2)=(labels(2)-point)/(labels(2)-labels(1)) ELSEIF(i_p == N_i+1) THEN t(1)=table(N_i-1) t(2)=table(N_i) IF(canExtrapolate) f(2)=(point-labels(N_i-1))/(labels(N_i)-labels(N_i-1)) ENDIF f(1)=1.0_SRK-f(2) ENDIF ENDSUBROUTINE Get_points_and_weights ! !------------------------------------------------------------------------------- !> @brief This routine gets the interval in labels that contains point. !> If the point lays outside of labels then the closest interval !> in labels is returned. !> @param labels axis labels at which data points in table are defined !> @param point point at which interpolation needs to be done to return a value !> @param i_p index of the interval around point in labels !> SUBROUTINE Get_interval(labels,point,i_p) REAL(SRK),INTENT(IN) :: labels(:) REAL(SRK),INTENT(IN) :: point INTEGER(SIK),INTENT(OUT) :: i_p(2) IF(labels(1) < labels(2)) THEN !Ascending order i_p(2)=getFirstGreater(labels,point) ELSE !Descending order i_p(2)=getFirstGreaterEqual(labels,point) ENDIF IF(i_p(2) == 1) i_p(2)=i_p(2)+1 IF(i_p(2) > SIZE(labels)) i_p(2)=i_p(2)-1 i_p(1)=i_p(2)-1 ENDSUBROUTINE Get_interval ! ENDMODULE InterpolatorsModule
src/ParameterLists.f90 +969 −516 File changed.Preview size limit exceeded, changes collapsed. Show changes
