module MathClass use fMathClass implicit none interface fstring module procedure fstring_Int, fstring_Real, fstring_Int_len, fstring_Real_len end interface fstring interface input module procedure input_Int,input_Real,input_IntVec,input_RealVec,input_IntArray,input_RealArray,input_String,input_logical end interface input interface zeroif module procedure zeroif_Int,zeroif_Real end interface zeroif interface removeWord module procedure removeWord_String end interface contains !######################################## function norm(vec) result(a) real(8),intent(in)::vec(:) integer :: n real(8) :: a n=size(vec) a=dsqrt(dot_product(vec,vec) ) end function !######################################## !######################################## function SearchNearestCoord(Array,x) result(id) real(8),intent(in) :: Array(:,:) real(8),intent(in) :: x(:) integer,allocatable::xr(:) integer :: i,id,n,m,norm,tr_norm n=size(Array,1) m=size(Array,2) if(m/=size(x) )then stop "ERROR :: SearchNearestCoord >> size(Array,2) should be =size(x)" endif allocate(xr(m) ) do i=1,n xr(:)=Array(i,:) tr_norm=dot_product(xr-x,xr-x) if(i==1)then norm=tr_norm id =i else if(norm > tr_norm)then norm=tr_norm id =i else cycle endif endif enddo end function !######################################## !################################################## function SearchIDIntVec(Vec,val) result(id_) integer,intent(in) :: Vec(:) integer,intent(in) :: val integer :: i,id_ do i=1,size(Vec) if(Vec(i)==val )then id_=i return endif enddo end function !################################################## subroutine heapsort(n,array) integer,intent(in) :: n integer,intent(inout) :: array(1:n) integer ::i,k,j,l integer :: t if(n.le.0)then write(6,*)"Error, at heapsort"; stop endif if(n.eq.1)return l=n/2+1 k=n do while(k.ne.1) if(l.gt.1)then l=l-1 t=array(L) else t=array(k) array(k)=array(1) k=k-1 if(k.eq.1) then array(1)=t exit endif endif i=l j=l+l do while(j.le.k) if(j.lt.k)then if(array(j).lt.array(j+1))j=j+1 endif if (t.lt.array(j))then array(i)=array(j) i=j j=j+j else j=k+1 endif enddo array(i)=t enddo return end subroutine heapsort !========================================================== !calculate cross product !--------------------------- function cross_product(a,b) result (c) real(8), intent(in) :: a(:),b(:) real(8), allocatable :: c(:) if(size(a) /= size(b)) then stop "wrong number on size a, b" endif allocate(c(size(a,1))) if(size(c,1)==3) then c(1) = a(2)*b(3) - a(3)*b(2) c(2) = a(3)*b(1) - a(1)*b(3) c(3) = a(1)*b(2) - a(2)*b(1) else stop "wrong number at cross_product" endif end function cross_product !========================================================= !calculate diadic !---------------------- function diadic(a,b) result(c) real(8), intent(in) :: a(:), b(:) real(8), allocatable :: c(:,:) integer n,i,j allocate(c(size(a),size(b) ) ) do i=1,size(a) do j=1,size(b) c(i,j)=a(i)*b(j) enddo enddo end function diadic !========================================================== !calculate gz !-------------- subroutine calcgz(x2,x11,x12,nod_coord,gzi) real(8), intent(in) :: nod_coord(:,:) real(8),intent(out) :: gzi integer,intent(in):: x2,x11,x12 real(8) l real(8),allocatable::avec(:) allocate(avec(2)) l = dot_product( nod_coord(x12,1:2) - nod_coord(x11,1:2), & nod_coord(x12,1:2) - nod_coord(x11,1:2) ) l=l**(1.0d0/2.0d0) avec(1:2) = ( nod_coord(x12,1:2) - nod_coord(x11,1:2) )/l if(l==0.0d0)then print *, "calcgz l=0" gzi=0.0d0 else gzi=1.0d0/l*dot_product( nod_coord(x2,1:2) -nod_coord(x11,1:2),avec(1:2) ) endif deallocate(avec) end subroutine calcgz !========================================================== subroutine eigen_2d(Amat,eigenvector) real(8),intent(in)::Amat(:,:) real(8),intent(inout)::eigenvector(:,:) real(8)::b,c,phy,eigenvalue(2) integer i,j eigenvalue(:)=0.0d0 eigenvector(:,:)=0.0d0 b=-1.0d0*(Amat(1,1)+Amat(2,2)) c=Amat(1,1)*Amat(2,2)-Amat(1,2)*Amat(1,2) if(Amat(1,2)/=Amat(2,1) )then stop "input matrice is not symmetric" endif do i=1,2 eigenvalue(i)=(-1.0d0*b+((-1.0d0)**dble(i))*(b*b-4.0d0*c)**(1.0d0/2.0d0))*(0.50d0) enddo do i=1,2 if(Amat(1,2)==0 )then cycle elseif(Amat(1,2)/=0 )then phy=atan( (eigenvalue(i)-Amat(1,1))/Amat(1,2) ) do j=1,2 eigenvector(i,1:2)=(/cos(phy),sin(phy)/) enddo endif enddo do i=1,2 eigenvector(i,:)=eigenvalue(i)*eigenvector(i,:) enddo end subroutine eigen_2d !========================================================== function signmm(a) result(b) real(8),intent(in)::a real(8) b if(a>0)then b=1.0d0 elseif(a<0)then b=-1.0d0 elseif(a==0)then b=0.0d0 else stop "ERROR: Invalid Real(8) in function_signm" endif end function signmm !========================================================== recursive function det_mat(a,n) result(det) integer, intent(in) :: n real(8), intent(in) :: a(n, n) real(8) det, b(n-1, n-1) integer i if (n > 1) then det = 0.0d0 do i = 1, n b(1:i-1, 1:n-1) = a(1:i-1, 2:n) b(i:n-1, 1:n-1) = a(i+1:n, 2:n) det = det + (-1.0d0) ** (i + 1) & * a(i, 1) * det_mat(b, n-1) enddo else det = a(1,1) endif end function det_mat !===================================================================================== subroutine trans_rank_2(A,A_T) real(8),intent(in)::A(:,:) real(8),allocatable,intent(out)::A_T(:,:) integer n,m,i,j n=size(A,1) m=size(A,2) if(.not. allocated(A_T) )allocate(A_T(m,n)) do i=1,n do j=1, m A_T(j,i)=A(i,j) enddo enddo end subroutine trans_rank_2 !================================================================================== function trans1(A) result(A_T) real(8),intent(in)::A(:) real(8),allocatable::A_T(:,:) integer n,m,i,j n=size(A) if(.not. allocated(A_T) )allocate(A_T(1,n)) do i=1,n A_T(1,i)=A(i) enddo end function trans1 !================================================================================== function trans2(A) result(A_T) real(8),intent(in)::A(:,:) real(8),allocatable::A_T(:,:) integer n,m,i,j n=size(A,1) m=size(A,2) if(.not. allocated(A_T) )allocate(A_T(m,n)) do i=1,n do j=1, m A_T(j,i)=A(i,j) enddo enddo end function trans2 !================================================================================== subroutine inverse_rank_2(A,A_inv) real(8),intent(in)::A(:,:) real(8),allocatable::A_inv(:,:) real(8) detA,detA_1 integer m,n m=size(A,1) n=size(A,2) if(.not. allocated(A_inv) )allocate(A_inv(m,n)) detA=det_mat(A,n) if(detA==0.0d0) stop "ERROR: inverse, detA=0" detA_1=1.0d0/detA if(n==2)then A_inv(1,1)=detA_1*A(2,2) A_inv(1,2)=-detA_1*A(1,2) A_inv(2,1)=-detA_1*A(2,1) A_inv(2,2)=detA_1*A(1,1) elseif(n==3)then A_inv(1,1)=detA_1*(A(2,2)*A(3,3)-A(2,3)*A(3,2)) A_inv(1,2)=detA_1*(A(1,3)*A(3,2)-A(1,2)*A(3,3)) A_inv(1,3)=detA_1*(A(1,2)*A(2,3)-A(1,3)*A(2,2)) A_inv(2,1)=detA_1*(A(2,3)*A(3,1)-A(2,1)*A(3,3)) A_inv(2,2)=detA_1*(A(1,1)*A(3,3)-A(1,3)*A(3,1)) A_inv(2,3)=detA_1*(A(1,3)*A(2,1)-A(1,1)*A(2,3)) A_inv(3,1)=detA_1*(A(2,1)*A(3,2)-A(2,2)*A(3,1)) A_inv(3,2)=detA_1*(A(1,2)*A(3,1)-A(1,1)*A(3,2)) A_inv(3,3)=detA_1*(A(1,1)*A(2,2)-A(1,2)*A(2,1)) else print *, "ERROR: Aij with i=j=",n,"/=2or3" endif end subroutine inverse_rank_2 !================================================================================== subroutine tensor_exponential(A,expA,TOL,itr_tol) real(8),intent(in)::A(:,:),TOL real(8),allocatable,intent(inout)::expA(:,:) integer, intent(in)::itr_tol real(8),allocatable::increA(:,:) real(8) increment,NN integer i,j,n if(.not. allocated(expA) )allocate(expA(size(A,1),size(A,2) )) allocate(increA(size(A,1),size(A,2) )) if(size(A,1)/=size(A,2)) stop "ERROR:tensor exp is not a square matrix" expA(:,:)=0.0d0 do n=1,size(expA,1) expA(n,n)=1.0d0 enddo NN=1.0d0 increA(:,:)=expA(:,:) do n=1,itr_tol if(n>1)then NN = NN*(NN+1.0d0) endif increA(:,:)=matmul(increA,A) expA(:,:)= expA(:,:)+1.0d0/NN*increA(:,:) increment=0.0d0 do i=1,size(A,1) do j=1,size(A,2) increment=increment+1.0d0/NN*increA(i,j)*increA(i,j) enddo enddo if(increment<=TOL)then exit else if(n>=itr_tol)then stop "tensor exponential is not converged" endif cycle endif enddo deallocate(increA) end subroutine tensor_exponential !================================================================================== subroutine tensor_expo_der(A,expA_A,TOL,itr_tol) real(8),intent(in)::A(:,:),TOL real(8),allocatable,intent(inout)::expA_A(:,:,:,:) integer, intent(in)::itr_tol real(8),allocatable::increA_1(:,:),increA_2(:,:),increA_3(:,:,:,:),I_ij(:,:),A_inv(:,:) real(8) increment,NN integer i,j,k,l,n,m,o if(.not. allocated(expA_A) )allocate(expA_A(size(A,1),size(A,1),size(A,1),size(A,1) )) allocate(I_ij(size(A,1),size(A,1) )) allocate(increA_1(size(A,1),size(A,1) )) allocate(increA_2(size(A,1),size(A,1) )) allocate(increA_3(size(A,1),size(A,1),size(A,1),size(A,1) ) ) if(size(A,1)/=size(A,2)) stop "ERROR:tensor exp is not a square matrix" call inverse_rank_2(A,A_inv) I_ij(:,:)=0.0d0 do n=1,size(expA_A,1) I_ij(n,n)=1.0d0 enddo NN=1.0d0 do i=1,size(A,1) do j=1,size(A,1) do k=1, size(A,1) do l=1, size(A,1) expA_A(i,j,k,l)=I_ij(i,k)*I_ij(l,j) enddo enddo enddo enddo increA_1(:,:)=I_ij(:,:) increA_2(:,:)=I_ij(:,:) do n=1,itr_tol if(n>2)then NN = NN*(NN+1.0d0) endif increA_1(:,:)=A_inv(:,:) increA_2(:,:)=matmul(increA_2,A) increA_3(:,:,:,:)=0.0d0 do m=1,n increA_1(:,:)=matmul(increA_1,A ) increA_2(:,:)=matmul(increA_2,A_inv) do i=1,size(A,1) do j=1,size(A,1) do k=1, size(A,1) do l=1, size(A,1) increA_3(i,j,k,l)=increA_3(i,j,k,l)+increA_1(i,k)*increA_2(l,j) expA_A(i,j,k,l)=expA_A(i,j,k,l)+1.0d0/NN*increA_3(i,j,k,l) enddo enddo enddo enddo enddo do i=1,size(A,1) do j=1,size(A,1) do k=1, size(A,1) do l=1, size(A,1) increment=increment+1.0d0/NN*increA_3(i,j,k,l)& *increA_3(i,j,k,l)& *increA_3(i,j,k,l)& *increA_3(i,j,k,l) enddo enddo enddo enddo if(increment<=TOL)then exit else if(n>=itr_tol)then stop "tensor exponential is not converged" endif cycle endif enddo deallocate(increA_1,increA_2,increA_3,I_ij,A_inv) end subroutine tensor_expo_der !================================================================================== function GetNormRe(a) result(b) real(8),intent(in)::a(:) real(8) :: b b=dot_product(a,a) end function !================================================================================== function GetNormMatRe(a) result(b) real(8),intent(in)::a(:,:) real(8) :: b integer :: i,j b=0 do i=1,size(a,1) do j=1,size(a,2) b=b+a(i,j)*a(i,j) enddo enddo end function !================================================================================== function trace(a) result(b) real(8),intent(in)::a(:,:) real(8) :: b integer :: i,j b=0 do i=1,size(a,1) b=b+a(i,i) enddo end function !================================================================================== !================================================================================== function pi(n) result(res) integer,intent(in)::n real(8) :: ptr real(8) :: an,bn,tn,pn real(8) :: atr,btr,ttr real(8) :: res integer :: i an=1.0d0 bn=1.0d0/sqrt(2.0d0) tn=0.250d0 pn=1.00d0 do i=1,n atr=0.50d0*(an+bn) btr=dsqrt(an*bn) ttr=tn-pn*(atr-an)*(atr-an) ptr=2.0d0*pn an=atr bn=btr tn=ttr pn=ptr res=(atr+btr)*(atr+btr)/4.0d0/ttr enddo end function !================================================================================== !================================================================================== function fstring_int(x) result(a) integer,intent(in) :: x character(len=20) :: a write(a,*) x end function !================================================================================== !================================================================================== function fstring_int_len(x,length) result(a) integer,intent(in) :: x integer,intent(in) :: length character(len=length) :: a write(a,*) x end function !================================================================================== !================================================================================== function fstring_real(x) result(a) real(8),intent(in) :: x character(len=20) :: a write(a,'(f0.8)') x end function !================================================================================== !================================================================================== function fstring_real_len(x,length) result(a) real(8),intent(in) :: x integer,intent(in) :: length character(len=60) :: a character*40 :: form write(a,'(f0.10)') x end function !================================================================================== !================================================================================== function fint(ch) result(a) character(*),intent(in) :: ch integer :: a read(ch,*) a end function !================================================================================== !================================================================================== function freal(ch) result(a) character(*),intent(in) :: ch real(8) :: a read(ch,*) a end function !================================================================================== !================================================================================== function input_Int(default,option) result(val) integer,intent(in) :: default integer,optional,intent(in)::option integer :: val if(present(option) )then val=option else val=default endif end function !================================================================================== !================================================================================== function input_Real(default,option) result(val) real(8),intent(in) :: default real(8),optional,intent(in)::option real(8) :: val if(present(option) )then val=option else val=default endif end function !================================================================================== !================================================================================== function input_IntVec(default,option) result(val) integer,intent(in) :: default(:) integer,optional,intent(in)::option(:) integer,allocatable :: val(:) integer :: n,m if(present(option) )then n=size(option,1) allocate(val(n) ) val(:)=option(:) else n=size(default,1) allocate(val(n) ) val(:)=default(:) endif end function !================================================================================== !================================================================================== function input_Realvec(default,option) result(val) real(8),intent(in) :: default(:) real(8),optional,intent(in)::option(:) real(8),allocatable :: val(:) integer :: n,m if(present(option) )then n=size(option,1) allocate(val(n) ) val(:)=option(:) else n=size(default,1) allocate(val(n) ) val(:)=default(:) endif end function !================================================================================== !================================================================================== function input_IntArray(default,option) result(val) integer,intent(in) :: default(:,:) integer,optional,intent(in)::option(:,:) integer,allocatable :: val(:,:) integer :: n,m if(present(option) )then n=size(option,1) m=size(option,2) allocate(val(n,m) ) val(:,:)=option(:,:) else n=size(default,1) m=size(default,2) allocate(val(n,m) ) val(:,:)=default(:,:) endif end function !================================================================================== !================================================================================== function input_RealArray(default,option) result(val) real(8),intent(in) :: default(:,:) real(8),optional,intent(in)::option(:,:) real(8),allocatable :: val(:,:) integer :: n,m if(present(option) )then n=size(option,1) m=size(option,2) allocate(val(n,m) ) val(:,:)=option(:,:) else n=size(default,1) m=size(default,2) allocate(val(n,m) ) val(:,:)=default(:,:) endif end function !================================================================================== !================================================================================== function input_String(default,option) result(val) character(*),intent(in) :: default character(*),optional,intent(in)::option character(len(default) ) :: val if(present(option) )then val=option else val=default endif end function !================================================================================== !================================================================================== function input_logical(default,option) result(val) logical,intent(in) :: default logical,optional,intent(in)::option logical :: val if(present(option) )then val=option else val=default endif end function !================================================================================== function zeroif_Int(val,negative,positive) result(retval) integer,intent(in)::val integer :: retval logical,optional,intent(in) :: negative,positive if(val/=val)then print *, "ERROR :: MAthClass >> zeroif_Int is invalid" endif retval=val if(present(negative) )then if(negative .eqv. .true.)then if(val<0)then retval=0 endif endif endif if(present(positive) )then if(positive .eqv. .true.)then if(val>0)then retval=0 endif endif endif end function function zeroif_Real(val,negative,positive) result(retval) real(8),intent(in)::val real(8) :: retval logical,optional,intent(in) :: negative,positive if(val/=val)then print *, "ERROR :: MAthClass >> zeroif_Int is invalid" endif retval=val if(present(negative) )then if(negative .eqv. .true.)then if(val<0.0d0)then retval=0.0d0 endif endif endif if(present(positive) )then if(positive .eqv. .true.)then if(val>0.0d0)then retval=0.0d0 endif endif endif end function ! ######################################################## subroutine removeWord_String(str,keyword,itr,Compare) character(*),intent(inout)::str character(*),intent(in )::keyword integer :: len_total,len_kw,i,j,n,itr_max integer,optional,intent(in)::itr logical,optional,intent(in)::Compare logical :: bk if(present(Compare))then if(Compare .eqv. .true.)then print *, "Before :: ",str endif endif itr_max=input(default=1,option=itr) bk=.false. len_total=len(str) len_kw =len(keyword) do i=1,itr_max n=index(str,keyword) do j=n,n+len_kw str(j:j)=" " enddo if(n==0)then exit endif enddo if(present(Compare))then if(Compare .eqv. .true.)then print *, "After :: ",str endif endif end subroutine ! ######################################################## end module MathClass