Switch to pardiso solver

Author: HarryMclean

src/heatflow/mod_evolve.f90

@@ -31,7 +31,7 @@ module evolution
   use boundary_vector, only: boundary
   use cattaneo, only: S_catS
 !   use tempdep, only: ChangeProp 
-  use sparse_solver, only: bicgstab
+  use sparse_solver, only: bicgstab, solve_pardiso
   implicit none
 
   private
@@ -46,8 +46,8 @@ module evolution
 !!!#################################################################################################
   subroutine simulate(itime)
     integer(int12), intent(in) :: itime
-    real(real12), dimension(NA) :: S, x, Q, Qdens, S_CAT, B
-    real(real12), dimension(:), allocatable :: x0
+    real(real12), dimension(NA) :: S, Q, Qdens, S_CAT, B
+    real(real12), dimension(:), allocatable :: x
     integer:: ncg, itol, itmax !, iss
     integer :: iter
     real(real12) :: e, err, tol
@@ -140,21 +140,23 @@ subroutine simulate(itime)
    ! iter:  Output - gives the number of the final iteration.
    ! err:   Output - records the error of the final iteration.
    ! iss:   Input - sets the Sparse Storage type (1=SRS, 2=SDS).
-    x=Temp_p+(Temp_p-Temp_pp) 
-    if (any(x-Temp_p .lt. TINY)) x=x+TINY !avoid nan solver issue
+   !  x=Temp_p+(Temp_p-Temp_pp)
+   !  if (any(x-Temp_p .lt. TINY)) x=x+TINY !avoid nan solver issue
     itol=1
     tol=1.e-32_real12
     itmax=50000
     ncg = 0
-    iter=ncg
+    iter= 0
     err=E
 
 
-   CALL bicgstab(acsr, ia, ja, S, itmax, x, x0, iter)
+   !  CALL bicgstab(acsr, ia, ja, S, itmax, x, x0, iter)
+
+   CALL solve_pardiso(acsr, S, ia, ja, x)
    !  CALL linbcg(S,x,itol=int(itol,I4B),tol=tol, itmax=int(itmax,I4B), iter=iter, &
-   !       err=E)
+         ! err=E)
 
-   x=x0
+   ! 
     if (any(isnan(x(:)))) then
        write(0,*) "fatal error: NAN in x tempurature vector"
        write(0,*) 'time step ', itime, "      T   ", sum(Temp_p)/size(Temp_p), E ,iter

src/heatflow/mod_global.f90

@@ -14,8 +14,8 @@ module globe_data
   TYPE(sprs2_dp) :: ra !Techniqually rH
   TYPE(diag_sprs_dp) :: da !Techniqually dH
   real(real12), dimension(:), allocatable :: acsr
-  integer (kind=8), dimension(:), allocatable :: ja
-  integer (kind=8), dimension(:), allocatable :: ia
+  integer, dimension(:), allocatable :: ja
+  integer, dimension(:), allocatable :: ia
 
 end module globe_data
 

src/heatflow/mod_setup.f90

@@ -81,7 +81,8 @@ subroutine set_global_variables()
          CALL sparse_Hmatrix()
          ! Allocate the arrays to hold the H matrix in CSR format
          allocate(acsr(ra%len), ja(ra%len), ia(ra%n+1))
-         CALL coo2csr(ra%n, ra%len, ra%val, ra%irow, ra%jcol, acsr, ja, ia)               
+         CALL coo2csr(ra%n, ra%len, ra%val, ra%irow, ra%jcol, acsr, ja, ia)
+         ! print*, ra%val               
       end if
       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

src/heatflow/mod_sparse_solver.f90

@@ -56,8 +56,8 @@ subroutine coo2csr( nrow, &
     integer (kind=8), dimension(nnz), intent(in) :: ir
     integer (kind=8), dimension(nnz), intent(in) :: jc   
     real(8), dimension(nnz), intent(out) :: acsr
-    integer (kind=8), dimension(nnz), intent(out) :: ja
-    integer (kind=8), dimension(nrow+1), intent(out) :: ia
+    integer, dimension(nnz), intent(out) :: ja
+    integer, dimension(nrow+1), intent(out) :: ia
 
     ! Local variables.
     integer (kind=8) :: i, iad, j, k, k0
@@ -151,16 +151,17 @@ subroutine bicgstab( acsr, &
     call mkl_dcsrgemv("N",n,acsr,ia,ja,x,temp1)
 
     r = b - temp1
-  
+
     call random_number(rst)
 
     p = r
+
     delta = dot_product(rst,r)
-    
+
     write(*,'(a,1x,f15.3)') "Starting delta: ", delta
 
     delta0 = delta
-    
+
     do i = 1, maxiter
 
        if ( norm2(r) /= norm2(r) ) then
@@ -187,7 +188,7 @@ subroutine bicgstab( acsr, &
        delta = dot_product(rst,r)
        beta = (delta/delta_old)*(alpha/omega)
        p = r + beta*(p - omega*temp1)
-       
+
        if(norm2(r) .lt. cc) then
           iter = i
           return
@@ -247,7 +248,9 @@ subroutine solve_pardiso( acsr, &
     mnum = 1
 
     if (.not.(allocated(x))) allocate(x(n))
-    
+
+
+
     allocate(iparm(64))		!set up pardiso control parameter
 
     do i=1,64