subroutine pardisoinit_sol(self)

type(pardiso_format), intent(inout) :: self

!self%pt = 0
self%mtype = 13
self%solver = 1
!self%iparm(1) = 0
! self%iparm(3) = 4

! call pardisoinit(self%pt, self%mtype, self%solver, self%iparm, self%dparm, self%error)

call pardisoinit(self%pt, self%mtype, self%iparm)

!if(self%error /= 0) then
!    write(*,*) 'Error happens in PARDISOINIT with digit',self%error
!    stop
!end if


end subroutine pardisoinit_sol

! -----------------------------------------------------------------------------------------
subroutine pardiso_sol(self, this)

type(pardiso_format), intent(inout) :: self
type(chebychev_solver), intent(inout) :: this
integer :: i,j,k,sum,h
integer :: reduced_matrix_size

!self%pt = 0
self%maxfct = 1
self%mnum = 1
!self%mtype = 13
self%phase = 13
self%msglvl = 1
self%nrhs = 1
!self%iparm(1) = 0

k=1
! compress sparse row format (CRS format)

self%ia(1) = 1

do i=1,this%n_order
    sum = 0
    do j=1,this%n_order
        if (uncompressed_matrix(i,j) /= (0.0,0.0)) then
            self%a(k) = uncompressed_matrix(i,j)
            self%ja(k) = j
            self%ia(k) = i
            k=k+1
            sum = sum + 1
        end if
    end do 
    
    self%ia(i+1) = self%ia(i) + sum
    
end do

reduced_matrix_size = k-1
h = this%n_order+1



do i=1,this%n_order
    self%b(i,self%nrhs) = this%x(i)
  !  self%x(i,self%nrhs) = this%result_vector(i)
end do

self%iparm(27) = 1

do i=1,reduced_matrix_size
    write(20,*)  self%a(i)
end do

write(21,*)  self%ia(1:h)
write(22,*)  self%ja(1:reduced_matrix_size) 

write(*,*) 'ia ja a pt',h,reduced_matrix_size,reduced_matrix_size,size(self%pt)

stop


self%phase = 11
call pardiso( self%pt, self%maxfct, self%mnum, self%mtype, self%phase, this%n_order, &
    self%a(1:reduced_matrix_size), self%ia(1:h), &
    self%ja(1:reduced_matrix_size), self%perm, self%nrhs, self%iparm, &
    self%msglvl, self%b(1:this%n_order,1), self%x(1:this%n_order,1), self%error)

self%phase = 22
call pardiso( self%pt, self%maxfct, self%mnum, self%mtype, self%phase, this%n_order, &
    self%a(1:reduced_matrix_size), self%ia(1:h), &
    self%ja(1:reduced_matrix_size), self%perm, self%nrhs, self%iparm, &
    self%msglvl, self%b(1:this%n_order,1), self%x(1:this%n_order,1), self%error)

self%phase = 33
call pardiso( self%pt, self%maxfct, self%mnum, self%mtype, self%phase, this%n_order, &
    self%a(1:reduced_matrix_size), self%ia(1:h), &
    self%ja(1:reduced_matrix_size), self%perm, self%nrhs, self%iparm, &
    self%msglvl, self%b(1:this%n_order,1), self%x(1:this%n_order,1), self%error)


do i=1,this%n_order
    this%result_vector(i) = self%x(i,self%nrhs)
end do

self%phase = -1
call pardiso( self%pt, self%maxfct, self%mnum, self%mtype, self%phase, this%n_order, &
    self%a(1:reduced_matrix_size), self%ia(1:h), &
    self%ja(1:reduced_matrix_size), self%perm, self%nrhs, self%iparm, &
    self%msglvl, self%b(1:this%n_order,1), self%x(1:this%n_order,1), self%error)

!-----------Original format-------------
!call pardiso( self%pt, self%maxfct, self%mnum, self%mtype, self%phase, this%n_order, &
!    self%a(1:reduced_matrix_size), self%ia(1:h), &
!    self%ja(1:reduced_matrix_size), self%perm, self%nrhs, self%iparm, &
!    self%msglvl, self%b(1:this%n_order,1), self%x(1:this%n_order,1), self%error)



!if(self%error /= 0) then
!    write(*,*) 'Error happens in PARDISO with error digit',self%error
!   ! stop
!end if


end subroutine pardiso_sol