text:
*> \brief \b CGETSQRHRT
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CGETSQRHRT + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
*      $                       LWORK, INFO )
*       IMPLICIT NONE
*
*       .. Scalar Arguments ..
*       INTEGER           INFO, LDA, LDT, LWORK, M, N, NB1, NB2, MB1
*       ..
*       .. Array Arguments ..
*       COMPLEX*16        A( LDA, * ), T( LDT, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CGETSQRHRT computes a NB2-sized column blocked QR-factorization
*> of a complex M-by-N matrix A with M >= N,
*>
*>    A = Q * R.
*>
*> The routine uses internally a NB1-sized column blocked and MB1-sized
*> row blocked TSQR-factorization and perfors the reconstruction
*> of the Householder vectors from the TSQR output. The routine also
*> converts the R_tsqr factor from the TSQR-factorization output into
*> the R factor that corresponds to the Householder QR-factorization,
*>
*>    A = Q_tsqr * R_tsqr = Q * R.
*>
*> The output Q and R factors are stored in the same format as in CGEQRT
*> (Q is in blocked compact WY-representation). See the documentation
*> of CGEQRT for more details on the format.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrix A. M >= N >= 0.
*> \endverbatim
*>
*> \param[in] MB1
*> \verbatim
*>          MB1 is INTEGER
*>          The row block size to be used in the blocked TSQR.
*>          MB1 > N.
*> \endverbatim
*>
*> \param[in] NB1
*> \verbatim
*>          NB1 is INTEGER
*>          The column block size to be used in the blocked TSQR.
*>          N >= NB1 >= 1.
*> \endverbatim
*>
*> \param[in] NB2
*> \verbatim
*>          NB2 is INTEGER
*>          The block size to be used in the blocked QR that is
*>          output. NB2 >= 1.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA,N)
*>
*>          On entry: an M-by-N matrix A.
*>
*>          On exit:
*>           a) the elements on and above the diagonal
*>              of the array contain the N-by-N upper-triangular
*>              matrix R corresponding to the Householder QR;
*>           b) the elements below the diagonal represent Q by
*>              the columns of blocked V (compact WY-representation).
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] T
*> \verbatim
*>          T is COMPLEX array, dimension (LDT,N))
*>          The upper triangular block reflectors stored in compact form
*>          as a sequence of upper triangular blocks.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*>          LDT is INTEGER
*>          The leading dimension of the array T.  LDT >= NB2.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          The dimension of the array WORK.
*>          LWORK >= MAX( LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) ),
*>          where
*>             NUM_ALL_ROW_BLOCKS = CEIL((M-N)/(MB1-N)),
*>             NB1LOCAL = MIN(NB1,N).
*>             LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL,
*>             LW1 = NB1LOCAL * N,
*>             LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) ),
*>          If LWORK = -1, then a workspace query is assumed.
*>          The routine only calculates the optimal size of the WORK
*>          array, returns this value as the first entry of the WORK
*>          array, and no error message related to LWORK is issued
*>          by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2020
*
*> \ingroup comlpexOTHERcomputational
*
*> \par Contributors:
*  ==================
*>
*> \verbatim
*>
*> November 2020, Igor Kozachenko,
*>                Computer Science Division,
*>                University of California, Berkeley
*>
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
     $                       LWORK, INFO )
      IMPLICIT NONE
*
*  -- LAPACK computational routine (version 3.10.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2020
*
*     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDT, LWORK, M, N, NB1, NB2, MB1
*     ..
*     .. Array Arguments ..
      COMPLEX           A( LDA, * ), T( LDT, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            CONE
      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            I, IINFO, J, LW1, LW2, LWT, LDWT, LWORKOPT,
     $                   NB1LOCAL, NB2LOCAL, NUM_ALL_ROW_BLOCKS
*     ..
*     .. External Subroutines ..
      EXTERNAL           CCOPY, CLATSQR, CUNGTSQR_ROW, CUNHR_COL,
     $                   XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          CEILING, REAL, CMPLX, MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      INFO = 0
      LQUERY  = LWORK.EQ.-1
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 .OR. M.LT.N ) THEN
         INFO = -2
      ELSE IF( MB1.LE.N ) THEN
         INFO = -3
      ELSE IF( NB1.LT.1 ) THEN
         INFO = -4
      ELSE IF( NB2.LT.1 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -7
      ELSE IF( LDT.LT.MAX( 1,  MIN( NB2, N ) ) ) THEN
         INFO = -9
      ELSE
*
*        Test the input LWORK for the dimension of the array WORK.
*        This workspace is used to store array:
*        a) Matrix T and WORK for CLATSQR;
*        b) N-by-N upper-triangular factor R_tsqr;
*        c) Matrix T and array WORK for CUNGTSQR_ROW;
*        d) Diagonal D for CUNHR_COL.
*
         IF( LWORK.LT.N*N+1 .AND. .NOT.LQUERY ) THEN
            INFO = -11
         ELSE
*
*           Set block size for column blocks
*
            NB1LOCAL = MIN( NB1, N )
*
            NUM_ALL_ROW_BLOCKS = MAX( 1,
     $                   CEILING( REAL( M - N ) / REAL( MB1 - N ) ) )
*
*           Length and leading dimension of WORK array to place
*           T array in TSQR.
*
            LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL

            LDWT = NB1LOCAL
*
*           Length of TSQR work array
*
            LW1 = NB1LOCAL * N
*
*           Length of CUNGTSQR_ROW work array.
*
            LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) )
*
            LWORKOPT = MAX( LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) )
*
            IF( ( LWORK.LT.MAX( 1, LWORKOPT ) ).AND.(.NOT.LQUERY) ) THEN
               INFO = -11
            END IF
*
         END IF
      END IF
*
*     Handle error in the input parameters and return workspace query.
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGETSQRHRT', -INFO )
         RETURN
      ELSE IF ( LQUERY ) THEN
         WORK( 1 ) = CMPLX( LWORKOPT )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( MIN( M, N ).EQ.0 ) THEN
         WORK( 1 ) = CMPLX( LWORKOPT )
         RETURN
      END IF
*
      NB2LOCAL = MIN( NB2, N )
*
*
*     (1) Perform TSQR-factorization of the M-by-N matrix A.
*
      CALL CLATSQR( M, N, MB1, NB1LOCAL, A, LDA, WORK, LDWT,
     $              WORK(LWT+1), LW1, IINFO )
*
*     (2) Copy the factor R_tsqr stored in the upper-triangular part
*         of A into the square matrix in the work array
*         WORK(LWT+1:LWT+N*N) column-by-column.
*
      DO J = 1, N
         CALL CCOPY( J, A( 1, J ), 1, WORK( LWT + N*(J-1)+1 ), 1 )
      END DO
*
*     (3) Generate a M-by-N matrix Q with orthonormal columns from
*     the result stored below the diagonal in the array A in place.
*

      CALL CUNGTSQR_ROW( M, N, MB1, NB1LOCAL, A, LDA, WORK, LDWT,
     $                   WORK( LWT+N*N+1 ), LW2, IINFO )
*
*     (4) Perform the reconstruction of Householder vectors from
*     the matrix Q (stored in A) in place.
*
      CALL CUNHR_COL( M, N, NB2LOCAL, A, LDA, T, LDT,
     $                WORK( LWT+N*N+1 ), IINFO )
*
*     (5) Copy the factor R_tsqr stored in the square matrix in the
*     work array WORK(LWT+1:LWT+N*N) into the upper-triangular
*     part of A.
*
*     (6) Compute from R_tsqr the factor R_hr corresponding to
*     the reconstructed Householder vectors, i.e. R_hr = S * R_tsqr.
*     This multiplication by the sign matrix S on the left means
*     changing the sign of I-th row of the matrix R_tsqr according
*     to sign of the I-th diagonal element DIAG(I) of the matrix S.
*     DIAG is stored in WORK( LWT+N*N+1 ) from the CUNHR_COL output.
*
*     (5) and (6) can be combined in a single loop, so the rows in A
*     are accessed only once.
*
      DO I = 1, N
         IF( WORK( LWT+N*N+I ).EQ.-CONE ) THEN
            DO J = I, N
               A( I, J ) = -CONE * WORK( LWT+N*(J-1)+I )
            END DO
         ELSE
            CALL CCOPY( N-I+1, WORK(LWT+N*(I-1)+I), N, A( I, I ), LDA )
         END IF
      END DO
*
      WORK( 1 ) = CMPLX( LWORKOPT )
      RETURN
*
*     End of CGETSQRHRT
*
      END


!
!     ------------------------------------------------------------------
!       S E T U P m
!     ------------------------------------------------------------------
!
      SUBROUTINE SETUPM(ISH, J1, J2, JA, JB, NA, NB)
      use medefn_C
      use ndims_C
      use non30_C
      use red_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  19:44:38  11/20/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: ISH
      INTEGER , INTENT(IN) :: J1
      INTEGER , INTENT(IN) :: J2
      INTEGER , INTENT(IN) :: JA
      INTEGER , INTENT(IN) :: JB
      INTEGER , INTENT(IN) :: NA
      INTEGER , INTENT(IN) :: NB
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I1, I2HSH, I, IC, JC, NCC, K, JD, I2
!-----------------------------------------------
!
!      COMMON/MEDEFN/IHSH,NJ(16),LJ(16),NOSH(16,2),J1QN(31,3,2),IJFUL(16)
!      POINTER(QNOC,NOCCSH(1)),(QNELCSH,NELCSH(8,1)),
!     :       (QNOCORB,NOCORB(8,1)),(QJ1,J1QNRD(15,1))
!      POINTER(QIAJCMP,IAJCMP(1)),(QLJCOMP,LJCOMP(1)),
!     :       (QNJCOMP,NJCOMP(1)),(QIAJCLD,IAJCLD(1)),
!     :       (QLJCLSD,LJCLSD(1))
!      COMMON /NDIMS/ QNOC,QNELCSH,QNOCORB,QJ1,NCFG
!      COMMON /NON30/ QIAJCMP,QNJCOMP,QLJCOMP,QIAJCLD,QLJCLSD,MAXORB
!      COMMON/RED/INDL(20),INDR(20),IORST(20),NCI1,KPL,KPR,NC,LIS(16),
!     : JIST,JFST,NCONTR
!
!     Inserts the current subshell i1 into the left and right
!     coupling tree at position ish.
!     j1, j2 : location in the configurations of the shell
!     na, nb : occupied or not: =1 not occupied, =2 occupied
!
      IF (NA == 2) THEN
         I1 = NOCORB(J1,JA)
      ELSE
         I1 = NOCORB(J2,JB)
      ENDIF
!
!     This code allows a maximum of 8 open shells but at
!     most two shells may differ in the present case.  IHSH is an
!     upper bound to the number of shells.
!
      I2HSH = IHSH + ISH - 1
!
! --- FIRST CONSIDER THE L.H.S. (I=1) OF THE MATRIX ELEMENT.  NCC=1
!     MEANS UNOCCUPIED, REPRESENTED BY A DUMMY SINGLET S SHELL, AND THE
!     ADDITIONAL SET OF COUPLING QUANTUM NUMBERS WILL BE THE SAME AS THE
!     LAST SET OF COUPLING QUANTUM NUMBERS ALREADY OBTAINED.
!     NCC=2 MEANS OCCUPIED.  THEN ALL THE NEW QUANTUM NUMBERS (BOTH FOR
!     THE SHELL AND FOR THE COUPLING OF THIS SHELL TO THE RESULTANT OF
!     THE PREVIOUS ONES) ARE DEFINED IN THE CORRESPONDING J1QNRD ARRAY.
!     NOSH - THE NUMBER OF ELECTRONS IN THIS SHELL, IS DEFINED BY THE
!     APPROPRIATE ENTRY IN NELCSH .  THE R.H.S. IS THEN CONSIDERED
!     SIMILARLY (I=2)
!
      I = 1
   25 CONTINUE
      IF (I == 1) THEN
         IC = J1
         JC = JA
         NCC = NA
      ELSE
         IC = J2
         JC = JB
         NCC = NB
      ENDIF
      IF (NCC == 1) THEN
!       .. shell is not occupied
         NOSH(ISH,I) = 0
         J1QN(ISH,1,I) = 0
         J1QN(ISH,2,I) = 1
         J1QN(ISH,3,I) = 1
         IF (ISH == 2) THEN
            J1QN(I2HSH,1,I) = 0
            J1QN(I2HSH,2,I) = J1QN(1,2,I)
            J1QN(I2HSH,3,I) = J1QN(1,3,I)
         ELSE IF (ISH > 2) THEN
            J1QN(I2HSH,:,I) = J1QN(I2HSH-1,:,I)
         ENDIF
      ELSE
!       .. shell is occupied
         NOSH(ISH,I) = NELCSH(IC,JC)
         JD = J1QNRD(IC,JC)
         J1QN(ISH,1,I) = MOD(JD,64)
         JD = JD/64
         J1QN(ISH,2,I) = MOD(JD,64)
         J1QN(ISH,3,I) = JD/64
!
         IF (ISH > 1) THEN
!
!         .. a resultant coupling is present
!
!     IS THIS THE FIRST OCCUPIED SHELL OF THIS CONFIGURATION, THOUGH NOT
!     THE FIRST OF THE OTHER CONFIGURATION.  IF SO, THE INTERMEDIATE
!     COUPLING FORMED HAS THE SAME  L,S  VALUES AS THIS OCCUPIED SHELL,
!     SINCE WE COUPLE THE SHELL TO A DUMMY SINGLET S.
!
            IF (IC == 1) THEN
               I2 = 1
            ELSE
               I2 = NOCCSH(JC) + IC - 1
            ENDIF
            JD = J1QNRD(I2,JC)
            IF (IC == 1) THEN
               J1QN(I2HSH,1,I) = 0
            ELSE
               J1QN(I2HSH,1,I) = MOD(JD,64)
            ENDIF
            JD = JD/64
            J1QN(I2HSH,2,I) = MOD(JD,64)
            J1QN(I2HSH,3,I) = JD/64
         ENDIF
      ENDIF
!
      I = I + 1
      IF (I == 2) GO TO 25
!
! --- SET THE NJ AND LJ VALUES OF THE OCCUPIED SHELLS
!
      NJ(ISH) = NJCOMP(I1)
      IJFUL(ISH) = I1
      LJ(ISH) = LJCOMP(I1)
      IF (NC /= 0) THEN
         IF (INDL(NCI1) == I1) KPL = ISH
         IF (INDR(NCI1) == I1) KPR = ISH
      ENDIF
      RETURN
!
      END SUBROUTINE SETUPM


text:
!##############################################################################
!# ****************************************************************************
!# <name> poisson2d_method2_collect </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module is a demonstration program how to solve a simple Poisson
!# problem with constant coefficients on a simple domain.
!#
!# The routine splits up the tasks of reading the domain, creating
!# triangulations, discretisation, solving, postprocessing and creanup into
!# different subroutines. The communication between these subroutines
!# is done using a collection structure that saves problem-dependent data.
!# </purpose>
!##############################################################################

module poisson2d_method2_collect

  use fsystem
  use genoutput
  use storage
  use linearsolver
  use boundary
  use derivatives
  use element
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use filtersupport
  use linearsystemscalar
  use linearsystemblock
  use matrixfilters
  use vectorfilters
  use discretebc
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use scalarpde
  use ucd
  use pprocerror

  use collection

  use poisson2d_callback

  implicit none

contains

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

!<subroutine>

  subroutine pm2_initParamTriang (ilv,rcollection)

!<description>
  ! This routine initialises the parametrisation and triangulation of the
  ! domain. The corresponding .prm/.tri files are read from disc and
  ! the triangulation is refined as described by the parameter ilv.
!</description>

!<inputoutput>
  ! The level up to where we refine the coarse mesh.
  ! If the number is too large, ilv is reduced to the maximum allowed value.
  integer, intent(inout) :: ilv

  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! An object for saving the domain:
    type(t_boundary), pointer :: p_rboundary

    ! An object for saving the triangulation on the domain
    type(t_triangulation), pointer :: p_rtriangulation

    ! Path to the mesh
    character(len=SYS_STRLEN) :: spredir

    ! Get the path $PREDIR from the environment, where to read .prm/.tri files
    ! from. If that does not exist, write to the directory "./pre".
    if (.not. sys_getenv_string("PREDIR", spredir)) spredir = "./pre"

    ! At first, read in the parametrisation of the boundary and save
    ! it to p_rboundary.
    allocate(p_rboundary)
    call boundary_read_prm(p_rboundary, trim(spredir)//"/QUAD.prm")

    ! Now read in the basic triangulation.
    allocate(p_rtriangulation)
    call tria_readTriFile2D (p_rtriangulation, trim(spredir)//"/QUAD.tri", p_rboundary)

    ! Refine it.
    call tria_quickRefine2LevelOrdering (ilv-1,p_rtriangulation,p_rboundary)

    ! And create information about adjacencies and everything one needs from
    ! a triangulation.
    call tria_initStandardMeshFromRaw (p_rtriangulation,p_rboundary)

    ! The TRIAS(,)-array is now part pf the triangulation structure,
    ! we do not need it anymore.
    !
    ! *What* we need later is the definition of the boundary and the
    ! triangulation. Save these to the collection.
    call collct_setvalue_bdry(rcollection,"DOMAIN",p_rboundary,.true.)
    call collct_setvalue_tria(rcollection,"TRIA",p_rtriangulation,.true.)

  end subroutine

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

!<subroutine>

  subroutine pm2_initDiscretisation (rcollection)

!<description>
  ! This routine initialises the discretisation structure of the underlying
  ! problem and saves it to the collection.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! An object for saving the domain:
    type(t_boundary), pointer :: p_rboundary

    ! An object for saving the triangulation on the domain
    type(t_triangulation), pointer :: p_rtriangulation

    ! An object for the spatial discretisation
    type(t_blockDiscretisation), pointer :: p_rdiscretisation

    ! A cubature information structure
    type(t_scalarCubatureInfo), pointer :: p_rcubatureInfo

    ! Ask the collection to give us the boundary and triangulation.
    ! We need it for the discretisation.
    p_rboundary => collct_getvalue_bdry(rcollection,"DOMAIN")
    p_rtriangulation => collct_getvalue_tria(rcollection,"TRIA")

    ! Now we can start to initialise the discretisation. At first, set up
    ! a block discretisation structure that specifies the blocks in the
    ! solution vector. In this simple problem, we only have one block.
    allocate(p_rdiscretisation)
    call spdiscr_initBlockDiscr (p_rdiscretisation,1,&
                                 p_rtriangulation, p_rboundary)

    ! p_rdiscretisation%Rdiscretisations is a list of scalar discretisation
    ! structures for every component of the solution vector.
    ! Initialise the first element of the list to specify the element
    ! for this solution component:
    call spdiscr_initDiscr_simple (p_rdiscretisation%RspatialDiscr(1), &
                                   EL_Q1,p_rtriangulation, p_rboundary)

    ! Add the discretisation structure to the collection so that
    ! we can use it later.
    call collct_setvalue_bldiscr(rcollection,"DISCR2D",p_rdiscretisation,.true.)

    ! Set up an cubature info structure to tell the code which cubature
    ! formula to use
    ! Create an assembly information structure which tells the code
    ! the cubature formula to use. Standard: Gauss 3x3.
    allocate(p_rcubatureInfo)
    call spdiscr_createDefCubStructure(&
        p_rdiscretisation%RspatialDiscr(1),p_rcubatureInfo,CUB_GEN_AUTO_G3)
    call collct_setvalue_cubinfo(rcollection,"CUBINFO2D",p_rcubatureInfo,.true.)

  end subroutine

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

!<subroutine>

  subroutine pm2_initMatVec (rcollection)

!<description>
  ! Calculates the system matrix and RHS vector of the linear system
  ! by discretising the problem with the default discretisation structure
  ! in the collection.
  ! Sets up a solution vector for the linear system.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! A bilinear and linear form describing the analytic problem to solve
    type(t_bilinearForm) :: rform
    type(t_linearForm) :: rlinform

    ! A matrix and vector. The vector accepts the RHS of the problem
    ! in scalar form. The matrix will receive the discrete Laplace operator.
    type(t_matrixBlock), pointer :: p_rmatrix
    type(t_vectorBlock), pointer :: p_rrhs,p_rvector

    ! A pointer to the discretisation structure with the data.
    type(t_blockDiscretisation), pointer :: p_rdiscretisation

    ! A cubature information structure
    type(t_scalarCubatureInfo), pointer :: p_rcubatureInfo

    ! Ask the collection to give us the discretisation structure
    p_rdiscretisation => collct_getvalue_bldiscr(rcollection,"DISCR2D")

    ! Get the cubature information structure
    p_rcubatureInfo => collct_getvalue_cubinfo(rcollection,"CUBINFO2D")

    ! Create the matrix and the vectors on the heap
    allocate(p_rmatrix)
    allocate(p_rrhs)
    allocate(p_rvector)

    ! Initialise the block matrix with default values based on
    ! the discretisation.
    call lsysbl_createMatrix (p_rdiscretisation,p_rmatrix)

    ! Save matrix and vectors to the collection.
    ! They maybe used later, expecially in nonlinear problems.
    call collct_setvalue_vec(rcollection,"RHS",p_rrhs,.true.)
    call collct_setvalue_vec(rcollection,"SOLUTION",p_rvector,.true.)
    call collct_setvalue_mat(rcollection,"LAPLACE",p_rmatrix,.true.)

    ! Now as the discretisation is set up, we can start to generate
    ! the structure of the system matrix which is to solve.
    ! We create that directly in the block (1,1) of the block matrix
    ! using the discretisation structure of the first block.
    call bilf_createMatrixStructure (p_rmatrix, 1, 1, LSYSSC_MATRIX9)

    ! And now to the entries of the matrix. For assembling of the entries,
    ! we need a bilinear form, which first has to be set up manually.
    ! We specify the bilinear form (grad Psi_j, grad Phi_i) for the
    ! scalar system matrix in 2D.

    rform%itermCount = 2
    rform%Idescriptors(1,1) = DER_DERIV_X
    rform%Idescriptors(2,1) = DER_DERIV_X
    rform%Idescriptors(1,2) = DER_DERIV_Y
    rform%Idescriptors(2,2) = DER_DERIV_Y

    ! In the standard case, we have constant coefficients:
    rform%ballCoeffConstant = .true.
    rform%Dcoefficients(1)  = 1.0
    rform%Dcoefficients(2)  = 1.0

    ! Now we can build the matrix entries.
    ! We specify the callback function coeff_Laplace for the coefficients.
    ! As long as we use constant coefficients, this routine is not used.
    ! By specifying ballCoeffConstant = .FALSE. above,
    ! the framework will call the callback routine to get analytical data.
    call bilf_buildMatrixScalar (rform,.true.,&
        p_rmatrix%RmatrixBlock(1,1),p_rcubatureInfo)

    ! Next step: Create a RHS vector and a solution vector and a temporary
    ! vector. All are filled with zero.
    call lsysbl_createVector (p_rdiscretisation,p_rrhs,.true.)
    call lsysbl_createVector (p_rdiscretisation,p_rvector,.true.)

    ! The vector structure is done but the entries are missing.
    ! So the next thing is to calculate the content of that vector.
    !
    ! At first set up the corresponding linear form (f,Phi_j):
    rlinform%itermCount = 1
    rlinform%Idescriptors(1) = DER_FUNC

    ! ... and then discretise the RHS to the first subvector of
    ! the block vector using the discretisation structure of the
    ! first block.
    call linf_buildVectorScalar (&
              rlinform,.true.,p_rrhs%RvectorBlock(1),p_rcubatureInfo,&
              coeff_RHS_2D)

  end subroutine

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

!<subroutine>

  subroutine pm2_initDiscreteBC (rcollection)

!<description>
  ! This calculates the discrete version of the boundary conditions and
  ! assigns it to the system matrix and RHS vector.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! A pointer to the system matrix and the RHS vector as well as
    ! the discretisation
    type(t_matrixBlock), pointer :: p_rmatrix
    type(t_vectorBlock), pointer :: p_rrhs,p_rvector
    type(t_blockDiscretisation), pointer :: p_rdiscretisation

    type(t_boundaryRegion) :: rboundaryRegion

    ! Pointer to structure for saving discrete BC`s:
    type(t_discreteBC), pointer :: p_rdiscreteBC

    ! A pointer to the domain
    type(t_boundary), pointer :: p_rboundary

    ! Ask the collection to give us the discretisation structure and
    p_rdiscretisation => collct_getvalue_bldiscr(rcollection,"DISCR2D")

    ! Get the domain from the discretisation
    p_rboundary => p_rdiscretisation%p_rboundary

    ! Get our matrix and right hand side from the collection.
    p_rrhs    => collct_getvalue_vec(rcollection,"RHS")
    p_rvector => collct_getvalue_vec(rcollection,"SOLUTION")
    p_rmatrix => collct_getvalue_mat(rcollection,"LAPLACE")

    ! Create a t_discreteBC structure where we store all discretised boundary
    ! conditions.
    allocate (p_rdiscreteBC)
    call bcasm_initDiscreteBC(p_rdiscreteBC)
    !
    ! We "know" already (from the problem definition) that we have four boundary
    ! segments in the domain. Each of these, we want to use for enforcing
    ! some kind of boundary condition.
    !
    ! We ask the boundary routines to create a "boundary region" - which is
    ! simply a part of the boundary corresponding to a boundary segment.
    ! A boundary region roughly contains the type, the min/max parameter value
    ! and whether the endpoints are inside the region or not.
    call boundary_createRegion(p_rboundary,1,1,rboundaryRegion)

    ! We use this boundary region and specify that we want to have Dirichlet
    ! boundary there. The following call does the following:
    ! - Create Dirichlet boundary conditions on the region rboundaryRegion.
    !   We specify icomponent="1" to indicate that we set up the
    !   Dirichlet BC`s for the first (here: one and only) component in the
    !   solution vector.
    ! - Discretise the boundary condition so that the BC`s can be applied
    !   to matrices and vectors
    ! - Add the calculated discrete BC`s to rdiscreteBC for later use.
    call bcasm_newDirichletBConRealBD (p_rdiscretisation,1,&
        rboundaryRegion,p_rdiscreteBC,&
        getBoundaryValues_2D,rcollection)

    ! Now to the edge 2 of boundary component 1 the domain.
    call boundary_createRegion(p_rboundary,1,2,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (p_rdiscretisation,1,&
        rboundaryRegion,p_rdiscreteBC,&
        getBoundaryValues_2D,rcollection)

    ! Edge 3 of boundary component 1.
    call boundary_createRegion(p_rboundary,1,3,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (p_rdiscretisation,1,&
        rboundaryRegion,p_rdiscreteBC,&
        getBoundaryValues_2D,rcollection)

    ! Edge 4 of boundary component 1. That is it.
    call boundary_createRegion(p_rboundary,1,4,rboundaryRegion)
    call bcasm_newDirichletBConRealBD (p_rdiscretisation,1,&
        rboundaryRegion,p_rdiscreteBC,&
        getBoundaryValues_2D,rcollection)

    ! Save the structures to the collection, so we can access them
    ! without having the matrix or the vector.
    call collct_setvalue_discbc(rcollection,"DISCBC",p_rdiscreteBC,.true.)

  end subroutine

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

!<subroutine>

  subroutine pm2_implementBC (rcollection)

!<description>
  ! Implements boundary conditions into the RHS and into a given solution vector.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! A pointer to the system matrix and the RHS vector as well as
    ! the discretisation
    type(t_matrixBlock), pointer :: p_rmatrix
    type(t_vectorBlock), pointer :: p_rrhs,p_rvector
    type(t_discreteBC), pointer :: p_rdiscreteBC

    ! Get our matrix and right hand side from the collection.
    p_rrhs    => collct_getvalue_vec(rcollection,"RHS")
    p_rvector => collct_getvalue_vec(rcollection,"SOLUTION")
    p_rmatrix => collct_getvalue_mat(rcollection,"LAPLACE")
    p_rdiscreteBC => collct_getvalue_discbc(rcollection,"DISCBC")

    ! Next step is to implement boundary conditions into the RHS,
    ! solution and matrix. This is done using a vector/matrix filter
    ! for discrete boundary conditions.
    call vecfil_discreteBCrhs (p_rrhs,p_rdiscreteBC)
    call vecfil_discreteBCsol (p_rvector,p_rdiscreteBC)
    call matfil_discreteBC (p_rmatrix,p_rdiscreteBC)

  end subroutine

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

!<subroutine>

  subroutine pm2_solve (rcollection)

!<description>
  ! Solves the given problem by applying a linear solver.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! Error indicator during initialisation of the solver
    integer :: ierror

    ! A filter chain to filter the vectors and the matrix during the
    ! solution process.
    type(t_filterChain), dimension(1), target :: RfilterChain

    ! Number of filters in the filter chain
    integer :: nfilters

    ! A pointer to the system matrix and the RHS vector as well as
    ! the discretisation
    type(t_matrixBlock), pointer :: p_rmatrix
    type(t_vectorBlock), pointer :: p_rrhs,p_rvector
    type(t_vectorBlock), target :: rvecTmp

    ! A solver node that accepts parameters for the linear solver
    type(t_linsolNode), pointer :: p_rsolverNode,p_rpreconditioner

    ! Pointer to structure for saving discrete BC`s:
    type(t_discreteBC), pointer :: p_rdiscreteBC

    ! Get our matrix and right hand side from the problem structure.
    p_rrhs    => collct_getvalue_vec(rcollection,"RHS")
    p_rvector => collct_getvalue_vec(rcollection,"SOLUTION")
    p_rmatrix => collct_getvalue_mat(rcollection,"LAPLACE")

    ! Create a temporary vector for the solver - it needs that.
    call lsysbl_createVecBlockIndirect (p_rrhs, rvecTmp, .false.)

    ! During the linear solver, the boundary conditions must
    ! frequently be imposed to the vectors. This is done using
    ! a filter chain. As the linear solver does not work with
    ! the actual solution vectors but with defect vectors instead,
    ! a filter for implementing the real boundary conditions
    ! would be wrong.
    ! Therefore, create a filter chain with one filter only,
    ! which implements Dirichlet-conditions into a defect vector.

    p_rdiscreteBC => collct_getvalue_discbc(rcollection,"DISCBC")

    call filter_initFilterChain (RfilterChain,nfilters)
    call filter_newFilterDiscBCDef (RfilterChain,nfilters,p_rdiscreteBC)

    ! Create a BiCGStab-solver. Attach the above filter chain
    ! to the solver, so that the solver automatically filters
    ! the vector during the solution process.
    nullify(p_rpreconditioner)
    call linsol_initBiCGStab (p_rsolverNode,p_rpreconditioner,RfilterChain)

    ! Set the output level of the solver to 2 for some output
    p_rsolverNode%ioutputLevel = 2

    ! Attach the system matrix to the solver.
    call linsol_setMatrix(p_RsolverNode,p_rmatrix)

    ! Initialise structure/data of the solver. This allows the
    ! solver to allocate memory / perform some precalculation
    ! to the problem.
    call linsol_initStructure (p_rsolverNode, ierror)

    if (ierror .ne. LINSOL_ERR_NOERROR) then
      call output_line("Matrix structure invalid!",OU_CLASS_ERROR)
      call sys_halt()
    end if

    call linsol_initData (p_rsolverNode, ierror)

    if (ierror .ne. LINSOL_ERR_NOERROR) then
      call output_line("Matrix singular!",OU_CLASS_ERROR)
      call sys_halt()
    end if

    ! Finally solve the system. As we want to solve Ax=b with
    ! b being the real RHS and x being the real solution vector,
    ! we use linsol_solveAdaptively. If b is a defect
    ! RHS and x a defect update to be added to a solution vector,
    ! we would have to use linsol_precondDefect instead.
    call linsol_solveAdaptively (p_rsolverNode,p_rvector,p_rrhs,rvecTmp)

    ! Release solver data and structure
    call linsol_doneData (p_rsolverNode)
    call linsol_doneStructure (p_rsolverNode)

    ! Release the solver node and all subnodes attached to it (if at all):
    call linsol_releaseSolver (p_rsolverNode)

    ! Release the filter chain
    call filter_doneFilterChain (RfilterChain,nfilters)

    ! Release the temporary vector
    call lsysbl_releaseVector (rvecTmp)

  end subroutine

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

!<subroutine>

  subroutine pm2_postprocessing (rcollection)

!<description>
  ! Writes the solution into a VTK file.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

  ! local variables

    ! Output block for UCD output to VTK file
    type(t_ucdExport) :: rexport
    character(len=SYS_STRLEN) :: sucddir

    ! A cubature information structure
    type(t_scalarCubatureInfo), pointer :: p_rcubatureInfo

    ! A pointer to the solution vector and to the triangulation.
    type(t_vectorBlock), pointer :: p_rvector
    type(t_triangulation), pointer :: p_rtriangulation

    ! Error of FE function to reference function
    real(DP) :: derror

    ! Get the solution vector from the collection.
    p_rvector => collct_getvalue_vec(rcollection,"SOLUTION")

    ! Get the cubature information structure
    p_rcubatureInfo => collct_getvalue_cubinfo(rcollection,"CUBINFO2D")

    ! From the attached discretisation, get the underlying triangulation
    p_rtriangulation => &
      p_rvector%RvectorBlock(1)%p_rspatialDiscr%p_rtriangulation

    ! p_rvector now contains our solution. We can now
    ! start the postprocessing.
    !
    ! Get the path for writing postprocessing files from the environment variable
    ! $UCDDIR. If that does not exist, write to the directory "./gmv".
    if (.not. sys_getenv_string("UCDDIR", sucddir)) sucddir = "./gmv"

    ! Start UCD export to VTK file:
    call ucd_startVTK (rexport,UCD_FLAG_STANDARD,p_rtriangulation,&
                       trim(sucddir)//"/u2d_2_collect.vtk")

    ! Add the solution to the UCD exporter
    call ucd_addVectorByVertex (rexport, "sol", UCD_VAR_STANDARD, &
        p_rvector%RvectorBlock(1))

    ! Write the file to disc, that is it.
    call ucd_write (rexport)
    call ucd_release (rexport)

    ! Calculate the error to the reference function.
    call pperr_scalar (PPERR_L2ERROR,derror,p_rvector%RvectorBlock(1),&
                       getReferenceFunction_2D, rcubatureInfo=p_rcubatureInfo)
    call output_line ("L2-error: " // sys_sdEL(derror,10) )

    call pperr_scalar (PPERR_H1ERROR,derror,p_rvector%RvectorBlock(1),&
                       getReferenceFunction_2D, rcubatureInfo=p_rcubatureInfo)
    call output_line ("H1-error: " // sys_sdEL(derror,10) )

  end subroutine

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

!<subroutine>

  subroutine pm2_doneMatVec (rcollection)

!<description>
  ! Releases system matrix and vectors.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

    ! local variables

    ! A pointer to the system matrix and the RHS vector as well as
    ! the discretisation
    type(t_matrixBlock), pointer :: p_rmatrix
    type(t_vectorBlock), pointer :: p_rrhs,p_rvector

    ! Get our matrix and right hand side from the collection.
    p_rrhs    => collct_getvalue_vec(rcollection,"RHS")
    p_rvector => collct_getvalue_vec(rcollection,"SOLUTION")
    p_rmatrix => collct_getvalue_mat(rcollection,"LAPLACE")

    ! Release them from memory
    call lsysbl_releaseVector (p_rvector)
    call lsysbl_releaseVector (p_rrhs)
    call lsysbl_releaseMatrix (p_rmatrix)
    deallocate(p_rmatrix)
    deallocate(p_rrhs)
    deallocate(p_rvector)

    ! Delete the variables from the collection.
    call collct_deletevalue (rcollection,"RHS")
    call collct_deletevalue (rcollection,"SOLUTION")
    call collct_deletevalue (rcollection,"LAPLACE")

  end subroutine

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

!<subroutine>

  subroutine pm2_doneBC (rcollection)

!<description>
  ! Releases discrete and analytic boundary conditions from the heap.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

    ! local variables

    ! Pointer to structure for saving discrete BC`s:
    type(t_discreteBC), pointer :: p_rdiscreteBC

    ! Get pointers to the discrete and analytic boundary conditions
    ! from the collection structure.
    p_rdiscreteBC => collct_getvalue_discbc(rcollection,"DISCBC")

    ! Release our discrete version of the boundary conditions
    call bcasm_releaseDiscreteBC (p_rdiscreteBC)
    deallocate(p_rdiscreteBC)

    ! Delete the variables from the collection.
    call collct_deletevalue (rcollection,"BDCOND")
    call collct_deletevalue (rcollection,"DISCBC")

  end subroutine

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

!<subroutine>

  subroutine pm2_doneDiscretisation (rcollection)

!<description>
  ! Releases the discretisation from the heap.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

    ! A pointer to the discretisation structure with the data.
    type(t_blockDiscretisation), pointer :: p_rdiscretisation

      ! A cubature information structure
    type(t_scalarCubatureInfo), pointer :: p_rcubatureInfo

    ! Ask the collection to give us the discretisation structure
    p_rdiscretisation => collct_getvalue_bldiscr(rcollection,"DISCR2D")

    ! Delete the discretisation...
    call spdiscr_releaseBlockDiscr(p_rdiscretisation)

    ! remove the allocated block discretisation structure
    deallocate(p_rdiscretisation)

    ! and remove it from the collection.
    call collct_deletevalue (rcollection,"DISCR2D")

    ! A pointer to the cubature information structure with the data.
    p_rcubatureInfo => collct_getvalue_cubinfo(rcollection,"CUBINFO2D")

    ! Delete the structure...
    call spdiscr_releaseCubStructure(p_rcubatureInfo)

    ! remove the allocated block discretisation structure
    deallocate(p_rcubatureInfo)

    ! and remove it from the collection.
    call collct_deletevalue (rcollection,"CUBINFO2D")

  end subroutine

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

!<subroutine>

  subroutine pm2_doneParamTriang (rcollection)

!<description>
  ! Releases the triangulation and parametrisation from the heap.
!</description>

!<inputoutput>
  ! A collection object for saving structural data and some problem-dependent
  ! information.
  type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!</subroutine>

    ! An object for saving the domain:
    type(t_boundary), pointer :: p_rboundary

    ! An object for saving the triangulation on the domain
    type(t_triangulation), pointer :: p_rtriangulation

    ! Ask the collection to give us the boundary and triangulation.
    ! We need it for the discretisation.
    p_rboundary => collct_getvalue_bdry(rcollection,"DOMAIN")
    p_rtriangulation => collct_getvalue_tria(rcollection,"TRIA")

    ! Release the triangulation.
    call tria_done (p_rtriangulation)
    deallocate(p_rtriangulation)

    ! Finally release the domain.
    call boundary_release (p_rboundary)
    deallocate(p_rboundary)

    ! Remove everything from the collection.
    call collct_deletevalue (rcollection,"TRIA")
    call collct_deletevalue (rcollection,"DOMAIN")

  end subroutine

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

!<subroutine>

  subroutine poisson2d_2_collect

!<description>
  ! This is a "separated" poisson solver for solving a Poisson
  ! problem. The different tasks of the problem are separated into
  ! subroutines. The problem uses a collection structure for the communication:
  ! All subroutines add their generated information to the collection, so that
  ! the other subroutines can work with them.
  !
  ! The following tasks are performed by the subroutines:
  !
  ! 1.) Read in parametrisation
  ! 2.) Read in triangulation
  ! 3.) Set up RHS
  ! 4.) Set up matrix
  ! 5.) Create solver structure
  ! 6.) Solve the problem
  ! 7.) Write solution to VTK file
  ! 8.) Release all variables, finish
!</description>

!</subroutine>

    ! NLMAX receives the level where we want to solve.
    integer :: NLMAX

    ! A collection structure for our problem
    type(t_collection) :: rcollection

    ! Ok, let us start.
    ! We want to solve our Poisson problem on level...

    NLMAX = 7

    ! Initialise the collection.
    call collct_init (rcollection)

    ! So now the different steps - one after the other.
    !
    ! Initialisation
    call pm2_initParamTriang (NLMAX,rcollection)
    call pm2_initDiscretisation (rcollection)
    call pm2_initMatVec (rcollection)
    call pm2_initDiscreteBC (rcollection)

    ! Implementation of boundary conditions
    call pm2_implementBC (rcollection)

    ! Solve the problem
    call pm2_solve (rcollection)

    ! Postprocessing
    call pm2_postprocessing (rcollection)

    ! Cleanup
    call pm2_doneMatVec (rcollection)
    call pm2_doneBC (rcollection)
    call pm2_doneDiscretisation (rcollection)
    call pm2_doneParamTriang (rcollection)

    ! Print some statistical data about the collection - anything forgotten?
    call output_lbrk ()
    call output_line ("Remaining collection statistics:")
    call output_line ("--------------------------------")
    call output_lbrk ()
    call collct_printStatistics (rcollection)

    ! Finally release the collection.
    call collct_done (rcollection)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
C    %W% %G%
        subroutine dcaux (dname,descr,icount,t,work)
c       SUBROUTINE DCAUX (DESCR,ICOUNT,T,WORK)
          character dname*(*),descr*(*)                                 !dem
          character*20 atem
          dimension t(*), work(*)
C * * *
C * * * THIS SUBROUTINE WRITES REQUESTED DC DATA TO THE AUXILIARY
C * * * OUTPUT FILE (FOR011).  IT IS CALLED BY NDCOUT AND CALLS
C * * * VARFMT.
C Revs:
c  Oct/05/92 - DEM: Expanded field sizes for alternate format; 6 dec
c                     places for time and 15 chars for parameters.
c  May/07/92 - DEM: Added DNAME argument so that every t-parm list could
c                   have branch name in its header for the new aux
c                   format.
c       DNAME . . Name of DC branch
C       DESCR . . Parm description
C       ICOUNT  . No of time steps
C       T . . . . Time step array
C       WORK  . . Array containing parm values
C       L11 . . . AUXILIARY OUTPUT FILE
      include 'tspinc/params.inc'                                              !dem
      include 'tspinc/blkcom1.inc'
      include 'tspinc/link56.inc'                                              !dem
c       -     begin     begin     begin     begin     begin     begin
      if (auxfmt .eq. 'STD') then                                       !dem
        write(l11,100) descr
100       format(a)
      else                                                              !dem
        write (l11,'(A)') '/CURVE'                                      !dem
        write (l11,'(2A)') '  LEGEND  ',dname                           !dem
        write (l11,'(2A)') '  ! ',descr                                 !dem
        write (l11,'(2A)') '  DATA'                                     !dem
      endif                                                             !dem
      do 105 j = 1, icount
c       -  New format has wider fields
        if (auxfmt .eq. 'STD') then                                     !dem
          call varfmt(work(j),atem,10)                                  !dem
          write(l11,'(F9.2,1X,A10)') t(j), atem                         !dem
        else                                                            !dem
          call varfmt(work(j),atem,15)                                  !dem
          write(l11,'(F13.6,1X,A15)') t(j), atem                        !dem
        endif                                                           !dem
c       CALL VARFMT(WORK(J),ATEM,10)
c       WRITE(L11,110) T(J), ATEM
c110      FORMAT(F9.2,1X,A10)
105   continue
      if (auxfmt .eq. 'STD') then                                       !dem
        write(l11,120)
120       format('-100.0')
      else                                                              !dem
        write (l11,'(A)') '/END'                                        !dem
      endif                                                             !dem
c
      return
      end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
! MIT License
!
! Copyright (c) 2019 YuyaM
!
! 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.
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !  MAIN PRORGRAM  phaseAverage                             !
      !                                                          !
      !                           2019/03/08                     !
      !                           WRITTEN BY YUYA MIKI           !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !<
      !! @param ME :: maximum element number
      !! @param MP :: maximum node number
      !! @param N  :: the number of node in a element
      !! @param NODE :: node number list
      !! @param X,Y,Z :: Position
      !! @param U,V,W :: Velocity
      !! @param P     :: Pressure at element
      !! @param PN    :: Pressure at node
      !! @param numPhase :: phase average interval
      !! @param NSUM  :: Number of Flow in FLOWS
      program phaseAverage
      implicit none
!
      include 'gf2.h'
      ! parameter for GFALL
      DATA IWRITE / 2 /
      DATA INAME  / 1 /

      integer(4) :: ME,MP
      integer(4),parameter :: N=8
!
      integer(4), allocatable :: NODE(:,:)
      real(4),allocatable,dimension(:)   :: X,Y,Z,U,V,W,P,PN
      real(4),allocatable,dimension(:,:) :: UA,VA,WA,PA,PNA
!
      integer(4)    :: ISTEP,NPCHK,NEPRS,NPPRS
      real(4)       :: TIMEP,INVNUM
      character(60) :: FILEMS,FILEAR,FILEFF,FILEAV,FILEAVE,FILEFS
      integer(4)    :: numPhase,ISUM,iPhase,NSUM,IFLOWS,NFLOWS
!     [work]
      integer(4)    :: NP,NE,NDUM,IP,IE
      character(4)  :: CNUM
      ! fortran input output
      integer(4),parameter :: IUT0  = 0
      integer(4),parameter :: IUT5  = 5
      integer(4),parameter :: IUT6  = 6
      integer(4),parameter :: IUTMS = 11
      integer(4),parameter :: IUTFF = 12
      integer(4),parameter :: IUTAV = 13
      integer(4) :: IERR,FLAGOK,FLAGGFSEP,FLAGDELIM
!
! read MESH DATA
!
      write(IUT6,*)
      write(IUT6,*) '** phaseAverage: Phase averaging FLOWS data **'
!
      do
          write(IUT6,*)
          write(IUT6,*) '[In]Specify filename for MESH data'
          read (IUT5,'(A60)') FILEMS
          write(IUT6,*) '[In]Specify filename for FLOWS data'
          read (IUT5,'(A60)') FILEFF
          write(IUT6,*) '[Out]Specify filename for AVE  data'
          read (IUT5,'(A60)') FILEAV
          write(IUT6,*) 'Specify interval number for phase-average'
          read (IUT5,*) numPhase
          write(IUT6,*) 'if FLOWS IS SEPARATED    :1'
          write(IUT6,*) 'if FLOWS IS NOT SEPARATED:0'
          read (IUT5,*) FLAGGFSEP
          write(IUT6,*) 'Specify number of flowsfor phase-average'
          read (IUT5,*) NFLOWS
          write(IUT6,*) 'if want to change delimiter .P****'
          read (IUT5,*) FLAGDELIM
          IF(FLAGDELIM.eq.1)THEN
            write(IUT6,*) 'USE .P**** :0'
            write(IUT6,*) 'USE .t**** :1'
            read (IUT5,*) FLAGDELIM
          ENDIF

          write(IUT6,'(A10,A60)') "FILEMS   =", FILEMS
          write(IUT6,'(A10,A60)') "FILEFF   =", FILEFF
          write(IUT6,'(A10,A60)') "FILEAV   =", FILEAV
          write(IUT6,'(A10,I6)')  "numPhase    =", numPhase
          write(IUT6,'(A10,I6)')  "FLAGGFSEP=", FLAGGFSEP
          write(IUT6,'(A10,I6)')  "NFLOWS   =", NFLOWS
          write(IUT6,*) 'These parameters are ok? Yes:1,No:0'
          read(IUT5,*) FLAGOK
          if(FLAGOK.EQ.1) exit
      enddo
!
! check data size using MESH file
!
      call CHKMSH(FILEMS,IUTMS,IUT0,IUT6,MP,ME,IERR)
      if (IERR.NE.0) STOP
!
! allocate variables
!
      write(IUT6,*)
      write(IUT6,*)
      write(IUT6,*) 'phaseAverage: ALLOCATING VARIABLES... '
      allocate(   NODE(N,ME), STAT=LERR(01))
      allocate(        X(MP), STAT=LERR(02))
      allocate(        Y(MP), STAT=LERR(03))
      allocate(        Z(MP), STAT=LERR(04))
      allocate(        U(MP), STAT=LERR(05))
      allocate(        V(MP), STAT=LERR(06))
      allocate(        W(MP), STAT=LERR(07))
      allocate(        P(ME), STAT=LERR(08))
      allocate(       PN(MP), STAT=LERR(09))
      allocate( UA(MP,numPhase), STAT=LERR(10))
      allocate( VA(MP,numPhase), STAT=LERR(11))
      allocate( WA(MP,numPhase), STAT=LERR(12))
      allocate( PA(ME,numPhase), STAT=LERR(13))
      allocate(PNA(MP,numPhase), STAT=LERR(14))
      call CHKALC(14,LERR,IUT6,IERR)
      if(IERR.NE.0) then
          write(IUT6,*) 'phaseAverage: allocating error       '
          STOP
      endif
      write(IUT6,*) 'phaseAverage: allocating finish       '

      write(IUT6,*)
      write(IUT6,*) 'reading MESH data...'
!
! read MESH file
!
      IACT = 1
      call GFALL(IUT0,IUT6,IUTMS,FILEMS,
     *           MCOM,NCOMFL,COMFLE,
     *           MCOM,NCOMST,COMSET,
     *           IACT,IWRITE,INAME,IRESV,
     *           ICAST,IDATA0,IALL,ISKIP,IERR,
     *           '*GRID_3D *NODE_3D !',
     *           NAME,MP,NP,X,Y,Z,
     *           NAME,ME,N,NE,NDUM,NODE,
     *           ICHECK)
      if(IERR.NE.0) STOP
!
! read FLOWS file
!
      IF(FLAGGFSEP.EQ.0)THEN
        IACT = 3
        ! Open FLOWS
        call GFALL(IUT0,IUT6,IUTFF,FILEFF,
     *             MCOM,NCOMFL,COMFLE,
     *             MCOM,NCOMST,COMSET,
     *             IACT,IWRITE,INAME,IRESV,
     *             ICAST,IDATA0,IALL,ISKIP,IERR,
     *             ' !',
     *             ICHECK)
        if(IERR.NE.0) STOP

        UA = 0.0E0
        VA = 0.0E0
        WA = 0.0E0
        PA = 0.0E0
        PNA = 0.0E0
        NSUM = 0
        NFLOWS = 0
        do
          IACT = 5
          ! Read FLOWS
          call GFALL(IUT0,IUT6,IUTFF,FILEFF,
     *               MCOM,NCOMFL,COMFLE,
     *               MCOM,NCOMST,COMSET,
     *               IACT,IWRITE,INAME,IRESV,
     *               ICAST,IDATA0,IALL,ISKIP,IERR,
     *               '*TIME_PS *STEP_PS *VELO_3D
     *                *PRES_3E *PRES_3D !',
     *               NAME,TIMEP,
     *               NAME,ISTEP,
     *               NAME,MP,NPCHK,U,V,W,
     *               NAME,ME,NEPRS,P,
     *               NAME,MP,NPPRS,PN,
     *               ICHECK)
          if(IERR.NE.0)   STOP
          ! NFLOWSはFLOWSの数
          NFLOWS = NFLOWS + 1
          ! ISUM = 0 - numPhase-1
          ! 1,2,3,..,0(numPhase),1,2,3..,0(numPhase)
          ISUM = mod(NFLOWS,numPhase)
          ! NSUMは平均した回数
          IF(ISUM == 0) NSUM = NSUM + 1
          ! ISUM+1 = 1 - numPhase
          do IP=1,NP
            UA (IP,ISUM+1) = UA (IP,ISUM+1) + U (IP)
            VA (IP,ISUM+1) = VA (IP,ISUM+1) + V (IP)
            WA (IP,ISUM+1) = WA (IP,ISUM+1) + W (IP)
            PNA(IP,ISUM+1) = PNA(IP,ISUM+1) + PN(IP)
          enddo
          do IE=1,NE
            PA (IE,ISUM+1) = PA (IE,ISUM+1) + P(IE)
          enddo
          write(IUT6,*) "averaged:when NFLOWS=",NFLOWS,
     *                                 ",NSUM=",NSUM,
     *                                 ",ISUM=",ISUM
          if(IACT .EQ. 7) exit
        end do

        IACT = 7
        ! Close FLOWS
        call GFALL(IUT0,IUT6,IUTFF,FILEFF,
     *             MCOM,NCOMFL,COMFLE,
     *             MCOM,NCOMST,COMSET,
     *             IACT,IWRITE,INAME,IRESV,
     *             ICAST,IDATA0,IALL,ISKIP,IERR,
     *             ' !',
     *             ICHECK)
        if(IERR.NE.0) STOP
      ELSEIF(FLAGGFSEP.EQ.1)THEN
        ! SEPARATED
        write(IUT6,*) "FLAGGFSEP",FLAGGFSEP
        write(IUT6,*) "NFLOWS",NFLOWS
        UA = 0.0E0
        VA = 0.0E0
        WA = 0.0E0
        PA = 0.0E0
        PNA = 0.0E0
        NSUM = 0
        do IFLOWS=1,NFLOWS
          write(CNUM,'(I4.4)') IFLOWS
          IF(FLAGDELIM.eq.0)THEN
            FILEFS=trim(FILEFF)//".P"//CNUM
          ELSEIF(FLAGDELIM.eq.1)THEN
            FILEFS=trim(FILEFF)//".t"//CNUM
          ENDIF
          ! Read FLOWS
          IACT = 1
          call GFALL(IUT0,IUT6,IUTFF,FILEFS,
     *               MCOM,NCOMFL,COMFLE,
     *               MCOM,NCOMST,COMSET,
     *               IACT,IWRITE,INAME,IRESV,
     *               ICAST,IDATA0,IALL,ISKIP,IERR,
     *               '*TIME_PS *STEP_PS *VELO_3D
     *                *PRES_3E *PRES_3D !',
     *               NAME,TIMEP,
     *               NAME,ISTEP,
     *               NAME,MP,NPCHK,U,V,W,
     *               NAME,ME,NEPRS,P,
     *               NAME,MP,NPPRS,PN,
     *               ICHECK)
          if(IERR.NE.0)   STOP
          ! ISUM = 0 - numPhase-1
          ! 1,2,3,..,0(numPhase),1,2,3..,0(numPhase)
          ISUM = mod(IFLOWS,numPhase)
          ! NSUMは平均した回数
          IF(ISUM == 0) NSUM = NSUM + 1
          ! ISUM+1 = 1 - numPhase
          do IP=1,NP
            UA (IP,ISUM+1) = UA (IP,ISUM+1) + U(IP)
            VA (IP,ISUM+1) = VA (IP,ISUM+1) + V(IP)
            WA (IP,ISUM+1) = WA (IP,ISUM+1) + W(IP)
            PNA(IP,ISUM+1) = PNA(IP,ISUM+1) + PN(IP)
          end do
          do IE=1,NE
            PA (IE,ISUM+1) = PA (IE,ISUM+1) + P(IE)
          end do
          write(IUT6,*) "averaged:when NFLOWS=",IFLOWS,",NSUM=",NSUM,
     *                                 ",ISUM=",ISUM
        end do
      ELSE
         write(IUT6,*) "FLAGGFSEP",FLAGGFSEP
         STOP
      ENDIF
!
! conduct average
!
      ! ゼロ割りを防ぐ
      if(NSUM < 1)then
          WRITE(IUT6,*) "ERROR IN FLAG1"
          STOP
      endif
      INVNUM = 1.0E0 / float(NSUM)
      do iPhase=1,numPhase
        do IP=1,NP
          UA (IP,iPhase) = UA (IP,iPhase) * INVNUM
          VA (IP,iPhase) = VA (IP,iPhase) * INVNUM
          WA (IP,iPhase) = WA (IP,iPhase) * INVNUM
          PNA(IP,iPhase) = PNA(IP,iPhase) * INVNUM
        end do
        do IE=1,NE
          PA (IE,iPhase) = PA (IE,iPhase) * INVNUM
        end do
      end do

!
! output AVE file
!
      do iPhase=1,numPhase
          TIMEP = 0.0E0
          ISTEP = NSUM
          do IP=1,NP
            U  (IP)= UA (IP,iPhase)
            V  (IP)= VA (IP,iPhase)
            W  (IP)= WA (IP,iPhase)
            PN (IP)= PNA(IP,iPhase)
          end do
          do IE=1,NE
            P(IE)  = PA (IE,iPhase)
          end do

          IACT = 2
          write(CNUM,'(I4.4)') iPhase
          FILEAVE = trim(FILEAV)//'.P'//CNUM
          call GFALL(IUT0,IUT6,IUTAV,FILEAVE,
     *               MCOM,NCOMFL,COMFLE,
     *               MCOM,NCOMST,COMSET,
     *               IACT,IWRITE,INAME,IRESV,
     *               ICAST,IDATA0,IALL,ISKIP,IERR,
     *               '*TIME_PS *STEP_PS *VELO_3D
     *                *PRES_3E *PRES_3D !',
     *               NAME,TIMEP,
     *               NAME,ISTEP,
     *               NAME,MP,NP,U,V,W,
     *               NAME,ME,NE,P,
     *               NAME,MP,NP,PN,
     *               ICHECK)
      end do

      STOP
      end program phaseAverage


text:
A Student Organizations student organization at UC Davis.  It has two divisions:

   Undergraduate Division  presumably premed students.
   Medical Division  for medical students.

They are a http://www.psr.org/spsr/ branch of http://www.psr.org/ Physicians for Social Responsibility.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE LINK5

! LINK5 takes the L-set displacements solved for in LINK3 (statics) or LINK4 (eigenvalues) and builds it back up to the G-set.
! See Appendix B to the MYSTRAN User's Reference Guide for an explanation of how this is done.

! In addition, for Craig-Bampton model generation (SOL = GEN CB MODEL or 31), array PHIXA is expanded to G-set size

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_BUG, WRT_ERR, WRT_LOG, ERR, F04, F06, L1H, L2A, L2E, L2F, L3A, L5A, L5B, SC1
      USE IOUNT1, ONLY                :  LINK1H, LINK2A, LINK2E, LINK2F, LINK3A, LINK5A, LINK5B
      USE IOUNT1, ONLY                :  L1H_MSG, L2A_MSG, L2E_MSG, L2F_MSG, L3A_MSG, L5A_MSG, L5B_MSG
      USE IOUNT1, ONLY                :  ERRSTAT, L1HSTAT, L2ESTAT, L2FSTAT, L3ASTAT
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, COMM, FATAL_ERR, LINKNO, MBUG, NDOFA, NDOFF, NDOFG, NDOFL, NDOFM,           &
                                         NDOFN, NDOFO, NDOFR, NDOFS, NDOFSE, NGRID, NSUB, NTERM_GMN, NTERM_GOA, NTERM_PO,          &
                                         NUM_CB_DOFS, NUM_EIGENS, NVEC, SOL_NAME, WARN_ERR
      USE TIMDAT, ONLY                :  YEAR, MONTH, DAY, HOUR, MINUTE, SEC, SFRAC, STIME, TSEC
      USE CONSTANTS_1, ONLY           :  ZERO, ONE
      USE PARAMS, ONLY                :  EIGNORM2, SUPINFO, SUPWARN
      USE NONLINEAR_PARAMS, ONLY      :  LOAD_ISTEP
      USE EIGEN_MATRICES_1 , ONLY     :  EIGEN_VAL, EIGEN_VEC, GEN_MASS, MODE_NUM
      USE FULL_MATRICES, ONLY         :  PHIZG_FULL
      USE SPARSE_MATRICES, ONLY       :  I_GMN, J_GMN, GMN, I_GOA, J_GOA, GOA
      USE OUTPUT4_MATRICES, ONLY      :  NUM_OU4_REQUESTS
      USE MISC_MATRICES, ONLY         :  UG_T123_MAT
      USE COL_VECS, ONLY              :  UG_COL, YSe, UO0_COL, UL_COL
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG
      USE DOF_TABLES, ONLY            :  TDOF, TDOFI
      USE MODEL_STUF, ONLY            :  GRID, GRID_ID, INV_GRID_SEQ, EIG_COMP, EIG_GRID, EIG_NORM, MAXMIJ, MIJ_COL, MIJ_ROW

      USE LINK5_USE_IFs

      IMPLICIT NONE

      LOGICAL                         :: VEC_SIGN_CHG(NDOFL) ! Indicators of whether user wants to change sign of an eigenvector

      CHARACTER, PARAMETER            :: CR13 = CHAR(13)   ! This causes a carriage return simulating the "+" action in a FORMAT
      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'LINK5'
      CHARACTER( 1*BYTE)              :: CLOSE_IT          ! Input to subr READ_MATRIX_i. 'Y'/'N' whether to close a file or not
      CHARACTER( 8*BYTE)              :: CLOSE_STAT        ! What to do with file when it is closed
      CHARACTER( 1*BYTE)              :: DO_IT             ! If 'Y' execute some code
      CHARACTER( 1*BYTE)              :: MIJ_COL_FOUND='N' ! 'Y' if MIJ_ROW is processed as a solution vector in this LINK
      CHARACTER( 1*BYTE)              :: MIJ_ROW_FOUND='N' ! 'Y' if MIJ_ROW is processed as a solution vector in this LINK
      CHARACTER(54*BYTE)              :: MODNAM            ! Name to write to screen to describe module being run
      CHARACTER( 1*BYTE)              :: READ_NTERM        ! 'Y' or 'N' Input to subr READ_MATRIX_1
      CHARACTER( 1*BYTE)              :: OPND              ! Input to subr READ_MATRIX_i. 'Y'/'N' whether to open  a file or not
      CHARACTER( 1*BYTE)              :: READ_UO0          ! If 'Y' then read UO0 data from file L2F

      INTEGER(LONG)                   :: COL_NUM           ! Arg passed to subr BUILD_A_LR
      INTEGER(LONG)                   :: EIG_NORM_GSET_DOF ! A-set DOF no. for EIG_GRID/EIG_COMP
      INTEGER(LONG)                   :: EIGNORM2_ERR        ! Error indicator for reading param EIGNORM2 data
      INTEGER(LONG)                   :: G_SET_COL         ! Col number in TDOF, TDOFI where G-set DOF's exist
      INTEGER(LONG)                   :: I,J,K,L           ! DO loop indices
      INTEGER(LONG)                   :: IERROR            ! Error count
      INTEGER(LONG)                   :: IGRID             ! Internal grid numbER for EIG_GRID
      INTEGER(LONG)                   :: IOCHK             ! IOSTAT error number when opening/reading a file
      INTEGER(LONG)                   :: NUM_COMPS         ! 6 if GRID_NUM is an physical grid, 1 if an SPOINT
      INTEGER(LONG)                   :: NUM_SOLNS    = 0  ! No. of solutions to process (e.g. NSUB for STATICS)
      INTEGER(LONG)                   :: OUNT(2)           ! File units to write messages to. Input to subr UNFORMATTED_OPEN
      INTEGER(LONG)                   :: P_LINKNO          ! Prior LINK no's that should have run before this LINK can execute
      INTEGER(LONG)                   :: REC_NO            ! Record number when reading a file

      REAL(DOUBLE)                    :: MIJ_COL_SCALE=ZERO! Scale fac for a col of gen mass matrix to renorm MAXMIJ from LINK4
      REAL(DOUBLE)                    :: MIJ_ROW_SCALE=ZERO! Scale fac for a col of gen mass matrix to renorm MAXMIJ from LINK4
      REAL(DOUBLE)                    :: PHI_SCALE_FAC     ! Scale factor that the eigenvector was renormalized to in subr RENORM

! **********************************************************************************************************************************
      LINKNO = 5

! Set time initializing parameters

      CALL TIME_INIT

! Initialize WRT_BUG

      DO I=0,MBUG-1
         WRT_BUG(I) = 0
      ENDDO

! Get date and time, write to screen

      CALL OURDAT
      CALL OURTIM
      WRITE(SC1,152) LINKNO

! Make units for writing errors the screen until we open output files

      OUNT(1) = SC1
      OUNT(2) = SC1

! Make units for writing errors the error file and output file

      OUNT(1) = ERR
      OUNT(2) = F06

! Write info to text files

      WRITE(F06,150) LINKNO
      IF (WRT_LOG > 0) THEN
         WRITE(F04,150) LINKNO
      ENDIF
      WRITE(ERR,150) LINKNO

! Read LINK1A file

      CALL READ_L1A ( 'KEEP', 'Y' )
! Check COMM for successful completion of prior LINKs

      IF (SOL_NAME(1:7) == 'STATICS') THEN

         P_LINKNO = 3
         IF (COMM(P_LINKNO) /= 'C') THEN
            WRITE(SC1,9998) P_LINKNO,P_LINKNO,LINKNO
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF

      ELSE IF (SOL_NAME(1:5) == 'MODES') THEN

         P_LINKNO = 4
         IF (COMM(P_LINKNO) /= 'C') THEN
            WRITE(ERR,9998) P_LINKNO,P_LINKNO,LINKNO
            WRITE(F06,9998) P_LINKNO,P_LINKNO,LINKNO
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF

      ELSE IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN

         P_LINKNO = 6
         IF (COMM(P_LINKNO) /= 'C') THEN
            WRITE(ERR,9998) P_LINKNO,P_LINKNO,LINKNO
            WRITE(F06,9998) P_LINKNO,P_LINKNO,LINKNO
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF

      ELSE IF (SOL_NAME(1:8) == 'BUCKLING') THEN

         IF      (LOAD_ISTEP == 1) THEN
            P_LINKNO = 3
            IF (COMM(P_LINKNO) /= 'C') THEN
               WRITE(ERR,9998) P_LINKNO,P_LINKNO,LINKNO
               WRITE(F06,9998) P_LINKNO,P_LINKNO,LINKNO
               FATAL_ERR = FATAL_ERR + 1
               CALL OUTA_HERE ( 'Y' )
            ENDIF
         ELSE IF (LOAD_ISTEP == 2) THEN
            P_LINKNO = 4
            IF (COMM(P_LINKNO) /= 'C') THEN
               WRITE(ERR,9998) P_LINKNO,P_LINKNO,LINKNO
               WRITE(F06,9998) P_LINKNO,P_LINKNO,LINKNO
               FATAL_ERR = FATAL_ERR + 1
               CALL OUTA_HERE ( 'Y' )
            ENDIF
         ELSE
            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,5002) SUBR_NAME, LOAD_ISTEP
            WRITE(F06,5002) SUBR_NAME, LOAD_ISTEP
            CALL OUTA_HERE ( 'Y' )
         ENDIF

      ELSE IF (SOL_NAME(1:8) == 'NLSTATIC') THEN

         P_LINKNO = 3
         IF (COMM(P_LINKNO) /= 'C') THEN
            WRITE(ERR,9998) P_LINKNO,P_LINKNO,LINKNO
            WRITE(F06,9998) P_LINKNO,P_LINKNO,LINKNO
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF

      ELSE

         WRITE(ERR,5001) SOL_NAME
         WRITE(F06,5001) SOL_NAME
         FATAL_ERR = FATAL_ERR + 1
         CALL OUTA_HERE ( 'Y' )

      ENDIF

! **********************************************************************************************************************************
! Allocate arrays for YSe, GMN, GOA

      CALL OURTIM
      MODNAM = 'ALLOCATE SEVERAL ARRAYS'
      WRITE(SC1,5092) LINKNO,MODNAM,HOUR,MINUTE,SEC,SFRAC
      CALL ALLOCATE_SPARSE_MAT ( 'GMN', NDOFM, NTERM_GMN, SUBR_NAME )
      CALL ALLOCATE_SPARSE_MAT ( 'GOA', NDOFO, NTERM_GOA, SUBR_NAME )
      CALL ALLOCATE_COL_VEC ( 'YSe' , NDOFS, SUBR_NAME )

! Read GMN matrix if there are MPC's

      IF (NTERM_GMN > 0) THEN

         CALL OURTIM
         MODNAM = 'READ GMN MATRIX'
         WRITE(SC1,5092) LINKNO,MODNAM,HOUR,MINUTE,SEC,SFRAC
         READ_NTERM = 'Y'
         OPND       = 'N'
         CLOSE_IT   = 'Y'
         CALL READ_MATRIX_1 ( LINK2A, L2A, OPND, CLOSE_IT, 'KEEP', L2A_MSG, 'GMN', NTERM_GMN, READ_NTERM, NDOFM                    &
                            , I_GMN, J_GMN, GMN )
      ENDIF

! Read GOA matrix if there are omitted DOFs

      IF (NTERM_GOA > 0) THEN
         CALL OURTIM
         MODNAM = 'READ GOA MATRIX'
         WRITE(SC1,5092) LINKNO,MODNAM,HOUR,MINUTE,SEC,SFRAC
         READ_NTERM = 'Y'
         OPND       = 'N'
         CLOSE_IT   = 'Y'
!        CLOSE_STAT = L2ESTAT
         CLOSE_STAT = 'KEEP'
         CALL READ_MATRIX_1 ( LINK2E, L2E, OPND, CLOSE_IT, CLOSE_STAT, L2E_MSG, 'GOA', NTERM_GOA, READ_NTERM, NDOFO                &
                            , I_GOA, J_GOA, GOA )
      ENDIF

! Read enforced displ's

      IF (NDOFSE > 0) THEN

         IF ((SOL_NAME(1:7) == 'STATICS') .OR. (SOL_NAME(1:8) == 'NLSTATIC') .OR.                                                  &
            ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 1))) THEN

            CALL FILE_OPEN ( L1H, LINK1H, OUNT, 'OLD', L1H_MSG, 'READ_STIME', 'UNFORMATTED', 'READ', 'REWIND', 'Y', 'N', 'Y' )

            CALL OURTIM
            MODNAM = 'READ YSe ENFORCED DISPLACEMENTS'
            WRITE(SC1,5092) LINKNO,MODNAM,HOUR,MINUTE,SEC,SFRAC

            IERROR = 0
            DO I=1,NDOFSE
               READ(L1H,IOSTAT=IOCHK) YSe(I)
               IF (IOCHK /= 0) THEN
                  IERROR = IERROR + 1
                  REC_NO = I
                  CALL READERR ( IOCHK, LINK1H, L1H_MSG, REC_NO, OUNT, 'Y' )
               ENDIF
            ENDDO
            IF (IERROR /= 0) THEN
               WRITE(ERR,9995) LINKNO,IERROR
               WRITE(F06,9995) LINKNO,IERROR
               CALL OUTA_HERE ( 'Y' )
            ENDIF

            CALL FILE_CLOSE ( L1H, LINK1H, L1HSTAT, 'Y' )

         ENDIF

      ENDIF

! Read eigenvalue data from file L1M if this is an eigenvalue problem. Check to see if user wants to change sign of any vector.
! EIGEN_VAL was not deallocated in LINK4 (see LINK4 comment 01/11/19) so we do not allocate it here anymore

      DO_IT = 'N'
      IF ((SOL_NAME(1:5) == 'MODES') .OR. (SOL_NAME(1:12) == 'GEN CB MODEL')) THEN
         DO_IT = 'Y'
      ENDIF
      IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN
         DO_IT = 'Y'
      ENDIF

      IF (DO_IT == 'Y') THEN
         IERROR = 0
         CALL ALLOCATE_EIGEN1_MAT ( 'MODE_NUM' , NUM_EIGENS, 1, SUBR_NAME )
!xx      CALL ALLOCATE_EIGEN1_MAT ( 'EIGEN_VAL', NUM_EIGENS, 1, SUBR_NAME )
         CALL ALLOCATE_EIGEN1_MAT ( 'GEN_MASS' , NUM_EIGENS, 1, SUBR_NAME )
         CALL READ_L1M ( IERROR )

         IF (IERROR /= 0) THEN
            WRITE(ERR,9995) LINKNO,IERROR
            WRITE(F06,9995) LINKNO,IERROR
            CALL OUTA_HERE ( 'Y' )
         ENDIF

      ENDIF

! Open file that has L-set displs, or eigenvectors ('MODES') or PHIZL ('GEN CB MODEL')

      IF (NDOFL > 0) THEN
         CALL FILE_OPEN ( L3A, LINK3A, OUNT, 'OLD', L3A_MSG, 'READ_STIME', 'UNFORMATTED', 'READ', 'REWIND', 'Y', 'N', 'Y' )
      ENDIF

! Open file for writing displs to.

      CALL FILE_CLOSE ( L5A, LINK5A, 'KEEP', 'Y' )
      CALL FILE_OPEN  ( L5A, LINK5A, OUNT, 'REPLACE', L5A_MSG, 'WRITE_STIME', 'UNFORMATTED', 'WRITE', 'REWIND', 'Y', 'N', 'Y' )

! Open file that has UO0

      IF (NTERM_PO > 0) THEN
         CALL FILE_OPEN ( L2F, LINK2F, OUNT, 'OLD', L2F_MSG, 'READ_STIME', 'UNFORMATTED', 'READ', 'REWIND', 'Y', 'N', 'Y' )
      ENDIF

! Set NUM_SOLNS for use in loop (below) to get outputs for each subcase/solution vector

      IF      ((SOL_NAME(1:7) == 'STATICS') .OR. (SOL_NAME(1:8) == 'NLSTATIC')) THEN
         NUM_SOLNS = NSUB
      ELSE IF (SOL_NAME(1:5) == 'MODES') THEN
         NUM_SOLNS = NVEC
      ELSE IF (SOL_NAME(1:8) == 'BUCKLING') THEN
         IF (LOAD_ISTEP == 1) THEN
            NUM_SOLNS = 1
         ELSE IF (LOAD_ISTEP == 2) THEN
            NUM_SOLNS = NVEC
         ENDIF
      ELSE IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
         NUM_SOLNS = NUM_CB_DOFS
      ENDIF

! Allocate memory to EIGEN_VEC so that we can write the G-set eigenvectors to it for possible OUTPUT4 later

      IF (SOL_NAME(1:5) == 'MODES') THEN
         CALL ALLOCATE_EIGEN1_MAT ( 'EIGEN_VEC', NDOFG, NUM_SOLNS, SUBR_NAME )
      ENDIF

! Allocate memory for CB matrix PHIZG_FULL

      IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
         CALL ALLOCATE_FULL_MAT ( 'PHIZG_FULL', NDOFG, NUM_CB_DOFS, SUBR_NAME )
      ENDIF

! Check for renorm = MAX or POINT.

      DO_IT = 'N'
      IF ((SOL_NAME(1:5) == 'MODES') .OR. (SOL_NAME(1:12) == 'GEN CB MODEL')) THEN
         DO_IT = 'Y'
      ENDIF
      IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN
         DO_IT = 'Y'
      ENDIF

      IF (DO_IT == 'Y') THEN
         EIG_NORM_GSET_DOF = 0
         IGRID = 0
         IF (EIG_NORM == 'POINT   ') THEN                  ! User requested to renormalize eigenvectors on POINT

            CALL TDOF_COL_NUM ( 'G ',  G_SET_COL )
i_do:       DO I=1,NDOFG
               IF (TDOF(I,1) == EIG_GRID) THEN
                  IGRID = TDOF(I,3)
                  EIG_NORM_GSET_DOF = TDOF(I,G_SET_COL) + EIG_COMP - 1
                  EXIT i_do
               ENDIF
            ENDDO i_do

            IF (IGRID > 0) THEN
               WRITE(ERR,5102) EIG_GRID, EIG_COMP
               IF (SUPINFO == 'N') THEN
                  WRITE(F06,5102) EIG_GRID, EIG_COMP
               ENDIF
            ELSE
               WARN_ERR = WARN_ERR + 1
               WRITE(ERR,5111) EIG_GRID
               IF (SUPWARN == 'N') THEN
                  WRITE(F06,5111) EIG_GRID
               ENDIF
               PHI_SCALE_FAC = ONE
            ENDIF

         ELSE IF (EIG_NORM == 'MAX     ') THEN             ! User requested to renormalize eigenvectors on MAX

            WRITE(ERR,5103)
            IF (SUPINFO == 'N') THEN
               WRITE(F06,5103)
            ENDIF

         ELSE                                              ! No renorm needed here (was done in LINK4 if not POINT or MASS)

            WRITE(ERR,5104)
            IF (SUPINFO == 'N') THEN
               WRITE(F06,5104)
            ENDIF

         ENDIF

      ENDIF

      IF ((EIGNORM2 == 'Y') .AND. (NVEC > 0)) THEN
         CALL READ_EIGNORM2
      ENDIF

! **********************************************************************************************************************************
      READ_UO0 = 'Y'
      IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN
         READ_UO0 = 'N'
      ENDIF

! Begin loop for reading L-set displs and building up to G-set displs one subcase/solution vector at a time

j_do: DO J = 1,NUM_SOLNS

         CALL ALLOCATE_COL_VEC ('UL_COL', NDOFL, SUBR_NAME)! Allocate array UL_COL

         CALL OURTIM                                       ! Read UL displs for the current subcase/vector from LINK3A
         IF     ((SOL_NAME(1: 7) == 'STATICS') .OR. (SOL_NAME(1:8) == 'NLSTATIC')) THEN
            MODNAM = 'READ  L-SET DISPLACEMENTS                      Subcase'
         ELSE IF (SOL_NAME(1: 5) == 'MODES') THEN
            MODNAM = 'READ  L-SET EIGENVECTORS                       Vector'
         ELSE IF (SOL_NAME(1: 8) == 'BUCKLING') THEN
            IF (LOAD_ISTEP == 1) THEN
               MODNAM = 'READ  L-SET DISPLACEMENTS                      Subcase'
            ELSE IF (LOAD_ISTEP == 2) THEN
               MODNAM = 'READ  L-SET EIGENVECTORS                       Vector'
            ENDIF
         ELSE IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
            MODNAM = 'READ  L-SET CB VECTORS (PHIZL)                 CB vec'
         ENDIF
         WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC

         REC_NO = 0
         IERROR = 0
         DO I=1,NDOFL
            REC_NO = REC_NO + 1
            READ(L3A,IOSTAT=IOCHK) UL_COL(I)               ! For CB, a col of PHIZL. So UG_COL, calc'd in this subr, is a PHIZG col
            IF (IOCHK /= 0) THEN
               CALL READERR ( IOCHK, LINK3A, L3A_MSG, REC_NO, OUNT, 'Y' )
               IERROR = IERROR + 1
            ENDIF
         ENDDO

         IF (IERROR /= 0) THEN
            WRITE(ERR,9995) LINKNO,IERROR
            WRITE(F06,9995) LINKNO,IERROR
            CALL OUTA_HERE ( 'Y' )
         ENDIF
                                                           ! Build UA from UL and UR
         CALL ALLOCATE_COL_VEC ( 'UA_COL', NDOFA, SUBR_NAME )
         CALL ALLOCATE_COL_VEC ( 'UR_COL', NDOFR, SUBR_NAME )
         CALL OURTIM
         MODNAM = 'BUILD UA DISPLS FROM UL, UR:                      "'
         WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC
         COL_NUM = 0
         IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
            IF ((J > NDOFR+nvec) .AND. (J <= NUM_CB_DOFS)) THEN
               COL_NUM = J
            ENDIF
         ENDIF
         CALL BUILD_A_LR ( COL_NUM )
         CALL DEALLOCATE_COL_VEC ( 'UL_COL' )
         CALL DEALLOCATE_COL_VEC ( 'UR_COL' )

                                                           ! Solve for UO and Build UF from UA and UO
         CALL ALLOCATE_COL_VEC ( 'UF_COL' , NDOFF, SUBR_NAME )
         CALL ALLOCATE_COL_VEC ( 'UO_COL' , NDOFO, SUBR_NAME )
         CALL ALLOCATE_COL_VEC ( 'UO0_COL', NDOFO, SUBR_NAME )
         CALL OURTIM
         MODNAM = 'BUILD UF DISPLS FROM UA, UO:                      "'
         WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC
         IF (READ_UO0 == 'Y') THEN
            IF (NDOFO > 0) THEN
               IF (NTERM_PO > 0) THEN
                  CALL OURTIM
                  MODNAM = '  READ UO0 DISPLS,                                "'
                  WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC

                  IERROR = 0
                  DO I=1,NDOFO
                     READ(L2F,IOSTAT=IOCHK) UO0_COL(I)
                     IF (IOCHK /= 0) THEN
                        REC_NO = I+1
                        CALL READERR ( IOCHK, LINK2F, L2F_MSG, REC_NO, OUNT, 'Y' )
                        IERROR = IERROR + 1
                     ENDIF
                  ENDDO
                  IF (IERROR /= 0) THEN
                     WRITE(ERR,9995) LINKNO,IERROR
                     WRITE(F06,9995) LINKNO,IERROR
                     CALL OUTA_HERE ( 'Y' )
                  ENDIF
               ELSE
                  DO I=1,NDOFO
                     UO0_COL(I) = ZERO
                  ENDDO
               ENDIF
            ENDIF
         ENDIF
         CALL BUILD_F_AO
         CALL DEALLOCATE_COL_VEC ( 'UA_COL' )
         CALL DEALLOCATE_COL_VEC ( 'UO_COL' )
         CALL DEALLOCATE_COL_VEC ( 'UO0_COL' )
                                                           ! Build UN from UF and US
         CALL ALLOCATE_COL_VEC ( 'UN_COL', NDOFN , SUBR_NAME)
         CALL ALLOCATE_COL_VEC ( 'US_COL', NDOFS, SUBR_NAME )
         CALL OURTIM
         MODNAM = 'BUILD UN DISPLS FROM UF, US:                      "'
         WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC
         CALL BUILD_N_FS
         CALL DEALLOCATE_COL_VEC ( 'UF_COL' )
         CALL DEALLOCATE_COL_VEC ( 'US_COL' )
                                                           ! Build UG from UN and UM
         CALL DEALLOCATE_COL_VEC ( 'UG_COL' )
         CALL ALLOCATE_COL_VEC ( 'UG_COL', NDOFG, SUBR_NAME )
         CALL ALLOCATE_COL_VEC ( 'UM_COL', NDOFM, SUBR_NAME )
         CALL OURTIM
         MODNAM = 'BUILD UG DISPLS FROM UN, UM:                      "'
         WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC
         CALL BUILD_G_NM
         CALL DEALLOCATE_COL_VEC ( 'UN_COL' )
         CALL DEALLOCATE_COL_VEC ( 'UM_COL' )

         IERROR = 0
!xx      IF ((SOL_NAME(1:8) == 'DIFFEREN') .OR. (SOL_NAME(1:8) == 'NLSTATIC')) THEN
            IF (ALLOCATED(UG_T123_MAT)) THEN
               CALL DEALLOCATE_MISC_MAT ( 'UG_T123_MAT' )
            ENDIF
            CALL ALLOCATE_MISC_MAT ( 'UG_T123_MAT', NGRID, 3, SUBR_NAME )
            CALL GET_UG_123_IN_GRD_ORD ( IERROR )
!xx      ENDIF
         IF (IERROR /= 0) THEN
            WRITE(ERR,9995) LINKNO,IERROR
            WRITE(F06,9995) LINKNO,IERROR
            CALL OUTA_HERE ( 'Y' )
         ENDIF

         IF      (SOL_NAME(1: 5) == 'MODES'       ) THEN   ! For modes, write UG_COL to jth col of EIGEN_VEC (for later OUTPUT4)
            DO I=1,NDOFG
               EIGEN_VEC(I,J)  = UG_COL(I)
            ENDDO
         ELSE IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN   ! For CB   , write UG_COL to jth col of PHiZG     (for later OUTPUT4)
            DO I=1,NDOFG
               PHIZG_FULL(I,J) = UG_COL(I)
            ENDDO
         ENDIF
                                                           ! Renorm eigenvecs, if requested
         IF (J <= NVEC) THEN

            DO_IT = 'N'
            IF ((SOL_NAME(1:5) == 'MODES') .OR. (SOL_NAME(1:12) == 'GEN CB MODEL')) THEN
               DO_IT = 'Y'
            ENDIF
            IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN
               DO_IT = 'Y'
            ENDIF

            IF (DO_IT == 'Y') THEN
               IF ((EIG_NORM == 'POINT   ') .OR. (EIG_NORM == 'MAX     ')) THEN
                  CALL RENORM (J, EIG_GRID, EIG_COMP, EIG_NORM, EIG_NORM_GSET_DOF, GEN_MASS(J), PHI_SCALE_FAC )
                  IF      (J == MIJ_ROW) THEN
                     MIJ_ROW_FOUND = 'Y'
                     MIJ_ROW_SCALE = PHI_SCALE_FAC
                  ELSE IF (J == MIJ_COL) THEN
                     MIJ_COL_FOUND = 'Y'
                     MIJ_COL_SCALE = PHI_SCALE_FAC
                  ENDIF
                  CALL WRITE_L1M                           ! Need to update GEN_MASS in L1M if vecs are renorm'd here
               ENDIF

            ENDIF

         ENDIF
                                                           ! See if user requested eigenvector sign changes
         IF ((EIGNORM2 == 'Y') .AND. (EIGNORM2_ERR == 0)) THEN
            IF (VEC_SIGN_CHG(J)) THEN
               DO I=1,NDOFG
                  UG_COL(I) = -UG_COL(I)
               ENDDO
            ENDIF
         ENDIF

         CALL OURTIM                                       ! Write UG displs for this subcase to file LINK5A
         MODNAM = 'WRITE UG DISPLS TO FILE,                          "'
         WRITE(SC1,5093) LINKNO,MODNAM,J,HOUR,MINUTE,SEC,SFRAC
         WRITE(SC1, * )                                    ! Separator between UG_COL calcs
         DO I=1,NDOFG
            WRITE(L5A) UG_COL(I)                           ! For CB this is a col of PHIZG (which is never processed as an array)
         ENDDO

      ENDDO j_do                                           ! End of loop on NUM_SOLNS

! If CB soln, expand PHIXA to G-set size and write to file unit L5B

     IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
                                                           ! Open file for writing cols of PHIXG
         CALL FILE_OPEN ( L5B, LINK5B, OUNT, 'REPLACE', L5B_MSG, 'WRITE_STIME', 'UNFORMATTED', 'WRITE', 'REWIND', 'Y', 'N', 'Y' )

         CALL DEALLOCATE_COL_VEC ( 'UG_COL' )
         CALL EXPAND_PHIXA_TO_PHIXG                        ! Expand PHIXA to PHIXG and write cols to file L5B
   !xx   WRITE(SC1, * )                                    ! Advance 1 line for screen messages
         WRITE(SC1,12345,ADVANCE='NO') '       Deallocate PHIXA', CR13   ;   CALL DEALLOCATE_SPARSE_MAT ( 'PHIXA' )

      ENDIF

! Code to write PHIZG in full format

      IF ((DEBUG(55) == 2) .OR. (DEBUG(55) == 3)) THEN
         IF (SOL_NAME == 'GEN CB MODEL') THEN
            WRITE(F06,*)
            WRITE(F06,99885)
            WRITE(F06,99886) (J,J=1,NUM_CB_DOFS)
            L = 0
            DO I=1,NGRID
               CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
               DO K=1,NUM_COMPS
                  L = L + 1
                  IF (K == 1) THEN
                     WRITE(F06,99887) TDOFI(L,1), K, (PHIZG_FULL(L,J),J=1,NUM_CB_DOFS)
                  ELSE
                     WRITE(F06,99888)             K, (PHIZG_FULL(L,J),J=1,NUM_CB_DOFS)
                  ENDIF
               ENDDO
               WRITE(F06,*)
            ENDDO
            WRITE(F06,*)
         ENDIF
      ENDIF

! **********************************************************************************************************************************
! Write eigen analysis summary if we have renormed vectors here in LINK5

      DO_IT = 'N'
      IF ((SOL_NAME(1:5) == 'MODES') .OR. (SOL_NAME(1:12) == 'GEN CB MODEL')) THEN
         DO_IT = 'Y'
      ENDIF
      IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN
         DO_IT = 'Y'
      ENDIF

      IF (DO_IT == 'Y') THEN
         IF ((EIG_NORM == 'POINT   ') .OR. (EIG_NORM == 'MAX     ')) THEN

            IF ((MIJ_ROW_FOUND == 'Y') .AND. (MIJ_COL_FOUND == 'Y')) THEN
               MAXMIJ = MAXMIJ/(MIJ_ROW_SCALE * MIJ_COL_SCALE)
            ELSE
               WARN_ERR = WARN_ERR + 1
               WRITE(ERR,5101) LINKNO
               IF (SUPWARN == 'N') THEN
                  WRITE(F06,5101) LINKNO
               ENDIF
            ENDIF

            CALL EIG_SUMMARY
         ENDIF
      ENDIF

! Close files

      IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 1)) THEN
         CLOSE_STAT = 'KEEP'
      ELSE
!xx      CLOSE_STAT = L2FSTAT
         CLOSE_STAT = 'KEEP'
      ENDIF
      CALL FILE_CLOSE ( L2F, LINK2F, CLOSE_STAT, 'Y' )

      IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
         CALL FILE_CLOSE ( L3A, LINK3A, 'KEEP', 'Y' )
      ELSE
         CALL FILE_CLOSE ( L3A, LINK3A, L3ASTAT, 'Y' )
      ENDIF

      CALL FILE_CLOSE ( L5A, LINK5A, 'KEEP', 'Y' )
      CALL FILE_CLOSE ( L5B, LINK5B, 'KEEP', 'Y' )

! Call OUTPUT4 processor to process output requests for OUTPUT4 matrices generated in this link

      IF (NUM_OU4_REQUESTS > 0) THEN
         CALL OURTIM
         MODNAM = 'WRITE OUTPUT4 NATRICES      '
         WRITE(SC1,5092) LINKNO,MODNAM,HOUR,MINUTE,SEC,SFRAC
         WRITE(F06,*)
         CALL OUTPUT4_PROC ( SUBR_NAME )
      ENDIF

! Deallocate arrays (except EIGEN_VAL, may be needed later)

!xx   WRITE(SC1, * ) '     DEALLOCATE SOME ARRAYS'
!xx   WRITE(SC1, * )                                       ! Advance 1 line for screen messages
      WRITE(SC1,12345,ADVANCE='NO') '       Deallocate GMN      ', CR13  ;   CALL DEALLOCATE_SPARSE_MAT ( 'GMN' )
      WRITE(SC1,12345,ADVANCE='NO') '       Deallocate GOA      ', CR13  ;   CALL DEALLOCATE_SPARSE_MAT ( 'GOA' )
      WRITE(SC1,12345,ADVANCE='NO') '       Deallocate MODE_NUM ', CR13  ;   CALL DEALLOCATE_EIGEN1_MAT ( 'MODE_NUM' )
!xx   WRITE(SC1,12345,ADVANCE='NO') '       Deallocate EIGEN_VAL', CR13  ;   CALL DEALLOCATE_EIGEN1_MAT ( 'EIGEN_VAL' )
      WRITE(SC1,12345,ADVANCE='NO') '       Deallocate GEN_MASS ', CR13  ;   CALL DEALLOCATE_EIGEN1_MAT ( 'GEN_MASS' )
      CALL DEALLOCATE_COL_VEC    ( 'YSe' )
      CALL DEALLOCATE_EIGEN1_MAT ( 'EIGEN_VEC' )

! Process is now complete so set COMM(LINKNO)

      COMM(LINKNO) = 'C'

! Write data to L1A
      CALL WRITE_L1A ( 'KEEP', 'Y', 'Y' )

! Check allocation status of allocatable arrays, if requested

      IF (DEBUG(100) > 0) THEN
         CALL CHK_ARRAY_ALLOC_STAT
         IF (DEBUG(100) > 1) THEN
            CALL WRITE_ALLOC_MEM_TABLE ( 'at the end of '//SUBR_NAME )
         ENDIF
      ENDIF

! Write LINK5 end to F04, F06

      CALL OURTIM
      IF (WRT_LOG > 0) THEN
         WRITE(F04,151) LINKNO
      ENDIF
      WRITE(F06,151) LINKNO

! Close files

      IF (( DEBUG(193) == 5) .OR. (DEBUG(193) == 999)) THEN
         CALL FILE_INQUIRE ( 'near end of LINK5' )
      ENDIF

! Write LINK5 end to screen
      WRITE(SC1,153) LINKNO

! **********************************************************************************************************************************
  101 format(32767(1es22.14))

  150 FORMAT(/,' >> LINK',I3,' BEGIN',/)

  151 FORMAT(/,' >> LINK',I3,' END',/)

  152 FORMAT(/,' >> LINK',I3,' BEGIN')

  153 FORMAT(  ' >> LINK',I3,' END')

 5001 FORMAT(' *ERROR  5001: INVALID SOLUTION NAME: SOL = ',A)

 5002 FORMAT(' *ERROR  5002: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                     ,/,14X,'VARIABLE LOAD_ISTEP MUST BE 1 OR 2 BUT VALUE IS = ',I8)

 5092 FORMAT(1X,I2,'/',A54,8X,2X,I2,':',I2,':',I2,'.',I3)

 5093 FORMAT(1X,I2,'/',A54,I8,2X,I2,':',I2,':',I2,'.',I3)

 5094 FORMAT(/,' >> LINK',I2,' END',19X,I2,':',I2,':',I2,'.',I3,/)

 5101 FORMAT(' *WARNING    : THE LARGEST OFF-DIAGONAL GENERALIZED MASS TERM REPORTED IN THE EIGENVALUE ANALYSIS SUMMARY CANNOT BE' &
                    ,/,14X,' RESCALED IN LINK ',I2,' BASED ON THE EIGENVALUE RENORMALIZATION REQUESTED.'                           &
                    ,/,14X,' IT IS RENORMALIZED BASED ON UNIT GENERALIZED MASS',/)

 5102 FORMAT(' *INFORMATION: ATTEMPTING TO RENORMALIZE EIGENVECTORS BASED ON GRID POINT-COMPONENT ',2I8                            &
                    ,/,14X,' THE RENORMALIZATION WILL BE DONE FOR ALL EIGENVECTORS WHOSE AMPLITUDE IS NONZERO AT THIS DOF',/)

 5103 FORMAT(' *INFORMATION: EIGENVECTORS WILL BE RENORMALIZED BASED ON MAX VALUE = 1.0',/)

 5104 FORMAT(' *INFORMATION: EIGENVECTORS WERE NORMALIZED TO MASS IN LINK4',/)

 5111 FORMAT(' *WARNING    : REQUEST TO NORMALIZE EIGENVECTORS BASED ON GRID POINT ',I8,', HOWEVER, THIS NOT A VALID GRID POINT.'  &
                    ,/,14X,' EIGENVECTORS WILL NOT BE RENORMALIZED',/)

 9995 FORMAT(/,' PROCESSING ENDED IN LINK ',I3,' DUE TO ABOVE ',I8,' ERRORS')

 9998 FORMAT(' *ERROR  9998: COMM ',I3,' INDICATES UNSUCCESSFUL LINK ',I2,' COMPLETION.'                                           &
                    ,/,14X,' FATAL ERROR - CANNOT START LINK ',I2)

98001 FORMAT(41X,'R E A L   E I G E N V A L U E S')

98002 FORMAT(82X,'BEFORE RENORMALIZATION OF EIGENVECTORS')

98003 FORMAT(3X,' MODE  EXTRACTION     EIGENVALUE            RADIANS              CYCLES            GENERALIZED         GENERALIZED&
&        ',/,3X,'NUMBER   ORDER                                                                        MASS              STIFFNESS'&
          ,/)

98004 FORMAT(1X,2I8,5(1ES20.6))

99885 FORMAT(82X,'MATRIX PHIZG',/,82X,'------------')

99886 FORMAT(5X,32676(I14))

99887 FORMAT(I8,'-',I1,32767(1ES14.6))

99888 FORMAT(8X,'-',I1,32767(1ES14.6))

12345 FORMAT(A,10X,A)

! ##################################################################################################################################

      CONTAINS

! ##################################################################################################################################

      SUBROUTINE READ_EIGNORM2

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  EIN, EINFIL, ERR, F06
      USE SCONTR, ONLY                :  IERRFL, JCARD_LEN

      IMPLICIT NONE

      LOGICAL                         :: FILE_EXIST

      CHARACTER(JCARD_LEN*BYTE)       :: DATA_FIELD
      CHARACTER(80*BYTE)              :: TITLE             ! First record in EINFIL

      INTEGER(LONG)                   :: IOCHK             ! Vaue of IOSTAT in file open
      INTEGER(LONG)                   :: NUM_CHANGED       ! Number of eigenvectors to have their sign changed
      INTEGER(LONG)                   :: REC_NUM           ! Number of the record read from file
      INTEGER(LONG)                   :: VEC_NUM           ! Number of a eigenvector read from EINFIL
      INTEGER(LONG)                   :: VECS_CHANGED(NVEC)! Numbers of eigenvectors to have their sign changed

! **********************************************************************************************************************************
      DO I=1,NDOFL
         VEC_SIGN_CHG(I) = .FALSE.
      ENDDO

      EIGNORM2_ERR = 0
      INQUIRE ( FILE=EINFIL, EXIST=FILE_EXIST )

      IF (FILE_EXIST) THEN

         EIGNORM2_ERR = 0

         OPEN (EIN, FILE=EINFIL, STATUS='OLD', IOSTAT=IOCHK)

         IF (IOCHK == 0) THEN                              ! File was opened successfully

            WRITE(ERR,*)
            WRITE(ERR,5996) EIGNORM2
            IF (SUPINFO == 'N') THEN
               WRITE(F06,*)
               WRITE(F06,5996) EIGNORM2
            ENDIF
            CALL WRITE_FILNAM ( EINFIL, F06, 15 )
            READ(EIN,'(A)') TITLE

            NUM_CHANGED = 0

            REC_NUM = 0
            DO

               READ(EIN,'(A)',IOSTAT=IOCHK) DATA_FIELD   ;   REC_NUM = REC_NUM + 1

               IF (IOCHK == 0) THEN

                  IF (DATA_FIELD(1:3) == 'END') EXIT

                  CALL I4FLD (DATA_FIELD(1:JCARD_LEN), 1, VEC_NUM )
                  IF (IERRFL(1) == 'N') THEN
                     IF ((VEC_NUM >= 1) .AND. (VEC_NUM <= NVEC)) THEN
                        VEC_SIGN_CHG(VEC_NUM) = .TRUE.
                        NUM_CHANGED = NUM_CHANGED + 1
                        VECS_CHANGED(NUM_CHANGED) = VEC_NUM
                     ELSE
                        WRITE(ERR,5992) REC_NUM, DATA_FIELD, NVEC
                        WRITE(F06,5992) REC_NUM, DATA_FIELD, NVEC
                     ENDIF
                  ELSE
                     WRITE(ERR,5995) REC_NUM, DATA_FIELD
                     WRITE(F06,5995) REC_NUM, DATA_FIELD
                     EXIT
                  ENDIF

                  CYCLE

               ELSE

                  WRITE(ERR,5991)
                  CALL WRITE_FILNAM ( EINFIL, ERR, 15 )
                  WRITE(ERR,59911)

                  WRITE(F06,5991)
                  CALL WRITE_FILNAM ( EINFIL, F06, 15 )
                  WRITE(F06,59911)

                  EIGNORM2_ERR = EIGNORM2_ERR + 1
                  EXIT

               ENDIF

               IF (REC_NUM > NDOFL) THEN
                  EIGNORM2_ERR = EIGNORM2_ERR + 1
                  WRITE(ERR,5998) EINFIL, REC_NUM, NVEC
                  WRITE(F06,5998) EINFIL, REC_NUM, NVEC
               ENDIF

            ENDDO

         ELSE                                              ! File not opened successfully, so no vecs will have signs changed

            WRITE(ERR,5997)
            CALL WRITE_FILNAM ( EINFIL, ERR, 15 )
            WRITE(ERR,59971)

            WRITE(F06,5997)
            CALL WRITE_FILNAM ( EINFIL, F06, 15 )
            WRITE(F06,59971)

            EIGNORM2_ERR = EIGNORM2_ERR + 1

         ENDIF

      ELSE

         WRITE(ERR,5990)
         CALL WRITE_FILNAM ( EINFIL, ERR, 15 )
         WRITE(ERR,59901)

         WRITE(F06,5990)
         CALL WRITE_FILNAM ( EINFIL, F06, 15 )
         WRITE(F06,59901)

         EIGNORM2_ERR = EIGNORM2_ERR + 1

      ENDIF

      IF (EIGNORM2_ERR == 0) THEN
         WRITE(ERR,*)
         WRITE(ERR,5993)
         WRITE(ERR,5994) (VECS_CHANGED(I),I=1,NUM_CHANGED)
         WRITE(ERR,*)
         IF (SUPINFO == 'N') THEN
            WRITE(F06,*)
            WRITE(F06,5993)
            WRITE(F06,5994) (VECS_CHANGED(I),I=1,NUM_CHANGED)
            WRITE(F06,*)
         ENDIF
      ELSE
         WRITE(F06,5999)
         CALL WRITE_FILNAM ( EINFIL, F06, 15 )
      ENDIF

      RETURN

! **********************************************************************************************************************************
 5990 FORMAT(' *WARNING    : PARAMETER EIGNORM2 REQUESTS THAT SOME EIGENVECTORS ARE TO HAVE THEIR SIGN CHANGED, BUT FILE:')

59901 FORMAT('               WHICH SHOULD HAVE THE NUMBERS OF THE VECTORS TO HAVE THEIR SIGN CHANGED CANNOT BE FOUND')

 5991 FORMAT(' *WARNING    : EOF/ERR READING FILE:')

59911 FORMAT('               ON RECORD NUMBER ',I8,'. FIELD THAT HAS ERROR = "',A,'"')

 5992 FORMAT(' *WARNING    : ON REC NUM ',I8,' EIGENVEC ',A,' FOR SIGN CHANGE WAS OUT OF RANGE.',                                  &
                           ' MUST BE 1 TO NVEC = ',I8,'. VALUE IGNORED')

 5993 FORMAT(' *INFORMATION: THE FOLLOWING EIGENVECTORS HAVE HAD THEIR SIGN CHANGED:')

 5994 FORMAT(10I8)

 5995 FORMAT(' *WARNING    : ERROR READING RECORD NUMBER ',I8,' = ',A,'. REMAINDER OF EIGENVECTOR SIGN CHANGE FILE IGNORED')

 5996 FORMAT(' *INFORMATION: BASED ON PARAMETER EIGNORM2 = ',A,' SOME EIGENVECS WILL HAVE THEIR SIGN CHANGED AS DEFINED IN FILE:')

 5997 FORMAT(' *WARNING    : PARAMETER EIGNORM2 REQUESTS THAT SOME EIGENVECTORS ARE TO HAVE THEIR SIGN CHANGED, BUT FILE:')

59971 FORMAT('               WHICH SHOULD HAVE THE NUMBERS OF THE VECTORS TO HAVE THEIR SIGN CHANGED COULD NOT BE OPENED')

 5998 FORMAT(' *WARNING    : ERROR READING FILE: ',A                                                                               &
                    ,/,14X,' ATTEMPTING TO READ RECORD NUMBER ',I8,' BUT THERE SHOULD ONLY BE NVEC = ',I8,' RECORDS')

 5999 FORMAT(' *WARNING    : NO EIGENVECTORS WILL HAVE THEIR SIGN CHANGED DUE TO ABOVE ERRORS READING FILE:')

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

      END SUBROUTINE READ_EIGNORM2

      END SUBROUTINE LINK5