text:
+incdir+../source/

//Source File (PKG, RTL, MACROS)
../source/alu.sv
../source/cpu.sv
../source/definitions.sv
../source/ram.sv
../source/rom.sv
../source/cpu_garage.sv

//Test Bench
../source/cpu_garage_tb.sv

ext: f

max_stars_repo_name: amichai-bd/CPUGarageChallenge

max_issues_repo_name: amichai-bd/CPUGarageChallenge

max_forks_repo_name: amichai-bd/CPUGarageChallenge


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2017 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
! Author: Miriam Hinzen

PROGRAM test

  use iso_c_binding

  interface
    subroutine chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt ) bind( c, name = 'dchase_' )
      use iso_c_binding
      real(c_double)                :: h(n,*), v(n,*)
      integer(c_int)                :: n, deg, nev, nex
      real(c_double)                :: ritzv(*), tol
      character(len=1,kind=c_char)  :: mode, opt
    end subroutine chase_r
  end interface

  real(c_double), ALLOCATABLE   :: h(:,:), v(:,:), ritzv(:)
  integer(c_int)                :: n, deg, nev, nex
  real(c_double)                :: tol
  character(len=1,kind=c_char)  :: mode, opt

  CALL chase_r( h, n, v, ritzv, nev, nex, deg, tol, mode, opt )

END

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
c
c Copyright (c) 2018, Lawrence Livermore National Security, LLC.
c Produced at the Lawrence Livermore National Laboratory
c
c Written by Jeffrey Banks banksj3@rpi.edu (Rensselaer Polytechnic Institute,
c Amos Eaton 301, 110 8th St., Troy, NY 12180); Jeffrey Hittinger
c hittinger1@llnl.gov, William Arrighi arrighi2@llnl.gov, Richard Berger
c berger5@llnl.gov, Thomas Chapman chapman29@llnl.gov (LLNL, P.O Box 808,
c Livermore, CA 94551); Stephan Brunner stephan.brunner@epfl.ch (Ecole
c Polytechnique Federale de Lausanne, EPFL SB SPC-TH, PPB 312, Station 13,
c CH-1015 Lausanne, Switzerland).
c CODE-744849
c
c All rights reserved.
c
c This file is part of Loki.  For details, see.
c
c Permission is hereby granted, free of charge, to any person obtaining a
c copy of this software and associated documentation files (the "Software"),
c to deal in the Software without restriction, including without limitation
c the rights to use, copy, modify, merge, publish, distribute, sublicense,
c and/or sell copies of the Software, and to permit persons to whom the
c Software is furnished to do so, subject to the following conditions:
c
c The above copyright notice and this permission notice shall be included in
c all copies or substantial portions of the Software.
c
c THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
c OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
c FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
c THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
c LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
c FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
c DEALINGS IN THE SOFTWARE.
c
c Fortran functions called by ShapedRampedCosineCurrentDriver.
c
      subroutine evaluateShapedRampedCurrentDriver(
     &     Jx,
     &     Jy,
     &     Jz,
     &     nd1a, nd1b, nd2a, nd2b,
     &     n1a, n1b, n2a, n2b,
     &     xlo, xhi, dx,
     &     t,
     &     e_ext_param)
c
c.. function to compute applied current
      implicit none
c
      integer nd1a, nd1b, nd2a, nd2b
      integer n1a, n1b, n2a, n2b
      real xlo(2), xhi(2), dx(2), t
      real Jx(nd1a:nd1b, nd2a:nd2b)
      real Jy(nd1a:nd1b, nd2a:nd2b)
      real Jz(nd1a:nd1b, nd2a:nd2b)
      real e_ext_param(*)
c
      integer i1_apply, i2
c
      real width, apply_dir, shape, omega, J_0, t0, t_ramp, t_off, x0
      real envel, ycoord, J_ext, pi, i1_exact
c
      real one, four
c
      one = 1.0
      four = 4.0
c
      pi = four*atan(one)
c
      width     = e_ext_param(1)
      apply_dir = e_ext_param(2)
      shape     = e_ext_param(3)
      omega     = e_ext_param(4)
      J_0       = e_ext_param(5)
      t0        = e_ext_param(6)
      t_ramp    = e_ext_param(7)
      t_off     = e_ext_param(8)
      x0        = e_ext_param(9)
c
      if ((t .lt. (t0+t_ramp+t_off)) .and. t .ge. t0) then
        i1_exact = (x0 - xlo(1)) / dx(1)
        if (i1_exact .lt. 0) then
           i1_exact = i1_exact - 0.5
        else
           i1_exact = i1_exact + 0.5
        end if
        i1_apply = i1_exact
        if (t .lt. t0+t_ramp) then
          envel = 0.5+0.5*tanh(4.0*(2.0*(t-t0)/t_ramp-1.0))
        else
          envel = 0.5-0.5*tanh(4.0*(2.0*(t-t0-t_ramp)/t_off-1.0))
        end if
        do i2 = n2a, n2b
c
          ycoord = xlo(2)+dx(2)*i2
          if (abs(ycoord) .lt. 0.5*width) then
            J_ext = J_0*envel*cos(omega*(t-t0))*
     *            (1.0-shape*(sin(pi*ycoord/width))**2)
          else
            J_ext = J_0*envel*cos(omega*(t-t0))*(1.0-shape)
          endif
c
          if (apply_dir .eq. 1.0) then
             Jx(i1_apply, i2) = Jx(i1_apply, i2) + J_ext
          else if (apply_dir .eq. 2.0) then
             Jy(i1_apply, i2) = Jy(i1_apply, i2) + J_ext
          else
             Jz(i1_apply, i2) = Jz(i1_apply, i2) + J_ext
          end if
        end do
      end if
c
      return
      end

max_stars_repo_name: LLNL/LOKI

max_issues_repo_name: LLNL/LOKI

max_forks_repo_name: LLNL/LOKI


      PROGRAM LA_SSTEV_EXAMPLE

!  -- LAPACK95 EXAMPLE DRIVER ROUTINE (VERSION 1.0) --
!     UNI-C, DENMARK
!     DECEMBER, 1999
!
!  .. "Use Statements" ..
      USE LA_PRECISION, ONLY: WP => SP
      USE F95_LAPACK, ONLY: LA_STEV, LA_STEVD
!  .. "Implicit Statement" ..
      IMPLICIT NONE
!  .. "Local Scalars" ..
      INTEGER :: I, J, INFO, N
!  .. "Local Arrays" ..
      REAL(WP), ALLOCATABLE :: D(:), E(:), Z(:,:)
!  .. "Executable Statements" ..
      WRITE (*,*) 'LA_STEV Example Program Results'
      N = 5;
      ALLOCATE( D(N), E(N), Z(N,N) )

      OPEN(UNIT=21,FILE='stevd.ma',STATUS='UNKNOWN')
      DO J=1,N
         READ(21,*) D(J)
      ENDDO
      CLOSE(21)

      OPEN(UNIT=21,FILE='steve.ma',STATUS='UNKNOWN')
      DO J=1,N
         READ(21,*) E(J)
      ENDDO
      CLOSE(21)

      WRITE(*,*)' Vector D : '
      DO I=1,N
      WRITE(*,"(5(I5,1X))") INT(D(I))
      ENDDO

      WRITE(*,*)' Vector E : '
      DO I=1,N
      WRITE(*,"(5(I5,1X))") INT(E(I))
      ENDDO

      WRITE(*,*) 'CALL LA_STEV( D, E )'
      CALL LA_STEV( D, E )

      WRITE(*,*) ' D on exit : '
      DO I=1,N
      WRITE(*,"(5(F12.5,1X))") D(I)
      ENDDO
      WRITE(*,*)
      WRITE(*,*) ' * EXAMPLE 2 * '

      OPEN(UNIT=21,FILE='stevd.ma',STATUS='UNKNOWN')
      DO J=1,N
         READ(21,*) D(J)
      ENDDO
      CLOSE(21)

      OPEN(UNIT=21,FILE='steve.ma',STATUS='UNKNOWN')
      DO J=1,N
         READ(21,*) E(J)
      ENDDO
      CLOSE(21)

      WRITE(*,*)' Vector D : '
      DO I=1,N
      WRITE(*,"(5(I5,1X))") INT(D(I))
      ENDDO

      WRITE(*,*)' Vector E : '
      DO I=1,N
      WRITE(*,"(5(I5,1X))") INT(E(I))
      ENDDO

      WRITE(*,*) "CALL LA_STEVD( D, E, Z, INFO )"
      CALL LA_STEVD( D, E, Z, INFO )

      WRITE(*,*)' Z on exit : '
      DO I=1,N;
         WRITE(*,"(5(F12.6,1X))") Z(I,:);
      ENDDO

      WRITE(*,*)' D on exit : '
      DO I=1,N
      WRITE(*,"(5(F12.5,1X))") D(I)
      ENDDO

      WRITE(*,*) ' INFO = ', INFO

      END PROGRAM LA_SSTEV_EXAMPLE


text:
      SUBROUTINE EXCSD(Q,Z,KEY,NSTEP,HEXTTL,AOPT)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
C    *****************************************************************
C    *                                                               *
C    *   THIS IS A CSC EXECUTIVE ROUTINE FOR FAST FOURIER TRANSFORM  *
C    *   ANALYSIS. MODIFIED TO CHANGE CALLING SEQUENCE SEPT. 1977.   *
C    *                                                               *
C    *****************************************************************
C
      COMPLEX*16 HEXTTL,AOPT
      COMPLEX*8    Z(1025),ZSTORE(1025)
      REAL*4 DELT, Q(2050),QSTORE(2050),DUMMY(7)
C     EQUIVALENCE (Q(1),Z(1))
      COMMON/RMAIN1/DELTAT,FACTOR,FREQ,TSTOP,DELMIT
      WRITE(6,937)
 937  FORMAT(1H1)
      WRITE(6,938)
 938  FORMAT(' ********************************     FAST    FOURIER
     1 TRANSFORM    ANALYSIS    ***********************')
C
      WRITE(6,600) KEY
  600 FORMAT(/,20X,'&&&&    FOR KEY =',I3,/)
C
      WRITE(6,6000) HEXTTL,AOPT
 6000 FORMAT('0',5X,10('*'),5X,'DATA RELATED TO ',2A8,2A8,' PLOT',5X,
     *       10('*'),/)
C
C
      WRITE(6,625) NSTEP
  625 FORMAT(/,20X,'NO. OF POINTS =',I4,/)
C
      NHALF=(NSTEP+1)/2
      DELT=FREQ
      DO 201 I=1,NHALF
  201 ZSTORE(I)=Z(I)
      DO 301 I=1,NSTEP
  301 QSTORE(I)=Q(I)
      CALL CSD(Z,Q,NSTEP,NHALF,IWIN,DELT,DUMMY)
      DO 202 I=1,NHALF
  202 Z(I)=ZSTORE(I)
      DO 302 I=1,NSTEP
  302 Q(I)=QSTORE(I)
   15 CONTINUE
      RETURN
      END

ext: for


text:
Describe Users/tombrokaw here.

20080323 16:28:51 nbsp Please consider the importance of using your RealName. Users/JasonAller
  For all we know, his name might coincidentally be Tom Brokaw. Its like if we get a new user whose name happened to be William Smith or Jack Chan or something...  Paul Amnuaypayoat

20080324 13:29:35 nbsp Isnt there a Walter Cronkite and a Dan Rather here too?  Amazing! Users/RichLindvall

20080324 18:16:21 nbsp Equally as amazing is the revelation that, according to Tom, his mom used to cook beef noodle soup for him just like they do at Hometown Chinese Food! ;) Users/DukeMcAdow

20080324 23:58:25 nbsp No Im not really tom brokaw.  Sorry about that.  Just a gaming handle I used.  Ill consider using my real name after I get more familiar with the place. Users/tombrokaw

20100207 01:00:52 nbsp Hey, many people didnt move here at at all... Users/NickSchmalenberger

20100801 23:45:12 nbsp Yeah bro, totally what Nick Said Users/StevenDaubert

20100802 12:28:46 nbsp ? Users/tombrokaw

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      MODULE qed_slfen_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.3E  06:33:54  12/28/06
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      SUBROUTINE qed_slfen (SLFINT)
      USE vast_kind_param, ONLY: DOUBLE
      USE parameter_def,   ONLY: NNNW
      REAL(DOUBLE), DIMENSION(NNNW), INTENT(OUT) :: SLFINT
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: compas/grasp-rci-qed

max_forks_repo_name: compas/grasp-rci-qed


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!>   This module provide common functions of 3D truss elements
module m_static_LIB_1d
  use hecmw, only : kint, kreal
  use elementInfo
  implicit none

contains
  !
  !=====================================================================*
  !  STF_C1
  !=====================================================================*
  !>  This subroutine calculate stiff matrix of 2-nodes truss element
  subroutine STF_C1( etype,nn,ecoord,area,gausses,stiff, u ,temperature )
    use mMechGauss
    integer(kind=kint), intent(in)  :: etype               !< element type
    integer(kind=kint), intent(in)  :: nn                  !< number of elemental nodes
    real(kind=kreal),   intent(in)  :: ecoord(3,nn)        !< coordinates of elemental nodes
    real(kind=kreal),   intent(in)  :: area                !< section area
    type(tGaussStatus), intent(in)  :: gausses(:)          !< status of qudrature points
    real(kind=kreal),   intent(out) :: stiff(:,:)          !< stiff matrix
    real(kind=kreal),   intent(in), optional :: u(:,:)     !< nodal displacemwent
    real(kind=kreal),   intent(in), optional :: temperature(nn)     !< temperature

    real(kind=kreal) DET,WG, llen, llen0, elem(3,nn)
    logical :: ierr
    real(kind=kreal) ina(1), outa(1), direc(3), direc0(3), coeff, strain
    integer(kind=kint) :: i,j

    ! we suppose the same material type in the element
    if( present(u) ) then
      elem(:,:) = ecoord(:,:) + u(:,:)
    else
      elem(:,:) = ecoord(:,:)
    endif

    direc = elem(:,2)-elem(:,1)
    llen = dsqrt( dot_product(direc, direc) )
    direc = direc/llen
    direc0 = ecoord(:,2)-ecoord(:,1)
    llen0 = dsqrt( dot_product(direc0, direc0) )

    if( present(temperature) ) then
      ina(1) = 0.5d0*(temperature(1)+temperature(2))
      call fetch_TableData( MC_ISOELASTIC, gausses(1)%pMaterial%dict, outa, ierr, ina )
    else
      call fetch_TableData( MC_ISOELASTIC, gausses(1)%pMaterial%dict, outa, ierr )
    endif
    if( ierr ) outa(1) = gausses(1)%pMaterial%variables(M_YOUNGS)
    coeff = outa(1)*area*llen0/(llen*llen)
    strain = gausses(1)%strain(1)

    stiff(:,:) = 0.d0
    do i=1,3
      stiff(i,i) = coeff*strain
      do j=1,3
        stiff(i,j) = stiff(i,j) + coeff*(1.d0-2.d0*strain)*direc(i)*direc(j)
      enddo
    enddo

    stiff(4:6,1:3) = -stiff(1:3,1:3)
    stiff(1:3,4:6) = transpose(stiff(4:6,1:3))
    stiff(4:6,4:6) = stiff(1:3,1:3)

  end subroutine STF_C1

  !
  !> Update strain and stress inside element
  !---------------------------------------------------------------------*
  subroutine UPDATE_C1( etype, nn, ecoord, area, u, du, qf ,gausses, TT, T0 )
    !---------------------------------------------------------------------*
    use m_fstr
    use mMechGauss
    ! I/F VARIAVLES
    integer(kind=kint), intent(in)     :: etype           !< \param [in] element type
    integer(kind=kint), intent(in)     :: nn              !< \param [in] number of elemental nodes
    real(kind=kreal),   intent(in)     :: ecoord(3,nn)    !< \param [in] coordinates of elemental nodes
    real(kind=kreal),   intent(in)     :: area            !< section area
    real(kind=kreal),   intent(in)     :: u(3,nn)         !< \param [in] nodal dislplacements
    real(kind=kreal),   intent(in)     :: du(3,nn)       !< \param [in] nodal displacement ( solutions of solver )
    real(kind=kreal),   intent(out)    :: qf(nn*3)        !< \param [out] Internal Force
    type(tGaussStatus), intent(inout)  :: gausses(:)      !< \param [out] status of qudrature points
    real(kind=kreal),   intent(in), optional :: TT(nn)    !< current temperature
    real(kind=kreal),   intent(in), optional :: T0(nn)    !< reference temperature

    ! LCOAL VARIAVLES
    real(kind=kreal)   :: direc(3), direc0(3)
    real(kind=kreal)   :: llen, llen0, ina(1), outa(1)
    real(kind=kreal)   :: elem(3,nn)
    real(kind=kreal)   :: young
    real(kind=kreal)   :: ttc, tt0, alp, alp0, epsth
    logical            :: ierr

    qf(:)              = 0.d0
    ! we suppose the same material type in the element
    elem(:,:) = ecoord(:,:) + u(:,:) + du(:,:)

    direc = elem(:,2)-elem(:,1)
    llen = dsqrt( dot_product(direc, direc) )
    direc = direc/llen
    direc0 = ecoord(:,2)-ecoord(:,1)
    llen0 = dsqrt( dot_product(direc0, direc0) )

    epsth = 0.d0
    if( present(tt) .and. present(t0) ) then
      ttc = 0.5d0*(TT(1)+TT(2))
      tt0 = 0.5d0*(T0(1)+T0(2))

      ina(1) = ttc
      call fetch_TableData( MC_ISOELASTIC, gausses(1)%pMaterial%dict, outa, ierr, ina )
      if( ierr ) outa(1) = gausses(1)%pMaterial%variables(M_YOUNGS)
      young = outa(1)

      call fetch_TableData( MC_THEMOEXP, gausses(1)%pMaterial%dict, outa(:), ierr, ina )
      if( ierr ) outa(1) = gausses(1)%pMaterial%variables(M_EXAPNSION)
      alp = outa(1)

      ina(1) = tt0
      call fetch_TableData( MC_THEMOEXP, gausses(1)%pMaterial%dict, outa(:), ierr, ina )
      if( ierr ) outa(1) = gausses(1)%pMaterial%variables(M_EXAPNSION)
      alp0 = outa(1)

      epsth=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
    else
      call fetch_TableData( MC_ISOELASTIC, gausses(1)%pMaterial%dict, outa, ierr )
      if( ierr ) outa(1) = gausses(1)%pMaterial%variables(M_YOUNGS)
      young = outa(1)
    endif

    gausses(1)%strain(1) = dlog(llen/llen0)
    gausses(1)%stress(1) = young*(gausses(1)%strain(1)-epsth)

    !set stress and strain for output
    gausses(1)%strain_out(1) = gausses(1)%strain(1)
    gausses(1)%stress_out(1) = gausses(1)%stress(1)

    qf(1) = gausses(1)%stress(1)*area*llen0/llen
    qf(1:3) = -qf(1)*direc
    qf(4:6) = -qf(1:3)

  end subroutine UPDATE_C1

  !----------------------------------------------------------------------*
  subroutine NodalStress_C1(ETYPE,NN,gausses,ndstrain,ndstress)
    !----------------------------------------------------------------------*
    !
    ! Calculate Strain and Stress increment of solid elements
    !
    use mMechGauss
    integer(kind=kint), intent(in) :: ETYPE,NN
    type(tGaussStatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out)  :: ndstrain(NN,6)
    real(kind=kreal), intent(out)  :: ndstress(NN,6)

    ndstrain(1,1:6) = gausses(1)%strain(1:6)
    ndstress(1,1:6) = gausses(1)%stress(1:6)
    ndstrain(2,1:6) = gausses(1)%strain(1:6)
    ndstress(2,1:6) = gausses(1)%stress(1:6)

  end subroutine
  !
  !
  !
  !----------------------------------------------------------------------*
  subroutine ElementStress_C1(ETYPE,gausses,strain,stress)
    !----------------------------------------------------------------------*
    !
    ! Calculate Strain and Stress increment of solid elements
    !
    use mMechGauss
    integer(kind=kint), intent(in) :: ETYPE
    type(tGaussStatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out)  :: strain(6)
    real(kind=kreal), intent(out)  :: stress(6)

    strain(:) = gausses(1)%strain(1:6)
    stress(:) = gausses(1)%stress(1:6)

  end subroutine

  !----------------------------------------------------------------------*
  subroutine DL_C1(etype, nn, xx, yy, zz, rho, thick, ltype, params, &
      vect, nsize)
    !----------------------------------------------------------------------*
    !**  SET DLOAD
    !   GRAV LTYPE=4  :GRAVITY FORCE
    ! I/F VARIABLES
    integer(kind = kint), intent(in)  :: etype, nn
    real(kind = kreal), intent(in)    :: xx(:), yy(:), zz(:)
    real(kind = kreal), intent(in)    :: params(0:6)
    real(kind = kreal), intent(inout) :: vect(:)
    real(kind = kreal) :: rho, thick
    integer(kind = kint) :: ltype, nsize, surtype
    ! LOCAL VARIABLES
    integer(kind = kint) :: ndof = 3
    integer(kind = kint) :: ivol, isuf, nsur, i
    integer(kind = kint) :: nod(nn)
    real(kind = kreal) :: vx, vy, vz, val, a, AA
    !--------------------------------------------------------------------
    val = params(0)
    !--------------------------------------------------------------

    ivol = 0
    isuf = 0
    if( ltype .LT. 10 ) then
      ivol = 1
    else if( ltype .GE. 10 ) then
      isuf = 1
    end if

    !--------------------------------------------------------------------
    nsize = nn*ndof
    !--------------------------------------------------------------------
    vect(1:nsize) = 0.0D0

    ! Volume force
    if( ivol .EQ. 1 ) then
      if( ltype .EQ. 4 ) then
        AA = dsqrt( ( xx(2)-xx(1) )*( xx(2)-xx(1) )   &
          +( yy(2)-yy(1) )*( yy(2)-yy(1) )   &
          +( zz(2)-zz(1) )*( zz(2)-zz(1) ) )

        a = thick
        vx = params(1)
        vy = params(2)
        vz = params(3)
        vx = vx/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
        vy = vy/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
        vz = vz/dsqrt( params(1)**2+params(2)**2+params(3)**2 )

        do i = 1, 2
          vect(3*i-2) = val*rho*a*0.5D0*AA*vx
          vect(3*i-1) = val*rho*a*0.5D0*AA*vy
          vect(3*i  ) = val*rho*a*0.5D0*AA*vz
        end do

      end if
    end if
    !--------------------------------------------------------------------

    return
  end subroutine

  !
  !
  !----------------------------------------------------------------------*
  subroutine truss_diag_modify(hecMAT,hecMESH)
    !----------------------------------------------------------------------*
    !
    use hecmw
    type (hecmwST_matrix)     :: hecMAT
    type (hecmwST_local_mesh) :: hecMESH
    integer(kind=kint) :: itype, is, iE, ic_type, icel, jS, j, n

    do itype = 1, hecMESH%n_elem_type
      ic_type = hecMESH%elem_type_item(itype)
      if(ic_type == 301)then
        is = hecMESH%elem_type_index(itype-1) + 1
        iE = hecMESH%elem_type_index(itype  )
        do icel = is, iE
          jS = hecMESH%elem_node_index(icel-1)
          do j=1,2
            n = hecMESH%elem_node_item(jS+j)
            if( hecMAT%D(9*n-8) == 0.0d0)then
              hecMAT%D(9*n-8) = 1.0d0
              !call search_diag_modify(n,1,hecMAT,hecMESH)
            endif
            if( hecMAT%D(9*n-4) == 0.0d0)then
              hecMAT%D(9*n-4) = 1.0d0
              !call search_diag_modify(n,2,hecMAT,hecMESH)
            endif
            if( hecMAT%D(9*n  ) == 0.0d0)then
              hecMAT%D(9*n  ) = 1.0d0
              !call search_diag_modify(n,3,hecMAT,hecMESH)
            endif
          enddo
        enddo
      endif
    enddo

  end subroutine

  !
  !
  !----------------------------------------------------------------------*
  subroutine search_diag_modify(n,nn,hecMAT,hecMESH)
    !----------------------------------------------------------------------*
    !
    use hecmw
    type (hecmwST_matrix)     :: hecMAT
    type (hecmwST_local_mesh) :: hecMESH
    integer :: n, nn, is, iE, i, j, in
    integer :: flagl, flagu, flagb, a

    if(nn == 1)then
      a = 0
      is = hecMAT%IndexL(n-1)+1
      iE = hecMAT%IndexL(n  )
      do i=is,iE
        if(hecMAT%AL(9*i-8) /= 0.0d0) a = 1
        if(hecMAT%AL(9*i-7) /= 0.0d0) a = 1
        if(hecMAT%AL(9*i-6) /= 0.0d0) a = 1
      enddo
      is = hecMAT%IndexU(n-1)+1
      iE = hecMAT%IndexU(n  )
      do i=is,iE
        if(hecMAT%AU(9*i-8) /= 0.0d0) a = 1
        if(hecMAT%AU(9*i-7) /= 0.0d0) a = 1
        if(hecMAT%AU(9*i-6) /= 0.0d0) a = 1
      enddo
      if(hecMAT%D(9*n-7) /= 0.0d0) a = 1
      if(hecMAT%D(9*n-6) /= 0.0d0) a = 1
      if(a == 0)then
        hecMAT%D(9*n-8) = 1.0d0
        !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:1"
      endif
    endif
    if(nn == 2)then
      a = 0
      is = hecMAT%IndexL(n-1)+1
      iE = hecMAT%IndexL(n  )
      do i=is,iE
        if(hecMAT%AL(9*i-5) /= 0.0d0) a = 1
        if(hecMAT%AL(9*i-4) /= 0.0d0) a = 1
        if(hecMAT%AL(9*i-3) /= 0.0d0) a = 1
      enddo
      is = hecMAT%IndexU(n-1)+1
      iE = hecMAT%IndexU(n  )
      do i=is,iE
        if(hecMAT%AU(9*i-5) /= 0.0d0) a = 1
        if(hecMAT%AU(9*i-4) /= 0.0d0) a = 1
        if(hecMAT%AU(9*i-3) /= 0.0d0) a = 1
      enddo
      if(hecMAT%D(9*n-5) /= 0.0d0) a = 1
      if(hecMAT%D(9*n-3) /= 0.0d0) a = 1
      if(a == 0)then
        hecMAT%D(9*n-4) = 1.0d0
        !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:2"
      endif
    endif
    if(nn == 3)then
      a = 0
      is = hecMAT%IndexL(n-1)+1
      iE = hecMAT%IndexL(n  )
      do i=is,iE
        if(hecMAT%AL(9*i-2) /= 0.0d0) a = 1
        if(hecMAT%AL(9*i-1) /= 0.0d0) a = 1
        if(hecMAT%AL(9*i  ) /= 0.0d0) a = 1
      enddo
      is = hecMAT%IndexU(n-1)+1
      iE = hecMAT%IndexU(n  )
      do i=is,iE
        if(hecMAT%AU(9*i-2) /= 0.0d0) a = 1
        if(hecMAT%AU(9*i-1) /= 0.0d0) a = 1
        if(hecMAT%AU(9*i  ) /= 0.0d0) a = 1
      enddo
      if(hecMAT%D(9*n-2) /= 0.0d0) a = 1
      if(hecMAT%D(9*n-1) /= 0.0d0) a = 1
      if(a == 0)then
        hecMAT%D(9*n  ) = 1.0d0
        !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:3"
      endif
    endif

  end subroutine

end module m_static_LIB_1d


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.02.03.01/09                        *
C  *    TEST TITLE : Character up vector                   *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

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

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

C text path
      INTEGER    PRIGHT,     PLEFT,     PUP,     PDOWN
      PARAMETER (PRIGHT = 0, PLEFT = 1, PUP = 2, PDOWN = 3)

C text precision
      INTEGER    PSTRP,     PCHARP,     PSTRKP
      PARAMETER (PSTRP = 0, PCHARP = 1, PSTRKP = 2)

      INTEGER    PICSTR, TXCI, NUMSTR, NGSTR

      CALL INITGL ('04.02.03.01/09')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C set-up of workstation and dialogue area
      PICSTR = 101
      TXCI = 1
      CALL SETDLG (PICSTR, 801,TXCI)

      CALL POPST (PICSTR)
C  by convention, view #1 is for picture
      CALL PSVWI (1)
C  use individual attributes
      CALL SETASF (PINDIV)
      CALL PSTXPR (PSTRKP)
      CALL PSCHH  (1.0)
      CALL PEXST (102)
      CALL PCLST
      CALL POPST (102)

C  use font #1 for horizontal path
      CALL TSCHUP (1, PRIGHT, 'Going up?', NUMSTR, NGSTR)
      CALL SETMSG ('16 27 31 66 67 69 74 75 80 81', 'A text ' //
     1             'primitive with STROKE precision and a '   //
     2             'horizontal text path should be rotated '  //
     3             'according to the character up vector.')
      CALL DCHPFV ('CHARACTER UP FOR HORIZONTAL TEXT PATH: Which ' //
     1             'text primitive is NOT aligned within the '     //
     2             'dotted text extent rectangle on the dashed '   //
     3             'baseline?', NUMSTR, NGSTR)
      CALL PEMST (102)

C  use font #2 for vertical path
      CALL TSCHUP (2, PDOWN, 'Up&@M', NUMSTR, NGSTR)
      CALL SETMSG ('16 27 31 66 67 68 74 75 80 81', 'A text '        //
     1             'primitive with STROKE precision and a vertical ' //
     2             'text path should be rotated according to the '   //
     3             'character up vector.')
      CALL DCHPFV ('CHARACTER UP FOR VERTICAL TEXT PATH: Which ' //
     1             'text primitive is NOT aligned within the '   //
     2             'dotted text extent rectangle on the dashed ' //
     3             'centerline?', NUMSTR, NGSTR)
      CALL PEMST (102)

666   CONTINUE
C  wrap it up.
      CALL ENDIT
      END


text:
C$Procedure      VSCL ( Vector scaling, 3 dimensions )

      SUBROUTINE VSCL ( S, V1, VOUT )

C$ Abstract
C
C     Multiply a scalar and a 3-dimensional double precision vector.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     VECTOR
C
C$ Declarations

      DOUBLE PRECISION   S
      DOUBLE PRECISION   V1   ( 3 )
      DOUBLE PRECISION   VOUT ( 3 )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     S          I   Scalar to multiply a vector.
C     V1         I   Vector to be multiplied.
C     VOUT       O   Product vector, S*V1.
C
C$ Detailed_Input
C
C     S    This is a double precision scalar used to multiply the
C          vector V1.
C
C     V1   This is a 3-dimensional, double precision vector which is
C          to be scaled by S.
C
C$ Detailed_Output
C
C     VOUT   This is a 3-dimensional, double precision vector which
C            is the scalar multiple of V1.  VOUT = S*V1.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     VSCL multiplies each component of V1 by S to form the respective
C     components of VOUT.  No error checking is performed.
C
C$ Examples
C
C     The following table shows the output VOUT as a function of the
C     the inputs V1, and S from the subroutine VSCL.
C
C        V1                   S         VOUT
C        ----------------------------------------------
C        (1D0, -2D0, 0D0)   -1D0       (-1D0, 2D0, 0D0)
C        (0D0, 0D0, 0D0)     5D0       (0D0, 0D0, 0D0)
C
C$ Restrictions
C
C     The user is responsible for insuring that no floating point
C     overflow occurs from multiplying S by any component of V1. No
C     error recovery or reporting scheme is incorporated in this
C     subroutine.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.M. Owen       (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.2, 23-APR-2010 (NJB)
C
C        Header correction: assertions that the output
C        can overwrite the input have been removed.
C
C-    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.0, 31-JAN-1990 (WMO)
C
C-&

C$ Index_Entries
C
C     3-dimensional vector scaling
C
C-&

      VOUT(1) = S * V1(1)
      VOUT(2) = S * V1(2)
      VOUT(3) = S * V1(3)

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


c    Subroutines for the program sdf_lister
c
c     Contents:
c
c             subroutine  assign_sectored_fluxes
c
c             subroutine assign_sectored_vse
c
c
c**********************************************************
c
      subroutine assign_sectored_fluxes
     1         (local_sectored_flag,
     1          local_tel_flag,
     1          flux_sect, flux_sect_sigma, flux_sect_flag)
c
c**********************************************************
c
c    calculates sectored fluxes and sigmas for 1 or both telescopes
c
c     11/10/95 by J. Dwyer
c             modified 1/8/96 by J. Dwyer divide sector exposures by 8.0
c
      include 'sdf_include.inc'   ! include type declarations
c           modified 4/22/96 corrected addition of two tel. when effic_flag = -1
c
      real    flux_tel(2,8), sigma_tel2(2,8) ! arrays containing the fluxes and sigmas for all the elements
      real    flux_sect(8), flux_sect_sigma(8)   ! arrays containing the fluxes and sigmas for all the elements
      integer   flux_flag_tel(2,8)   ! .true. when flux is measured and is not saturated, .false. otherwise
      integer   flux_sect_flag(8)   ! .true. when flux is measured and is not saturated, .false. otherwise
      integer  local_tel_flag   ! 1=use tel. 1 only, 2=use tel. 1 only,-1=use tel. 1 and tel. 2
      integer   Ntel
      logical  local_sectored_flag(41)  ! .true.=include in sectored flux calculation, .false.=don't include
c
      do 5 k=1,8
          flux_sect_flag(k) = 0
          flux_sect(k) = 0.0
          flux_sect_sigma(k) = 0.0
          do 5 j=1,2
              flux_flag_tel(j,k) = 0
              flux_tel(j,k) = 0.0
5            sigma_tel2(j,k) = 0.0
c
      do 30 j=1,2
         do 30 k=1,8
            do 20 i=1,41
                if (local_sectored_flag(i)) then
                   if (mdata_flag(i,j,k).eq.-1) then   ! don't include saturated counts
                      flux_flag_tel(j,k) = -1
                      flux_tel(j,k) = 0.0
                      goto 30
                   end if
                   if (mdata_flag(i,j,k).eq.-2) then   ! don't include saturated counts
                      flux_flag_tel(j,k) = -2
                      flux_tel(j,k) = 0.0
                   end if
                   if (mdata_flag(i,j,k).eq.1) then  ! include only values measured in this mf
                        flux_flag_tel(j,k) = 1
                        flux_tel(j,k)=flux_tel(j,k)+
     1                             mdata(i,j,k)/(exposure(i,j)/8.0)
                        sigma_tel2(j,k)= sigma_tel2(j,k)+
     1                   mdata(i,j,k)/(exposure(i,j)*exposure(i,j)/64.0)  ! variance
                   end if
                end if
20         continue
30        continue
c
c           calculate total fluxes (telescope 1, 2 or 1+2)
c        **************************************
      if (local_tel_flag.eq.1) then   ! telescope 1 only
         do 100 k=1,8
            flux_sect_flag(k) = flux_flag_tel(1,k)
            flux_sect(k) = flux_tel(1,k)
100      flux_sect_sigma(k) = sqrt(sigma_tel2(1,k))
      end if
       if (local_tel_flag.eq.2) then   ! telescope 2 only
         do 200 k=1,8
            flux_sect_flag(k) = flux_flag_tel(2,k)
            flux_sect(k) = flux_tel(2,k)
200      flux_sect_sigma(k) = sqrt(sigma_tel2(2,k))
      end if
      if (local_tel_flag.eq.-1) then   ! both telescopes
         do 300 k=1,8
           Ntel = 0
           flux_sect_flag(k) = 0
           flux_sect(k) = 0.0
           flux_sect_sigma(k) = 0.0
            if (flux_flag_tel(1,k).eq.1) then
                Ntel  = Ntel+1
                flux_sect_flag(k) = 1
                flux_sect(k) = flux_tel(1,k)
                flux_sect_sigma(k) = sigma_tel2(1,k)
            end if
             if (flux_flag_tel(2,k).eq.1) then
                Ntel  = Ntel+1
                flux_sect_flag(k) = 1
                flux_sect(k) = flux_sect(k)+flux_tel(2,k)
                flux_sect_sigma(k) = flux_sect_sigma(k)+
     1              sigma_tel2(2,k)
            end if
            if (Ntel.gt.0) then
               flux_sect(k) = flux_sect(k)/Ntel
               flux_sect_sigma(k) = sqrt(flux_sect_sigma(k))/Ntel
            end if
             if ((flux_flag_tel(1,k).eq.-2).and.
     1        (flux_flag_tel(2,k).eq.-2))  then
                flux_sect(k) = 0.0
                flux_sect_flag(k) = -2
            end if
            if ((flux_flag_tel(1,k).eq.-1).or.
     1        (flux_flag_tel(2,k).eq.-1))  then
                flux_sect(k) = -1.0
                flux_sect_flag(k) = -1
            end if
300     continue
      end if
c
      return
      end  ! end  assign_sectored_fluxes
c
c
c
c
c**********************************************************
c
      subroutine assign_sectored_vse
     1               (local_tel_flag, vse_sect,
     1                vse_sect_sigma, vse_sect_flag)
c
c**********************************************************
c
c    calculates sectored vse rates and sigmas for 1 or both telescopes
c
c     11/10/95 by J. Dwyer
c             modified 1/8/96 by J. Dwyer divide sector exposures by 8.0
c
      include 'sdf_include.inc'   ! include type declarations
c
      real    delta_tvse  ! delta time over which disc rate is calculated
      real    vse_sect(8), vse_sect_sigma(8)
      integer   vse_sect_flag(8)
      integer  local_tel_flag   ! 1=use tel. 1 only, 2=use tel. 1 only,-1=use tel. 1 and tel. 2
      logical  local_sectored_flag(41)  ! .true.=include in sectored flux calculation, .false.=don't include
c
       do 5 k=1,8
          vse_sect(k) = 0.0
          vse_sect_sigma(k) = 0.0
5      continue
c
      delta_tvse = time_interval   ! delta time over which disc rate is calculated
      if ((stateab_flag.eq.0).or.(stateab_flag.eq.1))   ! if selecting only state A or B
     1      delta_tvse = time_interval/2.0
c
c
      do 10 k=1,8
         if (local_tel_flag.eq.1) then   ! telescope 1 only
            vse_sect_flag(k) = vse_flag(1,k)
            vse_sect(k) = vse(1,k)/(delta_tvse/8.0)
            vse_sect_sigma(k) = sqrt(1.0*vse(1,k))/(delta_tvse/8.0)
         end if
          if (local_tel_flag.eq.2) then   ! telescope 2 only
            vse_sect_flag(k) = vse_flag(2,k)
            vse_sect(k) = vse(2,k)/(delta_tvse/8.0)
            vse_sect_sigma(k) = sqrt(1.0*vse(2,k))/(delta_tvse/8.0)
         end if
         if (local_tel_flag.eq.-1) then   ! both telescopes
            if ((vse_flag(1,k).eq.1).and.
     1        (vse_flag(2,k).eq.1)) then
                 vse_sect_flag(k) = 1
                 vse_sect(k) = (vse(1,k)+
     1                  vse(2,k))/(delta_tvse/8.0)
                 vse_sect_sigma(k) = sqrt(1.0*
     1                    (vse(1,k)+vse(2,k)))/(delta_tvse/8.0)
            else
                vse_sect_flag(k) = -1
            end if
         end if
10      continue
c
      return
      end   ! end  assign_sectored_vse
c
c
c


text: #redirect Cross Cultural Center


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! 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.

!>    @brief evaluate bayesian quality function for iterate p_k
!>    @param[in] p_k parameter (iteration k)
!>    @param[in] p_apr apriori parameter value
!>    @param[in] ismpl local sample index
!>    @return scalar quality function value
!>    @details
!> evaluate bayesian quality function for iterate p_k\n
!> corresponding to the definition of Tarantola (2004) p.68\n
!>    fkt = [(d -forw)^T *C_dd^-1 *(d -forw) \n
!>          +reg_func() *(p_k -p_apr)^T *C_pp^-1 *(p_k -p_apr)]\n
!>    [1] = {([ndata]-[ndata])^T *[ndata^2] *([ndata] -[ndata])\n
!>          +[1] *([mpara] -[mpara])^T *[mpara^2] *([mpara] -[mpara]}\n
!>    tmp1 = d -forw\n
!>    [ndata] = [ndata] -[ndata]\n
!>    tmp2 = C_dd^-1 *tmp1\n
!>    [ndata] = [ndata^2] *[ndata]\n
!>    tmp3 = p_k -p_apr\n
!>    [mpara] = [mpara] -[mpara]\n
!>    tmp4 = C_pp^-1 *tmp3\n
!>    [mpara] = [mpara^2] *[mpara]\n
!>    fkt = (tmp1^T *tmp2 +reg_func() *tmp3^T *tmp4)\n
!>    [1] = ([ndata]^T *[ndata] +[1] *[mpara]^T *[mpara])\n
      DOUBLE PRECISION FUNCTION qfunc(p_k,p_apr,ismpl)
        use arrays
        use mod_genrl
        use mod_genrlc
        use mod_inverse
        use mod_data
        use mod_linfos
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j

        DOUBLE PRECISION p_k(mpara), p_apr(mpara), t_d, t_p
        DOUBLE PRECISION, ALLOCATABLE :: tmp(:,:)

        DOUBLE PRECISION wp2_ps, wd2_ps, reg_func
!              [mpara^2], [ndata^2]
        EXTERNAL wp2_ps, wd2_ps, reg_func

        ALLOCATE(tmp(max(mpara,ndata),4))

!     tmp1 = d -forw
!     [ndata] = [ndata] -[ndata]
        DO i = 1, ndata
          tmp(i,1) = ddata(i,cdd_pv) - sdata(i,ismpl)
        END DO

!     tmp2 = C_dd^-1 *tmp1
!     [ndata] = [ndata^2] *[ndata]
        DO j = 1, ndata
          tmp(j,2) = 0.0D0
          DO i = 1, ndata
            tmp(j,2) = tmp(j,2) + wd2_ps(j,i,ismpl)*tmp(i,1)
          END DO
        END DO

!     tmp3 = p_k -p_apr
!     [mpara] = [mpara] -[mpara]
        DO i = 1, mpara
          tmp(i,3) = p_k(i) - p_apr(i)
        END DO

!     tmp4 = C_pp^-1 *tmp3
!     [mpara] = [mpara^2] *[mpara]
        DO j = 1, mpara
          tmp(j,4) = 0.0D0
          DO i = 1, mpara
            tmp(j,4) = tmp(j,4) + wp2_ps(j,i,ismpl)*tmp(i,3)
          END DO
        END DO

!     fkt = (tmp1^T *tmp2 + tmp3^T *tmp4)
!     [1] = ([ndata]^T *[ndata] +[mpara]^T *[mpara])
        CALL s_ddot(ndata,tmp(1,1),tmp(1,2),t_d)
        CALL s_ddot(mpara,tmp(1,3),tmp(1,4),t_p)
        qfunc = 0.5D0*(t_d+reg_func()*t_p)

        IF (linfos(2)>=2) WRITE(*,'(A,3(e15.9,A))') &
          '  objective function=', qfunc, ', ( d:', t_d, ', p:', t_p, ' )'

        rms_data = t_d
        rms_para = t_p
        DEALLOCATE(tmp)

        RETURN
      END


text:
Benjamin Catabas

 Relationship with the Davis Community
◙   UC Davis Student of International Relations, French, & Education

†   Active Member of the UCD Newman Center Newman Catholic Student Community and St. James Catholic Church, Davis

☺ UC Davis Campus Ambassador: UC Davis Tours Tour Guide & Visitors Service Desk Attendant

♠   Frequent Patron of local restaurants, grocery stores, & tea cafés.

20070609 14:30:22 nbsp hey ben, was the mural painted in Tercero D Building a map of the world? xie xie ni. Users/JessicaRockwell

20070609 15:25:56 nbsp why, yes it was... the one we painted... it was the map of the world with every countrys flag painted within each border. Users/BenjaminCatabas
   yay! i remembered! yeah, you all put so much time into that. friday nights in the stairwell, painting flags. good times.   )    Users/JessicaRockwell

20070614 21:14:20 nbsp Nice additions to the Urban Legends entry! Users/JabberWokky

20070819 18:35:39 nbsp While it probably would be a good idea to rewrite the Gunrock entry, stealing material that somebody else wrote isnt legal (and is kind of lazy).  Right under the edit screen it asks you not to do that.  I agree that the Gunrock entry is more pseudowit and less informative than it should be... so if you know a bit about the topic itself, please add it (in your own words). Users/JabberWokky

20070819 21:41:04 nbsp Sorry about this... actually, I did not make those changes to the gunrock page. I log into daviswiki using several different computers & I never log out. Looking at the IP address of that edit, I conclude that somebody has edited the gunrock content from a computer at work, where Im still logged in. I totally agree with the wrongness of stealing material from others.. thanks for the heads up! Users/BenjaminCatabas

20070820 03:16:23 nbsp No worries.  It was pretty clearly intended well, and swiping material from a public university site is likely the lowest level of theft.  Keeping everything free to use and open on the wiki is one of those things to keeping the wiki useful for everyone.  If your account was going to be used, at least it was used with good intent. Users/JabberWokky


text:
module util
    implicit none

    contains

    function arange(n) result(a)
    implicit none
    integer, intent(in) :: n
    integer :: i, a(n)
    do i = 1, n
        a(i) = i
    end do
    end function

    subroutine swap_int(a, b)
    implicit none
    integer, intent(inout) :: a, b
    integer :: c
    c = a
    a = b
    b = c
    end subroutine

    subroutine write_matrix(unit, a)
    implicit none
    integer, intent(in) :: unit, a(:, :)
    integer :: i
    character(len=7) :: fmt
    write(unit, '("shape",1x,2i6)') shape(a)
    write(fmt, '("(",i3,"i4)")') size(a, 2)
    do i = 1, size(a, 1)
        write(unit, fmt) a(i, :)
    end do
    end subroutine

end module

module combinations
    implicit none

    contains

    logical &
    function lex_gt(a, b)
    implicit none
    integer, intent(in) :: a(:), b(:)
    integer :: i
    lex_gt = .false.
    do i = 1, size(a)
        if (a(i) < b(i)) then
            lex_gt = .true.
            return
        else if (a(i) > b(i)) then
            return
        end if
    end do
    end function

    subroutine symmetries(a, c)
    implicit none
    integer, intent(in) :: a(:, :)
    integer, intent(out) :: c
    integer :: s(size(a, 1)), o(size(a, 1))
    integer :: m, n, k, sk, unit

    m = size(a, 1)
    n = size(a, 1)
    o = 0
    s = 0
    k = 1
    c = 0
    open(newunit=unit, file='symmetries.txt')

100 sk = s(k) + 1
    if (sk > n) goto 102
    s(k) = sk
    if (o(sk) > 0) goto 100
    if (any(a(1:k-1, k) /= a(s(1:k-1), s(k)))) goto 100
    o(sk) = 1
    if (k < m) then
        k = k + 1
        goto 100
    else
        c = c + 1
        write(unit, '(i12,6x,100i4)') c, s
    end if

101 o(s(k)) = 0
    goto 100

102 s(k) = 0
    k = k - 1
    if (k == 0) goto 103
    goto 101

103 close(unit)
    write(*, *) 'symmetries:', c
    end

    subroutine embeddings(t, a, isym)
    implicit none
    integer, intent(in) :: t(:), a(:, :)
    logical, intent(in), optional :: isym
    logical :: use_sym = .true.
    integer, allocatable :: sym(:, :)
    integer :: s(size(t)), o(size(a, 1))
    integer :: m, n, j, k, sk, c, nsym, unit

    m = size(t)
    n = size(a, 1)

    if (present(isym)) use_sym = isym

    if (use_sym) then
        call symmetries(a, nsym)
        allocate(sym(nsym, n))
        open(newunit=unit, file='symmetries.txt')
        do k = 1, nsym
            read(unit, *) j, sym(k, :)
        end do
        close(unit)
    end if

    o = 0
    s = 0
    k = 1
    c = 0
    open(newunit=unit, file='embeddings.txt')

100 sk = s(k) + 1
    if (sk > n) goto 102
    s(k) = sk
    if (o(sk) > 0) goto 100
    do j = 1, k - 1
       if (t(j) == t(k) .and. sk < s(j)) goto 100
    end do
    if (use_sym) then
        do j = 1, nsym
            if (lex_gt(sym(j, s(1:k)), s(1:k))) goto 100
        end do
    end if
    o(sk) = 1
    if (k < m) then
        k = k + 1
        goto 100
    else
        c = c + 1
        write(unit, '(i12,6x,100i4)') c, s
    end if

101 o(s(k)) = 0
    goto 100

102 s(k) = 0
    k = k - 1
    if (k == 0) goto 103
    goto 101

103 close(unit)
    write(*, *) 'embeddings:', c
    end

end module

module combinations_tests
    use util
    use combinations
    implicit none

    contains

    subroutine chain(n, a)
    implicit none
    integer, intent(in) :: n
    integer, intent(inout) :: a(n, n)
    integer :: i
    a = 0
    do i = 1, n-1
        a(i, i+1) = 1
        a(i+1, i) = 1
    end do
    end

    subroutine ring(n, a)
    implicit none
    integer, intent(in) :: n
    integer, intent(inout) :: a(n, n)
    call chain(n, a)
    a(1, n) = 1
    a(n, 1) = 1
    end

    subroutine test_symmetries
    use util, only: write_matrix
    use combinations, only: embeddings
    implicit none
    integer, parameter :: n=5
    integer :: a(n, n)
    integer :: nsym
    call chain(n, a)
    call write_matrix(0, a)
    call symmetries(a, nsym)
    end

    subroutine test_embeddings
    use combinations, only: embeddings
    implicit none
    integer, parameter :: n=4, m=2
    integer :: t(m), a(n, n), i
    t = arange(m)
    !t = 1
    call ring(n, a)
    call embeddings(t, a)
    end

    subroutine test_all
    call test_symmetries
    call test_embeddings
    end subroutine

end module

    program main
    use combinations
    use util
    implicit none
    integer :: n, i, j, k, l, unit, a, t, io
    allocatable :: a(:, :), t(:), io(:)
    open(newunit=unit, file='graph.txt')
    read(unit, *) n
    allocate(a(n, n), t(n), io(n))
    a = 0
    do i = 1, n
        read(unit, *) j, t(i), k, (io(l), l=1, k)
        if (i /= j) then
            print*, 'input error!'
            stop
        end if
        do l = 1, k
            a(i, io(l)) = 1
            a(io(l), i) = 1
        end do
    end do
    close(unit)
    !call write_matrix(0, a)
    !print*, t
    call embeddings(t, a)
    end program

max_stars_repo_name: amirhajibabaei/fortran-examples

max_issues_repo_name: amirhajibabaei/fortran-examples

max_forks_repo_name: amirhajibabaei/fortran-examples


! ##################################################################################################################################
! 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 BD_SLOAD0 ( CARD, NUM_PAIRS )

! Processes SLOAD Bulk Data Cards to count the number of SLOAD points/force mags on one entry

!   SETID, scalar point, load mag

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  WRT_LOG, F04
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, IERRFL, JCARD_LEN, JF
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_SLOAD0_BEGEND

      USE BD_SLOAD0_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BD_SLOAD0'
      CHARACTER(LEN=*),INTENT(IN)     :: CARD                ! A Bulk Data card
      CHARACTER(LEN=JCARD_LEN)        :: JCARD(10)           ! The 10 fields of characters making up CARD

      INTEGER(LONG), INTENT(OUT)      :: NUM_PAIRS           ! Number of pairs of SPOINT/force MAG on a SLOAD entry (can be up to 3)
      INTEGER(LONG)                   :: I                   ! DO loop index
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_SLOAD0_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! SLOAD Bulk Data Card routine

!   FIELD   ITEM           VARIABLE
!   -----   ------------   -------------
!    2      Load set ID    SLOAD_SIDS(I), SETID
!    3      Scalar point   GRID_NO
!    4      Load magnitude FMAG

! Fields (3,4) can be repeated in fields (5,6) and (7,8)

! Make JCARD from CARD

      CALL MKJCARD ( SUBR_NAME, CARD, JCARD )

      NUM_PAIRS = 0
      DO I=1,3
         IF (JCARD(3+2*(I-1))(1:) /= ' ') THEN
            NUM_PAIRS = NUM_PAIRS + 1
         ENDIF
      ENDDO

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

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

      END SUBROUTINE BD_SLOAD0