text:
! Copyright (c) 2017-2019 Etienne Descamps
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
!    this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation and/or
!    other materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be
!    used to endorse or promote products derived from this software without specific prior
!    written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
! INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
! EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Module mlf_utils
  Use ieee_arithmetic
  Use iso_c_binding
  Use mlf_cfuns
  IMPLICIT NONE
  PRIVATE

  Public :: Xchange, QSortIdx, QSort, Mean, PrintMatrix, mlf_countcl, mlf_reversecl, SetIf
  Public :: InitOrDefault, FirstEltVal, PresentAndTrue, OutsideBounds
  Public :: mlf_reverseid, mlf_cumulativefromvect, mlf_di_search, mlf_printmatrix_c, mlf_isSorted
  Public :: Median, AddCstr

  Interface FirstEltVal
    module procedure mlf_firstIntVal
    module procedure mlf_firstDoubleVal
  End Interface FirstEltVal

  Interface SetIf
    module procedure mlf_setIfInt
    module procedure mlf_setIfInt64
    module procedure mlf_setIfDouble
  End Interface SetIf
  Interface Xchange
    module procedure mlf_xchangeInt
    module procedure mlf_xchangeInt64
    module procedure mlf_xchangeDouble
    module procedure mlf_xchangeComplex4
    module procedure mlf_xchangeComplex8
    module procedure mlf_xchangeFloat
  End Interface Xchange
  Interface QSortIdx
    module procedure mlf_sortM
    module procedure mlf_sortV
  End Interface QSortIdx
  Interface QSort
    module procedure mlf_sortOutM
    module procedure mlf_sortOutV
  End Interface QSort
  Interface Mean
    module procedure mlf_MeanV
    module procedure mlf_MeanM
    module procedure mlf_MeanMV
  End Interface Mean
  Interface PrintMatrix
    module procedure mlf_printmatrix
    module procedure mlf_printmatrix2
  End Interface PrintMatrix
  Interface InitOrDefault
    module procedure mlf_initOrDefault_bool
    module procedure mlf_initOrDefault_int
    module procedure mlf_initOrDefault_real
    module procedure mlf_initOrDefault_double
    module procedure mlf_initOrDefault_int64
  End Interface InitOrDefault
  Interface Median
    module procedure mlf_median_real8
  End Interface Median
  real(c_double), public, parameter :: mlf_PI = ACOS(-1d0)
Contains
  Real(c_double) Function mlf_median_real8(X, idX) Result(Median)
    real(c_double), intent(in) :: X(:)
    integer, intent(in) :: idX(:)
    integer :: N
    N = SIZE(X)
    If(BTEST(N,0)) Then
      Median = X(idX((N+1)/2))
    Else
      Median = 0.5d0*(X(idX(N/2))+X(idX(N/2+1)))
    Endif
  End Function mlf_median_real8

  Elemental Real(c_double) Function OutsideBounds(X, XMin, XMax) Result(cstr)
    real(c_double), intent(in) :: X, XMin, XMax
    If(X < XMin) Then
      cstr = XMin - X
    Else If(X > XMax) Then
      cstr = XMax - X
    Else
      cstr = 0
    Endif
  End Function OutsideBounds

  Integer Function AddCstr(V, N, i) Result(k)
    integer, intent(inout) :: V(:)
    integer, intent(in) :: N, i
    If(ANY(V(1:N) == i)) Then
      k = N
    Else
      V(N+1) = i
      k = N+1
    Endif
  End Function AddCstr

  Logical Function PresentAndTrue(x) Result(y)
    logical, optional :: x
    y = .FALSE.
    If(PRESENT(x)) y = x
  End Function PresentAndTrue

  Pure Subroutine mlf_initOrDefault_bool(var, def, opt)
    logical, intent(out) :: var
    logical, intent(in) :: def
    logical, intent(in), optional :: opt
    If(PRESENT(opt)) Then
      var = opt
    Else
      var = def
    Endif
  End Subroutine mlf_initOrDefault_bool

  Pure Subroutine mlf_initOrDefault_int(var, def, opt)
    integer, intent(out) :: var
    integer, intent(in) :: def
    integer, intent(in), optional :: opt
    if(present(opt)) then
      var = opt
    else
      var = def
    endif
  End Subroutine mlf_initOrDefault_int

  Pure Subroutine mlf_initOrDefault_int64(var, def, opt)
    integer(8), intent(out) :: var
    integer(8), intent(in) :: def
    integer(8), intent(in), optional :: opt
    if(present(opt)) then
      var = opt
    else
      var = def
    endif
  End Subroutine mlf_initOrDefault_int64

  Pure Subroutine mlf_initOrDefault_real(var, def, opt)
    real, intent(out) :: var
    real, intent(in) :: def
    real, intent(in), optional :: opt
    if(present(opt)) then
      var = opt
    else
      var = def
    endif
  End Subroutine mlf_initOrDefault_real

  Pure Subroutine mlf_initOrDefault_double(var, def, opt)
    real(8), intent(out) :: var
    real(8), intent(in) :: def
    real(8), intent(in), optional :: opt
    if(present(opt)) then
      var = opt
    else
      var = def
    endif
  End Subroutine mlf_initOrDefault_double

  Pure Subroutine mlf_setIfInt(dest, def, val)
    integer(c_int), intent(out) :: dest
    integer(c_int), intent(in) :: def
    integer(c_int), intent(in), optional :: val
    if(present(val)) then
      dest = val
    else
      dest = def
    endif
  End Subroutine mlf_setIfInt

  Pure Subroutine mlf_setIfInt64(dest, def, val)
    integer(c_int64_t), intent(out) :: dest
    integer(c_int64_t), intent(in) :: def
    integer(c_int64_t), intent(in), optional :: val
    if(present(val)) then
      dest = val
    else
      dest = def
    endif
  End Subroutine mlf_setIfInt64

  Pure Subroutine mlf_setIfDouble(dest, def, val)
    real(c_double), intent(out) :: dest
    real(c_double), intent(in) :: def
    real(c_double), intent(in), optional :: val
    if(present(val)) then
      dest = val
    else
      dest = def
    endif
  End Subroutine mlf_setIfDouble

  Pure Integer Function mlf_firstIntVal(array, val) Result(k)
    integer(c_int32_t), intent(in) :: array(:), val
    integer :: i, N
    N = size(array)
    k = -1 ! Not found
    Do i = 1,N
      if(array(i) == val) then
        k = i
        RETURN
      endif
    End Do
  End Function mlf_firstIntVal

  Pure Integer Function mlf_firstDoubleVal(array, val) Result(k)
    real(c_double), intent(in) :: array(:), val
    integer :: i, N
    N = size(array)
    k = -1 ! Not found
    Do i = 1,N
      if(array(i) == val) then
        k = i
        RETURN
      endif
    End Do
  End Function mlf_firstDoubleVal

  ! Utility functions for exchanging to variables
  Pure Elemental Subroutine mlf_xchangeInt(a,b)
    integer(c_int32_t), intent(inout) :: a,b
    integer(c_int32_t) :: k
    k = a
    a = b
    b = k
  End Subroutine mlf_xchangeInt

  ! Utility functions for exchanging to variables
  Pure Elemental Subroutine mlf_xchangeInt64(a,b)
    integer(c_int64_t), intent(inout) :: a,b
    integer(c_int64_t) :: k
    k = a
    a = b
    b = k
  End Subroutine mlf_xchangeInt64

  Pure Elemental Subroutine mlf_xchangeDouble(a,b)
    real(c_double), intent(inout) :: a,b
    real(c_double) :: k
    k = a
    a = b
    b = k
  End Subroutine mlf_xchangeDouble

  Pure Elemental Subroutine mlf_xchangeFloat(a,b)
    real(c_float), intent(inout) :: a,b
    real(c_float) :: k
    k = a
    a = b
    b = k
  End Subroutine mlf_xchangeFloat

  Pure Elemental Subroutine mlf_xchangeComplex4(a,b)
    complex(kind = 4), intent(inout) :: a,b
    complex(kind = 4) :: k
    k = a
    a = b
    b = k
  End Subroutine mlf_xchangeComplex4

  Pure Elemental Subroutine mlf_xchangeComplex8(a,b)
    complex(kind = 8), intent(inout) :: a,b
    complex(kind = 8) :: k
    k = a
    a = b
    b = k
  End Subroutine mlf_xchangeComplex8

  ! Test function written to test QSort function
  Pure Logical Function mlf_isSorted(Y) result(T)
    real(c_double), intent(in) :: Y(:)
    integer :: i
    T=.FALSE.
    do i=1,size(Y)-1
      if(Y(i) > Y(i+1)) RETURN
    end do
    T=.TRUE.
  End Function mlf_isSorted

  ! Quick sort algorithm that uses FORTRAN90 array structures
  ! This method sort the index array idx (Y(idx) will generate a sorted Y array)
  ! Mu is the minimal number of sorted element in the begin of the array
  Subroutine mlf_sortV(Y, idx, mu)
    real(c_double), intent(in), contiguous :: Y(:)
    integer, intent(in), optional :: mu
    integer(c_int), intent(out) :: idx(:)
    integer(c_int) :: N, mu0
    N = size(Y,1)
    mu0 = N
    if(present(mu)) mu0 = mu
    call c_qsort(Y, idx, N, 1, 1, mu0)
  End Subroutine mlf_sortV
  Subroutine mlf_sortM(Y, idx, mu)
    real(c_double), intent(in), contiguous :: Y(:,:)
    integer, intent(in), optional :: mu
    integer(c_int), intent(out) :: idx(:)
    integer(c_int) :: ND, N, mu0
    N = size(Y,2); ND = size(Y,1)
    mu0 = N
    if(present(mu)) mu0 = mu
    call c_qsort(Y, idx, N, ND, ND, mu0)
  End Subroutine mlf_sortM

  Function mlf_sortOutV(Y, mu) result(Z)
    real(c_double), intent(in), contiguous :: Y(:)
    integer, intent(in), optional :: mu
    real(c_double), allocatable :: Z(:)
    integer(c_int), allocatable :: idx(:)
    integer(c_int) :: N, mu0
    N = size(Y,1)
    mu0 = N
    if(present(mu)) mu0 = mu
    allocate(idx(N))
    call c_qsort(Y, idx, N, 1, 1, mu0)
    Z = Y(idx)
  End Function mlf_sortOutV
  Function mlf_sortOutM(Y, mu) result(Z)
    real(c_double), intent(in), contiguous :: Y(:,:)
    integer, intent(in), optional :: mu
    real(c_double), allocatable :: Z(:,:)
    integer(c_int), allocatable :: idx(:)
    integer(c_int) :: ND, N, mu0
    N = size(Y,2); ND = size(Y,1)
    mu0 = N
    if(present(mu)) mu0 = mu
    allocate(idx(N))
    call c_qsort(Y, idx, N, ND, ND, mu0)
    Z = Y(:,idx)
  End Function mlf_sortOutM

  ! Useful functions for getting a mean value
  Pure Real(c_double) function mlf_MeanV(V)
    real(c_double), intent(in) :: V(:)
    mlf_MeanV = sum(V)/real(size(V),kind=8)
  end function mlf_MeanV

  Pure Real(c_double) Function mlf_MeanM(M)
    real(c_double), intent(in) :: M(:,:)
    mlf_MeanM = sum(M)/real(size(M), kind=8)
  End Function mlf_MeanM

  Pure Function mlf_MeanMV(M, dim) Result(V)
    real(c_double), allocatable :: V(:)
    real(c_double), intent(in) :: M(:,:)
    integer, intent(in) :: dim
    V = SUM(M, DIM=dim)/REAL(SIZE(M,dim), KIND=8)
  End Function mlf_MeanMV

  ! Count class apearance within an integer array
  Integer(c_int) Function mlf_countcl(id, cnt, l) result(N)
    integer(c_int), intent(in) :: id(:), l
    integer(c_int), intent(out) :: cnt(l:)
    integer :: u, i, k
    N = 0
    u = UBOUND(cnt,1)
    cnt = 0
    Do i=1,SIZE(id,1)
      k = id(i)
      If(k<l .OR. k>u) CYCLE ! Ignore elements outside the array
      cnt(k) = cnt(k) + 1
      N = N+1 ! Return the number of detected elements
    End Do
  End Function mlf_countcl

  ! Reverse class association array
  ! Use last appearance of the class
  ! Keep the value of cl if cl not present within id
  Integer(c_int) Function mlf_reversecl(id, cl, l) Result(N)
    integer(c_int), intent(in) :: id(:), l
    integer(c_int), intent(inout) :: cl(l:)
    integer :: u, i, k
    N = 0
    u = UBOUND(cl,1)
    Do i=1,SIZE(id,1)
      k = id(i)
      If(k<l .OR. k>u) CYCLE ! Ignore elements outside the array
      cl(k) = i
      N = N+1 ! Return the number of detected elements
    End Do
  End Function mlf_reversecl

  Function mlf_reverseid(id) Result(cl)
    integer(c_int), intent(in) :: id(:)
    integer(c_int), allocatable :: cl(:)
    integer :: info
    ALLOCATE(cl(SIZE(id)))
    info = mlf_reversecl(id,cl,1)
  End Function mlf_reverseid

  ! Utility function generating cumulative vector
  Pure Function mlf_cumulativefromvect(V) result(W)
    real(c_double), intent(in) :: V(:)
    real(c_double) :: W(size(V)), S
    integer :: i
    S=0
    Do i =1,size(V,1)
      S = S + V(i)
      W(i) = S
    End Do
  End Function mlf_cumulativefromvect

  ! Do a dichotomous search for a value
  Pure Integer(c_int) Function mlf_di_search(V, x) result(k)
    real(c_double), intent(in) :: V(:), x
    integer :: i, j
    i = 1; j = SIZE(V)
    Do While(j-i>1)
      k = (i+j)/2
      If(V(k) < x) Then
        i = k
      Else
        j = k
      Endif
    End Do
    If(x < V(j)) Then
      k = i
    Else
      k = j
    Endif
  End Function mlf_di_search
  ! Print matrix (useful for debug)
  Subroutine mlf_printmatrix(M)
    real(c_double), intent(in) :: M(:,:)
    integer :: i
    Do i=1,size(M,1)
      write(*, "(24f6.2)") M(i,:)
    End Do
  End Subroutine mlf_printmatrix
  ! Print matrix (useful for debug)
  Subroutine mlf_printmatrix2(M, P)
    real(c_double), intent(in) :: M(:,:), P(:,:)
    integer :: i, j
    Do i=1,size(M,1)
      write(*, "(24f6.2)") [([M(i,j), P(i,j)], j=1, size(M,2))]
    End Do
  End Subroutine mlf_printmatrix2

  Subroutine mlf_printmatrix_c(M, nL, nC) bind(C, name="mlf_printMat")
    integer(c_int), value :: nL, nC
    real(c_double), intent(in) :: M(nL, nC)
    call mlf_printmatrix(M)
  End Subroutine mlf_printmatrix_c
End Module mlf_utils

max_stars_repo_name: Etdescamps/MlFortran

max_issues_repo_name: Etdescamps/MlFortran

max_forks_repo_name: Etdescamps/MlFortran


text:
      program equiv14

C     Check offsets

      real x(100), y(100)
      common /foo/x
      equivalence (x(2), y), (y(50), z)
      common /foo/u

      do i = 2, 10
         y(i) = 0.
      enddo

      print *, i, j

      end

ext: f


text:
      SUBROUTINE PSGESVX( FACT, TRANS, N, NRHS, A, IA, JA, DESCA, AF,
     $                    IAF, JAF, DESCAF, IPIV, EQUED, R, C, B, IB,
     $                    JB, DESCB, X, IX, JX, DESCX, RCOND, FERR,
     $                    BERR, WORK, LWORK, IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     December 31, 1998
*
*     .. Scalar Arguments ..
      CHARACTER          EQUED, FACT, TRANS
      INTEGER            IA, IAF, IB, INFO, IX, JA, JAF, JB, JX, LIWORK,
     $                   LWORK, N, NRHS
      REAL               RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCAF( * ), DESCB( * ),
     $                   DESCX( * ), IPIV( * ), IWORK( * )
      REAL               A( * ), AF( * ), B( * ), BERR( * ), C( * ),
     $                   FERR( * ), R( * ), WORK( * ), X( * )
*     ..
*
*  Purpose
*  =======
*
*  PSGESVX uses the LU factorization to compute the solution to a real
*  system of linear equations
*
*        A(IA:IA+N-1,JA:JA+N-1) * X = B(IB:IB+N-1,JB:JB+NRHS-1),
*
*  where A(IA:IA+N-1,JA:JA+N-1) is an N-by-N matrix and X and
*  B(IB:IB+N-1,JB:JB+NRHS-1) are N-by-NRHS matrices.
*
*  Error bounds on the solution and a condition estimate are also
*  provided.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Description
*  ===========
*
*  In the following description, A denotes A(IA:IA+N-1,JA:JA+N-1),
*  B denotes B(IB:IB+N-1,JB:JB+NRHS-1) and X denotes
*  X(IX:IX+N-1,JX:JX+NRHS-1).
*
*  The following steps are performed:
*
*  1. If FACT = 'E', real scaling factors are computed to equilibrate
*     the system:
*        TRANS = 'N':  diag(R)*A*diag(C)     *inv(diag(C))*X = diag(R)*B
*        TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
*        TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
*     Whether or not the system will be equilibrated depends on the
*     scaling of the matrix A, but if equilibration is used, A is
*     overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
*     or diag(C)*B (if TRANS = 'T' or 'C').
*
*  2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
*     matrix A (after equilibration if FACT = 'E') as
*        A = P * L * U,
*     where P is a permutation matrix, L is a unit lower triangular
*     matrix, and U is upper triangular.
*
*  3. The factored form of A is used to estimate the condition number
*     of the matrix A.  If the reciprocal of the condition number is
*     less than machine precision, steps 4-6 are skipped.
*
*  4. The system of equations is solved for X using the factored form
*     of A.
*
*  5. Iterative refinement is applied to improve the computed solution
*     matrix and calculate error bounds and backward error estimates
*     for it.
*
*  6. If FACT = 'E' and equilibration was used, the matrix X is
*     premultiplied by diag(C) (if TRANS = 'N') or diag(R) (if
*     TRANS = 'T' or 'C') so that it solves the original system
*     before equilibration.
*
*  Arguments
*  =========
*
*  FACT    (global input) CHARACTER
*          Specifies whether or not the factored form of the matrix
*          A(IA:IA+N-1,JA:JA+N-1) is supplied on entry, and if not,
*          whether the matrix A(IA:IA+N-1,JA:JA+N-1) should be
*          equilibrated before it is factored.
*          = 'F':  On entry, AF(IAF:IAF+N-1,JAF:JAF+N-1) and IPIV con-
*                  tain the factored form of A(IA:IA+N-1,JA:JA+N-1).
*                  If EQUED is not 'N', the matrix
*                  A(IA:IA+N-1,JA:JA+N-1) has been equilibrated with
*                  scaling factors given by R and C.
*                  A(IA:IA+N-1,JA:JA+N-1), AF(IAF:IAF+N-1,JAF:JAF+N-1),
*                  and IPIV are not modified.
*          = 'N':  The matrix A(IA:IA+N-1,JA:JA+N-1) will be copied to
*                  AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
*          = 'E':  The matrix A(IA:IA+N-1,JA:JA+N-1) will be equili-
*                  brated if necessary, then copied to
*                  AF(IAF:IAF+N-1,JAF:JAF+N-1) and factored.
*
*  TRANS   (global input) CHARACTER
*          Specifies the form of the system of equations:
*          = 'N':  A(IA:IA+N-1,JA:JA+N-1) * X(IX:IX+N-1,JX:JX+NRHS-1)
*                = B(IB:IB+N-1,JB:JB+NRHS-1)     (No transpose)
*          = 'T':  A(IA:IA+N-1,JA:JA+N-1)**T * X(IX:IX+N-1,JX:JX+NRHS-1)
*                = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
*          = 'C':  A(IA:IA+N-1,JA:JA+N-1)**H * X(IX:IX+N-1,JX:JX+NRHS-1)
*                = B(IB:IB+N-1,JB:JB+NRHS-1)  (Transpose)
*
*  N       (global input) INTEGER
*          The number of rows and columns to be operated on, i.e. the
*          order of the distributed submatrix A(IA:IA+N-1,JA:JA+N-1).
*          N >= 0.
*
*  NRHS    (global input) INTEGER
*          The number of right-hand sides, i.e., the number of columns
*          of the distributed submatrices B(IB:IB+N-1,JB:JB+NRHS-1) and
*          X(IX:IX+N-1,JX:JX+NRHS-1).  NRHS >= 0.
*
*  A       (local input/local output) REAL pointer into
*          the local memory to an array of local dimension
*          (LLD_A,LOCc(JA+N-1)).  On entry, the N-by-N matrix
*          A(IA:IA+N-1,JA:JA+N-1).  If FACT = 'F' and EQUED is not 'N',
*          then A(IA:IA+N-1,JA:JA+N-1) must have been equilibrated by
*          the scaling factors in R and/or C.  A(IA:IA+N-1,JA:JA+N-1) is
*          not modified if FACT = 'F' or  'N', or if FACT = 'E' and
*          EQUED = 'N' on exit.
*
*          On exit, if EQUED .ne. 'N', A(IA:IA+N-1,JA:JA+N-1) is scaled
*          as follows:
*          EQUED = 'R':  A(IA:IA+N-1,JA:JA+N-1) :=
*                                      diag(R) * A(IA:IA+N-1,JA:JA+N-1)
*          EQUED = 'C':  A(IA:IA+N-1,JA:JA+N-1) :=
*                                      A(IA:IA+N-1,JA:JA+N-1) * diag(C)
*          EQUED = 'B':  A(IA:IA+N-1,JA:JA+N-1) :=
*                            diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  AF      (local input or local output) REAL pointer
*          into the local memory to an array of local dimension
*          (LLD_AF,LOCc(JA+N-1)).  If FACT = 'F', then
*          AF(IAF:IAF+N-1,JAF:JAF+N-1) is an input argument and on
*          entry contains the factors L and U from the factorization
*          A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by PSGETRF.
*          If EQUED .ne. 'N', then AF is the factored form of the
*          equilibrated matrix A(IA:IA+N-1,JA:JA+N-1).
*
*          If FACT = 'N', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
*          argument and on exit returns the factors L and U from the
*          factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original
*          matrix A(IA:IA+N-1,JA:JA+N-1).
*
*          If FACT = 'E', then AF(IAF:IAF+N-1,JAF:JAF+N-1) is an output
*          argument and on exit returns the factors L and U from the
*          factorization A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equili-
*          brated matrix A(IA:IA+N-1,JA:JA+N-1) (see the description of
*          A(IA:IA+N-1,JA:JA+N-1) for the form of the equilibrated
*          matrix).
*
*  IAF     (global input) INTEGER
*          The row index in the global array AF indicating the first
*          row of sub( AF ).
*
*  JAF     (global input) INTEGER
*          The column index in the global array AF indicating the
*          first column of sub( AF ).
*
*  DESCAF  (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix AF.
*
*  IPIV    (local input or local output) INTEGER array, dimension
*          LOCr(M_A)+MB_A. If FACT = 'F', then IPIV is an input argu-
*          ment and on entry contains the pivot indices from the fac-
*          torization A(IA:IA+N-1,JA:JA+N-1) = P*L*U as computed by
*          PSGETRF; IPIV(i) -> The global row local row i was
*          swapped with.  This array must be aligned with
*          A( IA:IA+N-1, * ).
*
*          If FACT = 'N', then IPIV is an output argument and on exit
*          contains the pivot indices from the factorization
*          A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the original matrix
*          A(IA:IA+N-1,JA:JA+N-1).
*
*          If FACT = 'E', then IPIV is an output argument and on exit
*          contains the pivot indices from the factorization
*          A(IA:IA+N-1,JA:JA+N-1) = P*L*U of the equilibrated matrix
*          A(IA:IA+N-1,JA:JA+N-1).
*
*  EQUED   (global input or global output) CHARACTER
*          Specifies the form of equilibration that was done.
*          = 'N':  No equilibration (always true if FACT = 'N').
*          = 'R':  Row equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1) has
*                  been premultiplied by diag(R).
*          = 'C':  Column equilibration, i.e., A(IA:IA+N-1,JA:JA+N-1)
*                  has been postmultiplied by diag(C).
*          = 'B':  Both row and column equilibration, i.e.,
*                  A(IA:IA+N-1,JA:JA+N-1) has been replaced by
*                  diag(R) * A(IA:IA+N-1,JA:JA+N-1) * diag(C).
*          EQUED is an input variable if FACT = 'F'; otherwise, it is an
*          output variable.
*
*  R       (local input or local output) REAL array,
*                                                  dimension LOCr(M_A).
*          The row scale factors for A(IA:IA+N-1,JA:JA+N-1).
*          If EQUED = 'R' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
*          on the left by diag(R); if EQUED='N' or 'C', R is not acces-
*          sed.  R is an input variable if FACT = 'F'; otherwise, R is
*          an output variable.
*          If FACT = 'F' and EQUED = 'R' or 'B', each element of R must
*          be positive.
*          R is replicated in every process column, and is aligned
*          with the distributed matrix A.
*
*  C       (local input or local output) REAL array,
*                                                  dimension LOCc(N_A).
*          The column scale factors for A(IA:IA+N-1,JA:JA+N-1).
*          If EQUED = 'C' or 'B', A(IA:IA+N-1,JA:JA+N-1) is multiplied
*          on the right by diag(C); if EQUED = 'N' or 'R', C is not
*          accessed.  C is an input variable if FACT = 'F'; otherwise,
*          C is an output variable.  If FACT = 'F' and EQUED = 'C' or
*          'B', each element of C must be positive.
*          C is replicated in every process row, and is aligned with
*          the distributed matrix A.
*
*  B       (local input/local output) REAL pointer
*          into the local memory to an array of local dimension
*          (LLD_B,LOCc(JB+NRHS-1) ).  On entry, the N-by-NRHS right-hand
*          side matrix B(IB:IB+N-1,JB:JB+NRHS-1). On exit, if
*          EQUED = 'N', B(IB:IB+N-1,JB:JB+NRHS-1) is not modified; if
*          TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
*          diag(R)*B(IB:IB+N-1,JB:JB+NRHS-1); if TRANS = 'T' or 'C'
*          and EQUED = 'C' or 'B', B(IB:IB+N-1,JB:JB+NRHS-1) is over-
*          written by diag(C)*B(IB:IB+N-1,JB:JB+NRHS-1).
*
*  IB      (global input) INTEGER
*          The row index in the global array B indicating the first
*          row of sub( B ).
*
*  JB      (global input) INTEGER
*          The column index in the global array B indicating the
*          first column of sub( B ).
*
*  DESCB   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix B.
*
*  X       (local input/local output) REAL pointer
*          into the local memory to an array of local dimension
*          (LLD_X, LOCc(JX+NRHS-1)).  If INFO = 0, the N-by-NRHS
*          solution matrix X(IX:IX+N-1,JX:JX+NRHS-1) to the original
*          system of equations.  Note that A(IA:IA+N-1,JA:JA+N-1) and
*          B(IB:IB+N-1,JB:JB+NRHS-1) are modified on exit if
*          EQUED .ne. 'N', and the solution to the equilibrated system
*          is inv(diag(C))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'N'
*          and EQUED = 'C' or 'B', or
*          inv(diag(R))*X(IX:IX+N-1,JX:JX+NRHS-1) if TRANS = 'T' or 'C'
*          and EQUED = 'R' or 'B'.
*
*  IX      (global input) INTEGER
*          The row index in the global array X indicating the first
*          row of sub( X ).
*
*  JX      (global input) INTEGER
*          The column index in the global array X indicating the
*          first column of sub( X ).
*
*  DESCX   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix X.
*
*  RCOND   (global output) REAL
*          The estimate of the reciprocal condition number of the matrix
*          A(IA:IA+N-1,JA:JA+N-1) after equilibration (if done).  If
*          RCOND is less than the machine precision (in particular, if
*          RCOND = 0), the matrix is singular to working precision.
*          This condition is indicated by a return code of INFO > 0.
*
*  FERR    (local output) REAL array, dimension LOCc(N_B)
*          The estimated forward error bounds for each solution vector
*          X(j) (the j-th column of the solution matrix
*          X(IX:IX+N-1,JX:JX+NRHS-1). If XTRUE is the true solution,
*          FERR(j) bounds the magnitude of the largest entry in
*          (X(j) - XTRUE) divided by the magnitude of the largest entry
*          in X(j).  The estimate is as reliable as the estimate for
*          RCOND, and is almost always a slight overestimate of the
*          true error.  FERR is replicated in every process row, and is
*          aligned with the matrices B and X.
*
*  BERR    (local output) REAL array, dimension LOCc(N_B).
*          The componentwise relative backward error of each solution
*          vector X(j) (i.e., the smallest relative change in
*          any entry of A(IA:IA+N-1,JA:JA+N-1) or
*          B(IB:IB+N-1,JB:JB+NRHS-1) that makes X(j) an exact solution).
*          BERR is replicated in every process row, and is aligned
*          with the matrices B and X.
*
*  WORK    (local workspace/local output) REAL array,
*                                                    dimension (LWORK)
*          On exit, WORK(1) returns the minimal and optimal LWORK.
*
*  LWORK   (local or global input) INTEGER
*          The dimension of the array WORK.
*          LWORK is local input and must be at least
*          LWORK = MAX( PSGECON( LWORK ), PSGERFS( LWORK ) )
*                  + LOCr( N_A ).
*
*          If LWORK = -1, then LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  IWORK   (local workspace/local output) INTEGER array,
*                                                  dimension (LIWORK)
*          On exit, IWORK(1) returns the minimal and optimal LIWORK.
*
*  LIWORK  (local or global input) INTEGER
*          The dimension of the array IWORK.
*          LIWORK is local input and must be at least
*          LIWORK = LOCr(N_A).
*
*          If LIWORK = -1, then LIWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, and i is
*                <= N:  U(IA+I-1,IA+I-1) is exactly zero.  The
*                       factorization has been completed, but the
*                       factor U is exactly singular, so the solution
*                       and error bounds could not be computed.
*                = N+1: RCOND is less than machine precision.  The
*                       factorization has been completed, but the
*                       matrix is singular to working precision, and
*                       the solution and error bounds have not been
*                       computed.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COLEQU, EQUIL, LQUERY, NOFACT, NOTRAN, ROWEQU
      CHARACTER          NORM
      INTEGER            CONWRK, I, IACOL, IAROW, IAFROW, IBROW, IBCOL,
     $                   ICOFFA, ICOFFB, ICOFFX, ICTXT, IDUMM,
     $                   IIA, IIB, IIX,
     $                   INFEQU, IROFFA, IROFFAF, IROFFB,
     $                   IROFFX, IXCOL, IXROW, J, JJA, JJB, JJX,
     $                   LCM, LCMQ,
     $                   LIWMIN, LWMIN, MYCOL, MYROW, NP, NPCOL, NPROW,
     $                   NQ, NQB, NRHSQ, RFSWRK
      REAL               AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
     $                   ROWCND, SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            CDESC( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DESCSET, PCHK2MAT,
     $                   INFOG2L, PSGECON, PSGEEQU, PSGERFS,
     $                   PSGETRF, PSGETRS, PSLACPY,
     $                   PSLAQGE, PSCOPY, PXERBLA, SGEBR2D,
     $                   SGEBS2D, SGAMN2D, SGAMX2D
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL, ILCM, INDXG2P, NUMROC
      REAL               PSLAMCH, PSLANGE
      EXTERNAL           ICEIL, ILCM, INDXG2P, LSAME, NUMROC, PSLANGE,
     $                   PSLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ICHAR, MAX, MIN, MOD, REAL
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -(800+CTXT_)
      ELSE
         CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 8, INFO )
         IF( LSAME( FACT, 'F' ) )
     $      CALL CHK1MAT( N, 3, N, 3, IAF, JAF, DESCAF, 12, INFO )
         CALL CHK1MAT( N, 3, NRHS, 4, IB, JB, DESCB, 20, INFO )
         CALL CHK1MAT( N, 3, NRHS, 4, IX, JX, DESCX, 24, INFO )
         NOFACT = LSAME( FACT, 'N' )
         EQUIL = LSAME( FACT, 'E' )
         NOTRAN = LSAME( TRANS, 'N' )
         IF( NOFACT .OR. EQUIL ) THEN
            EQUED = 'N'
            ROWEQU = .FALSE.
            COLEQU = .FALSE.
         ELSE
            ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
            COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
            SMLNUM = PSLAMCH( ICTXT, 'Safe minimum' )
            BIGNUM = ONE / SMLNUM
         END IF
         IF( INFO.EQ.0 ) THEN
            IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
     $                       NPROW )
            IAFROW = INDXG2P( IAF, DESCAF( MB_ ), MYROW,
     $                        DESCAF( RSRC_ ), NPROW )
            IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ),
     $                       NPROW )
            IXROW = INDXG2P( IX, DESCX( MB_ ), MYROW, DESCX( RSRC_ ),
     $                       NPROW )
            IROFFA = MOD( IA-1, DESCA( MB_ ) )
            IROFFAF = MOD( IAF-1, DESCAF( MB_ ) )
            ICOFFA = MOD( JA-1, DESCA( NB_ ) )
            IROFFB = MOD( IB-1, DESCB( MB_ ) )
            ICOFFB = MOD( JB-1, DESCB( NB_ ) )
            IROFFX = MOD( IX-1, DESCX( MB_ ) )
            ICOFFX = MOD( JX-1, DESCX( NB_ ) )
            CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW,
     $                    MYCOL, IIA, JJA, IAROW, IACOL )
            NP = NUMROC( N+IROFFA, DESCA( MB_ ), MYROW, IAROW,
     $                   NPROW )
            IF( MYROW.EQ.IAROW )
     $         NP = NP-IROFFA
            NQ = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL,
     $                   NPCOL )
            IF( MYCOL.EQ.IACOL )
     $         NQ = NQ-ICOFFA
            NQB = ICEIL( N+IROFFA, DESCA( NB_ )*NPCOL )
            LCM = ILCM( NPROW, NPCOL )
            LCMQ = LCM / NPCOL
            CONWRK = 2*NP + 2*NQ + MAX( 2, MAX( DESCA( NB_ )*
     $               MAX( 1, ICEIL( NPROW-1, NPCOL ) ), NQ +
     $               DESCA( NB_ )*
     $               MAX( 1, ICEIL( NPCOL-1, NPROW ) ) ) )
            RFSWRK = 3*NP
            IF( LSAME( TRANS, 'N' ) ) THEN
               RFSWRK = RFSWRK + NP + NQ +
     $                 ICEIL( NQB, LCMQ )*DESCA( NB_ )
            ELSE IF( LSAME( TRANS, 'T' ).OR.LSAME( TRANS, 'C' ) ) THEN
               RFSWRK = RFSWRK + NP + NQ
            END IF
            LWMIN = MAX( CONWRK, RFSWRK )
            WORK( 1 ) = REAL( LWMIN )
            LIWMIN = NP
            IWORK( 1 ) = LIWMIN
            IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND.
     $          .NOT.LSAME( FACT, 'F' ) ) THEN
               INFO = -1
            ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND.
     $         .NOT. LSAME( TRANS, 'C' ) ) THEN
               INFO = -2
            ELSE IF( IROFFA.NE.0 ) THEN
               INFO = -6
            ELSE IF( ICOFFA.NE.0 .OR. IROFFA.NE.ICOFFA ) THEN
               INFO = -7
            ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(800+NB_)
            ELSE IF( IAFROW.NE.IAROW ) THEN
               INFO = -10
            ELSE IF( IROFFAF.NE.0 ) THEN
               INFO = -10
            ELSE IF( ICTXT.NE.DESCAF( CTXT_ ) ) THEN
               INFO = -(1200+CTXT_)
            ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
     $         ( ROWEQU .OR. COLEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
               INFO = -13
            ELSE
               IF( ROWEQU ) THEN
                  RCMIN = BIGNUM
                  RCMAX = ZERO
                  DO 10 J = IIA, IIA + NP - 1
                     RCMIN = MIN( RCMIN, R( J ) )
                     RCMAX = MAX( RCMAX, R( J ) )
   10             CONTINUE
                  CALL SGAMN2D( ICTXT, 'Columnwise', ' ', 1, 1, RCMIN,
     $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
                  CALL SGAMX2D( ICTXT, 'Columnwise', ' ', 1, 1, RCMAX,
     $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
                  IF( RCMIN.LE.ZERO ) THEN
                     INFO = -14
                  ELSE IF( N.GT.0 ) THEN
                     ROWCND = MAX( RCMIN, SMLNUM ) /
     $                        MIN( RCMAX, BIGNUM )
                  ELSE
                     ROWCND = ONE
                  END IF
               END IF
               IF( COLEQU .AND. INFO.EQ.0 ) THEN
                  RCMIN = BIGNUM
                  RCMAX = ZERO
                  DO 20 J = JJA, JJA+NQ-1
                     RCMIN = MIN( RCMIN, C( J ) )
                     RCMAX = MAX( RCMAX, C( J ) )
   20             CONTINUE
                  CALL SGAMN2D( ICTXT, 'Rowwise', ' ', 1, 1, RCMIN,
     $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
                  CALL SGAMX2D( ICTXT, 'Rowwise', ' ', 1, 1, RCMAX,
     $                          1, IDUMM, IDUMM, -1, -1, MYCOL )
                  IF( RCMIN.LE.ZERO ) THEN
                     INFO = -15
                  ELSE IF( N.GT.0 ) THEN
                     COLCND = MAX( RCMIN, SMLNUM ) /
     $                        MIN( RCMAX, BIGNUM )
                  ELSE
                     COLCND = ONE
                  END IF
               END IF
            END IF
         END IF
*
         WORK( 1 ) = REAL( LWMIN )
         IWORK( 1 ) = LIWMIN
         LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
         IF( INFO.EQ.0 ) THEN
            IF( IBROW.NE.IAROW ) THEN
               INFO = -18
            ELSE IF( IXROW.NE.IBROW ) THEN
               INFO = -22
            ELSE IF( DESCB( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(2000+NB_)
            ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
               INFO = -(2000+CTXT_)
            ELSE IF( DESCX( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(2400+NB_)
            ELSE IF( ICTXT.NE.DESCX( CTXT_ ) ) THEN
               INFO = -(2400+CTXT_)
            ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
               INFO = -29
            ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
               INFO = -31
            END IF
            IDUM1( 1 ) = ICHAR( FACT )
            IDUM2( 1 ) = 1
            IDUM1( 2 ) = ICHAR( TRANS )
            IDUM2( 2 ) = 2
            IF( LSAME( FACT, 'F' ) ) THEN
               IDUM1( 3 ) = ICHAR( EQUED )
               IDUM2( 3 ) = 14
               IF( LWORK.EQ.-1 ) THEN
                  IDUM1( 4 ) = -1
               ELSE
                  IDUM1( 4 ) = 1
               END IF
               IDUM2( 4 ) = 29
               IF( LIWORK.EQ.-1 ) THEN
                  IDUM1( 5 ) = -1
               ELSE
                  IDUM1( 5 ) = 1
               END IF
               IDUM2( 5 ) = 31
               CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 8, N, 3,
     $                        NRHS, 4, IB, JB, DESCB, 20, 5, IDUM1,
     $                        IDUM2, INFO )
            ELSE
               IF( LWORK.EQ.-1 ) THEN
                  IDUM1( 3 ) = -1
               ELSE
                  IDUM1( 3 ) = 1
               END IF
               IDUM2( 3 ) = 29
               IF( LIWORK.EQ.-1 ) THEN
                  IDUM1( 4 ) = -1
               ELSE
                  IDUM1( 4 ) = 1
               END IF
               IDUM2( 4 ) = 31
               CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 8, N, 3,
     $                        NRHS, 4, IB, JB, DESCB, 20, 4, IDUM1,
     $                        IDUM2, INFO )
            END IF
         END IF
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PSGESVX', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
      IF( EQUIL ) THEN
*
*        Compute row and column scalings to equilibrate the matrix A.
*
         CALL PSGEEQU( N, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND,
     $                 AMAX, INFEQU )
         IF( INFEQU.EQ.0 ) THEN
*
*           Equilibrate the matrix.
*
            CALL PSLAQGE( N, N, A, IA, JA, DESCA, R, C, ROWCND, COLCND,
     $                    AMAX, EQUED )
            ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
            COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
         END IF
      END IF
*
*     Scale the right-hand side.
*
      CALL INFOG2L( IB, JB, DESCB, NPROW, NPCOL, MYROW, MYCOL, IIB,
     $              JJB, IBROW, IBCOL )
      NP = NUMROC( N+IROFFB, DESCB( MB_ ), MYROW, IBROW, NPROW )
      NRHSQ = NUMROC( NRHS+ICOFFB, DESCB( NB_ ), MYCOL, IBCOL, NPCOL )
      IF( MYROW.EQ.IBROW )
     $   NP = NP-IROFFB
      IF( MYCOL.EQ.IBCOL )
     $   NRHSQ = NRHSQ-ICOFFB
*
      IF( NOTRAN ) THEN
         IF( ROWEQU ) THEN
            DO 40 J = JJB, JJB+NRHSQ-1
               DO 30 I = IIB, IIB+NP-1
                  B( I+( J-1 )*DESCB( LLD_ ) ) = R( I )*
     $                                     B( I+( J-1 )*DESCB( LLD_ ) )
   30          CONTINUE
   40       CONTINUE
         END IF
      ELSE IF( COLEQU ) THEN
*
*        Transpose the Column scale factors
*
         CALL DESCSET( CDESC, 1, N+ICOFFA, 1, DESCA( NB_ ), MYROW,
     $                 IACOL, ICTXT, 1 )
         CALL PSCOPY( N, C, 1, JA, CDESC, CDESC( LLD_ ), WORK, IB, JB,
     $                DESCB, 1 )
         IF( MYCOL.EQ.IBCOL ) THEN
            CALL SGEBS2D( ICTXT, 'Rowwise', ' ', NP, 1, WORK( IIB ),
     $                    DESCB( LLD_ ) )
         ELSE
            CALL SGEBR2D( ICTXT, 'Rowwise', ' ', NP, 1, WORK( IIB ),
     $                    DESCB( LLD_ ), MYROW, IBCOL )
         END IF
         DO 60 J = JJB, JJB+NRHSQ-1
            DO 50 I = IIB, IIB+NP-1
               B( I+( J-1 )*DESCB( LLD_ ) ) = WORK( I )*
     $                                     B( I+( J-1 )*DESCB( LLD_ ) )
   50       CONTINUE
   60    CONTINUE
      END IF
*
      IF( NOFACT.OR.EQUIL ) THEN
*
*        Compute the LU factorization of A.
*
         CALL PSLACPY( 'Full', N, N, A, IA, JA, DESCA, AF, IAF, JAF,
     $                 DESCAF )
         CALL PSGETRF( N, N, AF, IAF, JAF, DESCAF, IPIV, INFO )
*
*        Return if INFO is non-zero.
*
         IF( INFO.NE.0 ) THEN
            IF( INFO.GT.0 )
     $         RCOND = ZERO
            RETURN
         END IF
      END IF
*
*     Compute the norm of the matrix A.
*
      IF( NOTRAN ) THEN
         NORM = '1'
      ELSE
         NORM = 'I'
      END IF
      ANORM = PSLANGE( NORM, N, N, A, IA, JA, DESCA, WORK )
*
*     Compute the reciprocal of the condition number of A.
*
      CALL PSGECON( NORM, N, AF, IAF, JAF, DESCAF, ANORM, RCOND, WORK,
     $              LWORK, IWORK, LIWORK, INFO )
*
*     Return if the matrix is singular to working precision.
*
      IF( RCOND.LT.PSLAMCH( ICTXT, 'Epsilon' ) ) THEN
         INFO = IA + N
         RETURN
      END IF
*
*     Compute the solution matrix X.
*
      CALL PSLACPY( 'Full', N, NRHS, B, IB, JB, DESCB, X, IX, JX,
     $              DESCX )
      CALL PSGETRS( TRANS, N, NRHS, AF, IAF, JAF, DESCAF, IPIV, X, IX,
     $              JX, DESCX, INFO )
*
*     Use iterative refinement to improve the computed solution and
*     compute error bounds and backward error estimates for it.
*
      CALL PSGERFS( TRANS, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF,
     $              DESCAF, IPIV, B, IB, JB, DESCB, X, IX, JX, DESCX,
     $              FERR, BERR, WORK, LWORK, IWORK, LIWORK, INFO )
*
*     Transform the solution matrix X to a solution of the original
*     system.
*
      CALL INFOG2L( IX, JX, DESCX, NPROW, NPCOL, MYROW, MYCOL, IIX,
     $              JJX, IXROW, IXCOL )
      NP = NUMROC( N+IROFFX, DESCX( MB_ ), MYROW, IXROW, NPROW )
      NRHSQ = NUMROC( NRHS+ICOFFX, DESCX( NB_ ), MYCOL, IXCOL, NPCOL )
      IF( MYROW.EQ.IBROW )
     $   NP = NP-IROFFX
      IF( MYCOL.EQ.IBCOL )
     $   NRHSQ = NRHSQ-ICOFFX
*
      IF( NOTRAN ) THEN
         IF( COLEQU ) THEN
*
*           Transpose the column scaling factors
*
            CALL DESCSET( CDESC, 1, N+ICOFFA, 1, DESCA( NB_ ), MYROW,
     $                    IACOL, ICTXT, 1 )
            CALL PSCOPY( N, C, 1, JA, CDESC, CDESC( LLD_ ), WORK, IX,
     $                   JX, DESCX, 1 )
            IF( MYCOL.EQ.IBCOL ) THEN
                CALL SGEBS2D( ICTXT, 'Rowwise', ' ', NP, 1,
     $                        WORK( IIX ), DESCX( LLD_ ) )
            ELSE
                CALL SGEBR2D( ICTXT, 'Rowwise', ' ', NP, 1,
     $                        WORK( IIX ), DESCX( LLD_ ), MYROW, IBCOL )
            END IF
*
            DO 80 J = JJX, JJX+NRHSQ-1
               DO 70 I = IIX, IIX+NP-1
                  X( I+( J-1 )*DESCX( LLD_ ) ) = WORK( I )*
     $                                      X( I+( J-1 )*DESCX( LLD_ ) )
   70          CONTINUE
   80       CONTINUE
            DO 90 J = JJX, JJX+NRHSQ-1
               FERR( J ) = FERR( J ) / COLCND
   90       CONTINUE
         END IF
      ELSE IF( ROWEQU ) THEN
         DO 110 J = JJX, JJX+NRHSQ-1
            DO 100 I = IIX, IIX+NP-1
               X( I+( J-1 )*DESCX( LLD_ ) ) = R( I )*
     $                                     X( I+( J-1 )*DESCX( LLD_ ) )
  100       CONTINUE
  110    CONTINUE
         DO 120 J = JJX, JJX+NRHSQ-1
            FERR( J ) = FERR( J ) / ROWCND
  120    CONTINUE
      END IF
*
      WORK( 1 ) = REAL( LWMIN )
      IWORK( 1 ) = LIWMIN
*
      RETURN
*
*     End of PSGESVX
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
      SUBROUTINE RDSINR (XNAME,XMIN,XMAX,X)
      IMPLICIT NONE

*     FORMAL_PARAMETERS:
      INTEGER X, XMIN, XMAX
      CHARACTER*(*) XNAME

**    local variables ; dummy set
      DOUBLE PRECISION D(1),DM
      REAL             R(1),RM
      INTEGER          I(1),IM
      CHARACTER*80     C(1),CM
      LOGICAL          L(1),LM

*     include files
      INCLUDE 'rdstainf.inc'

*     other
      INTEGER IL,IS,ILX
      CHARACTER Q*1,LXNAME*31
      SAVE

      DATA Q /''''/

*     read a single INTEGER number (value of IL is arbitrary) !!
      IL = 0
      CALL RDDATA (5,'RDSINR',0,0,' ',IS,XNAME,'I',D,R,I,C,L,
     $             1,1,1,1,1,IL,
     $             DM,RM,IM,CM,LM)

      X      = I(1)
      LXNAME = XNAME
      ILX    = LEN_TRIM (LXNAME)
      IF (XMAX.LT.XMIN) THEN
         IF (TOSCR) WRITE (*,'(1X,4A)')
     &     'ERROR in RDSINR: Range not valid of identifier ',
     &      Q,LXNAME(1:ILX),Q
         IF (TOLOG) WRITE (IULOG,'(1X,4A)')
     &     'ERROR in RDSINR: Range not valid of identifier ',
     &      Q,LXNAME(1:ILX),Q
         CALL FATALERR (' ',' ')
      ELSE IF ((X.LT.XMIN.OR.X.GT.XMAX).AND.X.NE.IM) THEN
         IF (TOSCR) WRITE (*,'(1X,5A,/,T19,A,I8,A,I8,A,I8,A)')
     &     'ERROR in RDSINR: Range error of identifier ',
     &      Q,LXNAME(1:ILX),Q,',',
     &     'value =',X,', range = [',XMIN,',',XMAX,']'
         IF (TOLOG) WRITE (IULOG,'(1X,5A,/,T19,A,I8,A,I8,A,I8,A)')
     &     'ERROR in RDSINR: Range error of identifier ',
     &      Q,LXNAME(1:ILX),Q,',',
     &     'value =',X,', range = [',XMIN,',',XMAX,']'
         CALL FATALERR (' ',' ')
      END IF

      RETURN
      END

ext: for

max_issues_repo_name: zhumiansanren/FIRI-WHS


text:
      subroutine LgetGL(id,idxx,ptr,n,itype)
      use maccbk
      use maccode
      implicit none

      character*(*) id,idw*(MAXPNAME)
      integer*8 ptr
      integer*4 idxx,n,itype
      integer*4 hsrchz
c
      idxx=hsrchz(id)
      if(idtype(idxx) .eq. icGLL) then
        ptr=idval(idxx)
        n= ilist(1,ptr)-1
        itype=ilist(2,ptr)
      else
        idw=id
        call errmsg('getglobal','type miss match . '//idw//'!',0,0)
        ptr=0
        n= 0
        itype=0
      endif
      return
c
      entry LsetGL(id,idxx,ptr,n,itype)
c
      idxx=hsrchz(id)
      if(idtype(idxx) .eq. icGLL) then
        idval(idxx)=ptr
        ilist(1,ptr)=n+1
        ilist(2,ptr)=itype
      else
        idw=id
        call errmsg('setglobal','type miss match . '//idw//'!',0,0)
      endif
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


    Program zgesdd_example

!     ZGESDD Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgesdd
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6, prerr = 0
!     .. Local Scalars ..
      Integer :: i, info, lda, ldu, ldvt, lwork, m, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), a_copy(:, :), b(:), u(:, :), &
        vt(:, :), work(:)
      Complex (Kind=dp) :: dummy(1, 1)
      Real (Kind=dp), Allocatable :: rwork(:), s(:)
      Integer, Allocatable :: iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: max, min, nint, real
!     .. Executable Statements ..
      Write (nout, *) 'ZGESDD Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      ldu = m
      ldvt = n
      Allocate (a(lda,n), a_copy(m,n), s(m), u(ldu,m), vt(ldvt,n), b(m), &
        rwork((5*m+7)*n), iwork(8*m))

!     Read the m by n matrix A from data file
      Read (nin, *)(a(i,1:n), i=1, m)

!     Read the right hand side of the linear system
      Read (nin, *) b(1:m)

      a_copy(1:m, 1:n) = a(1:m, 1:n)

!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call zgesdd('A', m, n, a, lda, s, u, ldu, vt, ldvt, dummy, lwork, rwork, &
        iwork, info)

!     Make sure that there is enough workspace for block size nb.
      lwork = max((2*m+2)*m+2*n+nb*(m+n), nint(real(dummy(1,1))))
      Allocate (work(lwork))

!     Compute the singular values and left and right singular vectors
!     of A.

      Call zgesdd('A', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, &
        iwork, info)

      If (info/=0) Then
        Write (nout, 100) 'Failure in ZGESDD. INFO =', info
100     Format (1X, A, I4)
        Go To 120
      End If

!     Print the significant singular values of A

      Write (nout, *) 'Singular values of A:'
      Write (nout, 110) s(1:min(m,n))
110   Format (1X, 4(3X,F11.4))

      If (prerr>0) Then
        Call compute_error_bounds(m, n, s)
      End If

      If (m<n) Then
!       Compute V*Inv(S)*U^T * b to get minimum norm solution.
        Call compute_minimum_norm(m, n, a_copy, m, u, ldu, vt, ldvt, s, b)
      End If

120   Continue

    Contains
      Subroutine compute_minimum_norm(m, n, a, lda, u, ldu, vt, ldvt, s, b)

!       .. Use Statements ..
        Use blas_interfaces, Only: dznrm2, zgemv
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In) :: lda, ldu, ldvt, m, n
!       .. Array Arguments ..
        Complex (Kind=dp), Intent (In) :: a(lda, n), u(ldu, m), vt(ldvt, n)
        Complex (Kind=dp), Intent (Inout) :: b(m)
        Real (Kind=dp), Intent (In) :: s(m)
!       .. Local Scalars ..
        Complex (Kind=dp) :: alpha, beta
        Real (Kind=dp) :: norm
!       .. Local Arrays ..
        Complex (Kind=dp), Allocatable :: x(:), y(:)
!       .. Intrinsic Procedures ..
        Intrinsic :: allocated, cmplx
!       .. Executable Statements ..
        Allocate (x(n), y(m))

!       Compute V*Inv(S)*U^H * b to get least squares solution.

!       y = U^H b
        alpha = cmplx(1.0_dp, 0.0_dp, kind=dp)
        beta = cmplx(0.0_dp, 0.0_dp, kind=dp)
        Call zgemv('C', m, m, alpha, u, ldu, b, 1, beta, y, 1)

        y(1:m) = y(1:m)/s(1:m)

!       x = V y
        Call zgemv('C', m, n, alpha, vt, ldvt, y, 1, beta, x, 1)

        Write (nout, *)
        Write (nout, *) 'Minimum norm solution:'
        Write (nout, 100) x(1:n)

        norm = dznrm2(n, x, 1)

        Write (nout, *)
        Write (nout, *) 'Norm of Solution:'
        Write (nout, 110) norm

!       Find norm of residual ||b-Ax||, should be zero.
        alpha = cmplx(-1.0_dp, 0.0_dp, kind=dp)
        beta = cmplx(1._dp, 0.0_dp, kind=dp)
        Call zgemv('N', m, n, alpha, a, lda, x, 1, beta, b, 1)

        norm = dznrm2(m, b, 1)

        Write (nout, *)
        Write (nout, *) 'Norm of Residual:'
        Write (nout, 110) norm

        If (allocated(x)) Then
          Deallocate (x)
        End If
        If (allocated(y)) Then
          Deallocate (y)
        End If

100     Format (4X, '(', F8.4, ',', F8.4, ')')
110     Format (4X, F11.4)

      End Subroutine

      Subroutine compute_error_bounds(m, n, s)

!       Error estimates for singular values and vectors is computed
!       and printed here.

!       .. Use Statements ..
        Use lapack_interfaces, Only: ddisna
        Use lapack_precision, Only: dp
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In) :: m, n
!       .. Array Arguments ..
        Real (Kind=dp), Intent (In) :: s(m)
!       .. Local Scalars ..
        Real (Kind=dp) :: eps, serrbd
        Integer :: i, info
!       .. Local Arrays ..
        Real (Kind=dp), Allocatable :: rcondu(:), rcondv(:), uerrbd(:), &
          verrbd(:)
!       .. Intrinsic Procedures ..
        Intrinsic :: epsilon
!       .. Executable Statements ..
        Allocate (rcondu(n), rcondv(n), uerrbd(n), verrbd(n))

!       Get the machine precision, EPS and compute the approximate
!       error bound for the computed singular values.  Note that for
!       the 2-norm, S(1) = norm(A)

        eps = epsilon(1.0E0_dp)
        serrbd = eps*s(1)

!       Call DDISNA to estimate reciprocal condition
!       numbers for the singular vectors

        Call ddisna('Left', m, n, s, rcondu, info)
        Call ddisna('Right', m, n, s, rcondv, info)

!       Compute the error estimates for the singular vectors

        Do i = 1, n
          uerrbd(i) = serrbd/rcondu(i)
          verrbd(i) = serrbd/rcondv(i)
        End Do

!       Print the approximate error bounds for the singular values
!       and vectors

        Write (nout, *)
        Write (nout, *) 'Error estimate for the singular values'
        Write (nout, 100) serrbd
        Write (nout, *)
        Write (nout, *) 'Error estimates for the left singular vectors'
        Write (nout, 100) uerrbd(1:n)
        Write (nout, *)
        Write (nout, *) 'Error estimates for the right singular vectors'
        Write (nout, 100) verrbd(1:n)

100     Format (4X, 1P, 6E11.1)

      End Subroutine

    End Program


text:
	SUBROUTINE JTB_TBLB ( np, parm, iscl, iref, ibit, punit, iret )
C************************************************************************
C* JTB_TBLB								*
C*									*
C* This subroutine scans information dumped by SUBROUTINE GETABDB to	*
C* get a list of the Table B entries.					*
C*									*
C* If NP = 0, there are no table B entries.				*
C*									*
C* JTB_TBLB ( NP, PARM, ISCL, IREF, IBIT, PUNIT, IRET )			*
C*									*
C* Input and output parameters						*
C*	NP		INTEGER		Max number of parm names	*
C*					OUTPUT:  actual # of parm names	*
C*									*
C* Output parameters:							*
C*	PARM   (NP)	CHAR*		Array of parameter names	*
C*	ISCL    (NP)	INTEGER		Power of 10 scaling		*
C*	IREF 	(NP)	INTEGER		Scaled reference value		*
C*	IBIT 	(NP)	INTEGER		Number of bits			*
C*	PUNIT	(NP)    CHAR*(*)	Parameter array			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -1 = Missing packing info	*
C*					 -2 = Too many Table B entries	*
C**									*
C* Log:									*
C* K. Brill/EMC		 9/98						*
C************************************************************************
	INCLUDE		'jtb.cmn'
C*
	CHARACTER*(*)	parm (*), punit (*)
	INTEGER		iscl (*), iref (*), ibit (*)
C*
	INTEGER		ipkp (3)
C-----------------------------------------------------------------------
	iret = 0
C
C*	Scan for Table B parameters.
C
	ip = 0
	ie = 0
	DO WHILE ( ie .lt. nntrys )
	    ie = ie + 1
	    IF ( entrys (ie) (1:1) .eq. 'B' ) THEN
		ip = ip + 1
		IF ( ip .gt. np ) THEN
		    iret = -2
		    RETURN
		END IF
		parm (ip) = entrys (ie) (3:10)
		punit (ip) = entrys (ie) (12:33)
		CALL ST_ILST ( entrys (ie) (35: ), ' ', -9999, 3,
     +			       ipkp, num, ier )
		IF ( ier .ne. 0 .or. num .ne. 3 ) THEN
		    iret = -1
		    RETURN
		END IF
		iscl (ip) = ipkp (1)
		iref (ip) = ipkp (2)
		ibit (ip) = ipkp (3)
	    END IF
	END DO
	np = ip
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine array

      integer a(3)
      integer equiv
      integer unknown

      equivalence (a(3), equiv)

c     check points
      integer a1, a2, a3 , a4 , a5, a6, a7, a8, a9, a10, a11

      a(1) = 1
      a(2) = 2
      if (a(3) .gt. 100) then
          a(2) =  2
          unknown = 1
      else
          a(2) = -2
          unknown = 3
      endif
      a1 = a(1)         ! Const (Int 1)
      a2 = a(2)         ! Bot
      a3 = a(3)         ! Top

c     equivalence check
      a(3) = 123
      a4 = equiv        ! Const(Int 123)

c     unknown index on rhs
      a5 = a(unknown)   ! Bot
c     unknown index on lhs
      a(unknown) = 100  ! this will set whole array to Bot
      a6 = a(1)         ! Bot
      a7 = a(2)         ! Bot
      a8 = a(3)         ! Bot
      a9 = equiv        ! Bot

c     update array after processing unknown index on the lhs

      a(1) = 456
      a10  = a(1)       ! Const (Int 456))
      a11 =  a(2)       ! Bot

      end

max_stars_repo_name: camfort/fortran-vars

max_issues_repo_name: camfort/fortran-vars

max_forks_repo_name: camfort/fortran-vars


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!- mod_namelist - (January - 2013)
!!
!!- Author: Nicolas.Grima@univ-brest.fr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!!
!!
!!
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
MODULE mod_namelist

  !===================!
  !- USE association -!
  !===================!
  USE mod_precision
  USE mod_cst
  USE mod_lun

  !===============================!
  !- DECLARATIONS of GLOBAL DATA -!
  !===============================!
  IMPLICIT NONE

  !----------!
  ! FSP item !
  !----------!
  !- FSP item pramaters -!
  INTEGER(kind = iprec)       :: ifsp = i_defn
  REAL   (kind = rprec)       :: rfsp = r_def
  LOGICAL                     :: lfsp = l_def
  CHARACTER(len = iSixtyFour) :: cfsp = c_def

  NAMELIST/FSP/ & !
       ifsp   , & ! an example of an  integer namelist parameter
       rfsp   , & ! an example of a      real namelist parameter
       lfsp   , & ! an example of a   logical namelist parameter
       cfsp       ! an example of a character namelist parameter

  !------------!
  ! PARAM item !
  !------------!
  !- PARAM item pramaters -!
  INTEGER(kind = iprec)       :: iparam = i_defn
  REAL   (kind = rprec)       :: rparam = r_def
  LOGICAL                     :: lparam = l_def
  CHARACTER(len = iSixtyFour) :: cparam = c_def

  NAMELIST/PARAM/ & !
       iparam   , & ! an example of an  integer namelist parameter
       rparam   , & ! an example of a      real namelist parameter
       lparam   , & ! an example of a   logical namelist parameter
       cparam       ! an example of a character namelist parameter

  !====================================!
  !- DECLARATIONS of INTERFACE BLOCKS -!
  !====================================!

CONTAINS

  !********************************************************************!
  !- sub_read_namelist - (January - 2013)
  !
  !- contact: Nicolas.Grima@univ-best.fr
  !********************************************************************!
  !
  !
  !
  !
  !
  !********************************************************************!

  SUBROUTINE sub_read_namelist( &
       verbose                 ) ! optional argument to select verbose mode or not

    !=============================!
    !- DECLARATIONS of ARGUMENTS -!
    !=============================!
    LOGICAL, OPTIONAL, INTENT(in) :: verbose

    !===================================!
    !- DECLARATIONS of LOCAL VARIABLES -!
    !===================================!
    LOGICAL :: print_info

    !=======================!
    != .oO) MAIN PART (Oo. =!
    !=======================!
    !!write(lun_dbg,*)' -db- Enter in sub_read_namelist'

    !- Test if the optional argument is present -!
    IF (PRESENT(verbose)) THEN
       IF (verbose) THEN
          print_info = .TRUE.
       ELSE
          print_info = .FALSE.
       ENDIF
    ELSE
       print_info = .FALSE.
    ENDIF

    !- print message in lun_standard -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)'============'
    WRITE(lun_standard,*)'= NAMELIST ='
    WRITE(lun_standard,*)'============'

    !-----------------------------!
    !- Opening the namelist file -!
    !-----------------------------!
    OPEN(unit=lun_nml, File='namelist', ACTION='READ')
    WRITE(lun_standard,*)' --- Successful Opening ---'

    !----------------------------------!
    !- Reading namelist file item FSP -!
    !----------------------------------!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' - Reading FSP item:'
    REWIND(unit=lun_nml)
    READ(unit=lun_nml, nml=FSP)

    !-------------------------------!
    !- Printing namelist parameter -!
    !-------------------------------!
    IF (print_info) THEN
       WRITE(lun_standard,*)'   - ifsp   =', ifsp
       WRITE(lun_standard,*)'   - rfsp   =', rfsp
       WRITE(lun_standard,*)'   - lfsp   =', lfsp
       WRITE(lun_standard,*)'   - cfsp   =', TRIM(cfsp)
    ENDIF

    !------------------------------------!
    !- Reading namelist file item PARAM -!
    !------------------------------------!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' - Reading PARAM item:'
    REWIND(unit=lun_nml)
    READ(unit=lun_nml, nml=PARAM)

    IF (print_info) THEN
       WRITE(lun_standard,*)'   - iparam =', iparam
       WRITE(lun_standard,*)'   - rparam =', rparam
       WRITE(lun_standard,*)'   - lparam =', lparam
       WRITE(lun_standard,*)'   - cparam =', TRIM(cparam)
    ENDIF

    !!write(lun_dbg,*)' -db- Exit sub_read_namelist'

  END SUBROUTINE sub_read_namelist

END MODULE mod_namelist

max_stars_repo_name: ngrigri/fsp

max_issues_repo_name: ngrigri/fsp

max_forks_repo_name: ngrigri/fsp


text:
*
*     -----------------------------------------------------------------
*      R M E W D 1 L S
*     -----------------------------------------------------------------
*
*     Written by G. Gaigalas,                                      *
*     Vilnius,  Lithuania                              June 1995   *
*
      SUBROUTINE RMEWD1LS(J,L,COEF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ID120N(136),ID120L(136)
      DATA ID120N/5*0,24696,3*0,-1008,3*0,-22050,18375,0,6720,0,630,
     :-23625,675,2*0,-8820,7*0,-15750,0,-300,4*0,-28224,3*0,26880,0,
     :-7056,3*0,-4032,14700,-34020,3*0,3528,3*0,-10080,0,-4200,3*0,
     :8820,5*0,31500,-900,0,14400,0,10080,4032,2376,8*0,-37800,9*0,
     :11880,5*0,-14580,11*0,-27720,-11088,-5544,0,15120,36036,11*0,
     :-26208,2*0,-3822,0/
      DATA ID120L/5*0,-226380,2*0,-25872,166782,4*0,-404250,2*0,
     :258720,-73920,0,4950,-161700,41580,4*0,-16500,0,36960,4*0,
     :118800,-14850,2*0,258720,223608,0,-110880,2*0,64680,2*0,-103488,
     :29568,0,110880,2*0,-316932,195657,7*0,92400,2*0,121275,5*0,
     :237600,2*0,374220,-10395,0,-9801,6*0,52800,66000,2*0,-226380,0,
     :363000,6*0,-232320,29040,2*0,-12705,0,-3300,-46875,8*0,-203280,
     :203280,0,-365904,2*0,66066,12*0,120120,-174720,0,95550/
      IF(L.EQ.0) THEN
        ID=ID120L(J)
        IV=2695
      ELSE
        ID=ID120N(J)
        IV=245
      ENDIF
      IF(ID.GE.0) THEN
        COEF=DSQRT(DBLE(ID))
      ELSE
        COEF=-DSQRT(-DBLE(ID))
      ENDIF
      COEF=COEF/DSQRT(DBLE(IV))
      RETURN
      END