!*****************
! File Description:  
!   
! Creator: Hoang-Trong Minh Tuan 
! 
! Rev: 0.1 
! Update: 
!    x.x: 
!*******************/
MODULE matrix_utility
  USE precisioninfo
  IMPLICIT NONE

  INTERFACE print_mat
     MODULE PROCEDURE print_mat_int, print_matdbl, print_matdbl_detail, print_matint_detail 
  END INTERFACE
  
  SAVE

CONTAINS

  !====================
  ! print out the matrix of size rows x cols
  !===================
  SUBROUTINE print_mat_int(matrix, rows, cols, filename)
    INTEGER, DIMENSION(rows,cols) :: matrix
    INTEGER, INTENT(IN) :: rows, cols
    CHARACTER(*), OPTIONAL :: filename
    !local
    INTEGER :: I, J
    LOGICAL:: UNITOK, UNITOP
    INTEGER :: fileID
    IF (PRESENT(filename)) THEN
       fileID = 20
       UNITOK = .FALSE.; unitop = .TRUE.
       DO WHILE  (.NOT. (UNITOK .AND. .NOT. UNITOP)) 
          INQUIRE (unit=fileID,exist=UNITOK,opened=UNITOP)
       END DO
       OPEN(fileID, FILE=TRIM(filename))
    ELSE
       fileID = 6 ! standard output
    ENDIF
       
    DO I = 1, rows 
       DO J = 1, cols
          WRITE(fileID, "(I5)", ADVANCE='NO') matrix(I,J)
       END DO
       WRITE(fileID,*)
    END DO
    IF (PRESENT(filename))  CLOSE(fileID)
       
  END SUBROUTINE print_mat_int
  
  SUBROUTINE print_mat_T(matrix, rows, cols)
    INTEGER, DIMENSION(rows,cols) :: matrix
    INTEGER, INTENT(IN) :: rows, cols
    !local
    INTEGER :: I, J
    DO J = 1, cols
       DO I = 1, rows 
          WRITE(*, "(I5)", ADVANCE='NO') matrix(I,J)
       END DO
       WRITE(*,*)
    END DO

  END SUBROUTINE print_mat_T

  SUBROUTINE print_matdbl(matrix, rows, cols)
    REAL(KIND=REAL8), DIMENSION(rows,cols) :: matrix
    INTEGER, INTENT(IN) :: rows, cols
    !local
    INTEGER :: I, J
    DO I = 1, rows 
       DO J = 1, cols
          WRITE(*, "(f15.7,',')", ADVANCE='NO') matrix(I,J)
       END DO
       WRITE(*,*)
    END DO

  END SUBROUTINE print_matdbl

  SUBROUTINE print_matdbl_T(matrix, rows, cols)
    REAL(KIND=REAL8), DIMENSION(rows,cols) :: matrix
    INTEGER, INTENT(IN) :: rows, cols
    !local
    INTEGER :: I, J
    DO J = 1, cols
       DO I = 1, rows 
          WRITE(*, "(f15.7,',')", ADVANCE='NO') matrix(I,J)
       END DO
       WRITE(*,*)
    END DO

  END SUBROUTINE print_matdbl_T
  
  !=========================
  ! set all matrix elements to zeros 
  !=========================
  SUBROUTINE mat_reset(matrix, rows, cols)
    INTEGER, INTENT(IN) :: rows, cols
    INTEGER, DIMENSION(:,:) :: matrix

    INTEGER :: I, J

    DO I = 1, rows
       DO J = 1, cols
          matrix(I,J) = 0
       END DO
    END DO
  END SUBROUTINE mat_reset


  ! Map from (i,j,k) to linear index
  ! MAXX = # of columns
  ! MAXY = # of rows
  FUNCTION map3D2linear(MAXX, MAXY, i, j, k) RESULT(res)
    INTEGER, INTENT(IN) :: MAXX, MAXY, i, j, k
    INTEGER :: res

    res = (k-1) * MAXX * MAXY + (j-1) * MAXX + i
  END FUNCTION map3D2linear

  SUBROUTINE print_matdbl_detail(matrix, rstart, nrows, cstart, ncols)
    REAL(KIND=REAL8), DIMENSION(rstart:rstart+nrows-1,cstart:cstart+ncols-1) :: matrix
    INTEGER, INTENT(IN) :: rstart, nrows, cstart, ncols
    !local
    INTEGER :: I, J
    DO I = rstart, rstart+nrows-1
       DO J = cstart, cstart+ncols-1
          WRITE(*, "(f15.7,',')", ADVANCE='NO') matrix(I,J)
       END DO
       WRITE(*,*)
    END DO

  END SUBROUTINE print_matdbl_detail

  SUBROUTINE print_matint_detail(matrix, rstart, nrows, cstart, ncols)
    INTEGER, DIMENSION(rstart:rstart+nrows-1,cstart:cstart+ncols-1) :: matrix
    INTEGER, INTENT(IN) :: rstart, nrows, cstart, ncols
    !local
    INTEGER :: I, J
    DO I = rstart, rstart+nrows-1
       DO J = cstart, cstart+ncols-1
          WRITE(*, "(I4,',')", ADVANCE='NO') matrix(I,J)
       END DO
       WRITE(*,*)
    END DO

  END SUBROUTINE print_matint_detail
  
END MODULE matrix_utility
