hecmw_matrix_dump.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 hecmw_matrix_dump
  use hecmw_util
  use hecmw_matrix_misc
  use m_hecmw_comm_f
 
  private
 
  public :: hecmw_mat_dump_type_none
  public :: hecmw_mat_dump_type_mm
  public :: hecmw_mat_dump_type_csr
  public :: hecmw_mat_dump_type_bsr
 
  public :: hecmw_mat_dump
  public :: hecmw_mat_dump_rhs
  public :: hecmw_mat_dump_solution
 
  integer(kind=kint), parameter :: hecmw_mat_dump_type_none = 0
  integer(kind=kint), parameter :: hecmw_mat_dump_type_mm   = 1
  integer(kind=kint), parameter :: hecmw_mat_dump_type_csr  = 2
  integer(kind=kint), parameter :: hecmw_mat_dump_type_bsr  = 3
 
  integer, save :: numcall = 0
 
contains
 
  subroutine hecmw_mat_dump( hecMAT, hecMESH )
    implicit none
    type(hecmwst_matrix) :: hecmat
    type(hecmwst_local_mesh) :: hecmesh
    numcall = numcall + 1
    select case( hecmw_mat_get_dump(hecmat) )
      case (hecmw_mat_dump_type_none)
        return
      case (hecmw_mat_dump_type_mm)
        call hecmw_mat_dump_mm(hecmat)
      case (hecmw_mat_dump_type_csr)
        call hecmw_mat_dump_csr(hecmat)
      case (hecmw_mat_dump_type_bsr)
        call hecmw_mat_dump_bsr(hecmat)
    end select
    call hecmw_mat_dump_rhs(hecmat)
    if (hecmw_mat_get_dump_exit(hecmat) /= 0) then
      call hecmw_barrier( hecmesh )
      stop "Exiting program after dumping matrix"
    end if
  end subroutine hecmw_mat_dump
 
  subroutine make_file_name(ext, fname)
    implicit none
    character(*) :: ext
    character(*) :: fname
    write(fname,"('dump_matrix_',I0,'_',I0,A)") &
      numcall, hecmw_comm_get_rank(), ext
  end subroutine make_file_name
 
  subroutine hecmw_mat_dump_mm( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecmat
    integer, parameter :: idump = 201
    character(len=64) :: fname
    integer :: i, j, i0, j0, idof, jdof, ii, jj
    integer :: idxl0, idxl, idxd, idxu0, idxu
    integer :: n, np, ndof, ndof2, nnz
    character(len=64), parameter :: lineformat = "(I0,' ',I0,' ',e20.12e3)"
    integer :: stat
    !n = hecMAT%N
    n = hecmat%NP
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof * ndof
    ! make fname
    call make_file_name('.mm', fname)
    ! open file
    open(idump, file=fname, status='replace', iostat=stat)
    if (stat /= 0) then
      write(*,*) 'WARNING: cannot open file ', fname, ' for matrix dump'
      return
    end if
    ! header
    write(idump,"(A)") '%%MatrixMarket matrix coordinate real general'
    nnz = ndof2 * (n + hecmat%indexL(n) + hecmat%indexU(n))
    write(idump,"(I0,' ',I0,' ',I0)") n*ndof, np*ndof, nnz
    idxd = 0
    do i = 1, n
      i0 = (i-1)*ndof
      do idof = 1, ndof
        ii = i0 + idof
        ! Lower
        do j = hecmat%indexL(i-1)+1,hecmat%indexL(i)
          j0 = (hecmat%itemL(j)-1)*ndof
          idxl0 = (j-1)*ndof2 + (idof-1)*ndof
          do jdof = 1, ndof
            jj = j0 + jdof
            idxl = idxl0 + jdof
            write(idump,lineformat) ii, jj, hecmat%AL(idxl)
          end do
        end do
        ! Diagonal
        j0 = i0
        do jdof = 1, ndof
          jj = j0 + jdof
          idxd = idxd + 1
          write(idump,lineformat) ii, jj, hecmat%D(idxd)
        end do
        ! Upper
        do j = hecmat%indexU(i-1)+1,hecmat%indexU(i)
          j0 = (hecmat%itemU(j)-1)*ndof
          idxu0 = (j-1)*ndof2 + (idof-1)*ndof
          do jdof = 1, ndof
            jj = j0 + jdof
            idxu = idxu0 + jdof
            write(idump,lineformat) ii, jj, hecmat%AU(idxu)
          end do
        end do
      end do
    end do
    ! close file
    close(idump)
  end subroutine hecmw_mat_dump_mm
 
  subroutine hecmw_mat_dump_csr( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecmat
    integer, parameter :: idump = 201
    character(len=64) :: fname
    integer :: i, j, i0, j0, idof, jdof, ii, jj
    integer :: idx, idxd, idxl, idxu, idxl0, idxu0
    integer :: n, np, ndof, ndof2, nnz, nnz1
    character(len=64), parameter :: lineformat = "(e20.12e3)"
    integer :: stat
    !n = hecMAT%N
    n = hecmat%NP
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof * ndof
    ! make fname
    call make_file_name('.csr', fname)
    ! open file
    open(idump, file=fname, status='replace', iostat=stat)
    if (stat /= 0) then
      write(*,*) 'WARNING: cannot open file ', fname, ' for matrix dump'
      return
    end if
    ! header
    write(idump,"(A)") '%%CSR matrix real general'
    nnz = ndof2 * (n + hecmat%indexL(n) + hecmat%indexU(n))
    write(idump,"(A)") '% nrow ncol nnonzero'
    write(idump,"(I0,' ',I0,' ',I0)") n*ndof, np*ndof, nnz
    ! index
    write(idump,"(A)") '% index(0:nrow)'
    idx = 0
    write(idump, "(I0)") idx
    do i = 1, n
      nnz1 = ndof * ((hecmat%indexL(i)-hecmat%indexL(i-1)) + &
        1 + (hecmat%indexU(i)-hecmat%indexU(i-1)))
      do idof = 1, ndof
        idx = idx + nnz1
        write(idump, "(I0)") idx
      end do
    end do
    ! item
    write(idump,"(A)") '% item(1:nnonzero)'
    do i = 1, n
      i0 = (i-1)*ndof
      do idof = 1, ndof
        ! Lower
        do j = hecmat%indexL(i-1)+1,hecmat%indexL(i)
          j0 = (hecmat%itemL(j)-1)*ndof
          do jdof = 1, ndof
            jj = j0 + jdof
            write(idump,"(I0)") jj
          end do
        end do
        ! Diagonal
        j0 = i0
        do jdof = 1, ndof
          jj = j0 + jdof
          write(idump,"(I0)") jj
        end do
        ! Upper
        do j = hecmat%indexU(i-1)+1,hecmat%indexU(i)
          j0 = (hecmat%itemU(j)-1)*ndof
          do jdof = 1, ndof
            jj = j0 + jdof
            write(idump,"(I0)") jj
          end do
        end do
      end do
    end do
    ! values
    write(idump,"(A)") '% value(1:nnonzero)'
    idxd = 0
    do i = 1, n
      i0 = (i-1)*ndof
      do idof = 1, ndof
        ii = i0 + idof
        ! Lower
        do j = hecmat%indexL(i-1)+1,hecmat%indexL(i)
          j0 = (hecmat%itemL(j)-1)*ndof
          idxl0 = (j-1)*ndof2 + (idof-1)*ndof
          do jdof = 1, ndof
            jj = j0 + jdof
            idxl = idxl0+jdof
            write(idump,lineformat) hecmat%AL(idxl)
          end do
        end do
        ! Diagonal
        j0 = i0
        do jdof = 1, ndof
          jj = j0 + jdof
          idxd = idxd + 1
          write(idump,lineformat) hecmat%D(idxd)
        end do
        ! Upper
        do j = hecmat%indexU(i-1)+1,hecmat%indexU(i)
          j0 = (hecmat%itemU(j)-1)*ndof
          idxu0 = (j-1)*ndof2 + (idof-1)*ndof
          do jdof = 1, ndof
            jj = j0 + jdof
            idxu = idxu0 + jdof
            write(idump,lineformat) hecmat%AU(idxu)
          end do
        end do
      end do
    end do
    ! close file
    close(idump)
  end subroutine hecmw_mat_dump_csr
 
  subroutine hecmw_mat_dump_bsr( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecmat
    integer, parameter :: idump = 201
    character(len=64) :: fname
    integer :: i, j
    integer :: idx, idxl0, idxd0, idxu0
    integer :: n, np, ndof, ndof2, nnz, nnz1
    character(len=64), parameter :: lineformat = "(e20.12e3)"
    integer :: stat
    !n = hecMAT%N
    n = hecmat%NP
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof * ndof
    ! make fname
    call make_file_name('.bsr', fname)
    ! open file
    open(idump, file=fname, status='replace', iostat=stat)
    if (stat /= 0) then
      write(*,*) 'WARNING: cannot open file ', fname, ' for matrix dump'
      return
    end if
    ! header
    write(idump,"(A)") '%%Block-CSR matrix real general'
    nnz = n + hecmat%indexL(n) + hecmat%indexU(n)
    write(idump,"(A)") '% nrow ncol nnonzero ndof'
    write(idump,"(I0,' ',I0,' ',I0,' ',I0)") n, np, nnz, ndof
    ! index
    write(idump,"(A)") '% index(0:nrow)'
    idx = 0
    write(idump, "(I0)") idx
    do i = 1, n
      nnz1 = (hecmat%indexL(i)-hecmat%indexL(i-1)) + &
        1 + (hecmat%indexU(i)-hecmat%indexU(i-1))
      idx = idx + nnz1
      write(idump, "(I0)") idx
    end do
    ! item
    write(idump,"(A)") '% item(1:nnonzero)'
    do i = 1, n
      ! Lower
      do j = hecmat%indexL(i-1)+1,hecmat%indexL(i)
        write(idump,"(I0)") hecmat%itemL(j)
      end do
      ! Diagonal
      write(idump,"(I0)") i
      ! Upper
      do j = hecmat%indexU(i-1)+1,hecmat%indexU(i)
        write(idump,"(I0)") hecmat%itemU(j)
      end do
    end do
    ! values
    write(idump,"(A)") '% value(1:nnonzero*ndof*ndof)'
    idxd0 = 0
    do i = 1, n
      ! Lower
      do j = hecmat%indexL(i-1)+1,hecmat%indexL(i)
        idxl0 = (j-1)*ndof2
        write(idump,lineformat) hecmat%AL(idxl0+1:idxl0+ndof2)
      end do
      ! Diagonal
      write(idump,lineformat) hecmat%D(idxd0+1:idxd0+ndof2)
      idxd0 = idxd0 + ndof2
      ! Upper
      do j = hecmat%indexU(i-1)+1,hecmat%indexU(i)
        idxu0 = (j-1)*ndof2
        write(idump,lineformat) hecmat%AU(idxu0+1:idxu0+ndof2)
      end do
    end do
    ! close file
    close(idump)
  end subroutine hecmw_mat_dump_bsr
 
  subroutine hecmw_mat_dump_rhs( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecmat
    integer, parameter :: idump = 201
    character(len=64) :: fname
    integer :: i
    integer :: n, np, ndof, ndof2
    character(len=64), parameter :: lineformat = "(e20.12e3)"
    integer :: stat
    if( hecmw_mat_get_dump(hecmat) == hecmw_mat_dump_type_none) return
    !n = hecMAT%N
    n = hecmat%NP
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof * ndof
    ! make fname
    call make_file_name('.rhs', fname)
    ! open file
    open(idump, file=fname, status='replace', iostat=stat)
    if (stat /= 0) then
      write(*,*) 'WARNING: cannot open file ', fname, ' for matrix dump'
      return
    end if
    do i = 1, np*ndof
      write(idump,lineformat) hecmat%B(i)
    end do
    ! close file
    close(idump)
  end subroutine hecmw_mat_dump_rhs
 
  subroutine hecmw_mat_dump_solution( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecmat
    integer, parameter :: idump = 201
    character(len=64) :: fname
    integer :: i
    integer :: n, np, ndof, ndof2
    character(len=64), parameter :: lineformat = "(e20.12e3)"
    integer :: stat
    if( hecmw_mat_get_dump(hecmat) == hecmw_mat_dump_type_none) return
    !n = hecMAT%N
    n = hecmat%NP
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof * ndof
    ! make fname
    call make_file_name('.sol', fname)
    ! open file
    open(idump, file=fname, status='replace', iostat=stat)
    if (stat /= 0) then
      write(*,*) 'WARNING: cannot open file ', fname, ' for matrix dump'
      return
    end if
    do i = 1, np*ndof
      write(idump,lineformat) hecmat%X(i)
    end do
    ! close file
    close(idump)
  end subroutine hecmw_mat_dump_solution
 
end module hecmw_matrix_dump