!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> A module to manage the process grid.
MODULE ProcessGridModule
  USE ErrorModule, ONLY : Error_t, ConstructError, SetGenericError
  USE LoggingModule, ONLY : ActivateLogger, EnterSubLog, ExitSubLog, &
       & WriteHeader, WriteListElement
  USE NTMPIModule

  IMPLICIT NONE
  PRIVATE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A datatype which stores a process grid and all its communicators.
  TYPE, PUBLIC :: ProcessGrid_t
     !! Describe the grid
     INTEGER, PUBLIC :: total_processors !< total processors in the grid.
     INTEGER, PUBLIC :: num_process_rows !< number of rows in the grid.
     INTEGER, PUBLIC :: num_process_columns !< number of columns in the grid.
     INTEGER, PUBLIC :: num_process_slices !< number of 2D slices in the grid.
     INTEGER, PUBLIC :: slice_size !< the size of a 2D slice.
     !! Identiy current process
     INTEGER, PUBLIC :: my_slice !< which slice is the current process in.
     INTEGER, PUBLIC :: my_row !< which row is the current process in.
     INTEGER, PUBLIC :: my_column !< which column is the current process in.
     !! Ranks for communication
     INTEGER, PUBLIC :: global_rank !< current process rank amongst processes.
     !> rank for within slice communication.
     INTEGER, PUBLIC :: within_slice_rank
     !> rank for between slice communication.
     INTEGER, PUBLIC :: between_slice_rank
     INTEGER, PUBLIC :: column_rank !< rank for within column communication.
     INTEGER, PUBLIC :: row_rank !< rank for within row communication.
     !! Communicators for communication
     INTEGER, PUBLIC :: global_comm !< communicator with every other process.
     INTEGER, PUBLIC :: row_comm !< communicator within a row.
     INTEGER, PUBLIC :: column_comm !< communicator within a column.
     INTEGER, PUBLIC :: within_slice_comm !< communicator within a slice.
     INTEGER, PUBLIC :: between_slice_comm !< communicator between slices.
     INTEGER, PUBLIC :: grid_error !< stores errors from MPI calls.
     INTEGER, PUBLIC :: RootID = 0 !< Which rank is root?
     !! Blocked communication
     INTEGER, PUBLIC :: block_multiplier !< Block scaling factor.
     INTEGER, PUBLIC :: number_of_blocks_columns !< number of column blocks.
     INTEGER, PUBLIC :: number_of_blocks_rows !< number of row blocks.
     !> blocked communicator within a row.
     INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: blocked_row_comm
     !> blocked communicator within a column.
     INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: blocked_column_comm
     !> blocked communicator within a slice.
     INTEGER, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: blocked_within_slice_comm
     !> blocked communicator between slices.
     INTEGER, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: blocked_between_slice_comm
  END TYPE ProcessGrid_t
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> The default process grid.
  TYPE(ProcessGrid_t), TARGET, PUBLIC :: global_grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  PUBLIC :: ConstructProcessGrid
  PUBLIC :: ConstructNewProcessGrid
  PUBLIC :: IsRoot
  PUBLIC :: SplitProcessGrid
  PUBLIC :: CopyProcessGrid
  PUBLIC :: DestructProcessGrid
  !! Accessors for grid information
  PUBLIC :: GetMySlice
  PUBLIC :: GetMyRow
  PUBLIC :: GetMyColumn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  INTERFACE ConstructProcessGrid
     MODULE PROCEDURE ConstructProcessGrid_full
     MODULE PROCEDURE ConstructProcessGrid_onlyslice
  END INTERFACE ConstructProcessGrid
  INTERFACE ConstructNewProcessGrid
     MODULE PROCEDURE ConstructNewProcessGrid_full
     MODULE PROCEDURE ConstructNewProcessGrid_onlyslice
  END INTERFACE ConstructNewProcessGrid
CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Setup the default process grid.
  SUBROUTINE ConstructProcessGrid_full(world_comm_, process_rows_, &
       & process_columns_, process_slices_, be_verbose_in)
    !> A communicator that every process in the grid is a part of.
    INTEGER, INTENT(IN) :: world_comm_
    !> The number of grid rows.
    INTEGER, INTENT(IN) :: process_rows_
    !> The number of grid columns.
    INTEGER, INTENT(IN) :: process_columns_
    !> The number of grid slices.
    INTEGER, INTENT(IN) :: process_slices_
    !> Set true to print process grid info.
    LOGICAL, INTENT(IN), OPTIONAL :: be_verbose_in
    !! Local Data
    LOGICAL :: be_verbose

    !! Process Optional Parameters
    IF (PRESENT(be_verbose_in)) THEN
       be_verbose = be_verbose_in
    ELSE
       be_verbose = .FALSE.
    END IF

    CALL ConstructNewProcessGrid(global_grid, world_comm_, process_rows_, &
         & process_columns_, process_slices_)

    !! Report
    IF (IsRoot(global_grid)) THEN
       CALL ActivateLogger
    END IF
    IF (be_verbose) THEN
       CALL WriteHeader("Process Grid")
       CALL EnterSubLog
       CALL WriteListElement("Process Rows", &
            & VALUE=global_grid%num_process_rows)
       CALL WriteListElement(key = "Process Columns", &
            & VALUE=global_grid%num_process_columns)
       CALL WriteListElement(key = "Process Slices", &
            & VALUE=global_grid%num_process_slices)
       CALL WriteListElement(key = "Column Blocks", &
            & VALUE=global_grid%number_of_blocks_columns)
       CALL WriteListElement(key = "Row Blocks", &
            & VALUE=global_grid%number_of_blocks_rows)
       CALL ExitSubLog
    END IF

  END SUBROUTINE ConstructProcessGrid_full
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Setup a process grid specifying only the slices
  SUBROUTINE ConstructProcessGrid_onlyslice(world_comm_, process_slices_in, &
       & be_verbose_in)
    !> A communicator that every process in the grid is a part of.
    INTEGER, INTENT(IN) :: world_comm_
    !> The number of grid slices.
    INTEGER, INTENT(IN), OPTIONAL :: process_slices_in
    !> Set true to print process grid info.
    LOGICAL, INTENT(IN), OPTIONAL :: be_verbose_in
    !! Local Data
    LOGICAL :: be_verbose
    INTEGER :: process_rows, process_columns, process_slices
    INTEGER :: total_processors
    INTEGER :: ierr

    !! Total processors
    CALL MPI_COMM_SIZE(world_comm_, total_processors, ierr)

    !! Process Optional Parameters
    IF (PRESENT(be_verbose_in)) THEN
       be_verbose = be_verbose_in
    ELSE
       be_verbose = .FALSE.
    END IF
    IF (PRESENT(process_slices_in)) THEN
       process_slices = process_slices_in
    ELSE
       CALL ComputeNumSlices(total_processors, process_slices)
    END IF

    !! Create a 3D grid
    CALL ComputeGridSize(total_processors, process_slices, process_rows, &
         & process_columns)

    !! Now call the full setup
    CALL ConstructProcessGrid(world_comm_, process_rows, process_columns, &
         & process_slices, be_verbose)
  END SUBROUTINE ConstructProcessGrid_onlyslice
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Construct a process grid.
  SUBROUTINE ConstructNewProcessGrid_full(grid, world_comm_, process_rows_, &
       & process_columns_, process_slices_)
    !> The grid to construct
    TYPE(ProcessGrid_t), INTENT(INOUT) :: grid
    !> A communicator that every process in the grid is a part of.
    INTEGER, INTENT(IN) :: world_comm_
    !> The number of grid rows.
    INTEGER, INTENT(IN) :: process_rows_
    !> The number of grid columns.
    INTEGER, INTENT(IN) :: process_columns_
    !> The number of grid slices.
    INTEGER, INTENT(IN) :: process_slices_
    !! Local Data
    INTEGER :: column_block_multiplier
    INTEGER :: row_block_multiplier
    INTEGER :: II, JJ

    INTEGER :: ierr
    TYPE(Error_t) :: err

    CALL ConstructError(err)
    CALL MPI_COMM_DUP(world_comm_, grid%global_comm, ierr)
    !! Grid Dimensions
    grid%num_process_rows = process_rows_
    grid%num_process_columns = process_columns_
    grid%num_process_slices = process_slices_
    CALL MPI_COMM_SIZE(grid%global_comm, grid%total_processors, ierr)
    grid%slice_size = grid%total_processors/grid%num_process_slices

    !! Do a sanity check
    IF (grid%num_process_rows*grid%num_process_columns*grid%num_process_slices &
         & .NE. grid%total_processors) THEN
       CALL SetGenericError(err, &
            & "you did not specify a consistent process grid size", .TRUE.)
    END IF
    IF (grid%num_process_slices .GT. 1) THEN
       IF (MOD(MAX(grid%num_process_rows, grid%num_process_columns), &
            & MIN(grid%num_process_rows, grid%num_process_columns)) &
            & .NE. 0) THEN
          CALL SetGenericError(err, &
               & "if slices >1, either rows or columns must be a multiple"//&
               & "of the other.", &
               & .TRUE.)
       END IF
    END IF

    !! Grid ID
    CALL MPI_COMM_RANK(grid%global_comm,grid%global_rank,ierr)
    grid%my_slice = grid%global_rank/grid%slice_size
    grid%my_row = MOD(grid%global_rank, grid%slice_size)/ &
         & grid%num_process_columns
    grid%my_column = MOD(grid%global_rank, grid%num_process_columns)

    !! Grid Communicators
    CALL MPI_COMM_SPLIT(grid%global_comm, grid%my_slice, grid%global_rank, &
         grid%within_slice_comm, ierr)
    CALL MPI_COMM_RANK(grid%within_slice_comm, grid%within_slice_rank, ierr)
    CALL MPI_COMM_SPLIT(grid%global_comm, grid%within_slice_rank, &
         & grid%global_rank, grid%between_slice_comm, ierr)
    CALL MPI_COMM_RANK(grid%between_slice_comm, grid%between_slice_rank, ierr)
    CALL MPI_COMM_SPLIT(grid%within_slice_comm, grid%my_row, grid%global_rank, &
         & grid%row_comm, ierr)
    CALL MPI_COMM_RANK(grid%row_comm, grid%row_rank, ierr)
    CALL MPI_COMM_SPLIT(grid%within_slice_comm, grid%my_column, &
         & grid%global_rank, grid%column_comm, ierr)
    CALL MPI_COMM_RANK(grid%column_comm, grid%column_rank, ierr)

    !! Blocking Information
    column_block_multiplier = (grid%num_process_rows/grid%num_process_columns)*&
         & grid%num_process_slices
    IF (column_block_multiplier .EQ. 0) THEN
       column_block_multiplier = 1*grid%num_process_slices
    END IF
    row_block_multiplier = (grid%num_process_columns/grid%num_process_rows)* &
         & grid%num_process_slices
    IF (row_block_multiplier .EQ. 0) THEN
       row_block_multiplier = 1*grid%num_process_slices
    END IF

    !! The rule right now seems to be to have at least half as many blocks as
    !! threads.
    grid%block_multiplier = 1

    grid%number_of_blocks_columns = &
         & column_block_multiplier*grid%block_multiplier
    grid%number_of_blocks_rows = &
         & row_block_multiplier*grid%block_multiplier

    !! Create Blocked Communicators
    ALLOCATE(grid%blocked_row_comm(grid%number_of_blocks_rows))
    ALLOCATE(grid%blocked_column_comm(grid%number_of_blocks_columns))
    ALLOCATE(grid%blocked_within_slice_comm(grid%number_of_blocks_rows, &
         & grid%number_of_blocks_columns))
    ALLOCATE(grid%blocked_between_slice_comm(grid%number_of_blocks_rows, &
         & grid%number_of_blocks_columns))

    DO JJ=1,grid%number_of_blocks_columns
       DO II=1,grid%number_of_blocks_rows
          CALL MPI_COMM_SPLIT(grid%global_comm, grid%my_slice, &
               & grid%global_rank, grid%blocked_within_slice_comm(II,JJ), &
               & ierr)
          CALL MPI_COMM_SPLIT(grid%global_comm, grid%within_slice_rank, &
               & grid%global_rank, grid%blocked_between_slice_comm(II,JJ), &
               & ierr)
       END DO
    END DO
    DO JJ=1,grid%number_of_blocks_columns
       CALL MPI_COMM_SPLIT(grid%within_slice_comm, grid%my_column, &
            & grid%global_rank, grid%blocked_column_comm(JJ), ierr)
    END DO
    DO II=1,grid%number_of_blocks_rows
       CALL MPI_COMM_SPLIT(grid%within_slice_comm, grid%my_row, &
            & grid%global_rank, grid%blocked_row_comm(II), ierr)
    END DO

  END SUBROUTINE ConstructNewProcessGrid_full
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Setup a process grid specifying only the slices
  SUBROUTINE ConstructNewProcessGrid_onlyslice(grid, world_comm_, &
       & process_slices_in)
    !> The grid to construct
    TYPE(ProcessGrid_t), INTENT(INOUT) :: grid
    !> A communicator that every process in the grid is a part of.
    INTEGER, INTENT(IN) :: world_comm_
    !> The number of grid slices.
    INTEGER, INTENT(IN), OPTIONAL :: process_slices_in
    !! Local Data
    INTEGER :: process_rows, process_columns, process_slices
    INTEGER :: total_processors
    INTEGER :: ierr

    !! Total processors
    CALL MPI_COMM_SIZE(world_comm_, total_processors, ierr)

    !! Process Optional Parameters
    IF (PRESENT(process_slices_in)) THEN
       process_slices = process_slices_in
    ELSE
       CALL ComputeNumSlices(total_processors, process_slices)
    END IF

    !! Create a 3D grid
    CALL ComputeGridSize(total_processors, process_slices, process_rows, &
         & process_columns)

    !! Now call the full setup
    CALL ConstructNewProcessGrid(grid, world_comm_, process_rows, &
         & process_columns, process_slices)
  END SUBROUTINE ConstructNewProcessGrid_onlyslice
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Copy a process grid.
  !> Note that this makes a complete and independent copy of the process grid.
  !> Which of course means that whatever is currently stored in new_grid will
  !> be destroyed, so do not leave any matrices pointing to it.
  SUBROUTINE CopyProcessGrid(old_grid, new_grid)
    !> The grid to copy.
    TYPE(ProcessGrid_t), INTENT(IN) :: old_grid
    !> New_grid = old_grid
    TYPE(ProcessGrid_t), INTENT(INOUT) :: new_grid
    INTEGER :: II, JJ, ierr

    !! Safe Copy
    CALL DestructProcessGrid(new_grid)

    !! Copy all the basic data about a process
    new_grid%total_processors = old_grid%total_processors
    new_grid%num_process_rows = old_grid%num_process_rows
    new_grid%num_process_columns = old_grid%num_process_columns
    new_grid%num_process_slices = old_grid%num_process_slices
    new_grid%slice_size = old_grid%slice_size
    new_grid%my_slice = old_grid%my_slice
    new_grid%my_row = old_grid%my_row
    new_grid%my_column = old_grid%my_column
    new_grid%global_rank = old_grid%global_rank
    new_grid%within_slice_rank = old_grid%within_slice_rank
    new_grid%between_slice_rank = old_grid%between_slice_rank
    new_grid%column_rank = old_grid%column_rank
    new_grid%row_rank = old_grid%row_rank
    new_grid%block_multiplier = old_grid%block_multiplier
    new_grid%number_of_blocks_columns = old_grid%number_of_blocks_columns
    new_grid%number_of_blocks_rows = old_grid%number_of_blocks_rows

    !! Allocate Blocks
    ALLOCATE(new_grid%blocked_row_comm(old_grid%number_of_blocks_rows))
    ALLOCATE(new_grid%blocked_column_comm(old_grid%number_of_blocks_columns))
    ALLOCATE(new_grid%blocked_within_slice_comm(&
         & old_grid%number_of_blocks_rows, old_grid%number_of_blocks_columns))
    ALLOCATE(new_grid%blocked_between_slice_comm( &
         & old_grid%number_of_blocks_rows, old_grid%number_of_blocks_columns))

    !! Copy the communicators
    DO II = 1, new_grid%number_of_blocks_rows
       CALL MPI_COMM_DUP(old_grid%blocked_row_comm(II), &
            & new_grid%blocked_row_comm(II), ierr)
    END DO
    DO JJ = 1, new_grid%number_of_blocks_columns
       CALL MPI_COMM_DUP(old_grid%blocked_column_comm(JJ), &
            & new_grid%blocked_column_comm(JJ), ierr)
    END DO
    DO JJ=1,new_grid%number_of_blocks_columns
       DO II=1,new_grid%number_of_blocks_rows
          CALL MPI_COMM_DUP(old_grid%blocked_within_slice_comm(II,JJ), &
               & new_grid%blocked_within_slice_comm(II,JJ), ierr)
       END DO
    END DO
    DO JJ=1,new_grid%number_of_blocks_columns
       DO II=1,new_grid%number_of_blocks_rows
          CALL MPI_COMM_DUP(old_grid%blocked_between_slice_comm(II,JJ), &
               & new_grid%blocked_between_slice_comm(II,JJ), ierr)
       END DO
    END DO

    CALL MPI_COMM_DUP(old_grid%global_comm, new_grid%global_comm, ierr)
    CALL MPI_COMM_DUP(old_grid%row_comm, new_grid%row_comm, ierr)
    CALL MPI_COMM_DUP(old_grid%column_comm, new_grid%column_comm, ierr)
    CALL MPI_COMM_DUP(old_grid%within_slice_comm, new_grid%within_slice_comm, &
         & ierr)
    CALL MPI_COMM_DUP(old_grid%between_slice_comm, &
         & new_grid%between_slice_comm, ierr)

  END SUBROUTINE CopyProcessGrid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Destruct a process grid.
  !> Be careful about doing this. Matrices have pointers to process grids. If
  !> you destruct a process grid without destructing the matrices pointing to
  !> it, they will become unusable.
  RECURSIVE SUBROUTINE DestructProcessGrid(grid_in)
    !> The grid to destruct. If none specified this destroys the global grid.
    TYPE(ProcessGrid_t), OPTIONAL, INTENT(INOUT) :: grid_in
    !! Counters
    INTEGER :: II, JJ
    INTEGER :: Ierr

    !! Handle optional parameters
    IF (.NOT. PRESENT(grid_in)) THEN
       CALL DestructProcessGrid(global_grid)
    ELSE !! Destruct
       IF (ALLOCATED(grid_in%blocked_row_comm)) THEN
          CALL MPI_COMM_FREE(grid_in%global_comm, ierr)
          CALL MPI_COMM_FREE(grid_in%row_comm, ierr)
          CALL MPI_COMM_FREE(grid_in%column_comm, ierr)
          CALL MPI_COMM_FREE(grid_in%within_slice_comm, ierr)
          CALL MPI_COMM_FREE(grid_in%between_slice_comm, ierr)
          DO II = 1, grid_in%number_of_blocks_rows
             CALL MPI_COMM_FREE(grid_in%blocked_row_comm(II), ierr)
          END DO
          DEALLOCATE(grid_in%blocked_row_comm)
       END IF

       IF (ALLOCATED(grid_in%blocked_column_comm)) THEN
          DO JJ = 1, grid_in%number_of_blocks_columns
             CALL MPI_COMM_FREE(grid_in%blocked_column_comm(JJ), ierr)
          END DO
          DEALLOCATE(grid_in%blocked_column_comm)
       END IF

       IF (ALLOCATED(grid_in%blocked_within_slice_comm)) THEN
          DO JJ=1,grid_in%number_of_blocks_columns
             DO II=1,grid_in%number_of_blocks_rows
                CALL MPI_COMM_FREE(grid_in%blocked_within_slice_comm(II,JJ), &
                     & ierr)
             END DO
          END DO
          DEALLOCATE(grid_in%blocked_within_slice_comm)
       END IF

       IF (ALLOCATED(grid_in%blocked_between_slice_comm)) THEN
          DO JJ=1,grid_in%number_of_blocks_columns
             DO II=1,grid_in%number_of_blocks_rows
                CALL MPI_COMM_FREE(grid_in%blocked_between_slice_comm(II,JJ), &
                     & ierr)
             END DO
          END DO
          DEALLOCATE(grid_in%blocked_between_slice_comm)
       END IF
    END IF

  END SUBROUTINE DestructProcessGrid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Given a process grid, this splits it into two grids of even size
  SUBROUTINE SplitProcessGrid(old_grid, new_grid, my_color, split_slice, &
       & between_grid_comm)
    !> The old grid to split
    TYPE(ProcessGrid_t), INTENT(INOUT) :: old_grid
    !> The new grid that we are creating
    TYPE(ProcessGrid_t), INTENT(INOUT) :: new_grid
    !> A color value indicating which set this process was split into
    INTEGER, INTENT(OUT) :: my_color
    !> True if we were able to split along slices.
    LOGICAL, INTENT(OUT) :: split_slice
    !> A communicator for sending messages between groups.
    INTEGER, INTENT(OUT) :: between_grid_comm
    !! Local Variables - new grid
    INTEGER :: new_comm
    INTEGER :: rows, cols, slices
    INTEGER :: midpoint
    !! For Between Comm
    INTEGER :: between_color, between_rank
    INTEGER :: left_grid_size
    !! Temporary
    INTEGER :: ierr

    split_slice = .FALSE.
    !! Handle base case
    IF (old_grid%total_processors .EQ. 1) THEN
       rows = 1
       cols = 1
       slices = 1
       my_color = 0
       between_rank = 0
       !! First preferentially try to split along slices
    ELSE IF (old_grid%num_process_slices .GT. 1) THEN
       midpoint = old_grid%num_process_slices/2
       cols = old_grid%num_process_columns
       rows = old_grid%num_process_rows
       IF (old_grid%my_slice .LT. midpoint) THEN
          my_color = 0
          slices = midpoint
       ELSE
          my_color = 1
          slices = old_grid%num_process_slices - midpoint
       END IF
       between_rank = old_grid%my_slice
       split_slice = .TRUE.
       left_grid_size = midpoint*cols*rows
       !! Next try to split the bigger direction
    ELSE IF (old_grid%num_process_rows .GT. old_grid%num_process_columns) THEN
       midpoint = old_grid%num_process_rows/2
       cols = old_grid%num_process_columns
       slices = 1
       IF (old_grid%my_row .LT. midpoint) THEN
          my_color = 0
          rows = midpoint
       ELSE
          my_color = 1
          rows = old_grid%num_process_rows - midpoint
       END IF
       between_rank = old_grid%my_row
       left_grid_size = midpoint*cols*slices
       !! Default Case
    ELSE
       midpoint = old_grid%num_process_columns/2
       slices = 1
       rows = old_grid%num_process_rows
       IF (old_grid%my_column .LT. midpoint) THEN
          my_color = 0
          cols = midpoint
       ELSE
          my_color = 1
          cols = old_grid%num_process_columns - midpoint
       END IF
       between_rank = old_grid%my_column
       left_grid_size = midpoint*slices*rows
    END IF

    !! Construct
    CALL MPI_COMM_SPLIT(old_grid%global_comm, my_color, old_grid%global_rank, &
         & new_comm, ierr)
    CALL ConstructNewProcessGrid(new_grid, new_comm, rows, cols, slices)

    !! For sending data between grids
    between_color = MOD(new_grid%global_rank, left_grid_size)
    CALL MPI_COMM_SPLIT(old_grid%global_comm, between_color, between_rank, &
         & between_grid_comm, ierr)

  END SUBROUTINE SplitProcessGrid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Check if the current process is the root process.
  FUNCTION IsRoot(grid) RESULT(is_root)
    !! Parameters
    !> The process grid.
    TYPE(ProcessGrid_t), INTENT(IN) :: grid
    !> True if the current process is root.
    LOGICAL :: is_root
    IF (grid%global_rank == 0) THEN
       is_root = .TRUE.
    ELSE
       is_root = .FALSE.
    END IF
  END FUNCTION IsRoot
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Get the slice of the current process.
  FUNCTION GetMySlice(grid) RESULT(return_val)
    !> The process grid.
    TYPE(ProcessGrid_t), INTENT(IN), OPTIONAL :: grid
    !> Slice number of the current process.
    INTEGER :: return_val

    IF (PRESENT(grid)) THEN
       return_val = grid%my_slice
    ELSE
       return_val = global_grid%my_slice
    END IF
  END FUNCTION GetMySlice
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Get the column of the current process.
  FUNCTION GetMyColumn(grid) RESULT(return_val)
    !> The process grid.
    TYPE(ProcessGrid_t), INTENT(IN), OPTIONAL :: grid
    !> The column number of the current process.
    INTEGER :: return_val

    IF (PRESENT(grid)) THEN
       return_val = grid%my_column
    ELSE
       return_val = global_grid%my_column
    END IF
  END FUNCTION GetMyColumn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Get the row of the current process.
  FUNCTION GetMyRow(grid) RESULT(return_val)
    !> The process grid.
    TYPE(ProcessGrid_t), INTENT(IN), OPTIONAL :: grid
    !> The row number of the current process.
    INTEGER :: return_val

    IF (PRESENT(grid)) THEN
       return_val = grid%my_row
    ELSE
       return_val = global_grid%my_row
    END IF
  END FUNCTION GetMyRow
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Sometimes we only want to specify for a process grid the number of slices
  !> and then automatically compute the right number of rows and columns.
  SUBROUTINE ComputeGridSize(total_processors, set_slices, rows, columns)
    !> Total processors in the grid
    INTEGER, INTENT(IN) :: total_processors
    !> Desired number of slices
    INTEGER, INTENT(IN) :: set_slices
    !> Computed number of rows
    INTEGER, INTENT(OUT) :: rows
    !> Computed number of columns
    INTEGER, INTENT(OUT) :: columns
    !! Local variables
    INTEGER :: slice_size, size_search
    INTEGER :: II

    rows   = 1
    columns = 1
    slice_size = total_processors/set_slices
    size_search = FLOOR(SQRT(REAL(slice_size)))
    DO II=size_search,1,-1
       IF (MOD(slice_size,II) .EQ. 0) THEN
          rows = II
          columns = slice_size/II
          EXIT
       END IF
    END DO

  END SUBROUTINE ComputeGridSize
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Pick an appropriate number of process slices for this calculation.
  !> This routine will focus on whether we can make a valid process grid with
  !> several slices.
  SUBROUTINE ComputeNumSlices(total_processors, slices)
    !> Total processors in the grid.
    INTEGER, INTENT(IN) :: total_processors
    !> Number of slices to use.
    INTEGER, INTENT(OUT) :: slices
    !! Local Variables
    INTEGER :: slice_size
    INTEGER :: slice_dim
    LOGICAL :: found

    !! Try manually values [4, 3, 2]. If they do not work, give up and use 1.
    found = .FALSE.
    DO slices = MIN(4, total_processors), 2, -1
       slice_size = total_processors / slices
       IF (slice_size * slices .NE. total_processors) CYCLE

       !! First try a square grid.
       slice_dim  = FLOOR(SQRT(REAL(slice_size)))
       IF (slice_dim*slice_dim .EQ. slice_size) THEN
          FOUND = .TRUE.
          EXIT
       END IF

       !! If not, we try a grid where the rows are twice the number of columns.
       slice_dim = FLOOR(SQRT(REAL(slice_size/2)))
       IF (slice_dim*slice_dim*2 .EQ. slice_size) THEN
          FOUND = .TRUE.
          EXIT
       END IF
    END DO
    IF (.NOT. FOUND) slices = 1

  END SUBROUTINE ComputeNumSlices
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE ProcessGridModule


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine modaldynamics(inpc,textpart,nmethod,tinc,tper,iexpl,
     &  istep,istat,n,iline,ipol,inl,ipoinp,inp,iperturb,isolver,
     &  cs,mcs,ipoinpc,idrct,ctrl,tmin,tmax,nforc,nload,nbody,iprestr,
     &  t0,t1,ithermal,nk,vold,veold,xmodal,set,nset,mi,cyclicsymmetry)
!
!     reading the input deck: *MODAL DYNAMIC
!
      implicit none
!
      logical steadystate,nodalset
!
      character*1 inpc(*)
      character*20 solver
      character*81 set(*),noset
      character*132 textpart(16)
!
      integer nmethod,istep,istat,n,key,iexpl,iline,ipol,inl,
     &  ipoinp(2,*),inp(3,*),iperturb(2),isolver,i,mcs,ipoinpc(0:*),
     &  idrct,nforc,nload,nbody,iprestr,ithermal,j,nk,ipos,nset,mi(*),
     &  cyclicsymmetry
!
      real*8 tinc,tper,cs(17,*),ctrl(*),tmin,tmax,t0(*),t1(*),
     &  vold(0:mi(2),*),veold(0:mi(2),*),xmodal(*)
!
      iexpl=0
      iperturb(1)=0
      iperturb(2)=0
      idrct=1
      tmin=0.d0
      tmax=0.d0
      steadystate=.false.
      if((mcs.ne.0).and.(cs(2,1).ge.0.d0)) then
         cyclicsymmetry=1
c      else
c         cyclicsymmetry=0
      endif
      nodalset=.false.
!
      if(istep.lt.1) then
         write(*,*) '*ERROR in modaldynamics: *MODAL DYNAMIC can only'
         write(*,*) '  be used within a STEP'
         call exit(201)
      endif
!
!     default solver
!
      solver='                    '
      if(isolver.eq.0) then
         solver(1:7)='SPOOLES'
      elseif(isolver.eq.2) then
         solver(1:16)='ITERATIVESCALING'
      elseif(isolver.eq.3) then
         solver(1:17)='ITERATIVECHOLESKY'
      elseif(isolver.eq.4) then
         solver(1:3)='SGI'
      elseif(isolver.eq.5) then
         solver(1:5)='TAUCS'
      elseif(isolver.eq.7) then
         solver(1:7)='PARDISO'
      endif
!
      do i=2,n
         if(textpart(i)(1:7).eq.'SOLVER=') then
            read(textpart(i)(8:27),'(a20)') solver
         elseif(textpart(i)(1:9).eq.'DIRECT=NO') then
            idrct=0
         elseif(textpart(i)(1:7).eq.'DELTMX=') then
            read(textpart(i)(8:27),'(f20.0)',iostat=istat) ctrl(27)
         elseif(textpart(i)(1:11).eq.'STEADYSTATE') then
            steadystate=.true.
c         elseif(textpart(i)(1:14).eq.'CYCLICSYMMETRY') then
c            cyclicsymmetry=.true.
c         elseif(textpart(i)(1:5).eq.'NSET=') then
c            nodalset=.true.
c            noset=textpart(i)(6:85)
c            noset(81:81)=' '
c            ipos=index(noset,' ')
c            noset(ipos:ipos)='N'
         else
            write(*,*)
     &        '*WARNING in modaldynamics: parameter not recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*MODAL DYNAMIC%")
         endif
      enddo
!
      if(solver(1:7).eq.'SPOOLES') then
         isolver=0
      elseif(solver(1:16).eq.'ITERATIVESCALING') then
         write(*,*) '*WARNING in modaldynamics: the iterative scaling'
         write(*,*) '         procedure is not available for modal'
         write(*,*) '         dynamic calculations; the default solver'
         write(*,*) '         is used'
      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
         write(*,*) '*WARNING in modaldynamics: the iterative scaling'
         write(*,*) '         procedure is not available for modal'
         write(*,*) '         dynamic calculations; the default solver'
         write(*,*) '         is used'
      elseif(solver(1:3).eq.'SGI') then
         isolver=4
      elseif(solver(1:5).eq.'TAUCS') then
         isolver=5
c      elseif(solver(1:13).eq.'MATRIXSTORAGE') then
c         isolver=6
      elseif(solver(1:7).eq.'PARDISO') then
         isolver=7
      else
         write(*,*) '*WARNING in modaldynamics: unknown solver;'
         write(*,*) '         the default solver is used'
      endif
!
      if((isolver.eq.2).or.(isolver.eq.3)) then
         write(*,*) '*ERROR in modaldynamics: the default solver ',
     & solver
         write(*,*) '       cannot be used for modal dynamic'
         write(*,*) '       calculations '
         call exit(201)
      endif
!
c      if(nodalset) then
c         do i=1,nset
c            if(set(i).eq.noset) exit
c         enddo
c         if(i.gt.nset) then
c            noset(ipos:ipos)=' '
c            write(*,*) '*ERROR in modaldynamics: node set ',noset
c            write(*,*) '  has not yet been defined.'
c            call exit(201)
c         endif
c         xmodal(10)=i+0.5d0
c      else
c         if(cyclicsymmetry) then
c            write(*,*) '*ERROR in modaldynamics: cyclic symmetric'
c            write(*,*) '       structure, yet no node set defined'
c            call exit(201)
c         endif
c      endif
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &     ipoinp,inp,ipoinpc)
      if((istat.lt.0).or.(key.eq.1)) then
         write(*,*) '*ERROR in modaldynamics: definition not complete'
         write(*,*) '       '
         call inputerror(inpc,ipoinpc,iline,
     &"*MODAL DYNAMIC%")
         call exit(201)
      endif
      read(textpart(1)(1:20),'(f20.0)',iostat=istat)tinc
      if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*MODAL DYNAMIC%")
      read(textpart(2)(1:20),'(f20.0)',iostat=istat)tper
      if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*MODAL DYNAMIC%")
      read(textpart(3)(1:20),'(f20.0)',iostat=istat) tmin
      if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*MODAL DYNAMIC%")
      read(textpart(4)(1:20),'(f20.0)',iostat=istat) tmax
      if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*MODAL DYNAMIC%")
!
      if(steadystate) then
!
!        modal dynamics calculation till steady state
!
         if(tper.le.0.d0) then
            write(*,*) '*ERROR in modaldynamics: relative error'
            write(*,*) '       is nonpositive'
            call exit(201)
         endif
         tper=-tper
         if(tinc.le.0.d0) then
            write(*,*) '*ERROR in modaldynamics: initial increment'
            write(*,*) '       size is nonpositive'
            call exit(201)
         endif
         if(tmin.lt.0.d0) then
            tmin=1.d-10
         endif
         if(tmax.lt.1.d-10) then
            tmax=1.d+30
         endif
      else
!
!        transient modal dynamics calculation
!
         if(tper.lt.0.d0) then
            write(*,*) '*ERROR in modaldynamics: step size is negative'
            call exit(201)
         elseif(tper.le.0.d0) then
            tper=1.d0
         endif
         if(tinc.lt.0.d0) then
            write(*,*) '*ERROR in modaldynamics: initial increment size
     &is negative'
            call exit(201)
         elseif(tinc.le.0.d0) then
            tinc=tper
         endif
         if(tinc.gt.tper) then
            write(*,*) '*ERROR in modaldynamics: initial increment size
     &exceeds step size'
            call exit(201)
         endif
!
         if(idrct.ne.1) then
            if(tmin.lt.1.d-10*tper) then
               tmin=min(tinc,1.d-10*tper)
            endif
            if(tmax.lt.1.d-10) then
               tmax=1.d+30
            endif
         endif
      endif
!
!     removing the present loading
!
      nforc=0
      nload=0
      nbody=0
      iprestr=0
      if((ithermal.eq.1).or.(ithermal.eq.3)) then
         do j=1,nk
            t1(j)=t0(j)
         enddo
      endif
!
!     resetting fields vold and veold after a frequency or
!     buckling step
!
      if((nmethod.eq.2).or.(nmethod.eq.3)) then
         do i=1,nk
            do j=1,3
               vold(j,i)=0.d0
               veold(j,i)=0.d0
            enddo
         enddo
      endif
!
      nmethod=4
!
!     correction for cyclic symmetric structures:
!     if the present step was not preceded by a frequency step
!     no nodal diameter has been selected. To make sure that
!     mastructcs is called instead of mastruct a fictitious
!     minimum nodal diameter is stored
!
      if((cyclicsymmetry.eq.1).and.(mcs.ne.0).and.(cs(2,1)<0.d0))
     &       cs(2,1)=0.d0
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &     ipoinp,inp,ipoinpc)
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      SUBROUTINE E01001 (USEPRM)

C  E01001 tests the handling of error number 1: function requires
C  state (PHCL,WSCL,STCL,ARCL).

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

      INTEGER   USEPRM
      LOGICAL   STREQ

      CURCON = 'PHIGS is open'
      CALL SETVS ('1', EXPERR,EXPSIZ)
      ERRSRS = '1'
      CALL ESETUP (USEPRM)

      CALL ROPPH (ERRFIL, MEMUN)
      CALL TSTIGN (STREQ('O***'))

      CALL ENDERR

      END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a
!!!!   copy of this software and associated documentation files (the "Software"),
!!!!   to deal in the Software without restriction, including without limitation
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!!!!   and/or sell copies of the Software, and to permit persons to whom the
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/master/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module Test_Matrix_mod

    use Matrix_mod
    use Err_mod, only: Err_type
    use Test_mod, only: Test_type
    implicit none

    private
    public :: test_Matrix

    type(Test_type) :: Test

!***********************************************************************************************************************************
!***********************************************************************************************************************************

contains

!***********************************************************************************************************************************
!***********************************************************************************************************************************

    subroutine test_Matrix()

        implicit none

        Test = Test_type(moduleName=MODULE_NAME)

        call test_sortPosDefMat()
        call test_getRegresCoef()
        call Test%finalize()
#ifdef CAF_ENABLED
        sync all
#endif

    end subroutine test_Matrix

!***********************************************************************************************************************************
!***********************************************************************************************************************************

    subroutine test_sortPosDefMat()

        use Constants_mod, only: RK, IK

        implicit none
        integer                 :: rank
        real(RK), allocatable   :: PosDefMat(:,:), OutPosDefMat(:,:), OutPosDefMatRef(:,:)
        integer                 :: i,j

        if (Test%Image%isFirst) call Test%testing("sortPosDefMat")

        rank = 5
        allocate(PosDefMat(rank,rank))
        do j = 1,rank
            do i = 1,j
                PosDefMat(i,j) = i*10 + j
            end do
        end do

        ! switch the variable 4 with 5, such that the output remains a positive-definite matrix.
        OutPosDefMat = sortPosDefMat(rank,PosDefMat,1_IK,[4],[5])

        if (Test%isDebugMode .and. Test%Image%isFirst) then
            write(Test%outputUnit,"(*(g0))")
            write(Test%outputUnit,"(*(g0))") "Original Matrix:"
            do i = 1,rank
                write(Test%outputUnit,"(*(F7.1))") (PosDefMat(i,j),j=1,rank)
            end do
            write(Test%outputUnit,"(*(g0))")
            write(Test%outputUnit,"(*(g0))") "Output Positive-Definite Matrix with variables 4 and 5 swapped:"
            do i = 1,rank
                write(Test%outputUnit,"(*(F7.1))") (OutPosDefMat(i,j),j=1,rank)
            end do
            write(Test%outputUnit,"(*(g0))")
        end if

        OutPosDefMatRef = reshape(  [ 11._RK, 12._RK, 13._RK, 15._RK, 14._RK &
                                    ,  0._RK, 22._RK, 23._RK, 25._RK, 24._RK &
                                    ,  0._RK,  0._RK, 33._RK, 35._RK, 34._RK &
                                    ,  0._RK,  0._RK,  0._RK, 55._RK, 45._RK &
                                    ,  0._RK,  0._RK,  0._RK,  0._RK, 44._RK ] &
                                    , shape = [rank,rank] )
        OutPosDefMatRef = transpose(OutPosDefMatRef)
        do j = 1, rank
            do i = 1, j
                Test%assertion = OutPosDefMatRef(i,j) == OutPosDefMat(i,j)
               !write(*,*) OutPosDefMatRef(i,j) , OutPosDefMat(i,j)
                call Test%verify()
            end do
        end do
        !call Test%skipping()

    end subroutine test_sortPosDefMat

!***********************************************************************************************************************************
!***********************************************************************************************************************************

    subroutine test_getRegresCoef()

        use Constants_mod, only: RK, IK

        implicit none
        real(RK)                ::  normalizedDifference
        integer , parameter     ::  rankPDM = 4, rankS11 = 3, rankS22 = 1
        real(RK), parameter     ::  PosDefMat(rankPDM,rankPDM) = reshape( &
                                    [ 4.414182620515998_RK, 1.173760167060120_RK, 0.757607629189287_RK, 5.075277296976230_RK &
                                    , 1.173760167060120_RK, 0.866956750091570_RK, 0.310654936099342_RK, 1.621274787164182_RK &
                                    , 0.757607629189287_RK, 0.310654936099342_RK, 0.955157221699132_RK, 1.254186231887444_RK &
                                    , 5.075277296976230_RK, 1.621274787164182_RK, 1.254186231887444_RK, 6.407791808961157_RK ] &
                                    , shape=shape(PosDefMat) )
        real(RK), parameter     ::  RegresCoefMatRef(rankS11,rankS22) = reshape( &
                                    [ 0.792047783119072 &
                                    , 0.253016145889269 &
                                    , 0.195728305363088 ] &
                                    , shape=shape(RegresCoefMatRef) )
        real(RK), parameter     ::  SchurComplementRef(rankS11,rankS11) = reshape( &
                                    [ 0.3943204887314210_RK, -0.110366933940115_RK, -0.235767795395625_RK &
                                    , -0.110366933940115_RK,  0.456748052015843_RK, -0.006674430520204_RK &
                                    , -0.235767795395625_RK, -0.006674430520204_RK,  0.709677475922085_RK ] &
                                    , shape=shape(SchurComplementRef) )
        real(RK)                :: SchurComplement(rankS11,rankS11)
        real(RK)                :: RegresCoefMat(rankS11,rankS22)
        integer                 :: i,j

        if (Test%Image%isFirst) call Test%testing("getRegresCoef")

        call getRegresCoef  ( rankPDM           = rankPDM           &
                            , rankS11           = rankS11           &
                            , rankS22           = rankS22           &
                            , PosDefMat         = PosDefMat         &
                            , RegresCoefMat     = RegresCoefMat     &
                            , SchurComplement   = SchurComplement   &
                            )

        if (Test%isDebugMode .and. Test%Image%isFirst) then
            write(Test%outputUnit,"(*(g0))")
            write(Test%outputUnit,"(*(g0))") "Original Covariance Matrix:"
            do i = 1,rankPDM
                write(Test%outputUnit,"(*(F22.15))") (PosDefMat(i,j),j=1,rankPDM)
            end do
        end if

        if (Test%isDebugMode .and. Test%Image%isFirst) then
            write(Test%outputUnit,"(*(g0))"); write(Test%outputUnit,"(*(g0))") "RegresCoefMat, RegresCoefMatRef, difference, normalizedDifference:"
        end if
        do i = 1,rankS11
            !write(Test%outputUnit,"(*(F22.15))") (RegresCoefMat(i,j),j=1,rankS22), (RegresCoefMatRef(i,j),j=1,rankS22)
            do j = 1,rankS22
                normalizedDifference = abs(RegresCoefMatRef(i,j) - RegresCoefMat(i,j)) / (RegresCoefMatRef(i,j) + RegresCoefMat(i,j))
                if (Test%isDebugMode .and. Test%Image%isFirst) then
                    write(Test%outputUnit,"(*(F22.15))" ) RegresCoefMat(i,j) &
                                                        , RegresCoefMatRef(i,j) &
                                                        , abs(RegresCoefMat(i,j)-RegresCoefMatRef(i,j)) &
                                                        , normalizedDifference
                end if
                Test%assertion = normalizedDifference < 1.e-5_RK
                call Test%verify()
            end do
        end do

        if (Test%isDebugMode .and. Test%Image%isFirst) then
            write(Test%outputUnit,"(*(g0))"); write(Test%outputUnit,"(*(g0))") "SchurComplement, SchurComplementRef, difference, normalizedDifference:"
        end if
        do i = 1,rankS11
            do j = 1,rankS11
                normalizedDifference = abs(SchurComplementRef(i,j) - SchurComplement(i,j)) / (SchurComplementRef(i,j) + SchurComplement(i,j))
                if (Test%isDebugMode .and. Test%Image%isFirst) then
                    write(Test%outputUnit,"(*(F22.15))" ) SchurComplement(i,j) &
                                                        , SchurComplementRef(i,j) &
                                                        , abs(SchurComplement(i,j)-SchurComplementRef(i,j)) &
                                                        , normalizedDifference
                end if
                Test%assertion = normalizedDifference < 1.e-6_RK
                call Test%verify()
            end do
        end do

        write(Test%outputUnit,"(*(g0))")

        !call Test%skipping()

    end subroutine test_getRegresCoef

!***********************************************************************************************************************************
!***********************************************************************************************************************************

end module Test_Matrix_mod


text:
JasmineMassoumi can describe himself or herself here.

20100719 22:27:28 nbsp Hello and welcome to the Wiki!  Users/hankim

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Abaqus utility modules
include 'abaqus_utils_mod.f90'          ! Do not include when running Abaqus
!include 'abaqus_utils_dummy_mod.f90'    ! Include when running Abaqus
include 'includes.f90'

program test_usub
use find_mod
implicit none
    double precision, allocatable   :: tmp(:)
    integer, allocatable            :: tmp_int(:)
    double precision, allocatable   :: utmp(:)
    integer                         :: a2d(3,5)
    integer                         :: a1d(5)
    integer                         :: indices(2)

    a1d = [11, 12, 13, 14, 15]
    a2d(1,:) = [11, 12, 13, 14, 15]
    a2d(2,:) = [21, 22, 23, 24, 25]
    a2d(3,:) = [31, 32, 33, 34, 35]

    indices = find(a2d, 23)

    write(*,*) find(a1d, 13)
    write(*,*) find(a2d, 23)

    write(*,*) find(a2d, 0)

end program test_usub

max_stars_repo_name: KnutAM/AbaqusRolloverSimulation

max_issues_repo_name: KnutAM/AbaqusRolloverSimulation

max_forks_repo_name: KnutAM/AbaqusRolloverSimulation


module mod_vtrans
!$$$   module documentation block
!                .      .    .                                       .
! module:    mod_vtrans      contains all routines for vertical modes
!
!   prgmmr: parrish          org: np23                date: 2006-06-26
!
! abstract:  contains routines and variables for vertical mode
!            generation and use in vertical transforms
!
! program history log:
!   2006-06-26
!   2007-05-08   kleist - finish vertical coordinate generalization
!   2011-07-04  todling  - fixes to run either single or double precision
!   2012-11-19  parrish - compute vertical mode information on pe 0 only, and
!                          then broadcast results to other processors.
!                          This fixes a problem with reproducibility encountered on tide.
!   2012-12-18  parrish - problem with reproducibility reappeared.  Even with all vertical mode
!                           computations done on pe 0, the library routine dgeev supplied on the WCOSS
!                           machine uses threading internally in ways that prevent reproducibility,
!                           even when computation is nominally restricted to one processor.  To fix this,
!                           available local subroutines within gsi src directory have been used in place of
!                           dgeev.  Subroutine eigen, located in raflib.f90, computes eigenvalues and
!                           eigenvectors of symmetric matrices.  The matrix qmat is not symmetric, but a
!                           rescaling produces an almost symmetric matrix.  The symmetric average of the
!                           rescaled qmat has eigenvalues and eigenvectors that differ from those of qmat
!                           by only a few percent.  Using these as starting values, inverse iteration can
!                           retrieve the desired eigenvalues and eigenvectors of qmat.  To reproduce the
!                           eigenvalues and eigenvectors of qmat obtained from dgeev to 15 digits, the
!                           inverse iteration refinement of the starting eigenvalues and eigenvectors is
!                           done in quad precision.  This requires a quad precision matrix inversion
!                           routine.  A copy of subroutine iminv in get_semimp_mats.f90 has been copied
!                           on the end of this module and given the name iminv_quad.  This is also used
!                           to get the inverse of qmat, instead of the previous method which was to
!                           construct the inverse of qmat from its eigenvalues and eigenvectors.  This
!                           required all of the eigenvalues and eigenvectors, which dgeev provided.
!                           However, by using iminv_quad, now only the first nvmodes_keep eigenvectors and
!                           eigenvalues need be computed.  The complete computation if the eigenvectors
!                           and eigenvalues without using dgeev uses less than 2 seconds for 8
!                           eigenvalue/vectors and all related computations.
!   2018-02-15  wu       - add code for fv3_regional option
!
! subroutines included:
!   sub init_vtrans              - initialize vertical mode related variables
!   sub create_vtrans            - allocate and load vert mode variables
!   sub destroy_vtrans           - deallocate vert mode variables
!   sub getabc                   -
!   sub vtrans_inv               - physical space u,v,T,p --> vert transformed u,v,phi
!   sub vtrans_inv_ad            - adjoint of vtrans_inv
!   sub vtrans                   - vert transformed u,v,phi --> physical space u,v,T,p
!   sub vtrans_ad                - adjoint of vtrans
!
! Variable Definitions:
!   def nvmodes_keep             - number of vertical modes to keep ( <= nsig )
!   def speeds                   - phase speeds of vertical modes
!   def vmodes                   - vertical modes
!   def dualmodes                - dual vertical modes
!   def phihat2t                 - matrix operator to convert phihat to T
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds,only: r_kind,i_kind
  implicit none

! set default to private
  private
! set subroutines to public
  public :: init_vtrans
  public :: create_vtrans
  public :: destroy_vtrans
  public :: getabc
  public :: vtrans_inv
  public :: vtrans_inv_ad
  public :: vtrans
  public :: vtrans_ad
! set passed variables to public
  public :: nvmodes_keep,depths,speeds

  integer(i_kind) nvmodes_keep
  real(r_kind),dimension(:),allocatable:: depths,speeds
  real(r_kind),dimension(:,:),allocatable:: vmodes,phihat2t,dualmodes
  real(r_kind),dimension(:,:),allocatable:: t2phihat
  real(r_kind),dimension(:),allocatable:: p2phihat,phihat2p

  character(len=*),parameter::myname='mod_vtrans'

contains

  subroutine init_vtrans
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_vtrans     set default values for vtrans variables
!   prgmmr: parrish         org: np23                date: 2006-06-26
!
! abstract: set default values for vtrans variables
!
! program history log:
!   2006-06-26  parrish
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
    implicit none

    nvmodes_keep=0

  end subroutine init_vtrans

  subroutine create_vtrans(mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    create_vtrans   get vertical functions for dynamic constraint
!   prgmmr: parrish        org: np23         date:  2006-06-26
!
! abstract:  using linearization of dynamics about rest state, derive
!             coupling matrix and obtain eigenvectors/values
!             (linearization follows Juang, 2005, NCEP Office Note 445)
!            The conversion from mass variable to T,p follows the
!            Machenhauer-Phillips approach, pointed out by R. Errico
!            in a conversation on 7-11-2006.
!
!           Briefly, the linear equations used to define the vertical
!            modes are
!
!             dD/dt = -laplacian ( H*p + A*T )                          (1)
!
!             dT/dt = -B*D                                              (2)
!
!             dp/dt = -S*D                                              (3)
!
!                 where D is divergence, p is surface pressure, T is virtual temperature,
!
!                    and the matrices H, A, B, S are as defined in ON 445.
!
!           Taking the time derivative of (1) and substituting from (2) and (3)
!           yields an equation for just divergence,
!
!              d2D/dt2 = -Q * laplacian(D)
!
!                 Q = H*S + A*B
!
!           The vertical modes are the right eigenvectors of Q and the
!           scale geopotential values for each vertical mode are the eigenvalues of Q.
!
!                 Q = U*E*V(transpose)
!
!           To transform from physical space to vertical modes, first form
!           the mass variable
!
!                 M = H*p + A*T
!
!           Then the transform variables are
!
!            (Mhat,uhat,vhat) = V(transpose)*(M,u,v)
!
!           To return from mode space to physical space, we have
!
!             (M,u,v) = U*(Mhat,uhat,vhat)
!
!           Finally, to get T,p from M using the Machenhauer-Phillips approach,
!
!              T = B*Q**(-1)*M
!
!              p = S*Q**(-1)*M
!
!           The above is only strictly valid for T and p small perturbations in gravity modes
!           only, but that is the application for which this code is intended.
!
!
! program history log:
!   2006-06-26  parrish
!   2007-02-26  yang    - replace IBM subroutine of dgeev by GSI dgeev
!   2010-04-01  treadon - move strip to gridmod
!   2012-11-19  parrish - compute vertical mode information on pe 0 only, and
!                          then broadcast results to other processors.
!                          This fixes a problem with reproducibility encountered on tide.
!   2012-12-18  parrish - substantial rewrite of this code to allow use of all local code.
!                          r_quad precision has been used in eigenvalue/eigenvector computation
!                          to get agreement with eigenvalues/eigenvectors from library routine dgeev,
!                          which has been replaced in this version.  It appears that reproducibility
!                          has been restored.  (See more extensive comment at top of this module.)
!   2013-10-19  todling - metguess now holds background
!
!
! usage:
!   input argument list:
!       mype     - current processor number
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
    use kinds, only: r_quad
    use constants,only: zero,one,one_tenth,ten
    use gridmod,only: lat2,lon2,nsig,nlat,nlon
    use mpimod,only: mpi_rtype,mpi_comm_world,ierror,mpi_integer
    use guess_grids, only: ntguessig
    use general_sub2grid_mod, only: general_sub2grid
    use general_commvars_mod, only: g1
    use gsi_metguess_mod, only: gsi_metguess_bundle
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use mpeu_util, only: die
    use gsi_io, only: verbose
    implicit none

!   Declare passed variables
    integer(i_kind),intent(in   ) :: mype

!   Declare local variables
    character(len=*),parameter::myname_=myname//'*create_vtrans'
    integer(i_kind) i,j,k,n
    real(r_kind) count,psbar
    real(r_kind) sum
    real(r_kind),dimension(1,lat2,lon2,1):: worksub
    real(r_kind),allocatable,dimension(:,:,:,:):: hwork
    real(r_kind),dimension(nsig+1)::pbar
    real(r_kind),dimension(nsig)::tbar
    real(r_kind),dimension(nsig+1)::ahat,bhat,chat
    real(r_kind),dimension(nsig)::hmat,smat,sqmatinv
    real(r_kind),dimension(nsig,nsig)::amat,bmat,qmat,qmatinv,bqmatinv
    real(r_kind),allocatable,dimension(:):: swww,swwwd
    real(r_kind),allocatable,dimension(:,:):: szzz,szzzd
    integer(i_kind),allocatable:: numlevs(:)
    real(r_kind),dimension(:,:  ),pointer:: ges_ps_nt=>NULL()
    real(r_kind),dimension(:,:,:),pointer:: ges_tv_nt=>NULL()
    integer(i_kind) workpe,ier,istatus
    integer(i_kind) lpivot(nsig),mpivot(nsig)
    real(r_quad) qmatinv_quad(nsig,nsig),detqmat_quad
    logical print_verbose

!   get work pe:

    print_verbose=.false.
    allocate(numlevs(0:g1%npe-1))
    numlevs(0:g1%npe-1)=g1%kend(0:g1%npe-1)-g1%kbegin(0:g1%npe-1)+1
    if(g1%mype==0) then
        workpe=-1
        do i=0,g1%npe-1
   !       write(6,*)' i,numlevs(i)=',i,numlevs(i)
           if(numlevs(i)==1) workpe=i
        end do
   !    write(6,*)' workpe=',workpe
    end if
    call mpi_bcast(workpe,1,mpi_integer,0,mpi_comm_world,ierror)
   !write(6,*)' mype,workpe=',mype,workpe
    if(verbose .and. g1%mype==workpe) print_verbose=.true.

!    obtain vertical coordinate constants ahat,bhat,chat
    if(g1%mype==workpe) call getabc(ahat,bhat,chat)

!   get global mean T and ps

    allocate(hwork(g1%inner_vars,g1%nlat,g1%nlon,g1%kbegin_loc:g1%kend_alloc))

!   count:
!  Not clear if area weighting would be better.
    count=one/float(nlat*nlon)

    ier=0
    call gsi_bundlegetpointer (gsi_metguess_bundle(ntguessig),'ps',ges_ps_nt, istatus)
    ier=ier+istatus
    call gsi_bundlegetpointer (gsi_metguess_bundle(ntguessig),'tv',ges_tv_nt, istatus)
    ier=ier+istatus
    if(ier/=0) call die(myname_,'missing require variables, ier=', ier)

!   psbar:

    do j=1,lon2
       do i=1,lat2
          worksub(1,i,j,1)=ges_ps_nt(i,j)
       end do
    end do
    call general_sub2grid(g1,worksub,hwork)
    if(g1%mype==workpe) then
       psbar=zero
       do j=1,nlon
          do i=1,nlat
             psbar=psbar+hwork(1,i,j,1)
          end do
       end do
       psbar=ten*count*psbar
       do k=1,nsig+1
          pbar(k)=ahat(k)+bhat(k)*psbar     !  + chat(k)*(T/T0)**(1/kappa)   --- add later
       end do
    end if

!   tbar:

    do k=1,nsig
       do j=1,lon2
          do i=1,lat2
             worksub(1,i,j,1)=ges_tv_nt(i,j,k)
          end do
       end do
       call general_sub2grid(g1,worksub,hwork)
       if(g1%mype==workpe) then
          tbar(k)=zero
          do j=1,nlon
             do i=1,nlat
                tbar(k)=tbar(k)+hwork(1,i,j,1)
             end do
          end do
          tbar(k)=count*tbar(k)
       end if
    end do

    if(print_verbose) then
       do k=1,nsig
          write(6,'(" k,pbar,tbar = ",i5,2f18.9)')k,pbar(k),tbar(k)
       end do
       k=nsig+1
       write(6,'(" k,pbar      = ",i5,f18.9)')k,pbar(k)
    end if

!    allocate variables used in vertical mode transformations:

    allocate(depths(nvmodes_keep),speeds(nvmodes_keep))
    allocate(vmodes(nsig,nvmodes_keep),dualmodes(nsig,nvmodes_keep))
    allocate(phihat2t(nsig,nvmodes_keep),phihat2p(nvmodes_keep))
    allocate(p2phihat(nvmodes_keep),t2phihat(nsig,nvmodes_keep))

if(mype==workpe) then    ! BEGIN MYPE=workpe SECTION !!!!!!!!!!!!!

    hmat=zero ; smat=zero ; amat=zero ; bmat=zero

! Get matrices for variable transforms/vertical modes
    call get_semimp_mats(tbar,pbar,bhat,chat,amat,bmat,hmat,smat)

!   qmat = hmat*smat + amat*bmat

    do j=1,nsig
       do i=1,nsig
          qmat(i,j)=hmat(i)*smat(j)
          do k=1,nsig
             qmat(i,j)=qmat(i,j)+amat(i,k)*bmat(k,j)
          end do
       end do
    end do

!  get inverse of Q using iminv_quad

    qmatinv_quad=qmat
    call iminv_quad(qmatinv_quad,nsig,detqmat_quad,lpivot,mpivot)
    qmatinv=qmatinv_quad

!   check inverse

 !errormax=zero
 !do j=1,nsig
 !   do i=1,nsig
 !      sum=zero
 !      if(i==j) sum=-one
 !      do k=1,nsig
 !         sum=sum+qmat(i,k)*qmatinv(k,j)
 !      end do
 !      errormax=max(abs(sum),errormax)
 !   end do
 !end do
 !write(6,*)' error in qmatinv =',errormax

!     next get eigenvalues and eigenvectors.

   allocate(swww(nvmodes_keep),swwwd(nvmodes_keep))
   allocate(szzz(nsig,nvmodes_keep),szzzd(nsig,nvmodes_keep))
   call special_eigvv(qmat,hmat,smat,nsig,swww,szzz,swwwd,szzzd,nvmodes_keep)

   do k=1,nvmodes_keep
      depths(k)=swww(k)
      speeds(k)=sqrt(depths(k))
      do j=1,nsig
         vmodes(j,k)=szzz(j,k)
         dualmodes(j,k)=szzzd(j,k)
      end do
   end do

!   next compute p2phihat and t2phihat

  t2phihat=zero
  p2phihat=zero
  do n=1,nvmodes_keep
     do k=1,nsig
        do j=1,nsig
           t2phihat(k,n)=t2phihat(k,n)+szzzd(j,n)*amat(j,k)
        end do
        p2phihat(n)=p2phihat(n)+szzzd(k,n)*hmat(k)
     end do
  end do
  p2phihat=ten*p2phihat ! in this code, p is in units of mb, but in gsi, p is in cb -- change later

!   finally compute phihat2p, phihat2t

  do j=1,nsig
     do i=1,nsig
        bqmatinv(i,j)=zero
        do k=1,nsig
           bqmatinv(i,j)=bqmatinv(i,j)+bmat(i,k)*qmatinv(k,j)
        end do
     end do
  end do
  do j=1,nsig
     sqmatinv(j)=zero
     do i=1,nsig
        sqmatinv(j)=sqmatinv(j)+smat(i)*qmatinv(i,j)
     end do
  end do

  do j=1,nvmodes_keep
     sum=zero
     do k=1,nsig
        sum=sum+sqmatinv(k)*szzz(k,j)
     end do
     phihat2p(j)=sum
     do i=1,nsig
        sum=zero
        do k=1,nsig
           sum=sum+bqmatinv(i,k)*szzz(k,j)
        end do
        phihat2t(i,j)=sum
     end do
  end do
  phihat2p=one_tenth*phihat2p ! local units are mb, but gsi units are cb--fix later
end if  ! END MYPE=workpe SECTION !!!!!!!!!!!!!

!  BROADCAST RESULTS FROM ABOVE SECTION TO ALL PES

    call mpi_bcast(depths,nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(speeds,nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(vmodes,nsig*nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(phihat2t,nsig*nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(dualmodes,nsig*nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(t2phihat,nsig*nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(p2phihat,nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)
    call mpi_bcast(phihat2p,nvmodes_keep,mpi_rtype,workpe,mpi_comm_world,ierror)

  end subroutine create_vtrans

  subroutine destroy_vtrans
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_vtrans   remove space used by vtrans routines
!   prgmmr: parrish        org: np23         date:  2006-06-26
!
! abstract:  remove space used in vertical mode transformations
!
! program history log:
!   2006-06-26  parrish
!
!
! usage:
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
    implicit none

!    deallocate variables used in vertical mode transformations:

    deallocate(depths,speeds,vmodes,dualmodes)
    deallocate(phihat2t,phihat2p,p2phihat,t2phihat)

  end subroutine destroy_vtrans

  subroutine getabc(ahat,bhat,chat)
!$$$  subprogram documentation block
!
! subprogram:    getabc       get pressure constants
!
!   prgmmr: parrish          org: np22                date: 2006-05-04
!
! abstract: return constants used to get 3d pressure field at interfaces based on
!            generalized vertical coordinate
!
! program history log:
!   2006-05-04  kleist
!   2010-12-17  pagowski - add cmaq
!
! usage:
!   input argument list:
!
!   output argument list:
!     ahat       -   p(i,j,k)  = ahat(k) + bhat(k)*psfc(i,j)+chat(k)*(T(i,j,k)/T0(k))**(1/kappa)
!     bhat       -
!     chat       -
!
!$$$ end documentation block
    use constants,only: zero,ten
    use gridmod,only: nsig,ak5,bk5,ck5,fv3_regional
    use gridmod,only: wrf_nmm_regional,nems_nmmb_regional,eta1_ll,eta2_ll,pdtop_ll,pt_ll,cmaq_regional
    implicit none

!   Declare passed variables
    real(r_kind),dimension(nsig+1),intent(  out) :: ahat,bhat,chat

!   Declare local variables
    integer(i_kind) k

    if(wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional) then
       do k=1,nsig+1
          ahat(k)=eta1_ll(k)*pdtop_ll-eta2_ll(k)*(pdtop_ll+pt_ll)+pt_ll
          bhat(k)=eta2_ll(k)
          chat(k)=zero
       end do
    elseif(fv3_regional) then
       do k=1,nsig+1
          ahat(k)=eta1_ll(k)*0.01_r_kind
          bhat(k)=eta2_ll(k)
          chat(k)=zero
       end do
    else
       do k=1,nsig+1
          ahat(k)=ak5(k)*ten
          bhat(k)=bk5(k)
          chat(k)=ck5(k)*ten
       end do
    end if

    return
  end subroutine getabc

  subroutine vtrans_inv(uhat,vhat,phihat,u,v,t,p)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    vtrans_inv
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-05-29  safford -- add subprogram doc block
!
!   input argument list:
!     uhat,vhat,phihat
!
!   output argument list:
!     u,v,t
!     p
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

    use gridmod,only: lat2,lon2,nsig
    use constants,only: zero
    implicit none

    real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(in   ) :: uhat,vhat,phihat
    real(r_kind),dimension(lat2,lon2,nsig)        ,intent(  out) :: u,v,t
    real(r_kind),dimension(lat2,lon2)             ,intent(  out) :: p

    integer(i_kind) i,j,k,n

!$omp parallel do schedule(dynamic,1) private(n,i,j,k)
    do k=1,nsig
       do n=1,nvmodes_keep
          do j=1,lon2
             do i=1,lat2
                u(i,j,k)=u(i,j,k)+vmodes(k,n)*uhat(i,j,n)
                v(i,j,k)=v(i,j,k)+vmodes(k,n)*vhat(i,j,n)
                t(i,j,k)=t(i,j,k)+phihat2t(k,n)*phihat(i,j,n)
             end do
          end do
          if( k == 1)then
             do j=1,lon2
                do i=1,lat2
                   p(i,j)=p(i,j)+phihat2p(n)*phihat(i,j,n)
                end do
             end do
          end if
       end do
    end do

  end subroutine vtrans_inv

  subroutine vtrans_inv_ad(uhat,vhat,phihat,u,v,t,p)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    vtrans_inv_ad
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-05-29  safford -- add subprogram doc block
!
!   input argument list:
!     uhat,vhat,phihat
!     u,v,t
!     p
!
!   output argument list:
!     uhat,vhat,phihat
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block
    use gridmod,only: lat2,lon2,nsig
    use constants,only: zero
    implicit none

    real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(inout) :: uhat,vhat,phihat
    real(r_kind),dimension(lat2,lon2,nsig)        ,intent(in   ) :: u,v,t
    real(r_kind),dimension(lat2,lon2)             ,intent(in   ) :: p

    integer(i_kind) i,j,k,n

!$omp parallel do schedule(dynamic,1) private(n,i,j,k)
    do n=1,nvmodes_keep
       do k=1,nsig
          do j=1,lon2
             do i=1,lat2
                uhat(i,j,n)=uhat(i,j,n)+vmodes(k,n)*u(i,j,k)
                vhat(i,j,n)=vhat(i,j,n)+vmodes(k,n)*v(i,j,k)
                phihat(i,j,n)=phihat(i,j,n)+phihat2t(k,n)*t(i,j,k)
             end do
          end do
       end do
       do j=1,lon2
          do i=1,lat2
             phihat(i,j,n)=phihat(i,j,n)+phihat2p(n)*p(i,j)
          end do
       end do
    end do

  end subroutine vtrans_inv_ad

  subroutine vtrans(u,v,t,p,uhat,vhat,phihat)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    vtrans
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-05-29  safford -- add subprogram doc block
!
!   input argument list:
!     u,v,t
!     p
!
!   output argument list:
!     uhat,vhat,phihat
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

    use gridmod,only: lat2,lon2,nsig
    use constants,only: zero
    implicit none

    real(r_kind),dimension(lat2,lon2,nsig)        ,intent(in   ) :: u,v,t
    real(r_kind),dimension(lat2,lon2)             ,intent(in   ) :: p
    real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(  out) :: uhat,vhat,phihat

    integer(i_kind) i,j,k,n

!$omp parallel do schedule(dynamic,1) private(n,i,j,k)
    do n=1,nvmodes_keep
       do j=1,lon2
          do i=1,lat2
             uhat(i,j,n)=zero
             vhat(i,j,n)=zero
             phihat(i,j,n)=p2phihat(n)*p(i,j)
          enddo
       enddo
       do k=1,nsig
          do j=1,lon2
             do i=1,lat2
                uhat(i,j,n)=uhat(i,j,n)+dualmodes(k,n)*u(i,j,k)
                vhat(i,j,n)=vhat(i,j,n)+dualmodes(k,n)*v(i,j,k)
                phihat(i,j,n)=phihat(i,j,n)+t2phihat(k,n)*t(i,j,k)
             end do
          end do
       end do
    end do

  end subroutine vtrans

  subroutine vtrans_ad(u,v,t,p,uhat,vhat,phihat)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    vtrans_ad
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-05-29  safford -- add subprogram doc block
!
!   input argument list:
!     u,v,t
!     p
!     uhat,vhat,phihat
!
!   output argument list:
!     u,v,t
!     p
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

    use gridmod,only: lat2,lon2,nsig
    implicit none

    real(r_kind),dimension(lat2,lon2,nsig)        ,intent(inout) :: u,v,t
    real(r_kind),dimension(lat2,lon2)             ,intent(inout) :: p
    real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(in   ) :: uhat,vhat,phihat

    integer(i_kind) i,j,k,n

!$omp parallel do schedule(dynamic,1) private(n,i,j,k)
    do k=1,nsig
       do n=1,nvmodes_keep
          do j=1,lon2
             do i=1,lat2
                u(i,j,k)=u(i,j,k)+dualmodes(k,n)*uhat(i,j,n)
                v(i,j,k)=v(i,j,k)+dualmodes(k,n)*vhat(i,j,n)
                t(i,j,k)=t(i,j,k)+t2phihat(k,n)*phihat(i,j,n)
             end do
          end do
          if(k == 1)then
             do j=1,lon2
               do i=1,lat2
                  p(i,j)=p(i,j)+p2phihat(n)*phihat(i,j,n)
               end do
             end do
          end if
       end do
    end do

  end subroutine vtrans_ad

end module mod_vtrans

subroutine special_eigvv(qmat0,hmat0,smat0,nmat,swww0,szzz0,swwwd0,szzzd0,nvmodes_keep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    special_eigvv   obtain eigenvalues and eigenvecgtors of qmat0
!   prgmmr: parrish        org: np22         date:  2012-12-18
!
! abstract:  Compute eigenvalues and eigenvectors of matrix qmat0.  Input and output
!            variables are r_kind, but internal calculations are done in r_quad to eliminate
!            possible failures in the inverse iteration algorithm used to refine
!            the accuracy of the eigenvectors and eigenvalues so as to agree
!            with those produced previously by library routine dgeev to 15
!            digits.  Because it is only necessary that the initial estimate of
!            eigenvectors and eigenvalues be near (but not too near) the desired
!            final values, the initial values are obtained using r_kind
!            precision only (subroutine eigen, located in raflib.F90).
!            Extensive testing with different initial eigenvalue/vector
!            estimates confirmed that identical r_kind precision
!            eigenvalues and eigenvectors are always obtained from the quad precision
!            inverse iteration algorithm.  This translates to exact
!            reproducibility on different numbers of processors for the full
!            analysis.
!
!
! program history log:
!   2012-12-18  parrish
!
!
! usage:
!   input argument list:
!     qmat0,hmat0,smat0,nmat,nvmodes_keep
!
!   output argument list:
!     swww0,szzz0,swwwd0,szzzd0
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_quad
  use gsi_io, only: verbose

  implicit none

  integer(i_kind),intent(in)::nmat,nvmodes_keep
  real(r_kind),intent(in)::qmat0(nmat,nmat),hmat0(nmat),smat0(nmat)
  real(r_kind),intent(inout):: swww0(nvmodes_keep),swwwd0(nvmodes_keep)
  real(r_kind),intent(inout):: szzz0(nmat,nvmodes_keep),szzzd0(nmat,nvmodes_keep)

  real(r_quad) qmat(nmat,nmat),hmat(nmat),smat(nmat)
  real(r_quad) swww(nvmodes_keep),swwwd(nvmodes_keep)
  real(r_quad) szzz(nmat,nvmodes_keep),szzzd(nmat,nvmodes_keep)
  real(r_quad) rmat(nmat),qtildemat(nmat,nmat),atemp(nmat*nmat),btemp(nmat,nmat)
  real(r_quad) eigvals(nmat)
  integer(i_kind) i,j,k,ia,mv,iret,istop
  real(r_quad) sum
! real(r_quad) orthoerror
  real(r_quad) aminv(nmat,nmat),aminvt(nmat,nmat)
  real(r_quad) eigval_this,eigval_next
  real(r_quad) zero_quad,half_quad,one_quad
  real(r_quad) dist_to_closest_eigval,perturb_eigval
  real(r_quad) err_aminv
  real(r_kind) atemp8(nmat*nmat),btemp8(nmat,nmat)
  logical noskip,print_verbose

  print_verbose = .false.
  if(verbose) print_verbose=.true.
  zero_quad=0.0_r_quad
  half_quad=0.5_r_quad
  one_quad =1.0_r_quad
  qmat=qmat0
  hmat=hmat0
  smat=smat0
  do i=1,nmat
     rmat(i)=sqrt(hmat(i)*smat(i))
  end do
  do j=1,nmat
     do i=1,nmat
        qtildemat(i,j)=rmat(i)*rmat(j)*qmat(i,j)/(hmat(i)*smat(j))
     end do
  end do

!  get eigenvalues, eigenvectors of symmetrized version
!    of qtildemat (gives most accurate eigenvalue estimates):
  atemp=0.0_r_quad
  do i=1,nmat
     do j=i,nmat
        ia=i+(j*j-j)/2
        atemp(ia)=half_quad*(qtildemat(j,i)+qtildemat(i,j))
     end do
  end do
  mv=0
  atemp8=atemp
  call eigen(atemp8,btemp8,nmat,mv)
  btemp=btemp8
  atemp=atemp8
  do i=1,nvmodes_keep
     j=i
     ia=i+(j*j-j)/2
     eigvals(i)=atemp(ia)
  !  write(6,*)' i,ave eigvals(i)=',i,eigvals(i)
  end do

!      iterative improvement to get eigvals, eigvects of qtildemat

  do i=1,nvmodes_keep
     szzz(:,i)=btemp(:,i)
     szzzd(:,i)=btemp(:,i)
!  renormalize szzz and szzzd
     sum=zero_quad
     do j=1,nmat
        sum=sum+szzz(j,i)**2
     end do
     sum=one_quad/sqrt(sum)
     if(szzz(1,i)<zero_quad) sum=-sum
     do j=1,nmat
        szzz(j,i)=sum*szzz(j,i)
     end do
     sum=zero_quad
     do j=1,nmat
        sum=sum+szzzd(j,i)**2
     end do
     sum=one_quad/sqrt(sum)
     if(szzzd(1,i)<zero_quad) sum=-sum
     do j=1,nmat
        szzzd(j,i)=sum*szzzd(j,i)
     end do
     eigval_this=eigvals(i)
     iret=0
     dist_to_closest_eigval=min( abs(eigvals(i)-eigvals(max(i-1,1           ))), &
                                 abs(eigvals(i)-eigvals(min(i+1,nvmodes_keep))) )
     if(i==1) then
        dist_to_closest_eigval=abs(eigvals(i+1)-eigvals(i))
     elseif(i==nvmodes_keep) then
        dist_to_closest_eigval=abs(eigvals(i)-eigvals(i-1))
     else
        dist_to_closest_eigval=min( abs(eigvals(i)-eigvals(max(i-1,1           ))), &
                                    abs(eigvals(i)-eigvals(min(i+1,nvmodes_keep))) )
     end if

     !write(6,*)' i,eigvals(i),dist_to_closest_eigval=',i,eigvals(i),dist_to_closest_eigval
     perturb_eigval=0.000001_r_quad*dist_to_closest_eigval
     !write(6,*)' i,perturb_eigval/eigvals(i)=',i,perturb_eigval/eigvals(i)
     !write(6,*)' precomputation of aminv with check on err_aminv'
     call iterative_improvement0(qtildemat,eigval_this,aminv,aminvt,nmat,iret,err_aminv)
     noskip=.false.
     if(err_aminv>10._r_quad**(-20).or.iret==1) then
        noskip=.true.
        eigval_this=eigval_this+perturb_eigval
     end if
     iret=0
     do j=1,10
        if(noskip)  call iterative_improvement0(qtildemat,eigval_this,aminv,aminvt,nmat,iret,err_aminv)
        noskip=.true.
        if(iret==1) then
           if(print_verbose)write(6,*)' det=0 in iterative_improvement0, eigenvalue converged, it = ',j
           exit
        end if
        call iterative_improvement(eigval_this,eigval_next,aminv,aminvt,szzz(:,i),szzzd(:,i),nmat,istop)
        if(eigval_this==eigval_next) then
           if(print_verbose)write(6,*)' no change in eigenvalue, convergence to machine precision achieved, it = ',j
           exit
        end if
        eigval_this=eigval_next
        if(istop==1) then
           if(print_verbose)write(6,*)' eigval relative change less than 10**(-24), no further iteration necessary, it = ',j
           exit
        end if
     end do
     swww(i)=eigval_next
     swwwd(i)=eigval_next
     !write(6,*)' i,j,swww(i)=',i,j,swww(i)
  end do

!  compute unscaled left and right eigenvectors of "corrected" matrix qmat_c

  do k=1,nvmodes_keep
     do i=1,nmat
        szzz(i,k)=hmat(i)*szzz(i,k)/rmat(i)
        szzzd(i,k)=smat(i)*szzzd(i,k)/rmat(i)
     end do
  end do

!  renormalize szzz and szzzd

  do i=1,nvmodes_keep
     sum=zero_quad
     do j=1,nmat
        sum=sum+szzz(j,i)**2
     end do
     sum=one_quad/sqrt(sum)
     if(szzz(1,i)<zero_quad) sum=-sum
     do j=1,nmat
        szzz(j,i)=sum*szzz(j,i)
     end do
     sum=zero_quad
     do j=1,nmat
        sum=sum+szzzd(j,i)**2
     end do
     sum=one_quad/sqrt(sum)
     if(szzzd(1,i)<zero_quad) sum=-sum
     do j=1,nmat
        szzzd(j,i)=sum*szzzd(j,i)
     end do
  end do

!  joint normalization:

  do i=1,nvmodes_keep
     sum=zero_quad
     do j=1,nmat
        sum=sum+szzz(j,i)*szzzd(j,i)
     end do
     sum=one_quad/sqrt(abs(sum))
     do j=1,nmat
        szzz(j,i)=sum*szzz(j,i)
        szzzd(j,i)=sum*szzzd(j,i)
     end do
  end do

 !orthoerror=zero_quad
 !do i=1,nvmodes_keep
 !   do j=1,nvmodes_keep
 !      sum=zero_quad
 !      if(i==j) sum=-one_quad
 !      do k=1,nmat
 !         sum=sum+szzz(k,i)*szzzd(k,j)
 !      end do
 !      orthoerror=max(abs(sum),orthoerror)
 !   end do
 !end do
 !write(6,*)' orthoerror for szzz,szzzd=',orthoerror

! check error in qmat*szzz - swww*szzz  and qmat_trans*szzzd - swww*szzzd

 !do i=1,nvmodes_keep
 !   sum=zero_quad
 !   rnorm=zero_quad
 !   sumd=zero_quad
 !   rnormd=zero_quad
 !   do j=1,nmat
 !      term=-swww(i)*szzz(j,i)
 !      term2=term
 !      do k=1,nmat
 !         term=term+qmat(j,k)*szzz(k,i)
 !      end do
 !      sum=sum+term**2
 !      rnorm=rnorm+term2**2
 !      term=-swww(i)*szzzd(j,i)
 !      term2=term
 !      do k=1,nmat
 !         term=term+qmat(k,j)*szzzd(k,i)
 !      end do
 !      sumd=sumd+term**2
 !      rnormd=rnormd+term2**2
 !   end do
 !   write(6,*)' i,eigval/vect error check: vec,dual errs=',i,sqrt(sum/rnorm),sqrt(sumd/rnormd)
 !end do
  swww0 =swww
  swwwd0=swwwd
  szzz0 =szzz
  szzzd0=szzzd
  !do i=1,nvmodes_keep
  !   write(6,*)' i,swww0(i),swwwd0(i)=',i,swww0(i),swwwd0(i)
  !   do j=1,nmat,nmat-1
  !      write(6,*)' i,j,szzz0(j,i),szzzd0(j,i)=',i,j,szzz0(j,i),szzzd0(j,i)
  !   end do
  !end do

end subroutine special_eigvv

subroutine iterative_improvement0(a,mu,aminv,aminvt,na,iret,errormax)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    iterative_improvement0  compute inverse of a - mu*I
!   prgmmr: parrish        org: np22         date:  2012-12-18
!
! abstract:  Compute inverse of matrix a - mu*I, where mu is current estimate of eigenvalue
!            in inverse iteration algorithm.
!
!
! program history log:
!   2012-12-18  parrish
!
!
! usage:
!   input argument list:
!     a,mu,na
!
!   output argument list:
!     aminv,aminvt,iret
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_quad,i_kind
  use gsi_io, only: verbose
  implicit none

  integer(i_kind),intent(in)::na
  real(r_quad),intent(in)::a(na,na),mu
  real(r_quad),intent(inout)::aminv(na,na)
  real(r_quad),intent(inout)::aminvt(na,na)
  real(r_quad),intent(out):: errormax
  integer(i_kind),intent(out)::iret

  real(r_quad) am(na,na)
  real(r_quad) sum,detam,errlimit
  integer(i_kind) i,j,k,lpivot(na),mpivot(na)
  real(r_quad) zero_quad,one_quad
  logical print_verbose

  print_verbose=.false.
  if(verbose)print_verbose=.true.
  errlimit=1000._r_quad
  zero_quad=0._r_quad
  one_quad=1._r_quad

!  compute (A-mu*I)**(-1)

  iret=0

  am=zero_quad
  do j=1,na
     am(j,j)=-mu
  end do
  do j=1,na
     do i=1,na
        am(i,j)=am(i,j)+a(i,j)
     end do
  end do
  aminv=am
  call iminv_quad(aminv,na,detam,lpivot,mpivot)
  if(detam==zero_quad) then
     iret=1
     return
  end if
  do j=1,na
     do i=1,na
        aminvt(i,j)=aminv(j,i)
     end do
  end do

!   check inverse

  errormax=zero_quad
  do j=1,na
     do i=1,na
        sum=zero_quad
        if(i==j) sum=-one_quad
        do k=1,na
           sum=sum+am(i,k)*aminv(k,j)
        end do
        errormax=max(abs(sum),errormax)
     end do
  end do
  if(print_verbose)write(6,*)' delta in aminv =',errormax
  if(errormax>errlimit) iret=-1

end subroutine iterative_improvement0

subroutine iterative_improvement(mu,mu_next,aminv,aminvt,bv,bw,na,istop)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    iterative_improvement  complete one iteration of inverse iteration
!   prgmmr: parrish        org: np22         date:  2012-12-18
!
! abstract:  Apply one iteration of inverse iteration, producing a new estimate for eigenvalue
!            and eigenvector.
!
!
! program history log:
!   2012-12-18  parrish
!
!
! usage:
!   input argument list:
!     mu,aminv,aminvt,bv,bw,na
!
!   output argument list:
!     mu_next,bv,bw,istop
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_quad,i_kind
  implicit none

  integer(i_kind),intent(in)::na
  integer(i_kind),intent(inout)::istop
  real(r_quad),intent(in)::mu
  real(r_quad),intent(in)::aminv(na,na),aminvt(na,na)
  real(r_quad),intent(inout):: mu_next
  real(r_quad),intent(inout)::bv(na),bw(na)

  real(r_quad) bnextv(na),bnextw(na)
  real(r_quad) factorv,factorw,sumv,sumw
  real(r_quad) zero_quad,half_quad,one_quad,two_quad,tol_quad
  integer(i_kind) i,j

  istop=0
  zero_quad=0.0_r_quad
  half_quad=0.5_r_quad
  one_quad=1.0_r_quad
  two_quad=2.0_r_quad
  tol_quad=10.0_r_quad**(-24)
  bnextv=zero_quad
  bnextw=zero_quad
  do i=1,na
     do j=1,na
        bnextv(i)=bnextv(i)+aminv(i,j)*bv(j)
        bnextw(i)=bnextw(i)+aminvt(i,j)*bw(j)
     end do
  end do

!  compute dot prod of bnext with b

  sumv=zero_quad
  sumw=zero_quad
  do i=1,na
     sumv=sumv+bnextv(i)*bv(i)
     sumw=sumw+bnextw(i)*bw(i)
  end do
  mu_next=mu+half_quad/sumv+half_quad/sumw
  if(two_quad*abs(mu_next-mu)/abs(mu+mu_next) < tol_quad) istop=1

!   normalize bnextv and bnextw

  factorv=zero_quad
  factorw=zero_quad
  do i=1,na
     factorv=factorv+bnextv(i)**2
     factorw=factorw+bnextw(i)**2
  end do
  factorv=one_quad/sqrt(factorv)
  factorw=one_quad/sqrt(factorw)
  if(bnextv(1)<zero_quad) factorv=-factorv
  if(bnextw(1)<zero_quad) factorw=-factorw
  do i=1,na
     bv(i)=factorv*bnextv(i)
     bw(i)=factorw*bnextw(i)
  end do

end subroutine iterative_improvement

      subroutine iminv_quad (a,n,d,l,m)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    iminv    invert a matrix
!
!   prgrmmr:
!
! abstract:      the standard gauss-jordan method is used. the determinant
!                is also calculated. a determinant of zero indicates that
!                the matrix is singular.
!
! remarks        matrix a must be a general matrix
!
! program history log:
!   2008-06-04  safford -- add subprogram doc block
!   2012-12-18  parrish -- create r_quad precision version
!   2013-01-10  parrish -- change so d is absolute value of determinant and
!                            is bounded by 1e50.  The variable d is only used
!                            in this application if it is exactly zero, so there
!                            is no harm in bounding by a large number.  This was
!                            done to prevent overflow trap in debug mode.
!
!   input argument list:
!     a - input matrix, destroyed in computation and replaced by resultant inverse
!     n - order of matrix a
!     d - resultant determinant
!     l - work vector of length n
!     m - work vector of length n
!
!   output argument list:
!     a - input matrix, destroyed in computation and replaced by resultant inverse
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

        use kinds,only: r_quad,i_kind
        implicit none

        integer(i_kind)             ,intent(in   ) :: n
        integer(i_kind),dimension(n),intent(inout) :: l,m
        real(r_quad)                ,intent(inout) :: d
        real(r_quad),dimension(n*n) ,intent(inout) :: a

        integer(i_kind):: nk,k,j,iz,i,ij
        integer(i_kind):: kj,ik,jr,jq,jk,ki,kk,jp,ji

        real(r_quad):: biga,hold
        real(r_quad):: zero_quad,one_quad
        real(r_quad) dmax
!
!        if a double precision version of this routine is desired, the
!        ! in column 1 should be removed from the double precision
!        statement which follows.
!
!     double precision a, d, biga, hold
!
!        the ! must also be removed from double precision statements
!        appearing in other routines used in conjunction with this
!        routine.
!
!        the double precision version of this sr........ must also
!        contain double precision fortran functions.  abs in statemen
!        10 must be changed to dabs  .
!
!        ...............................................................

         zero_quad=0.0_r_quad
         one_quad=1.0_r_quad
         dmax=(10.0_r_quad)**50
!
!        search for largest element
!
         d=one_quad
         nk=-n
         do 80 k=1,n
            nk=nk+n
            l(k)=k
            m(k)=k
            kk=nk+k
            biga=a(kk)
            do 20 j=k,n
               iz=n*(j-1)
               do 20 i=k,n
                  ij=iz+i
   10             if(abs(biga)-abs(a(ij))) 15,20,20
   15             biga=a(ij)
                  l(k)=i
                  m(k)=j
   20       continue
!
!        interchange rows
!
            j=l(k)
            if(j-k) 35,35,25
   25       ki=k-n
            do 30 i=1,n
               ki=ki+n
               hold=-a(ki)
               ji=ki-k+j
               a(ki)=a(ji)
   30          a(ji) =hold
!
!        interchange columns
!
   35          i=m(k)
               if(i-k) 45,45,38
   38          jp=n*(i-1)
               do 40 j=1,n
                  jk=nk+j
                  ji=jp+j
                  hold=-a(jk)
                  a(jk)=a(ji)
   40             a(ji) =hold
!
!        divide column by minus pivot (value of pivot element is
!        contained in biga)
!
   45             if(biga) 48,46,48
   46             d=zero_quad
                  return
   48             do 55 i=1,n
                     if(i-k) 50,55,50
   50                ik=nk+i
                     a(ik)=a(ik)/(-biga)
   55             continue
!
!        reduce matrix
!
                  do 65 i=1,n
                     ik=nk+i
                     ij=i-n
                     do 65 j=1,n
                        ij=ij+n
                        if(i-k) 60,65,60
   60                   if(j-k) 62,65,62
   62                   kj=ij-i+k
                        a(ij)=a(ik)*a(kj)+a(ij)
   65             continue
!
!        divide row by pivot
!
                  kj=k-n
                  do 75 j=1,n
                     kj=kj+n
                     if(j-k) 70,75,70
   70                a(kj)=a(kj)/biga
   75             continue
!
!        product of pivots
!
                  d=min(dmax,abs(d*biga))
!
!        replace pivot by reciprocal
!
                  a(kk)=one_quad/biga
   80    continue
!
!        final row and column interchange
!
         k=n
         do
  100       k=(k-1)
            if(k) 150,150,105
  105       i=l(k)
            if(i-k) 120,120,108
  108       jq=n*(k-1)
            jr=n*(i-1)
            do 110 j=1,n
               jk=jq+j
               hold=a(jk)
               ji=jr+j
               a(jk)=-a(ji)
               a(ji) =hold
  110       continue
  120       j=m(k)
            if(j-k) 100,100,125
  125       ki=k-n
            do 130 i=1,n
               ki=ki+n
               hold=a(ki)
               ji=ki-k+j
               a(ki)=-a(ji)
               a(ji) =hold
  130       continue
         end do
  150    return
      end subroutine iminv_quad


text:
      SUBROUTINE MF_QUARTC
     &          (b4,four b3,six b2,four b1,one b0,q)

c     Solution for the roots of a fourth-order polynomial (quartic
c     equation) taken from Burnside and Panton (1904), The Theory of
c     Equations.  A summary of pertinent equations is given in Radio
c     Science Vol. 3, Aug. 1968, pp.792-795.

      character*200 error_msg
      real          mag pos,mag neg
      complex       b4,four b3,six b2,four b1,one b0,q,
     &              h,i,g,h prime,g prime,
     &              sq root,p pos,p,log p,cub rt0,cub rt1,cub rt2,
     &              omega1,omega2,root p,root q,root r,twlv b3,twlv b2,
     &              f,dfdq,del q,fb4
      complex  * 16 b3,b2,b1,b0,b3 sq

      dimension     p ri(2),q(4)

      equivalence  (p ri(1),p)

      data          omega1/(-0.5, 0.8660254038)/,
     &              omega2/(-0.5,-0.8660254038)/

      b3=four b3/(4.*b4)
      b2=six b2/(6.*b4)
      b1=four b1/(4.*b4)
      b0=one b0/b4

      b3 sq=b3**2
      h=b2-b3 sq
      i=b0-4.*b3*b1+3.*b2**2
      g=b1+b3*(-3.*b2+2.*b3 sq)
      h prime=-i/12.
      g prime=-g**2*.25-h*(h**2+3.*h prime)

      sq root=SQRT(g prime**2+4.*h prime**3)
      p=(-g prime+sq root)*0.5
      mag pos=ABS(p ri(1))+ABS(p ri(2))
      p pos=p
      p=(-g prime-sq root)*0.5
      mag neg=ABS(p ri(1))+ABS(p ri(2))
      if (mag pos .gt. mag neg) p=p pos
      log p=LOG(p)
      cub rt0=EXP(log p/3.)
      cub rt1=omega1*cub rt0
      cub rt2=omega2*cub rt0

      root p=SQRT(cub rt0-h prime/cub rt0-h)
      root q=SQRT(cub rt1-h prime/cub rt1-h)
      root r=SQRT(cub rt2-h prime/cub rt2-h)
      if (ABS(g) .ge. 1.e-20) then
         sign=-root p*root q*root r*2./g
         if (sign .lt. 0.) root r=-root r
      end if
      q(1)=+root p+root q+root r-b3
      q(2)=+root p-root q-root r-b3
      q(3)=-root p+root q-root r-b3
      q(4)=-root p-root q+root r-b3

c     Newton-Raphson (first order) iterative improvement.
      fb4=4.*b4
      twlv b3=3.*four b3
      twlv b2=2.*six b2
      do j=1,4
         dlqmin=9.9e9
         do jj=1,4
            if (jj .ne. j) then
               dlq=ABS(q(j)-q(jj))
               if (dlq .lt. dlqmin) dlqmin=dlq
             end if
         end do
         dlqmax=dlqmin/3.

         last=0
         ncount=1
         do while (last .eq. 0 .and. ncount .le. 10)
            f=(((b4*q(j)+four b3)*q(j)+six b2)*q(j)+four b1)*q(j)+one b0
            dfdq=((fb4*q(j)+twlv b3)*q(j)+twlv b2)*q(j)+four b1
            del q=-f/dfdq
            if (ABS(delq) .gt. dlqmax) delq=delq*dlqmax/ABS(delq)
            q(j)=q(j)+del q
            if (ABS(del q/q(j)) .lt. 1.e-4) last=1
            ncount=ncount+1
         end do
         if (ABS(del q/q(j)) .gt. 1.e-2) then
             write(error_msg,
     &           '(''[MF_QUARTC]: q='',4(1pe13.3,e11.3))') q
             call LWPC_ERROR ('WARNING',error_msg)
             write(error_msg,
     &            '(''[MF_QUARTC]: q fails to converge'')')
             call LWPC_ERROR ('ERROR',error_msg)
         end if
      end do
      RETURN
      END      ! MF_QUARTC

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC