D:/demona/2.0/manual/doc/2-31set/processor/mf.f90

Go to the documentation of this file.

module mf


use global

use model


contains



function matVec(v)

implicit none

real(8) v(:)

real(8) matVec(size(v))

real(8) myEps


matVec=0.d0

if (dot_product(v,v).le.epsilon(0.d0)) return


! have to decide on eps!!!
myEps=dsqrt(epsilon(0.d0))





!myEps=epsilon(0.0)

!tempy(setVec)=y(setVec)  !<---------------???????????

tempy=y  !<---------------???????????

!print*,myEps


!print*,tempy

!stop


myEps=epsValue(tempy(setVec),v)

!print*,myeps

!stop


!tempy(myModel%physics(currentPhysics)%domain(currentDomain)%setVec)=tempy(myModel%physics(currentPhysics)%domain(currentDomain)%setVec)+myeps*v
tempy(setVec)=tempy(setVec)+myeps*v



matVec= ( discrete(1) - discrete(0) ) /myEps


tempy(setVec)=y(setVec)  !<---------------???????????
tempy=y  !<---------------???????????





end function matVec



function epsValue(x,v) result(eps)
implicit none

integer i,n

real(8) x(:),v(:)

real(8),allocatable :: typVec(:)

real(8) eps,myeps1,myeps2

real(8) normX,normV

real(8) typSize

real(8) dotXv

real(8) sum

real(8) sigma

integer myCase

real(8) myValue


myCase=2
!myCase=3

n=size(v)



myeps1=dsqrt(epsilon(0.d0))

myeps2=epsilon(0.0)

sigma=1e-14

!eps=1e-10

!print*,dsqrt(epsilon(0.d0))
!print*,epsilon(0.0)
!stop

normX=dsqrt(dot_product(X,X))
!normV=dot_product(v,v)
normV=dsqrt(dot_product(v,v))





dotXv=dot_product(X,v)

  
        select case(myCase)
 
 
        case (1) ! keyes
  
                !print*,1e-6

                myValue=0.000001d0

                !print*,myValue

                !stop

                allocate(typVec(n))

                call typicalSize(X,n,typSize)

                typVec=typSize

                !eps=1e-6*max(dabs(dotXv),dot_product(typVec,dabs(v)))/normV                    !A

                eps=myValue*max(dabs(dotXv),dot_product(typVec,dabs(v)))/normV                  !B

                ! A is better than B in poissons equation

                                                                                                        
                ! 1e-6 because of "Globalized Newton-Krylov-Schwarz algorithms and software for parallel implicit CFD"

                eps=dabs(eps)

                !print*,eps

                !if (isnan(eps)) then
                !       print*,normV
                !       !print*,'Nan!!'
                !       pause
                !       print*,normV,dotXv,normX,dot_product(typVec,dabs(v)),eps
                !       print*,v
                !       stop
                !end if
    

        case(2) ! saad & brown

                allocate(typVec(n))
                
                call typicalSize(X,n,typSize)

                typVec=typSize

                !eps=myeps1*max(dabs(dotXv),dot_product(typVec,dabs(v)))* (sign(dotXv,1.d0)/dotXv)/normV
                                                                                                                                        ! this part is sign function - returns 1 or -1
                eps=myeps1*max(dabs(dotXv),dot_product(typVec,dabs(v)))/normV
                !eps=myeps2*max(dabs(dotXv),dot_product(typVec,dabs(v)))/normV
                                                                                                        

                eps=dabs(eps)

                if (isnan(eps)) then

                        print*,'isNan'
                        stop
                        print*,normV,dotXv,normX,dot_product(typVec,dabs(v)),eps
                        print*,v
                        stop

                end if

                !if (eps.lt.0.d0) then

                !!print*,dotXv
                !!print*,(sign(dotXv,1.d0)/dotXv) 
                !!print*,dot_product(typVec,v)
                !!print*,max(dabs(dotXv),dot_product(typVec,dabs(v)))
                !!print*,typSize
                !!print*,eps
                !!stop

                !eps=dabs(eps)

                !end if

                !neg eps?

        case(3)

                eps=myeps1


        case(4)

                eps=myeps2

        case(5) !quin

                eps=dsqrt(sigma)/normV

        case(6) !kelley

                eps=dsqrt(sigma)*normX/normV


        case(7) !McHugh & Knoll

                sum=0.d0

                do i=1,n
                        
                        sum=X(i)+sum

                end do

                sum=sum/(n*1.d0)

                sum=sum+1.d0

                eps=dsqrt(sigma)*sum


        case(8) ! pernice & walker
                

        end select


!write(5555,'(E18.10)')eps


end function epsValue



subroutine typicalSize(X,n,typSize)
implicit none
integer n
real(8) X(n)
real(8) minTyp
real(8) maxTyp
real(8) typSize

typSize=0.d0


        if (minval(dabs(X)).eq.0.d0) then
                minTyp=0.d0
        else
                minTyp=dlog10(minval(dabs(X)))
        end if

        if (maxval(dabs(X)).eq.0.d0) then
                maxTyp=0.d0
        else
                maxTyp=dlog10(maxval(dabs(X)))
        end if

        !print*,'mintyp-----------------'
        !print*,mintyp
        !print*,'maxtyp-----------------'
        !print*,maxtyp

        typSize=10.d0**((minTyp+maxTyp)*0.5d0)
        typSize=dabs(typSize)

end subroutine typicalSize

!******************************************************************

subroutine colorJacobianCRS(val,col_ind,row_ptr,dof,symbolic,valSize,indSize,diagVecInd,returnNumberOfColors,numberOfColors)

implicit none

integer dof

real(8):: myF(dof)

real(8) val(:)

integer col_ind(:),row_ptr(:)

integer symbolic,valSize,indSize,ptrSize

integer  colorVector(dof,3)
real(8):: eps
real(8), allocatable ::  Fcolors(:,:)
integer:: i, j, k, z

integer unk

integer stencil
integer numberOfColors
integer,allocatable :: colors(:)

integer diagonalColor
integer aa



integer,pointer :: direction_stencil(:)
integer,allocatable :: newcolors(:)
!integer,allocatable :: direction_stencil(:),newcolors(:)



integer minvalue,minindex,mymod,column



integer crsColSum,TotalcrsColSum

real(8),allocatable :: crsColVal(:)
integer,allocatable :: crsColInd(:),sortedCrsColInd(:)
integer startIndex,endIndex
integer nonZeroIndex
integer, allocatable :: nonZeroVec(:)

integer diagVecInd(dof),diagVecIndDesignator,diagVecIndVal




!###################


! demona related defitions

integer, pointer :: myDofmatrix(:,:,:,:)

integer diagI,diagJ

integer dofOfmyDofmatrix

real(8) myEps


!integer startIndex,endIndex

integer,pointer :: nx,ny

integer myP,myD

integer,pointer :: localSetVec(:)

integer myZ

integer colorCorrector

integer returnNumberOfColors

integer zDof

integer,pointer :: minI,minJ,maxI,maxJ

eps=dsqrt(epsilon(0.d0))
!eps=epsilon(0.0)




stencil=5



diagVecIndDesignator=0
diagVecIndVal=0



!stop




colorVector=0

pVec=>myModel%pVec

endIndex=0
colorFunctionIndex=0
colorValMax=0

do myP=1,size(pVec)

        currentPhysics=pVec(myP)

        !print*,'---------',currentPhysics

        dVec=>myModel%physics(currentPhysics)%dVec

        do myD=1,size(dVec)

                currentDomain=dVec(myD)

                !print*,currentDomain

                myModel%physics(currentPhysics)%domain(currentDomain)%colorMin=colorValMax+1

                nx=>myModel%physics(currentPhysics)%domain(currentDomain)%nx
                ny=>myModel%physics(currentPhysics)%domain(currentDomain)%ny

                minI=>myModel%physics(currentPhysics)%domain(currentDomain)%minI
                maxI=>myModel%physics(currentPhysics)%domain(currentDomain)%maxI
                minJ=>myModel%physics(currentPhysics)%domain(currentDomain)%minJ
                maxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%maxJ

                localSetVec=>myModel%physics(currentPhysics)%domain(currentDomain)%setVec

                !unk=myModel%physics(currentPhysics)%dof
                unk=myModel%physics(currentPhysics)%domain(currentDomain)%setDof

                startIndex=endIndex+1

                endIndex=endIndex+size(localSetVec)


                !print*,startIndex,endIndex

                !print*,unk


                myModel%physics(currentPhysics)%domain(currentDomain)%colorStartIndex=startIndex
                myModel%physics(currentPhysics)%domain(currentDomain)%colorEndIndex=endIndex

                
                !print*,size(colorVector(setVec(startIndex:endIndex),:)),size(ColoringFunction(nx,ny,stencil,size(localSetVec),unk))
                !print*,size(ColoringFunction(nx,ny,stencil,size(localSetVec),unk))

                !print*,setVec(startIndex:endIndex)

                !colorVector(setVec(startIndex:endIndex),:)=ColoringFunction(nx,ny,stencil,size(localSetVec),unk)
                !colorVector(startIndex:endIndex,:)=ColoringFunction(nx,ny,stencil,size(localSetVec),unk)

                colorVector(startIndex:endIndex,:)=ColoringFunction(maxI-minI+1,maxJ-minJ+1,stencil,size(localSetVec),unk)


                tempy=y

                colorFunctionIndex=endIndex

                colorValMax=maxVal(colorVector(:,1))

                myModel%physics(currentPhysics)%domain(currentDomain)%colorMax=colorValMax
                
        end do

end do



numberOfColors=colorValMax









!print*,numberOfColors

if (returnNumberOfColors.eq.1) return




!************************
! crs relared definition
allocate(crsColVal(numberOfColors),crsColInd(numberOfColors))
crsColVal=0.d0
crsColInd=0
TotalcrsColSum=0
!************************





!allocate(Fcolors(dof,numberOfColors),colors(numberOfColors))
allocate(Fcolors(dof,numberOfColors),colors(stencil))
Fcolors=0.d0



!do i=1,numberOfColors
do i=1,stencil
        colors(i)=i
end do



!< evaluate Fcolors
Fcolors=ColoringEvaluation(colorVector,numberOfColors,dof);


!< calculate F(x)
tempy=y
myF=discrete(1)

tempy=y

if (test.eq.1) then

!print*,Fcolors(2,:)-myF(2)
!stop

end if



allocate(direction_stencil(stencil),newcolors(stencil))


direction_stencil=0
newcolors=0


do myP=1,size(pVec)

        currentPhysics=pVec(myP)

        !print*,'##############'
        !print*,'currentPhysics=',currentPhysics
        !print*,'##############'

        dVec=>myModel%physics(currentPhysics)%dVec

        do myD=1,size(dVec)

                currentDomain=dVec(myD)

                !print*,'**************'
                !print*,'currentDomain=',currentDomain
                !print*,'**************'

                
                myDofmatrix=>myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix

                startIndex=myModel%physics(currentPhysics)%domain(currentDomain)%colorStartIndex

                endIndex=myModel%physics(currentPhysics)%domain(currentDomain)%colorEndIndex

                colorMin=>myModel%physics(currentPhysics)%domain(currentDomain)%colorMin

                colorMax=>myModel%physics(currentPhysics)%domain(currentDomain)%colorMax


                !print*,startIndex,endIndex,colorMin,colorMax




                do  i=startIndex,endIndex
                
                        !if(test.eq.1) print*,'---------------------------',i


                        diagonalColor=colorVector(i,1);


    

                        diagI=colorVector(i,2)
                        diagJ=colorVector(i,3)

                        


                        !print*,diagonalColor,diagi,diagj
                        !stop

    

                        ! instead of defining direction_stecil beforehand
                        ! it should be defined in the z loop

                        
                        !!************************
                        !! crs relared definition
                        crsColVal=0.d0
                        crsColInd=0
                        crsColSum=1
                        !!************************


                                  
                        newcolors(1:stencil)=cshift(colors(1:stencil),(diagonalColor-1));

                        colorCorrector=colorMin-1
               
        
                        !do z=1,numberOfColors
                        do z=colorMin,colorMax


    
                                !if Fcolors(i,z)==0   % is zero or Fcolors(i,z)-myF(i) == 0   ?
                                !if ((Fcolors(i,z)-myF(i)) .eq. 0) then  ! is zero or Fcolors(i,z)-myF(i) == 0
                                if (dabs(Fcolors(i,z)-myF(i)) .le. epsilon(0.d0)) then  ! is zero or Fcolors(i,z)-myF(i) == 0
                                        cycle
                                end if


                                zDof=(colorMax-colorMin+1)/stencil

                                !if (test.eq.1) print*,'     z=',z


                                !row=i;
        
                                !stop
             
       
                                ! reorder color permutation for the diagonal color
                                ! here it does not matter if diagonalColor > stencil then the
                                ! circular shift occurs more than once like saying pi=3*pi=5*pi

                                ! puts the color of the diagonal and continues with colors in this order: CENSW
                                ! example, if the color of diag is 4 then 4-5-1-2-3
        
                                !do k=1,1
    
                                        !newcolors(1+(k-1)*stencil:stencil*k)=cshift(colors(1+(k-1)*stencil:stencil*k),-(diagonalColor-1));
                                !       newcolors(1+(k-1)*stencil:stencil*k)=cshift(colors(1+(k-1)*stencil:stencil*k),(diagonalColor-1));
                 
                                !end do

                                !do k=1,numberOfColors/stencil

                                        
                                !       newcolors(newcolors(1+(k-1)*stencil:stencil*k))=

                                !end do
                        
                                myZ=mod(z,stencil)

                                if(myZ.eq.0) myZ=stencil


        
                                !---------------
                                !newcolors=newcolors-z;
                                !newcolors=abs(newcolors);
        
                                !if (i.eq.5) print*,'     z=',z

        
        
                                ! j detection
        
                                ![minvalue minindex]=min(abs(newcolors-z));
        
                                !minvalue=minval(iabs(newcolors-z))
                                minvalue=minval(iabs(newcolors-myZ))


                                !gives the index of the current color z, on the newcolors vector
                                !minindex=minloc(iabs(newcolors-z),dim=1)
                                minindex=minloc(iabs(newcolors-myZ),dim=1)




                                minindex=minIndex-colorCorrector
                                
                                minindex=minindex+stencil*(z-myZ)/stencil

      
                                mymod=mod(minindex,stencil)

                                ! mymod=1 is C
                                ! mymod=2 is E
                                ! mymod=3 is N
                                ! mymod=4 is S
                                ! mymod=0 is W



                                dofOfmyDofmatrix=((minindex-mymod)/stencil+min(mymod,1))

                                !print*,minindex,dofOfmyDofmatrix,giveDirections(mymod),mymod


                                select case (mymod)

                                        case(1)

                                                column=myDofmatrix(diagI,diagJ,dofOfmyDofmatrix,relative)

                                        case(2)

                                                column=myDofmatrix(diagI+1,diagJ,dofOfmyDofmatrix,relative)

                                        case(3)

                                                column=myDofmatrix(diagI,diagJ+1,dofOfmyDofmatrix,relative)

                                        case(4)

                                                column=myDofmatrix(diagI,diagJ-1,dofOfmyDofmatrix,relative)

                                        case(0)

                                                column=myDofmatrix(diagI-1,diagJ,dofOfmyDofmatrix,relative)


                                end select



                if (column.eq.0) then           
                !print*,minindex,minValue,dofOfmyDofmatrix,giveDirections(mymod),z
                !print*,i,diagi,diagj,mymod


                !print*,'-----------------------------------------'
                
                !stop



                end if




                                crsColVal(crsColSum)=(Fcolors(i,z)-myF(i))/eps

                                crsColInd(crsColSum)=column



                                if (i.eq.column) then

                                        diagVecIndDesignator=crsColSum !+diagVecIndDesignator

                                        !print*,i,column,diagVecIndDesignator,crsColSum

                                        !pause

                                end if

                                !if (i.eq.5) print*,'column=',column


                                crsColSum=crsColSum+1

                                !if (test.eq.1) print*,'column=',column

       
                        end do


                        nonZeroIndex=0

                        ! since zero jacobian entries are skipped with the cycle command
                        ! at the beginning of the main z  loop, zero values  are located
                        ! at the end of on the crsColInd vector
                        !do z=1,numberOfColors
                        !       if (crsColInd(z).eq.0) then
                        !               exit
                        !       else
                        !               nonZeroIndex=nonZeroIndex+1
                        !       end if
                        !end do

                        !print*,nonZeroIndex,crsColSum-1,nonZeroIndex-(crsColSum-1)
                        !pause
                        nonZeroIndex=crsColSum-1
                        

                        allocate(nonZeroVec(nonZeroIndex),sortedCrsColInd(nonZeroIndex))
                        nonZeroVec=0
                        sortedCrsColInd=0

                        nonZeroVec(1:nonZeroIndex)=crsColInd(1:nonZeroIndex)


                        ! is sortedCrsColInd same as the newColors vector evaluated for z=1 or z of the center
                        !call qsortd(crsColInd,sortedCrsColInd,numberOfColors)
                        call qsortd(nonZeroVec,sortedCrsColInd,nonZeroIndex)


                        startIndex=TotalcrsColSum+1

                        TotalcrsColSum=crsColSum-1+TotalcrsColSum


                        ! this two lines are useful to detect the number of entries in the i'th row of the jacobian
                        !print*,i,crsColSum-1
                        !pause


                        endIndex=TotalcrsColSum

                        
                        diagVecIndDesignator=minloc(abs(sortedCrsColInd-diagVecIndDesignator),dim=1)


                        
                        !diagVecIndVal=sortedCrsColInd(diagVecIndDesignator)+diagVecIndVal
                        !diagVecIndVal=diagVecIndDesignator+diagVecIndVal
                        diagVecIndVal=diagVecIndDesignator+startIndex-1



                        !if (symbolic.eq.1) diagVecInd(i)=diagVecIndVal  ! if not symbolic - do not calculate
                        diagVecInd(i)=diagVecIndVal  ! if not symbolic - do not calculate


                
                




                        do z=1,nonZeroIndex

                                if ((startIndex+z-1).gt.size(val)) then

                                        print*,'woink'
                                        print*,startIndex+z-1,i,z,size(val)
                                        stop

                                end if


                        
                                val(startIndex+z-1)=crsColVal(sortedCrsColInd(z))
                                col_ind(startIndex+z-1)=crsColInd(sortedCrsColInd(z))  !maybe sortedCrsColInd should be used only.



                        end do

                        !pause

                        row_ptr(i)=startIndex



                        deallocate(nonZeroVec,sortedCrsColInd)



                end do
        !stop   
        end do ! dVec loop ends here


end do !pVec loop ends here





row_ptr(dof+1)=TotalcrsColSum+1



if (symbolic.eq.1) then

        valSize=TotalcrsColSum
        indSize=TotalcrsColSum

end if

!print*,valsize
!pause




end subroutine colorJacobianCRS

!##########################################################################################




!##########################################################################################


function ColoringFunction(nx,ny,stencil,dof,unk) 
implicit none

integer nx,ny,unk,stencil,dof

integer ColoringFunction(dof,3)

integer AA(Ny,Nx),Hcolors(stencil),Vcolors(stencil)

integer i,j,k,t,Vindex,Hindex,row,myColor,index

integer iStart,iEnd
integer jStart,jEnd
!integer kStart,kEnd 3D coloring is not yet available

integer,pointer :: myDofMatrix(:,:,:,:)

integer myGhostCells

integer,pointer :: dofVec(:)

integer myT

ColoringFunction=0


!A=zeros(Ny,Nx);
AA=0






!Hcolors=1:stencil;
do i=1,stencil
        Hcolors(i)=i
end do


select case(stencil)
    
    case(5)
        
        !%        K
        !%        |
        !%   G----R----Y
        !%        |
        !%        B
        !% K:blacK:1
        !% B:Blue:2
        !% G:Green:3
        !% R:Red:4 
        !% Y:Yellow:5
        
        !Vcolors=[1 3 5 2 4];
                Vcolors(1)=1
                Vcolors(2)=3
                Vcolors(3)=5
                Vcolors(4)=2
                Vcolors(5)=4

        do i=1,Ny
            
            if (mod(i,stencil).eq.0) then
                Vindex=stencil
            else
                Vindex=mod(i,stencil)
            end if

            row=Ny-i+1
            
            myColor=Vcolors(Vindex)
            
            AA(row,1)=myColor

                        AA(row,1)=AA(row,1)+colorValMax
            
            do j=2,Nx

                Hindex=myColor-1+j;
                
                if (mod(Hindex,stencil).eq.0) then
                    Hindex=stencil
                else
                    Hindex=mod(Hindex,stencil)
                end if

                AA(row,j)=Hcolors(Hindex)+colorValMax

            end do
            
        end do

   case default
      print*,'!-----------------------------------------------'
      print*,'no coloring for this stencil is defined!!!!!'
          stop
end select










! change the color matrix into a colorbased and dof based vector.


!Here the colors are defined for the full matrix
!below, the matrix turned into a vector. however, only values
! with index > 0 should be written.




!print*,currentPhysics,currentDomain
!print*,myModel%physics(currentPhysics)%pdof
!print*,myModel%physics(currentPhysics)%domain(currentDomain)%ddof
!stop

!print*,myModel%physics(currentPhysics)%domain(currentDomain)%dofmatrix(:,1,1,relative,1),dof

!print*,discType,unk

!stop


if (dFX.ne.dFV)  then
! finite difference setup

        istart=1
        jStart=1
        
        iEnd=Nx
        jend=Ny

        myGhostCells=0


else
!finite volume ! do not forget to consider the number of ghost cells options when available

!currently assuming that two ghost cells are used

        istart=1-ghostCells !-1
        jstart=1-ghostCells !-1
        
        iend=Nx-ghostCells !Nx-2
        jend=Ny-ghostCells !Ny-2


        myGhostCells=ghostCells




end if


myDofmatrix=>myModel%physics(currentPhysics)%domain(currentDomain)%dofmatrix

dofVec=>myModel%physics(currentPhysics)%domain(currentDomain)%dofVec

setDof=>myModel%physics(currentPhysics)%domain(currentDomain)%setDof

do j=jStart,jEnd

        do i=iStart,iEnd

                do myT=1,setDof

                        t=dofVec(myT)

                        if (myDofmatrix(i,j,t,relative).ne.0 ) then

                                ColoringFunction(myDofmatrix(i,j,t,relative)-colorFunctionIndex,1)=AA(Ny-myGhostCells+1-j,i+myGhostCells)+(t-1)*stencil
                                !ColoringFunction(myDofmatrix(i,j,t,relative)-colorFunctionIndex,1)=AA(Ny-myGhostCells+1-j,i+myGhostCells)+(myT-1)*stencil
                                ColoringFunction(myDofmatrix(i,j,t,relative)-colorFunctionIndex,2)=i
                                ColoringFunction(myDofmatrix(i,j,t,relative)-colorFunctionIndex,3)=j

                        end if

                end do

        end do

end do






end function ColoringFunction

!##########################################################################################

function ColoringEvaluation(colorVector,numberOfColors,dof);
implicit none

integer i,j,index

integer nx,ny,unk,stencil,dof,numberOfColors

integer colorVector(dof)

real(8) eps

real(8) ColoringEvaluation(dof,numberOfColors)

real(8) myEps

real(8) typSize,myX(1)

ColoringEvaluation=0.d0



eps=epsilon(0.0)
eps=dsqrt(epsilon(0.d0))


!% first fill the FC with original x
!% here we avoid storing colored x's 
!% in new memory locations.

do i=1,numberOfColors
    ColoringEvaluation(:,i)=tempy(setVec)
end do




!% this part should be index independent - i.e. no mesh related data should be given.
!% maybe it is a better idea to have the grid in vector rather than a matrix



!eps=dsqrt(epsilon(0.d0));


!do i=1,numberOfColors
!    ColoringEvaluation(:,i)=x
!end


!myX(1)=-1.d0

!call typicalSize(myX,1,typSize)

!print*,typsize


!stop


!% x's are modified with epsilon

!stop

do i=1,dof
    
   index=colorVector(i);

!   myEps=epsvalue(,) !<#########################################################<<<<<<<
   
   ColoringEvaluation(i,index)= ColoringEvaluation(i,index)+eps;

   !print*,index,ColoringEvaluation(i,index)
   !stop
       
end do

!print*,ColoringEvaluation(1,:)
!stop
!% x's are evaluated


do i=1,numberOfColors

        tempy(setVec)=ColoringEvaluation(:,i)

        ColoringEvaluation(:,i)=discrete(1)

        tempy=y


end do


!PRINT*,ColoringEvaluation(1,:)

!STOP






end function ColoringEvaluation

!##########################################################################################


function crsMatVecMult(val,col_ind,row_ptr,vec)

implicit none

real(8) val(:),vec(:)

integer col_ind(:),row_ptr(:)

real(8) crsMatVecMult(size(Vec))

integer i,j

crsMatVecMult=0.d0

do i=1,size(Vec)

        !print*,i
        
        do j=row_ptr(i),row_ptr(i+1)-1

                !print*,row_ptr(i),row_ptr(i+1)-1

                crsMatVecMult(i)=crsMatVecMult(i)+val(j)*vec(col_ind(j))

        end do

end do



end function crsMatVecMult

!##########################################################################################

function LUyX(val,col_ind,row_ptr,dof,diag_ptr,valSize,Vec)

implicit none

integer i,j

integer dof,valsize

real(8) Vec(dof),pivots(dof),z(dof)

real(8) val(valSize)

real(8) LUyX(dof)

integer col_ind(valSize),diag_ptr(dof),row_ptr(dof+1)

real(8) sum

LUyX=0.d0

z=0.d0

do i=1,dof
        
        pivots(i)=val(diag_ptr(i))
        pivots(i)=1.d0/pivots(i)

end do




! L solve
do i=1,dof
        
        sum=0.d0
        
        do j=row_ptr(i),diag_ptr(i)-1


                sum=sum+val(j)*z(col_ind(j))

        end do

        !z(i)=pivots(i)*(vec(i)-sum)
        z(i)=vec(i)-sum


end do

!print*,z


!stop

! U solve
do i=dof,1,-1
        
        sum=0.d0
        
        do j=diag_ptr(i)+1,row_ptr(i+1)-1

                sum=sum+val(j)*LUyX(col_ind(j))

        end do

        !LUyX(i)=z(i) - pivots(i)*sum
        LUyX(i)=(z(i) - sum)*pivots(i)


end do




end function LUyX



function giveDirections(value)
implicit none
character giveDirections
integer value

select case(value)

        case(1)

                givedirections='C'

        case(2)

                givedirections='E'
        
        case(3)

                givedirections='N'

        case(4)

                givedirections='S'

        case(0)

                givedirections='W'

end select

end function



!##########################################################################################
!##########################################################################################







!##########################################################################################

SUBROUTINE qsortd(x,ind,n)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2002-12-18  Time: 11:55:47

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER, INTENT(IN)    :: n

!REAL (dp), INTENT(IN)  :: x(n)


integer, INTENT(IN)  :: x(n)
INTEGER, INTENT(OUT)   :: ind(n)


!***************************************************************************

!                                                         ROBERT RENKA
!                                                 OAK RIDGE NATL. LAB.

!   THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO SORT A REAL (dp)
! ARRAY X INTO INCREASING ORDER.  THE ALGORITHM IS AS FOLLOWS.  IND IS
! INITIALIZED TO THE ORDERED SEQUENCE OF INDICES 1,...,N, AND ALL INTERCHANGES
! ARE APPLIED TO IND.  X IS DIVIDED INTO TWO PORTIONS BY PICKING A CENTRAL
! ELEMENT T.  THE FIRST AND LAST ELEMENTS ARE COMPARED WITH T, AND
! INTERCHANGES ARE APPLIED AS NECESSARY SO THAT THE THREE VALUES ARE IN
! ASCENDING ORDER.  INTERCHANGES ARE THEN APPLIED SO THAT ALL ELEMENTS
! GREATER THAN T ARE IN THE UPPER PORTION OF THE ARRAY AND ALL ELEMENTS
! LESS THAN T ARE IN THE LOWER PORTION.  THE UPPER AND LOWER INDICES OF ONE
! OF THE PORTIONS ARE SAVED IN LOCAL ARRAYS, AND THE PROCESS IS REPEATED
! ITERATIVELY ON THE OTHER PORTION.  WHEN A PORTION IS COMPLETELY SORTED,
! THE PROCESS BEGINS AGAIN BY RETRIEVING THE INDICES BOUNDING ANOTHER
! UNSORTED PORTION.

! INPUT PARAMETERS -   N - LENGTH OF THE ARRAY X.

!                      X - VECTOR OF LENGTH N TO BE SORTED.

!                    IND - VECTOR OF LENGTH >= N.

! N AND X ARE NOT ALTERED BY THIS ROUTINE.

! OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N PERMUTED IN THE SAME
!                          FASHION AS X WOULD BE.  THUS, THE ORDERING ON
!                          X IS DEFINED BY Y(I) = X(IND(I)).

!*********************************************************************

! NOTE -- IU AND IL MUST BE DIMENSIONED >= LOG(N) WHERE LOG HAS BASE 2.

!*********************************************************************

INTEGER   :: iu(21), il(21)
INTEGER   :: m, i, j, k, l, ij, it, itt, indx
REAL      :: r
REAL (dp) :: t

! LOCAL PARAMETERS -

! IU,IL =  TEMPORARY STORAGE FOR THE UPPER AND LOWER
!            INDICES OF PORTIONS OF THE ARRAY X
! M =      INDEX FOR IU AND IL
! I,J =    LOWER AND UPPER INDICES OF A PORTION OF X
! K,L =    INDICES IN THE RANGE I,...,J
! IJ =     RANDOMLY CHOSEN INDEX BETWEEN I AND J
! IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND
! INDX =   TEMPORARY INDEX FOR X
! R =      PSEUDO RANDOM NUMBER FOR GENERATING IJ
! T =      CENTRAL ELEMENT OF X

IF (n <= 0) RETURN

! INITIALIZE IND, M, I, J, AND R



DO  i = 1, n

  ind(i) = i
END DO


m = 1
i = 1
j = n
r = .375

! TOP OF LOOP

20 IF (i >= j) GO TO 70
IF (r <= .5898437) THEN
  r = r + .0390625
ELSE
  r = r - .21875
END IF

! INITIALIZE K

30 k = i

! SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T

ij = i + r*(j-i)
it = ind(ij)
t = x(it)

! IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T,
!   INTERCHANGE IT WITH T

indx = ind(i)
IF (x(indx) > t) THEN
  ind(ij) = indx
  ind(i) = it
  it = indx
  t = x(it)
END IF

! INITIALIZE L

l = j

! IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T,
!   INTERCHANGE IT WITH T

indx = ind(j)
IF (x(indx) >= t) GO TO 50
ind(ij) = indx
ind(j) = it
it = indx
t = x(it)

! IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T,
!   INTERCHANGE IT WITH T

indx = ind(i)
IF (x(indx) <= t) GO TO 50
ind(ij) = indx
ind(i) = it
it = indx
t = x(it)
GO TO 50

! INTERCHANGE ELEMENTS K AND L

40 itt = ind(l)
ind(l) = ind(k)
ind(k) = itt

! FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS
!   NOT LARGER THAN T

50 l = l - 1
indx = ind(l)
IF (x(indx) > t) GO TO 50

! FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS NOT SMALLER THAN T

60 k = k + 1
indx = ind(k)
IF (x(indx) < t) GO TO 60

! IF K <= L, INTERCHANGE ELEMENTS K AND L

IF (k <= l) GO TO 40

! SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE
!   ARRAY YET TO BE SORTED

IF (l-i > j-k) THEN
  il(m) = i
  iu(m) = l
  i = k
  m = m + 1
  GO TO 80
END IF

il(m) = k
iu(m) = j
j = l
m = m + 1
GO TO 80

! BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY

70 m = m - 1
IF (m == 0) RETURN
i = il(m)
j = iu(m)

80 IF (j-i >= 11) GO TO 30
IF (i == 1) GO TO 20
i = i - 1

! SORT ELEMENTS I+1,...,J.  NOTE THAT 1 <= I < J AND J-I < 11.

90 i = i + 1
IF (i == j) GO TO 70
indx = ind(i+1)
t = x(indx)
it = indx
indx = ind(i)
IF (x(indx) <= t) GO TO 90
k = i

100 ind(k+1) = ind(k)
k = k - 1
indx = ind(k)
IF (t < x(indx)) GO TO 100

ind(k+1) = it
GO TO 90
END SUBROUTINE qsortd



!##########################################################################################











end module mf

---