set_arrays_DirectSolver.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module m_set_arrays_directsolver_contact
  use hecmw_util
 
  implicit none
 
  integer (kind=kint)               :: numnon0
  integer (kind=kint), allocatable :: pointers(:)
  integer (kind=kint), allocatable :: indices(:)
  real (kind=kreal)   , allocatable :: values(:)          
 
  logical :: symmetricmatrixstruc
 
contains
 
 
  subroutine set_pointersandindices_directsolver(hecMAT,hecLagMAT,is_sym)
 
    type(hecmwst_matrix)                     :: hecMAT
    type (hecmwST_matrix_lagrange)           :: hecLagMAT
    logical :: is_sym
 
    integer (kind=kint)     :: np
    integer (kind=kint)     :: ndof
    integer (kind=kint)     :: num_lagrange
    integer (kind=kint)     :: nn
    integer (kind=kint)     :: ierr
    integer (kind=kint)     :: i, j, k, l, countNon0
 
    symmetricmatrixstruc = is_sym
 
    np = hecmat%NP ; ndof = hecmat%NDOF ; num_lagrange = heclagmat%num_lagrange
    nn = np*ndof + num_lagrange + 1
 
    if( symmetricmatrixstruc )then
      numnon0 = hecmat%NPU*ndof**2+hecmat%NP*ndof*(ndof+1)/2 &
        + (heclagmat%numU_lagrange)*ndof+heclagmat%num_lagrange
    else
      numnon0 = (hecmat%NPL+hecmat%NPU+hecmat%NP)*ndof**2 &
        + (heclagmat%numL_lagrange+heclagmat%numU_lagrange)*ndof
    endif
 
    if(allocated(pointers))deallocate(pointers)
    allocate(pointers(nn), stat=ierr)
    if( ierr /= 0 ) stop " Allocation error, mkl%pointers "
    pointers = 0
 
    if(allocated(indices))deallocate(indices)
    allocate(indices(numnon0), stat=ierr)
    if( ierr /= 0 ) stop " Allocation error, mkl%indices "
    indices = 0
 
    pointers(1) = 1
    countnon0 = 1
 
    do i = 1, np
      do j = 1, ndof
        if( .not. symmetricmatrixstruc )then
          do l = hecmat%indexL(i-1)+1, hecmat%indexL(i)
            do k = 1, ndof
              indices(countnon0) = (hecmat%itemL(l)-1)*ndof + k
              countnon0 = countnon0 + 1
            enddo
          enddo
          do k = 1, j-1
            indices(countnon0) = (i-1)*ndof + k
            countnon0 = countnon0 + 1
          enddo
        endif
        do k = j, ndof
          indices(countnon0) = (i-1)*ndof + k
          countnon0 = countnon0 + 1
        enddo
        do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
          do k = 1, ndof
            indices(countnon0) = (hecmat%itemU(l)-1)*ndof + k
            countnon0 = countnon0 + 1
          enddo
        enddo
        if( num_lagrange > 0 )then
          do l = heclagmat%indexU_lagrange(i-1)+1, heclagmat%indexU_lagrange(i)
            indices(countnon0) = np*ndof + heclagmat%itemU_lagrange(l)
            countnon0 = countnon0 + 1
          enddo
        endif
        pointers((i-1)*ndof+j+1) = countnon0
      enddo
    enddo
 
    if( num_lagrange > 0 )then
      do i = 1, num_lagrange
        if( symmetricmatrixstruc )then
          indices(countnon0) = np*ndof + i
          countnon0 = countnon0 + 1
        else
          do l = heclagmat%indexL_lagrange(i-1)+1, heclagmat%indexL_lagrange(i)
            do k = 1, ndof
              indices(countnon0) = (heclagmat%itemL_lagrange(l)-1)*ndof + k
              countnon0 = countnon0 + 1
            enddo
          enddo
        endif
        pointers(np*ndof+i+1) = countnon0
      enddo
    endif
 
  end subroutine set_pointersandindices_directsolver
 
 
  subroutine set_values_directsolver(hecMAT,hecLagMAT)
 
    type(hecmwst_matrix)                    :: hecMAT
    type (hecmwST_matrix_lagrange)          :: hecLagMAT
 
    integer (kind=kint)     :: np
    integer (kind=kint)     :: ndof
    integer (kind=kint)     :: num_lagrange
    integer (kind=kint)     :: ierr
    integer (kind=kint)     :: i, j, k, l
    integer (kind=kint)     :: countNon0, locINal, locINd, locINau, locINal_lag, locINau_lag
 
    np = hecmat%NP ; ndof = hecmat%NDOF ; num_lagrange = heclagmat%num_lagrange
 
    if(allocated(values))deallocate(values)
    allocate(values(numnon0), stat=ierr)
    if( ierr /= 0 ) stop " Allocation error, mkl%values "
    values = 0.0d0
 
    countnon0 = 1
    do i = 1, np
      do j = 1, ndof
        if( .not. symmetricmatrixstruc )then
          do l = hecmat%indexL(i-1)+1, hecmat%indexL(i)
            do k = 1, ndof
              locinal = ((l-1)*ndof+j-1)*ndof + k
              values(countnon0) = hecmat%AL(locinal)
              countnon0 = countnon0 + 1
            enddo
          enddo
          do k = 1, j-1
            locind = ((i-1)*ndof+j-1)*ndof + k
            values(countnon0) = hecmat%D(locind)
            countnon0 = countnon0 + 1
          enddo
        endif
        do k = j, ndof
          locind = ((i-1)*ndof+j-1)*ndof + k
          values(countnon0) = hecmat%D(locind)
          countnon0 = countnon0 + 1
        enddo
        do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
          do k = 1, ndof
            locinau = ((l-1)*ndof+j-1)*ndof + k
            values(countnon0) = hecmat%AU(locinau)
            countnon0 = countnon0 + 1
          enddo
        enddo
        if( num_lagrange > 0 )then
          do l = heclagmat%indexU_lagrange(i-1)+1, heclagmat%indexU_lagrange(i)
            locinau_lag = (l-1)*ndof + j
            values(countnon0) = heclagmat%AU_lagrange(locinau_lag)
            countnon0 = countnon0 + 1
          enddo
        endif
      enddo
    enddo
 
    if( .not.symmetricmatrixstruc .and. num_lagrange > 0 )then
      do i = 1, num_lagrange
        do l = heclagmat%indexL_lagrange(i-1)+1, heclagmat%indexL_lagrange(i)
          do k = 1, ndof
            locinal_lag = (l-1)*ndof + k
            values(countnon0) = heclagmat%AL_lagrange(locinal_lag)
            countnon0 = countnon0 + 1
          enddo
        enddo
      enddo
    endif
 
  end subroutine set_values_directsolver
 
  subroutine getapproximateb(ntdf,x,y)
 
    integer(kind=kint)     :: ntdf
    integer(kind=kint)     :: i, j, k
    real(kind=kreal)        :: x(ntdf)                    
    real(kind=kreal)        :: y(ntdf)                    
 
    y = 0.0d0
    do i = 1, ntdf
      do j = pointers(i), pointers(i+1)-1
        k = indices(j)
        y(i) = y(i) + values(j)*x(k)
        if( symmetricmatrixstruc .and. k/=i )&
          y(k)=y(k)+values(j)*x(i)
      enddo
    enddo
 
  end subroutine getapproximateb
 
 
  subroutine checkresidual(hecMAT,hecLagMAT)
    type(hecmwst_matrix)                    :: hecmat
    type (hecmwst_matrix_lagrange)          :: hecLagMAT
    integer(kind=kint)                       :: ntdf
    real(kind=kreal), allocatable            :: y(:)           
    real(kind=kreal)                          :: residual_max   
 
    ntdf = hecmat%NP*hecmat%NDOF + heclagmat%num_lagrange
 
    allocate(y(size(hecmat%B)))
    y = 0.0d0
    call getapproximateb(ntdf,hecmat%X,y)
    residual_max=maxval(dabs(y-hecmat%B))
    write(*,*)' maximum residual = ',residual_max
    if( dabs(residual_max) >= 1.0d-8) then
      write(*,*) ' ###Maximum residual exceeded 1.0d-8---Direct Solver### '
      !        stop
    endif
    deallocate(y)
 
  end subroutine checkresidual
 
 
end module m_set_arrays_directsolver_contact