text:
C$Id: testMemSizef.f 466 2005-10-03 10:37:50Z vmp $
      program testmemsize
      external countsetbits
      integer countsetbits
      integer i, j, ti(10), ci0(10), ci1(10)
      integer*8 i8, ti8(10), ci80(10), ci81(10)
      real r, tr(10), cr0(10), cr1(10)
      real*4 r4, tr4(10), cr40(10), cr41(10)
      real*8 r8, tr8(10), cr80(10), cr81(10)
c      real*16 r16, tr16(10), cr160(10), cr161(10)
      double precision d, td(10), cd0(10), cd1(10)
      complex cx, tcx(10), ccx0(10), ccx1(10)
      complex*8 cx8, tcx8(10), ccx80(10), ccx81(10)
      complex*16 cx16, tcx16(10), ccx160(10), ccx161(10)
      double complex cxd, tcxd(10), ccxd0(10), ccxd1(10)
      logical l, tl(10), cl0(10), cl1(10)
      character c, tc(10), cc0(10), cc1(10)
c
      integer tci0(10000)
      common /cbuf0/ tci0
      equivalence (tci0(1),ci0(1))
      equivalence (tci0(1),ci80(1))
      equivalence (tci0(1),cr0(1))
      equivalence (tci0(1),cr40(1))
      equivalence (tci0(1),cr80(1))
c      equivalence (tci0(1),cr160(1))
      equivalence (tci0(1),cd0(1))
      equivalence (tci0(1),ccx0(1))
      equivalence (tci0(1),ccx80(1))
      equivalence (tci0(1),ccx160(1))
      equivalence (tci0(1),ccxd0(1))
      equivalence (tci0(1),cl0(1))
      equivalence (tci0(1),cc0(1))
c
      integer tci1(10000)
      common /cbuf1/ tci1
      equivalence (tci1(1),ci1(1))
      equivalence (tci1(1),ci81(1))
      equivalence (tci1(1),cr1(1))
      equivalence (tci1(1),cr41(1))
      equivalence (tci1(1),cr81(1))
c      equivalence (tci1(1),cr161(1))
      equivalence (tci1(1),cd1(1))
      equivalence (tci1(1),ccx1(1))
      equivalence (tci1(1),ccx81(1))
      equivalence (tci1(1),ccx161(1))
      equivalence (tci1(1),ccxd1(1))
      equivalence (tci1(1),cl1(1))
      equivalence (tci1(1),cc1(1))
c
      integer nb
c
      do i = 1,10
         ti(i) = (i-1)*i*i*i*i*i*i*i*i*i
         ti8(i) = (i-1)*i*i*i*i*i*i*i*i*i
         tr(i) = (i-1)*i*exp(2.0*i+0.001)
         tr4(i) = (i-1)*i*exp(2.0*i+0.001)
         tr8(i) = (i-1)*i*exp(2.0*i+0.001)
c         tr16(i) = (i-1)*i*exp(2.0*i+0.001)
         td(i) = (i-1)*i*exp(2.d0*i+1.d-3)
         tcx(i) = cmplx(tr(i),tr(i)+12.3)
         tcx8(i) = cmplx(tr8(i),tr(i)+12.3)
         tcx16(i) = cmplx(tr8(i),tr(i)+12.3)
         tcxd(i) = cmplx(td(i),td(i)+12.3d0)
         tl(i) = (((i/3)*3).eq.i)
      enddo
      tc(1) = 'a'
      tc(2) = 'b'
      tc(3) = 'c'
      tc(4) = 'd'
      tc(5) = 'e'
      tc(6) = 'f'
      tc(7) = 'g'
      tc(8) = 'h'
      tc(9) = 'i'
      tc(10) = 'j'

      i = 12345
      i8 = 12345
      r = sqrt(12345.0)
      r4 = sqrt(12345.0)
      r8 = sqrt(12345.0)
c      r16 = sqrt(12345.0)
      d = sqrt(12345.d0)
      cx = cmplx(-12.34,56.78)
      cx8 = cmplx(-12.34,56.78)
      cx16 = cmplx(-12.34,56.78)
      cxd = cmplx(-12.34d0,56.78d0)
      l = .true.
      c = 'z'
c
      call allzero(tci0(1))
      call allones(tci1(1))
c      call displaybits(tci0(1), 3)
c      call displaybits(tci1(1), 3)
      print *, ' '
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         ci0(j) = ti(j+2)
         ci1(j) = ti(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'INTEGER:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         ci80(j) = ti8(j+2)
         ci81(j) = ti8(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'INTEGER*8:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         cr0(j) = tr(j+2)
         cr1(j) = tr(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'REAL:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         cr40(j) = tr4(j+2)
         cr41(j) = tr4(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'REAL*4:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         cr80(j) = tr8(j+2)
         cr81(j) = tr8(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'REAL*8:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         cd0(j) = td(j+2)
         cd1(j) = td(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'DOUBLE PRECISION:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         ccx0(j) = tcx(j+2)
         ccx1(j) = tcx(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'COMPLEX:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         ccx80(j) = tcx8(j+2)
         ccx81(j) = tcx8(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'COMPLEX*8:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         ccx160(j) = tcx16(j+2)
         ccx161(j) = tcx16(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'COMPLEX*16:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         ccxd0(j) = tcxd(j+2)
         ccxd1(j) = tcxd(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'DOUBLE COMPLEX:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         cl0(j) = tl(j+2)
         cl1(j) = tl(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'LOGICAL:', nb, ' bytes'
c
      call allzero(tci0(1))
      call allones(tci1(1))
      do j=1,8
         cc0(j) = tc(j+2)
         cc1(j) = tc(j+2)
      enddo
      nb = countsetbits(tci0(1), tci1(1), 8, 100)
      print 100, 'CHARACTER:', nb, ' bytes'
c
 100  format(a18,i3,a6)
      end

max_stars_repo_name: MateusAraujoBorges/civl

max_issues_repo_name: MateusAraujoBorges/civl

max_forks_repo_name: MateusAraujoBorges/civl


REAL*8 FUNCTION SPHBRU(N,X)
  !IMPLICIT REAL*8(A-H,O-Z)
  IMPLICIT NONE
  INTEGER, intent(in) :: n
  REAL*8,  intent(in) :: x
  ! locals
  REAL*8, PARAMETER :: TOL = 1.D-14
  REAL*8  :: BB, DBB, FAC, FACD
  INTEGER :: i, j
  I=0
  BB=1.D0
  DO
     I=I+1
     FAC=1.D0
     DO J=1,I
        FAC=FAC/DFLOAT(J)
     ENDDO
     FACD=1.D0
     DO J=1,I
        FACD=FACD/(2.D0*DFLOAT(N)+3.D0+2.D0*DFLOAT(J-1))
     ENDDO
     DBB=(X*X/2.D0)**I*FAC*FACD
     BB=BB+DBB*(-1.D0)**I
     IF(DBB.LT.TOL) EXIT
  ENDDO
  SPHBRU=BB
  RETURN
END FUNCTION SPHBRU


text:
      SUBROUTINE DT_RNDMST(Na1,Na2,Na3,Nb1)

      IMPLICIT NONE
      INTEGER ii1 , ii2 , ma1 , ma2 , ma3 , mat , mb1 , Na1 , Na2 ,
     &        Na3 , Nb1
      DOUBLE PRECISION s , t
      SAVE

C random number generator
      INCLUDE 'inc/dtrand'

      ma1 = Na1
      ma2 = Na2
      ma3 = Na3
      mb1 = Nb1
      I = 97
      J = 33
      DO ii2 = 1 , 97
         s = 0
         t = 0.5D0
         DO ii1 = 1 , 24
            mat = MOD(MOD(ma1*ma2,179)*ma3,179)
            ma1 = ma2
            ma2 = ma3
            ma3 = mat
            mb1 = MOD(53*mb1+1,169)
            IF ( MOD(mb1*mat,64).GE.32 ) s = s + t
            t = 0.5D0*t
         END DO
         U(ii2) = s
      END DO
      C = 362436.0D0/16777216.0D0
      CD = 7654321.0D0/16777216.0D0
      CM = 16777213.0D0/16777216.0D0
      END SUBROUTINE

ext: f

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


program entropy
    USE yepLibrary, ONLY: yepLibrary_Init, yepLibrary_Release, yepLibrary_GetTimerFrequency, yepLibrary_GetTimerTicks
    USE yepCore, ONLY: yepCore_DotProduct_V64fV64f_S64f
    USE yepMath, ONLY: yepMath_Log_V64f_V64f
    USE ISO_C_BINDING, ONLY : C_SIZE_T, C_INT64_T, C_INT, C_DOUBLE
    implicit none
    INTEGER(C_SIZE_T), PARAMETER :: n = 1024*1024*16
    INTEGER(C_INT64_T) :: ticksStart, ticksEnd, frequency
    real(C_DOUBLE) :: p(n), logP(n), H
    integer(C_INT) :: s
    s = yepLibrary_Init()
    s = yepLibrary_GetTimerFrequency(frequency)
    call RANDOM_NUMBER(p)

    s = yepLibrary_GetTimerTicks(ticksStart)
    s = yepMath_Log_V64f_V64f(p, logP, n)
    s = yepCore_DotProduct_V64fV64f_S64f(p, logP, H, n)
    H = -H
    s = yepLibrary_GetTimerTicks(ticksEnd)
    print*, "Yeppp!"
    print*, "    Entropy =", H
    print*, "    Time =", REAL(ticksEnd - ticksStart) / REAL(frequency)

    s = yepLibrary_GetTimerTicks(ticksStart)
    H = -DOT_PRODUCT(p, LOG(p))
    s = yepLibrary_GetTimerTicks(ticksEnd)
    print*, "Naive"
    print*, "    Entropy =", H
    print*, "    Time =", REAL(ticksEnd - ticksStart) / REAL(frequency)
    s = yepLibrary_Release()
end


text:
type FwEventIdType
type FwChanIdType
type FwOpcodeType
type FwPrmIdType
type NATIVE_INT_TYPE
type NATIVE_UINT_TYPE
type POINTER_CAST

ext: fpp


text:
MODULE symmetry_module

  use parallel_module, np_2d => node_partition

  implicit none

  PRIVATE
  PUBLIC :: init_symmetry, read_symmetry, prep_symmetry, sym_rho, sym_force &
           ,nsym,isymmetry,rgb,construct_matrix_symmetry

  integer :: isymmetry
  character(30) :: file_symdat

  integer :: nsym, nnp
  integer,allocatable :: rga(:,:,:), pga(:,:)
  real(8),allocatable :: rgb(:,:,:), pgb(:,:)

  integer :: a1b,b1b,a2b,b2b,a3b,b3b
  integer :: ML,ML1,ML2,ML3
  real(8) :: dV
  real(8) :: aa(3,3), bb(3,3)

CONTAINS

  SUBROUTINE read_symmetry( rank, unit )
    implicit none
    integer,intent(IN) :: rank,unit
    integer :: i,ierr
    character(4) :: cbuf,ckey

    isymmetry   = 0
    file_symdat = ""

    if ( rank == 0 ) then
       rewind unit
       do i=1,10000
          read(unit,*,END=999) cbuf
          call convert_capital(cbuf,ckey)
          if ( ckey == "ISYM" ) then
             backspace(unit)
             read(unit,*) cbuf, isymmetry, file_symdat
             exit
          end if
       end do
999    continue
       write(*,*) "isymmetry=",isymmetry
       write(*,*) "file_symdat=",file_symdat
    end if

    call MPI_BCAST(isymmetry,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(file_symdat,30,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)

  END SUBROUTINE read_symmetry

  SUBROUTINE init_symmetry(Ngrid,dV_in,aa_in,bb_in,MI,Kion,asi)
    implicit none
    integer,intent(IN) :: Ngrid(0:3)
    real(8),intent(IN) :: dV_in,aa_in(3,3),bb_in(3,3)
    integer,intent(IN) :: MI,Kion(MI)
    real(8),intent(IN) :: asi(3,MI)

    if ( isymmetry == 0 ) return

    call write_border( 0," init_symmetry(start)")

    ML  = Ngrid(0)
    ML1 = Ngrid(1)
    ML2 = Ngrid(2)
    ML3 = Ngrid(3)
    dV  = dV_in
    aa  = aa_in
    bb  = bb_in

    call input_symdat( MI, Kion, asi )

    call input_symdat_2( MI, Kion, asi )
!    call input_symdat_3( MI, Kion, asi )
!    call input_symdat_4( MI, Kion, asi )

    call write_border( 0," init_symmetry(end)")

  END SUBROUTINE init_symmetry

!--------1---------2---------3---------4---------5---------6---------7--
!
! "isymmetry>0" assume the atomic coordinates are given within [0,1]
! "isymmetry<0" assume the atomic coordinates are given within [-0.5,0.5]

  SUBROUTINE input_symdat( MI, Kion, asi )
    implicit none
    integer,intent(IN) :: MI
    integer,intent(IN) :: Kion(MI)
    real(8),intent(IN) :: asi(3,MI)
    integer,parameter :: u=2
    integer :: i,i1,i2,j,k,n,ierr,flag
    integer :: A11,A12,A13,A21,A22,A23,A31,A32,A33,gcounter
    integer :: AAA(9,19683)
    real(8) :: RR0(3),RR1(3),RR2(3),RRR
    real(8) :: R1q,R2q,R3q,R1p,R2p,R3p,R1s,R2s,R3s
    real(8) :: tmp2(3,3,2)
    real(8) :: c1
    real(8),allocatable :: XX(:,:), pga_tmp(:,:)
    character(12) :: label_sym

    if ( .not.(abs(isymmetry)==1 .or. abs(isymmetry)==2 .or. &
               abs(isymmetry)==3) )  return

    call write_border( 0, " input_symdat(start)" )

    if ( disp_switch_parallel ) then
       write(*,*) "isymmetry=",isymmetry
    end if

    nnp  = 1
    nsym = 0

!
! --- Read symmetry operation matrix ---
!
    if ( myrank == 0 ) then

       if ( abs(isymmetry) <= 2 ) open(u,file=file_symdat,status='old')

       select case( isymmetry )
!(1)
       case( 1, -1 )

          read(u,*) nsym, nnp
          allocate( rga(3,3,nsym),pga(3,nsym) ) ; rga=0 ; pga=0
          allocate( rgb(3,3,nsym),pgb(3,nsym) ) ; rgb=0 ; pgb=0
          do n=1,nsym
             read(u,*) ( (rga(i,j,n),j=1,3),i=1,3 ),pga(:,n)
          end do
!(2)
       case( 2, -2 )

          read(u,*) label_sym, nsym
          write(*,*) "label_sym= ",label_sym
          write(*,*) "nsym     = ",nsym
          allocate( rga(3,3,nsym),pga(3,nsym) ) ; rga=0 ; pga=0
          allocate( rgb(3,3,nsym),pgb(3,nsym) ) ; rgb=0 ; pgb=0
          do i=1,nsym,2
             j=min(i+1,nsym)
             n=j-i+1
             read(u,*) (((tmp2(i2,i1,j),i1=1,3),i2=1,3),j=1,n)
             do j=1,n
                rga(1:3,1:3,i+j-1)=nint( tmp2(1:3,1:3,j) )
             end do
          end do

          allocate( pga_tmp(3,nsym) ) ; pga_tmp=0.d0
          do i=1,nsym,3
             j=min(i+2,nsym)
             read(u,*) pga_tmp(1:3,i:j)
          end do
          c1=1.0d10
          do i=1,nsym
             do j=1,3
                if ( abs(pga_tmp(j,i)) > 1.0d-10 ) then
                   c1=min( c1, abs(pga_tmp(j,i)) )
                end if
             end do
          end do
          pga(1:3,1:nsym)=nint( pga_tmp(1:3,1:nsym)/c1 )
          nnp=nint( 1.0d0/c1 )
          deallocate( pga_tmp )
!(3)
       case( 3, -3 )

          allocate( XX(3,MI) ) ; XX=0.0d0

          do i=1,MI
             XX(1:3,i)=asi(1:3,i)
             do j=1,3
                do  while ( XX(j,i) >  0.5d0 )
                   XX(j,i)=XX(j,i)-1.0d0
                end do
                do  while ( XX(j,i) <= -0.5d0 )
                   XX(j,i)=XX(j,i)+1.0d0
                end do
             end do
          end do ! i

          do A11=-1,1
          do A12=-1,1
          do A13=-1,1
          do A21=-1,1
          do A22=-1,1
          do A23=-1,1
          do A31=-1,1
          do A32=-1,1
          do A33=-1,1

             gcounter=0

             do i=1,MI

                RR0(1:3) = XX(1:3,i)

                flag = 0

                do j=1,MI

                   if ( Kion(i) /= Kion(j) ) cycle

                   RR1(1:3) = XX(1:3,j)

                   RR2(1) = A11*RR1(1) + A12*RR1(2) + A13*RR1(3)
                   RR2(2) = A21*RR1(1) + A22*RR1(2) + A23*RR1(3)
                   RR2(3) = A31*RR1(1) + A32*RR1(2) + A33*RR1(3)

                   do k=1,3
                      do  while ( RR2(k) >  0.5d0 )
                         RR2(k)=RR2(k)-1.0d0
                      end do
                      do while ( RR2(k) <= -0.5d0 )
                         RR2(k)=RR2(k)+1.0d0
                      end do
                   end do ! k

                   RRR = sqrt( (RR2(1)-RR0(1))**2 &
                             + (RR2(2)-RR0(2))**2 + (RR2(3)-RR0(3))**2 )
                   if ( RRR < 1.d-6 ) then
                      flag=1
                   end if

                end do ! j

                gcounter = gcounter + flag

             end do ! i

             if ( gcounter == MI ) then

                R1q = sqrt( AA(1,1)**2 + AA(2,1)**2 + AA(3,1)**2 )
                R2q = sqrt( AA(1,2)**2 + AA(2,2)**2 + AA(3,2)**2 )
                R3q = sqrt( AA(1,3)**2 + AA(2,3)**2 + AA(3,3)**2 )

                R1s = real(A11)*1 + real(A12)*0 + real(A13)*0
                R2s = real(A21)*1 + real(A22)*0 + real(A23)*0
                R3s = real(A31)*1 + real(A32)*0 + real(A33)*0
                R1p = sqrt( (AA(1,1)*R1s + AA(1,2)*R2s + AA(1,3)*R3s)**2 &
                          + (AA(2,1)*R1s + AA(2,2)*R2s + AA(2,3)*R3s)**2 &
                          + (AA(3,1)*R1s + AA(3,2)*R2s + AA(3,3)*R3s)**2   )

                R1s = real(A11)*0 + real(A12)*1 + real(A13)*0
                R2s = real(A21)*0 + real(A22)*1 + real(A23)*0
                R3s = real(A31)*0 + real(A32)*1 + real(A33)*0
                R2p = sqrt( (AA(1,1)*R1s + AA(1,2)*R2s + AA(1,3)*R3s)**2 &
                          + (AA(2,1)*R1s + AA(2,2)*R2s + AA(2,3)*R3s)**2 &
                          + (AA(3,1)*R1s + AA(3,2)*R2s + AA(3,3)*R3s)**2   )

                R1s = real(A11)*0 + real(A12)*0 + real(A13)*1
                R2s = real(A21)*0 + real(A22)*0 + real(A23)*1
                R3s = real(A31)*0 + real(A32)*0 + real(A33)*1
                R3p = sqrt( (AA(1,1)*R1s + AA(1,2)*R2s + AA(1,3)*R3s)**2 &
                          + (AA(2,1)*R1s + AA(2,2)*R2s + AA(2,3)*R3s)**2 &
                          + (AA(3,1)*R1s + AA(3,2)*R2s + AA(3,3)*R3s)**2   )

                if ( abs(R1q-R1p) <= 1.d-6 .and. &
                     abs(R2q-R2p) <= 1.d-6 .and. &
                     abs(R3q-R3p) <= 1.d-6        ) then
                   nsym=nsym+1
                   AAA(1,nsym)=A11
                   AAA(2,nsym)=A12
                   AAA(3,nsym)=A13
                   AAA(4,nsym)=A21
                   AAA(5,nsym)=A22
                   AAA(6,nsym)=A23
                   AAA(7,nsym)=A31
                   AAA(8,nsym)=A32
                   AAA(9,nsym)=A33
                end if

             end if

          end do
          end do
          end do
          end do
          end do
          end do
          end do
          end do
          end do

          deallocate( XX )

          write(*,*) "------------------------------------------------------"
          write(*,'(I0,42X,I0)')  nsym, 1
          do i=1,nsym
             write(*,'(3(3I3,3X),6X,3I3)')  AAA(1:9,i),0,0,0
          end do !i
          write(*,*) "------------------------------------------------------"
          write(*,*) "BE CAREFUL!! Bubun heishin & Rasen are NOT considered."
          write(*,*) 'Intentional STOP, isymmetry=', isymmetry
          stop "stop@input_symdat"
!
!---
!
       end select

       close(u)

    end if

    call mpi_bcast(nsym,1,mpi_integer,0,mpi_comm_world,ierr)
    call mpi_bcast(nnp ,1,mpi_integer,0,mpi_comm_world,ierr)

    if ( myrank /= 0 ) then
       allocate( rga(3,3,nsym),pga(3,nsym) ) ; rga=0 ; pga=0
       allocate( rgb(3,3,nsym),pgb(3,nsym) ) ; rgb=0 ; pgb=0
    end if

    call mpi_bcast(rga,9*nsym,mpi_integer,0,mpi_comm_world,ierr)
    call mpi_bcast(pga,3*nsym,mpi_integer,0,mpi_comm_world,ierr)

    if ( disp_switch_parallel ) then
       write(*,*) "nsym,nnp=",nsym,nnp
    end if

    call write_border( 0, " input_symdat(end)" )

  END SUBROUTINE input_symdat

  SUBROUTINE chk_grp( rga, aa, bb )
    implicit none
    integer,intent(IN) :: rga(:,:,:)
    real(8),intent(IN) :: aa(3,3), bb(3,3)
    real(8) :: aa_inv(3,3),tmp0(3,3),tmp1(3,3)
    real(8),allocatable :: rtmp(:,:,:)
    integer :: n,ns
    ns=size(rga,3)
    allocate( rtmp(3,3,ns) ) ; rtmp=0.0d0
    aa_inv(:,:) = transpose( bb )/(2.0d0*acos(-1.0d0))
    do n=1,ns
       tmp0(:,:) = rga(:,:,n)
       tmp1(:,:) = matmul( tmp0, aa_inv )
       rtmp(:,:,n) = matmul( aa, tmp1 )
    end do
    do n=1,ns
       tmp0(:,:) = transpose( rtmp(:,:,n) )
       tmp1(:,:) = matmul( rtmp(:,:,n), tmp0 )
       write(*,*) n,sum(abs(tmp1))
    end do
    deallocate( rtmp )
    stop "stop@chk_grp"
  END SUBROUTINE chk_grp

  SUBROUTINE input_symdat_2( MI, Kion, asi )
    implicit none
    integer,intent(IN) :: MI, Kion(MI)
    real(8),intent(IN) :: asi(3,MI)
    real(8) :: pi2,c1,r0
    real(8) :: tmp(3,3),tmp1(3,3),tmp2(3,3,2),tmp0(3,3),tmp4(3)
    real(8) :: bb_inv(3,3),aa_inv(3,3)
    real(8),allocatable :: Rra(:,:)
    integer :: mtmp1(3,3),mtmp2(3,3),itmp1(3),itmp2(3)
    integer :: i,i1,i2,i3,isym,j,loop,m,m1,m2,n
    integer,allocatable :: itmp(:)

    call write_border( 0, " input_symdat_2(start)" )

    pi2 = 2.0d0*acos(-1.0d0)

! --- make rgb (symmetry operation matrix in bb-representation) ---
!
    bb_inv(1:3,1:3) = transpose( aa(1:3,1:3) )/pi2
    aa_inv(1:3,1:3) = transpose( bb(1:3,1:3) )/pi2
    tmp2(1:3,1:3,1) = matmul( bb_inv(1:3,1:3), aa(1:3,1:3) )
    tmp2(1:3,1:3,2) = matmul( aa_inv(1:3,1:3), bb(1:3,1:3) )
    do n=1,nsym
       tmp0(1:3,1:3)  = matmul( tmp2(1:3,1:3,1), rga(1:3,1:3,n) )
       tmp1(1:3,1:3)  = matmul( tmp0(1:3,1:3), tmp2(1:3,1:3,2) )
       rgb(1:3,1:3,n) = tmp1(1:3,1:3)
    end do

!
! --- Check the symmetry operations ---
!

! tr(rgb)=(rga)^-1
!
    do n=1,nsym
       tmp0(1:3,1:3) = transpose( rgb(1:3,1:3,n) )
       tmp1(1:3,1:3) = matmul( tmp0(1:3,1:3), rga(1:3,1:3,n) )
       itmp1(1)=nint( tmp1(1,1) )
       itmp1(2)=nint( tmp1(2,2) )
       itmp1(3)=nint( tmp1(3,3) )
       i1=sum(abs(nint(tmp1)))
       if ( .not.( all(itmp1==1) .and. i1==3 ) ) then
          write(*,*) "tr(rgb) is not the inverse of rga!"
          write(*,*) n,i1
          write(*,*) tmp1
          goto 900
       end if
    end do

! check duplication
!
    allocate( itmp(nsym) )
    do n=1,nsym
       j=0
       itmp=0
       do m=1,nsym
          mtmp1(1:3,1:3)=rga(1:3,1:3,n)-rga(1:3,1:3,m)
          itmp1(1:3)=pga(1:3,n)-pga(1:3,m)
          if ( all(mtmp1==0) .and. all(itmp1==0) ) then
             j=j+1
             itmp(j)=m
          end if
       end do
       if ( j /= 1 ) then
          write(*,*) "symmetry error (duplication)",j,n
          write(*,'(1x,i5,3x,9i3,3x,3i3)') n,rga(1:3,1:3,n),pga(1:3,n)
          do i=1,j
             write(*,'(1x,i5,3x,9i3,3x,3i3)') &
                  itmp(i),rga(1:3,1:3,itmp(i)),pga(1:3,itmp(i))
          end do
          goto 900
       end if
    end do
    deallocate( itmp )

! existence of the unit operator
!
    m1=0
    do n=1,nsym
       itmp1(1)=rga(1,1,n)
       itmp1(2)=rga(2,2,n)
       itmp1(3)=rga(3,3,n)
       i=sum(abs(rga(1:3,1:3,n)))
       if ( i==3 .and. all(itmp1==1) .and. all(pga(1:3,n)==0) ) then
          m1=m1+1
       end if
    end do
    if ( m1 /= 1 ) then
       write(*,*) "symmetry error (no unit operator)",m1
       goto 900
    end if

! (rga_1)(rga_2)=(rga)
!
    do n=1,nsym
       do m=1,nsym
          mtmp1(1:3,1:3)=matmul( rga(1:3,1:3,m),rga(1:3,1:3,n) )
          itmp1(1:3)=matmul( rga(1:3,1:3,m),pga(1:3,n) )+pga(1:3,m)
          do j=1,3
             do loop=1,10000
                if ( itmp1(j)<0 ) then
                   itmp1(j)=itmp1(j)+nnp
                else if ( itmp1(j)>=nnp ) then
                   itmp1(j)=itmp1(j)-nnp
                else
                   exit
                end if
             end do
          end do
          m2=0
          do i=1,nsym
             mtmp2(1:3,1:3)=mtmp1(1:3,1:3)-rga(1:3,1:3,i)
             itmp2(1:3)=itmp1(1:3)-pga(1:3,i)
             if ( all(mtmp2==0) .and. all(itmp2==0) ) then
                m2=m2+1
             end if
          end do ! i
          if ( m2 /= 1 ) then
             write(*,*) "symmetry error (this is not group)",m2,m,n
             write(*,'(1x,3(3i3,1x),2x,3i3)') mtmp1,itmp1
             goto 900
          end if
       end do ! m
    end do ! n

!
! --- Check the consistency of symmetry operations with atomic coordinates ---
!
    m1 = 1
    m2 = MI

    allocate( Rra(3,m1:m2) ) ; Rra=0.0d0

    c1=1.d0/dble(nnp)

    do isym=1,nsym

       do i=m1,m2

          Rra(1:3,i)=matmul( rga(1:3,1:3,isym),asi(1:3,i) )+c1*pga(1:3,isym)

          do j=1,3
             if ( abs(Rra(j,i)) >= 1.0d0 ) then
                n=int(Rra(j,i))
                Rra(j,i)=Rra(j,i)-n
             end if
          end do

          loop_j : do j=1,MI
             if( Kion(i) == Kion(j) )then
                do i3=-1,1
                do i2=-1,1
                do i1=-1,1
                   tmp4(1)=Rra(1,i)+i1
                   tmp4(2)=Rra(2,i)+i2
                   tmp4(3)=Rra(3,i)+i3
                   r0=sum( (tmp4(1:3)-asi(1:3,j))**2 )
                   if ( r0 < 1.d-10 ) then
                      Rra(1:3,i)=tmp4(1:3)
                      exit loop_j
                   end if
                end do
                end do
                end do
             end if
          end do loop_j
          if ( j > MI ) then
             write(*,*) "symmetry error (atomic coordinates)",isym,i
             goto 900
          end if

       end do ! i

    end do ! isym

    deallocate( Rra )

    call write_border( 0, " input_symdat_2(end)" )

    return

900 stop "stop@input_symdat_2"

  END SUBROUTINE input_symdat_2

  SUBROUTINE check_duplication(rga,pga)
    implicit none
    integer :: rga(:,:,:), pga(:,:)
    integer,allocatable :: itmp(:),ichk(:)
    integer :: n,j,m,i,mtmp1(3,3),itmp1(3)
    allocate( itmp(nsym) ) ; itmp=0
    allocate( ichk(nsym) ) ; ichk=0
    if ( myrank == 0 ) rewind 11
    do n=1,nsym
       ichk(n)=ichk(n)+1
       j=0
       itmp(:)=0
       do m=1,nsym
          if ( m == n ) cycle
          mtmp1(1:3,1:3)=rga(1:3,1:3,n)-rga(1:3,1:3,m)
          itmp1(1:3)=pga(1:3,n)-pga(1:3,m)
          if ( all(mtmp1==0) .and. all(itmp1==0) ) then
             j=j+1
             itmp(j)=m
             ichk(m)=ichk(m)+1
          end if
       end do ! m
       if ( j /= 0 ) then
          write(*,*) "symmetry error (duplication)",j,n
          write(*,'(1x,i5,3x,9i3,3x,3i3)') n,rga(1:3,1:3,n),pga(1:3,n)
          do i=1,j
             write(*,'(1x,i5,3x,9i3,3x,3i3)') &
                  itmp(i),rga(1:3,1:3,itmp(i)),pga(1:3,itmp(i))
          end do
       end if
       if ( ichk(n) == 1 ) then
          if ( myrank == 0 ) then
             write(11,'(1x,3(3i3,1x),2x,3i3)') &
                  ((rga(i,j,n),j=1,3),i=1,3),pga(1:3,n)
          end if
       end if
    end do ! n
    deallocate( ichk )
    deallocate( itmp )
  END SUBROUTINE check_duplication

  SUBROUTINE input_symdat_4( MI, Kion, asi )
    use aa_module, only: ax
    use bb_module, only: calc_bb
    implicit none
    integer,intent(IN) :: MI, Kion(MI)
    real(8),intent(IN) :: asi(3,MI)
    real(8) :: aa0(3,3),bb0(3,3),aa_inv(3,3),aa0_inv(3,3)
    real(8) :: pi2, tmp1(3,3),tmp2(3,3),tmp3(3,3),tmp4(3,3)
    integer :: isym, u, i, j

    pi2 = 2.0d0*acos(-1.0d0)

!    aa0(1:3,1) = (/ 0.0d0, 0.5d0, 0.5d0 /)
!    aa0(1:3,2) = (/ 0.5d0, 0.0d0, 0.5d0 /)
!    aa0(1:3,3) = (/ 0.5d0, 0.5d0, 0.0d0 /)
    aa0(1:3,1) = (/ 1.0d0, 0.0d0, 0.0d0 /)
    aa0(1:3,2) = (/ 0.0d0, 1.0d0, 0.0d0 /)
    aa0(1:3,3) = (/ 0.0d0, 0.0d0, 1.0d0 /)
    aa0=aa0*ax
    call calc_bb( aa0, bb0 )
    aa0_inv = transpose( bb0 )/pi2

    aa_inv = transpose( bb )/pi2

    u=10
    rewind u

    write(u,*) nsym,nnp

    do isym=1,nsym

       tmp1 = matmul( rga(:,:,isym), aa0_inv )
       tmp2 = matmul( aa0, tmp1 )

       tmp3 = matmul( aa_inv, tmp2 )
       tmp4 = matmul( tmp3, aa )

       write(*,*) "------------ isym=",isym
       write(*,*) tmp4,pga(:,isym)

       rga(:,:,isym) = nint( tmp4 )

       write(u,'(1x,3(3i3,1x),2x,3i3)') &
            ((rga(i,j,isym),j=1,3),i=1,3),pga(1:3,isym)

    end do ! isym
stop
  END SUBROUTINE input_symdat_4

  SUBROUTINE translation1(Rra)
    implicit none
    real(8),intent(INOUT) :: Rra(3)
    integer :: j,k
    do j=1,3
       if ( Rra(j) >= 1.0d0 ) then
          do k=1,1000
             Rra(j) = Rra(j) - k
             if ( 0.0d0 <= Rra(j) .and. Rra(j) < 1.0d0 ) exit
          end do
       else if ( Rra(j) < 0.0d0 ) then
          do k=1,1000
             Rra(j) = Rra(j) + k
             if ( 0.0d0 <= Rra(j) .and. Rra(j) < 1.0d0 ) exit
          end do
       end if
    end do
  END SUBROUTINE translation1

  SUBROUTINE translation2(Rra)
    implicit none
    real(8),intent(INOUT) :: Rra(3)
    integer :: j
    do j=1,3
       if ( Rra(j) >= 0.5d0 ) then
          Rra(j)=Rra(j)-1.0d0
       else if ( Rra(j) < -0.5d0 ) then
          Rra(j)=Rra(j)+1.0d0
       end if
    end do
  END SUBROUTINE translation2

  SUBROUTINE input_symdat_3( MI, Kion, asi )
    implicit none
    integer,intent(IN) :: MI, Kion(MI)
    real(8),intent(IN) :: asi(3,MI)
    real(8) :: tmp1(3), shift(3)
    real(8) :: fac,c1,r0
    real(8),allocatable :: Rra(:,:)
    integer :: i,isym,j,k,n,u,ierr,ierr2
    integer,allocatable :: ichk(:)

!
! --- make symmetry operation matrix ---
!

    allocate( ichk(MI) ) ; ichk=0

    fac = 4.0d0

    u = 20
    rewind u

    if ( myrank == 0 ) then
       write(*,*) "make symmetry opration matrix"
       rewind 10
       write(10,*) nsym,nint(fac)
    end if

    c1=1.d0/nnp
    allocate( Rra(3,MI) ) ; Rra=0.d0

    do isym=1,nsym

       ichk(:)=0
       ierr   =0
       ierr2  =0
       do i=1,MI

          Rra(1:3,i)=matmul( rga(1:3,1:3,isym), asi(1:3,i) )
          if ( isymmetry > 0 ) then
             call translation1( Rra(1,i) )
          else if ( isymmetry < 0 ) then
             call translation2( Rra(1,i) )
          end if

          if ( myrank == 0 ) then
             write(*,'(1x,2i6,2(3f10.5,2x))') isym,i, Rra(1:3,i), asi(1:3,i)
          end if

          do j=1,MI
             if ( Kion(i) /= Kion(j) ) cycle
             r0=sum( (Rra(:,i)-asi(:,j))**2 )
             if ( r0 < 1.d-10 ) then
                ichk(j)=ichk(j)+1
                exit
             end if
          end do ! j
          if ( j > MI ) then
             ierr=ierr+1
          end if
          if ( any( ichk >= 2 ) ) then
             ierr2=ierr2+1
          end if

       end do ! i

       if ( ierr2 /= 0 ) then
          do i=1,MI
             write(*,*) i,ichk(i)
          end do
          stop "stop@input_symda_3(1)"
       end if

       if ( ierr == 0 ) then
          if ( myrank == 0 ) then
             write(u,*) i,"okok"
             write(*,*) i,"okok"
             write(10,'(1x,3(3i3,1x),2x,3i3)') &
                  ((rga(i,j,isym),j=1,3),i=1,3),nint(fac*tmp1(1:3))
             write(*,*) "isym/nsym=",isym,nsym
          end if
       else
          write(*,*) "symmetry error (atomic coordinates)"
          write(*,*) "isym=",isym
          write(*,*) "ierr=",ierr
          stop "input_symdat_3"
          tmp1(1:3)=-0.25d0
          do i=1,MI
             Rra(1:3,i)=Rra(1:3,i)+tmp1(1:3)
             do j=1,3
                do k=1,1000
                   if ( 0.0d0 <= Rra(j,i) .and. Rra(j,i) < 1.0d0 ) exit
                   Rra(j,i) = Rra(j,i) - k
                end do
             end do
             write(*,'(1x,i6,2(3f10.5,2x))') i, Rra(1:3,i), asi(1:3,i)
             ierr=0
             do j=1,MI
                if( Kion(i) == Kion(j) )then
                   r0=sum( (Rra(:,i)-asi(:,j))**2 )
                   if( r0 < 1.d-10 )exit
                end if
             end do
             if ( j > MI ) then
                ierr=ierr+1
             end if
          end do ! i
          if ( ierr == 0 ) then
             write(u,*) i,"ok"
             write(*,*) i,"ok"
             if ( myrank == 0 ) then
                write(10,'(1x,3(3i3,1x),2x,3i3)') &
                     ((rga(i,j,isym),j=1,3),i=1,3),nint(fac*tmp1(1:3))
                write(*,*) "isym/nsym=",isym,nsym
             end if
          else
             write(u,*) i,"xx"
             write(*,*) i,"xx"
          end if
       end if

       pga(1:3,isym) = nint( fac*tmp1(1:3) )

    end do ! isym

    deallocate( Rra  )
    deallocate( ichk )

    call check_duplication(rga,pga)

    if ( myrank == 0 ) write(*,*) "see fort.10"
    stop "stop@input_symdat_3"

  END SUBROUTINE input_symdat_3

!--------1---------2---------3---------4---------5---------6---------7--

  SUBROUTINE prep_symmetry( Igrid )

    implicit none

    integer,intent(IN) :: Igrid(2,0:3)
    integer :: i,j,n,m,i1,i2,i3,isym,ierr,n1,n2,ML0
    integer :: j1,j2,j3,ii1,ii2,ii3,jj1,jj2,jj3,k1,k2,k3
    integer :: m1,m2,m3,ktmp(3,2),k
    integer,allocatable :: itmp(:),la(:,:),Rla(:,:),mtmp(:,:,:)
    integer,allocatable :: LLL2(:,:,:)
    real(8),allocatable :: Ra(:,:),Rra(:,:)
    real(8) :: r0,r1,r2,r3,tmp(3,3),c1

    if ( isymmetry == 0 ) return

    call write_border( 0, " prep_symmetry(start)" )

    n1  = idisp(myrank)+1
    n2  = idisp(myrank)+ircnt(myrank)
    ML0 = n2-n1+1
    c1  = 1.d0/dble(nnp)

    a1b = Igrid(1,1)
    b1b = Igrid(2,1)
    a2b = Igrid(1,2)
    b2b = Igrid(2,2)
    a3b = Igrid(1,3)
    b3b = Igrid(2,3)

    allocate( itmp(nsym) ) ; itmp=0
    allocate( LLL2(0:ML1-1,0:ML2-1,0:ML3-1) ) ; LLL2=0

    allocate( la(3,n1:n2)  ) ; la=0
    allocate( Rla(3,n1:n2) ) ; Rla=0

    n=-1
    i=0
    do j3=0,np_2d(3)-1
    do j2=0,np_2d(2)-1
    do j1=0,np_2d(1)-1
       n=n+1
       ii1=pinfo_grid(1,n) ; jj1=ii1+pinfo_grid(2,n)-1
       ii2=pinfo_grid(3,n) ; jj2=ii2+pinfo_grid(4,n)-1
       ii3=pinfo_grid(5,n) ; jj3=ii3+pinfo_grid(6,n)-1
       do i3=ii3,jj3
       do i2=ii2,jj2
       do i1=ii1,jj1
          i=i+1
          LLL2(i1,i2,i3)=i
       end do
       end do
       end do
    end do
    end do
    end do

    i=n1-1
    do i3=a3b,b3b
    do i2=a2b,b2b
    do i1=a1b,b1b
       i=i+1
       la(1,i)=i1
       la(2,i)=i2
       la(3,i)=i3
    end do
    end do
    end do

    if ( mod(ML1,nnp)/=0 .or. mod(ML2,nnp)/=0 .or. mod(ML2,nnp)/=0 ) then
       write(*,*) "ML1,ML2,ML3,nnp=",ML1,ML2,ML3,nnp
       write(*,*) "The grid is inconsistent with the symmetry"
       goto 900
    end if

    m1 = ML1/nnp
    m2 = ML2/nnp
    m3 = ML3/nnp

    do i=n1,n2

       itmp(1:nsym) = 0

       do isym=1,nsym

          Rla(1,i)=sum( rga(1,1:3,isym)*la(1:3,i) ) + m1*pga(1,isym)
          Rla(2,i)=sum( rga(2,1:3,isym)*la(1:3,i) ) + m2*pga(2,isym)
          Rla(3,i)=sum( rga(3,1:3,isym)*la(1:3,i) ) + m3*pga(3,isym)

          k1=Rla(1,i)/ML1 ; if ( Rla(1,i)<0 ) k1=(Rla(1,i)+1)/ML1-1
          k2=Rla(2,i)/ML2 ; if ( Rla(2,i)<0 ) k2=(Rla(2,i)+1)/ML2-1
          k3=Rla(3,i)/ML3 ; if ( Rla(3,i)<0 ) k3=(Rla(3,i)+1)/ML3-1
          Rla(1,i)=Rla(1,i)-k1*ML1
          Rla(2,i)=Rla(2,i)-k2*ML2
          Rla(3,i)=Rla(3,i)-k3*ML3
          i1=Rla(1,i)
          i2=Rla(2,i)
          i3=Rla(3,i)
          if ( i1<0 .or. i2<0 .or. i3<0 .or. &
               i1>=ML1 .or. i2>=ML2 .or. i3>=ML3 ) then
             write(*,*) i1,i2,i3
             goto 900
          end if

          itmp(isym)=LLL2(i1,i2,i3)

       end do ! isym

    end do ! i

    deallocate( Rla )
    deallocate( la )
    deallocate( LLL2 )
    deallocate( itmp )

    call write_border( 0, " prep_symmetry(end)" )

    return

900 stop "stop@prep_symmetry"

  END SUBROUTINE prep_symmetry

!--------1---------2---------3---------4---------5---------6---------7--

  SUBROUTINE sym_rho( n1,n2,n3,MSP_0,MSP_1,rho )
    implicit none
    integer,intent(IN) :: n1,n2,n3,MSP_0,MSP_1
    real(8),intent(INOUT) :: rho(n1:n2,n3)
    integer :: i,j,m,n,ierr,ispin,m1,m2,m3
    integer :: k1,k2,k3,i1,i2,i3,j1,j2,j3,l1,l2,l3,isym
    real(8) :: rho_i,c,c0,fac
    real(8),allocatable :: rho_tmp(:),rho3(:,:,:),work3(:,:,:)

    if ( isymmetry == 0 ) return

    call write_border( 1, " sym_rho(start)" )

    allocate( rho3(0:ML1-1,0:ML2-1,0:ML3-1)  ) ; rho3=0.0d0
    allocate( work3(0:ML1-1,0:ML2-1,0:ML3-1) ) ; work3=0.0d0

    m1  = ML1/nnp
    m2  = ML2/nnp
    m3  = ML3/nnp
    fac = 1.0d0/nsym

    do ispin=MSP_0,MSP_1

       work3(:,:,:)=0.0d0

       i=n1-1
       do i3=a3b,b3b
       do i2=a2b,b2b
       do i1=a1b,b1b
          i=i+1
          work3(i1,i2,i3)=rho(i,ispin)
       end do
       end do
       end do

       call mpi_allreduce(work3,rho3,ML,mpi_real8,mpi_sum,comm_grid,ierr)

       rho(n1:n2,ispin)=0.0d0

       i=n1-1
       do i3=a3b,b3b
       do i2=a2b,b2b
       do i1=a1b,b1b
          i=i+1
          do isym=1,nsym
             j1=i1-m1*pga(1,isym)
             j2=i2-m2*pga(2,isym)
             j3=i3-m3*pga(3,isym)
             l1=nint(rgb(1,1,isym)*j1+rgb(2,1,isym)*j2+rgb(3,1,isym)*j3)
             l2=nint(rgb(1,2,isym)*j1+rgb(2,2,isym)*j2+rgb(3,2,isym)*j3)
             l3=nint(rgb(1,3,isym)*j1+rgb(2,3,isym)*j2+rgb(3,3,isym)*j3)
             k1=l1/ML1 ; if ( l1<0 ) k1=(l1+1)/ML1-1
             k2=l2/ML2 ; if ( l2<0 ) k2=(l2+1)/ML2-1
             k3=l3/ML3 ; if ( l3<0 ) k3=(l3+1)/ML3-1
             l1=l1-k1*ML1
             l2=l2-k2*ML2
             l3=l3-k3*ML3
             rho(i,ispin) = rho(i,ispin) + fac*rho3(l1,l2,l3)
          end do
       end do
       end do
       end do

       c0=sum(rho(n1:n2,ispin))*dV
       call mpi_allreduce(c0,c,1,mpi_real8,mpi_sum,comm_grid,ierr)
!       if ( disp_switch_parallel ) write(*,*) "check sum(rho)=",c

    end do ! ispin

    deallocate( work3 )
    deallocate( rho3  )

    call write_border( 1, " sym_rho(end)" )

    return

  END SUBROUTINE sym_rho

!--------1---------2---------3---------4---------5---------6---------7--

  SUBROUTINE sym_force( MI, Kion, asi, force )
    implicit none
    integer,intent(IN) :: MI, Kion(MI)
    real(8),intent(IN) :: asi(3,MI)
    real(8),intent(INOUT) :: force(3,MI)
    integer :: a,b,isym,j,m,i1,i2,i3,n
    real(8) :: atmp(3),c1,c2,c3,r,r0,faa(3),tmp1(3)
    real(8),allocatable :: ftmp(:,:)

    if ( isymmetry /= 1 ) return

    call write_border( 1, " sym_force(start)" )

    c1 = 1.0d0/dble(nnp)
    c2 = 1.0d0/dble(nsym)
    c3 = 1.0d0/( 2.0d0*acos(-1.0d0) )

    allocate( ftmp(3,MI) )

    ftmp(:,:) = force(:,:)

    do a=1,MI
       force(1:3,a)=matmul( ftmp(1:3,a),bb(1:3,1:3) )*c3
    end do

    ftmp(:,:) = 0.0d0

    do a=1,MI

       do isym=1,nsym

          atmp(:) = matmul( rga(:,:,isym),asi(:,a) )+c1*pga(:,isym)

          do j=1,3
             if ( abs(atmp(j)) >= 1.0d0 ) then
                n=int(atmp(j))
                atmp(j)=atmp(j)-n
             end if
          end do

          loop_b : do b=1,MI

             if ( Kion(a) == Kion(b) ) then

                do i3=-1,1
                do i2=-1,1
                do i1=-1,1
                   tmp1(1)=atmp(1)+i1
                   tmp1(2)=atmp(2)+i2
                   tmp1(3)=atmp(3)+i3
                   r0=sum( (tmp1(1:3)-asi(1:3,b))**2 )
                   if ( r0 < 1.d-10 ) then
                      ftmp(1:3,b) = ftmp(1:3,b) &
                           + matmul( rga(:,:,isym),force(1:3,a) )
                      exit loop_b
                   end if
                end do
                end do
                end do
             end if

          end do loop_b

          if ( b > MI ) then
             write(*,*) "symmetry error (sym_force)",isym,a
             stop "stop@sym_force"
          end if

       end do ! isym

    end do ! a

    do a=1,MI
       force(1:3,a) = matmul( aa(1:3,1:3), ftmp(1:3,a) )*c2
    end do

    deallocate( ftmp )

    call write_border( 1, " sym_force(end)" )

    return
  END SUBROUTINE sym_force

!--------1---------2---------3---------4---------5---------6---------7--

  SUBROUTINE chk_sym_mat(ng,np,mat,vec)

    implicit none

    integer,intent(IN) :: ng,np
    integer,intent(IN) :: mat(3,3,ng),vec(3,ng)
    integer :: n,j,m,loop,i,m1,m2
    integer :: mtmp1(3,3),itmp1(3),mtmp2(3,3),itmp2(3)
    integer,allocatable :: itmp(:)

! check duplication
!
    allocate( itmp(ng) )

    do n=1,ng
       j=0
       itmp=0
       do m=1,ng
          mtmp1(1:3,1:3)=mat(1:3,1:3,n)-mat(1:3,1:3,m)
          itmp1(1:3)=vec(1:3,n)-vec(1:3,m)
          if ( all(mtmp1==0) .and. all(itmp1==0) ) then
             j=j+1
             itmp(j)=m
          end if
       end do
       if ( j/=1 ) then
          write(*,*) "symmetry error (duplication)",j,n
          write(*,'(1x,i5,3x,9i3,3x,3i3)') n,mat(1:3,1:3,n),vec(1:3,n)
          do i=1,j
             write(*,'(1x,i5,3x,9i3,3x,3i3)') &
                  itmp(i),mat(1:3,1:3,itmp(i)),vec(1:3,itmp(i))
          end do
          goto 900
       end if
    end do
    deallocate( itmp )

! existence of the unit operator
!
    m1=0
    do n=1,ng
       itmp1(1)=mat(1,1,n)
       itmp1(2)=mat(2,2,n)
       itmp1(3)=mat(3,3,n)
       i=sum(abs(mat(1:3,1:3,n)))
       if ( i==3 .and. all(itmp1==1) .and. all(vec(1:3,n)==0) ) then
          m1=m1+1
       end if
    end do
    if ( m1/=1 ) then
       write(*,*) "symmetry error (no unit operator)",m1
       goto 900
    end if

! (mat_1)(mat_2)=(mat)
!
    do n=1,ng
       do m=1,ng
          mtmp1(1:3,1:3)=matmul( mat(1:3,1:3,m),mat(1:3,1:3,n) )
          itmp1(1:3)=matmul( mat(1:3,1:3,m),vec(1:3,n) )+vec(1:3,m)
          do j=1,3
             do loop=1,10000
                if ( itmp1(j)<0 ) then
                   itmp1(j)=itmp1(j)+np
                else if ( itmp1(j)>=np ) then
                   itmp1(j)=itmp1(j)-np
                else
                   exit
                end if
             end do
          end do
          m2=0
          do i=1,ng
             mtmp2(1:3,1:3)=mtmp1(1:3,1:3)-mat(1:3,1:3,i)
             itmp2(1:3)=itmp1(1:3)-vec(1:3,i)
             if ( all(mtmp2==0) .and. all(itmp2==0) ) then
                m2=m2+1
             end if
          end do
          if ( m2/=1 ) then
             write(*,*) "symmetry error (this is not group)",m2,m,n
             goto 900
          end if
       end do ! m
    end do ! n

    return

900 stop "stop@chk_sym_mat"

  END SUBROUTINE chk_sym_mat

  SUBROUTINE construct_matrix_symmetry( isym, ML, SymMat )
    implicit none
    integer,intent(IN)  :: isym, ML
    real(8),intent(OUT) :: SymMat(ML,ML)
    integer :: i,j,n,m,i1,i2,i3,ierr,n1,n2,ML0
    integer :: j1,j2,j3,ii1,ii2,ii3,jj1,jj2,jj3,k1,k2,k3
    integer :: m1,m2,m3,ktmp(3,2),k
    integer,allocatable :: LLL2(:,:,:)
    real(8) :: la(3),Rla(3)
    real(8) :: r0,r1,r2,r3,tmp(3,3),c1

    if ( isymmetry == 0 ) return

    SymMat(:,:) = 0.0d0

    n1  = idisp(myrank)+1
    n2  = idisp(myrank)+ircnt(myrank)
    ML0 = n2-n1+1
    c1  = 1.d0/dble(nnp)

    allocate( LLL2(0:ML1-1,0:ML2-1,0:ML3-1) ) ; LLL2=0

    n=-1
    i=0
    do j3=0,np_2d(3)-1
    do j2=0,np_2d(2)-1
    do j1=0,np_2d(1)-1
       n=n+1
       ii1=pinfo_grid(1,n) ; jj1=ii1+pinfo_grid(2,n)-1
       ii2=pinfo_grid(3,n) ; jj2=ii2+pinfo_grid(4,n)-1
       ii3=pinfo_grid(5,n) ; jj3=ii3+pinfo_grid(6,n)-1
       do i3=ii3,jj3
       do i2=ii2,jj2
       do i1=ii1,jj1
          i=i+1
          LLL2(i1,i2,i3)=i
       end do
       end do
       end do
    end do
    end do
    end do

    m1 = ML1/nnp
    m2 = ML2/nnp
    m3 = ML3/nnp

    do i3=0,ML3-1
    do i2=0,ML2-1
    do i1=0,ML1-1

       i = 1 + i1 + i2*ML1 + i3*ML1*ML2

       la(1) = i1
       la(2) = i2
       la(3) = i3

       Rla(1)=sum( rga(1,1:3,isym)*la(1:3) ) + m1*pga(1,isym)
       Rla(2)=sum( rga(2,1:3,isym)*la(1:3) ) + m2*pga(2,isym)
       Rla(3)=sum( rga(3,1:3,isym)*la(1:3) ) + m3*pga(3,isym)

       k1=Rla(1)/ML1 ; if ( Rla(1)<0 ) k1=(Rla(1)+1)/ML1-1
       k2=Rla(2)/ML2 ; if ( Rla(2)<0 ) k2=(Rla(2)+1)/ML2-1
       k3=Rla(3)/ML3 ; if ( Rla(3)<0 ) k3=(Rla(3)+1)/ML3-1
       Rla(1)=Rla(1)-k1*ML1
       Rla(2)=Rla(2)-k2*ML2
       Rla(3)=Rla(3)-k3*ML3
       j1=Rla(1)
       j2=Rla(2)
       j3=Rla(3)

       j = LLL2(j1,j2,j3)

       SymMat(j,i) = 1.0d0

    end do ! i1
    end do ! i2
    end do ! i3

    deallocate( LLL2 )

    return

  END SUBROUTINE construct_matrix_symmetry

END MODULE symmetry_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
      subroutine rpt2(ixy,maxm,meqn,mwaves,mbc,mx,
     &      ql,qr,aux1,aux2,aux3,
     &      imp,asdq,bmasdq,bpasdq)
c     =====================================================
      implicit double precision (a-h,o-z)
c
c     # Riemann solver in the transverse direction for the shallow water
c     # equations  on a quadrilateral grid.
c
c     # Split asdq (= A^* \Delta q, where * = + or -)
c     # into down-going flux difference bmasdq (= B^- A^* \Delta q)
c     #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)
c
c     # Uses Roe averages and other quantities which were
c     # computed in rpn2sh and stored in the common block comroe.
c
      dimension     ql(1-mbc:maxm+mbc, meqn)
      dimension     qr(1-mbc:maxm+mbc, meqn)
      dimension   asdq(1-mbc:maxm+mbc, meqn)
      dimension bmasdq(1-mbc:maxm+mbc, meqn)
      dimension bpasdq(1-mbc:maxm+mbc, meqn)
      dimension   aux1(1-mbc:maxm+mbc, *)
      dimension   aux2(1-mbc:maxm+mbc, *)
      dimension   aux3(1-mbc:maxm+mbc, *)
c
      parameter (maxm2 = 1002)          !# assumes at most 1000x1000 grid with mbc=2
      dimension u(-1:maxm2),v(-1:maxm2),a(-1:maxm2),h(-1:maxm2)
      dimension wave(-1:maxm2, 4, 3)
      dimension    s(-1:maxm2, 3)
      dimension enx(-1:maxm2), eny(-1:maxm2), enz(-1:maxm2)
      dimension etx(-1:maxm2), ety(-1:maxm2), etz(-1:maxm2)
      dimension gamma(-1:maxm2)

      dimension delta(4)
      common /sw/  g
      common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
c

      if (-1.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then
         write(6,*) 'need to increase maxm2 in rpB'
         stop
      endif

      if(ixy.eq.1) then
         dx = dxcom
      else
         dx = dycom
      endif
c
c
      if(ixy.eq.1) then
         ioff = 7
      else
         ioff = 1
      endif

c
c     # imp is used to flag whether wave is going to left or right,
c     # since states and grid are different on each side
c
      if (imp.eq.1) then
c        # asdq = amdq, moving to left
         ix1 = 2-mbc
         ixm1 = mx+mbc
      else
c        # asdq = apdq, moving to right
         ix1 = 1-mbc
         ixm1 = mx+mbc
      endif
c
c     --------------
c     # up-going:
c     --------------
c

c     # determine rotation matrix for interface above cell, using
c     aux3
c     [ alf  beta ]
c     [-beta  alf ]
c
      do i=ix1,ixm1
c
         if (imp.eq.1) then
            i1 = i-1
         else
            i1 = i
         endif
c
         enx(i) =   aux3(i1,ioff+1)
         eny(i) =   aux3(i1,ioff+2)
         enz(i) =   aux3(i1,ioff+3)
         etx(i) =   aux3(i1,ioff+4)
         ety(i) =   aux3(i1,ioff+5)
         etz(i) =   aux3(i1,ioff+6)
         gamma(i) = dsqrt(etx(i)**2 + ety(i)**2 + etz(i)**2)
         etx(i) =   etx(i) / gamma(i)
         ety(i) =   ety(i) / gamma(i)
         etz(i) =   etz(i) / gamma(i)

         h(i) = ql(i1,1)
         u(i) = (enx(i)*ql(i1,2)+eny(i)*ql(i1,3)+enz(i)*ql(i1,4))
     &         / h(i)
         v(i) = (etx(i)*ql(i1,2)+ety(i)*ql(i1,3)+etz(i)*ql(i1,4))
     &         / h(i)
         a(i) = dsqrt(g*h(i))
      enddo
c
c
c
c     # now split asdq into waves:
c
c     # find a1 thru a3, the coefficients of the 3 eigenvectors:
      do 20 i = ix1,ixm1
         delta(1) = asdq(i,1)
         delta(2) = enx(i)*asdq(i,2)+eny(i)*asdq(i,3)+enz(i)*asdq(i,4)
         delta(3) = etx(i)*asdq(i,2)+ety(i)*asdq(i,3)+etz(i)*asdq(i,4)
         a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i))
         a2 = -v(i)*delta(1) + delta(3)
         a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i))
c
c        # Compute the waves.
c
         wave(i,1,1) = a1
         wave(i,2,1) = enx(i)*a1*(u(i)-a(i)) + etx(i)*a1*v(i)
         wave(i,3,1) = eny(i)*a1*(u(i)-a(i)) + ety(i)*a1*v(i)
         wave(i,4,1) = enz(i)*a1*(u(i)-a(i)) + etz(i)*a1*v(i)
         s(i,1) = (u(i)-a(i))*gamma(i)/dx
c
         wave(i,1,2) = 0.0d0
         wave(i,2,2) = etx(i)*a2
         wave(i,3,2) = ety(i)*a2
         wave(i,4,2) = etz(i)*a2
         s(i,2) = u(i) * gamma(i)/dx
c
         wave(i,1,3) = a3
         wave(i,2,3) = enx(i)*a3*(u(i)+a(i)) + etx(i)*a3*v(i)
         wave(i,3,3) = eny(i)*a3*(u(i)+a(i)) + ety(i)*a3*v(i)
         wave(i,4,3) = enz(i)*a3*(u(i)+a(i)) + etz(i)*a3*v(i)
         s(i,3) = (u(i)+a(i)) * gamma(i)/dx
   20 continue
c
c
c     # compute flux difference bpasdq
c     --------------------------------
c
      do 40 m=1,meqn
         do 40 i=ix1,ixm1
	    bpasdq(i,m) = 0.d0
	    do 30 mw=1,mwaves
	       bpasdq(i,m) = bpasdq(i,m)
     &               + dmax1(s(i,mw),0.d0)*wave(i,m,mw)
   30       continue
   40    continue
c
c        # project momentum component of bpasdq to tangent plane:
         do i=ix1, ixm1
            if (imp.eq.1) then
               i1 = i-1
            else
               i1 = i
            endif
            erx = aux3(i1,14)
            ery = aux3(i1,15)
            erz = aux3(i1,16)
            bn = erx*bpasdq(i,2)+ery*bpasdq(i,3)+erz*bpasdq(i,4)
            bpasdq(i,2) = bpasdq(i,2) - bn*erx
            bpasdq(i,3) = bpasdq(i,3) - bn*ery
            bpasdq(i,4) = bpasdq(i,4) - bn*erz

         enddo

c
c
c        --------------
c        # down-going:
c        --------------
c

c        # determine rotation matrix for interface below cell, using
c        aux2
c        [ alf  beta ]
c        [-beta  alf ]
c
         do i=ix1,ixm1
c
            if (imp.eq.1) then
               i1 = i-1
            else
               i1 = i
            endif
c
            enx(i) =   aux2(i1,ioff+1)
            eny(i) =   aux2(i1,ioff+2)
            enz(i) =   aux2(i1,ioff+3)
            etx(i) =   aux2(i1,ioff+4)
            ety(i) =   aux2(i1,ioff+5)
            etz(i) =   aux2(i1,ioff+6)
            gamma(i) = dsqrt(etx(i)**2 + ety(i)**2 + etz(i)**2)
            etx(i) =   etx(i) / gamma(i)
            ety(i) =   ety(i) / gamma(i)
            etz(i) =   etz(i) / gamma(i)
            u(i) = (enx(i)*ql(i1,2)+eny(i)*ql(i1,3)+enz(i)*ql(i1,4))
     &            / h(i)
            v(i) = (etx(i)*ql(i1,2)+ety(i)*ql(i1,3)+etz(i)*ql(i1,4))
     &            / h(i)
         enddo
c
c
c
c        # now split asdq into waves:
c
c        # find a1 thru a3, the coefficients of the 3 eigenvectors:
         do 80 i = ix1,ixm1
            delta(1) = asdq(i,1)
            delta(2) = enx(i)*asdq(i,2)+eny(i)*asdq(i,3)+enz(i)*asdq(i
     &            ,4)
            delta(3) = etx(i)*asdq(i,2)+ety(i)*asdq(i,3)+etz(i)*asdq(i
     &            ,4)

            a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i))
            a2 = -v(i)*delta(1) + delta(3)
            a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i))
c
c           # Compute the waves.
c
            wave(i,1,1) = a1
            wave(i,2,1) = enx(i)*a1*(u(i)-a(i)) + etx(i)*a1*v(i)
            wave(i,3,1) = eny(i)*a1*(u(i)-a(i)) + ety(i)*a1*v(i)
            wave(i,4,1) = enz(i)*a1*(u(i)-a(i)) + etz(i)*a1*v(i)
            s(i,1) = (u(i)-a(i)) * gamma(i)/dx
c
            wave(i,1,2) = 0.0d0
            wave(i,2,2) = etx(i)*a2
            wave(i,3,2) = ety(i)*a2
            wave(i,4,2) = etz(i)*a2
            s(i,2) = u(i) * gamma(i)/dx
c
            wave(i,1,3) = a3
            wave(i,2,3) = enx(i)*a3*(u(i)+a(i)) + etx(i)*a3*v(i)
            wave(i,3,3) = eny(i)*a3*(u(i)+a(i)) + ety(i)*a3*v(i)
            wave(i,4,3) = enz(i)*a3*(u(i)+a(i)) + etz(i)*a3*v(i)
            s(i,3) = (u(i)+a(i)) * gamma(i)/dx
   80    continue
c
c
c        # compute flux difference bmasdq
c        --------------------------------
c
         do 100 m=1,meqn
            do 100 i=ix1,ixm1
               bmasdq(i,m) = 0.d0
               do 90 mw=1,mwaves
                  bmasdq(i,m) = bmasdq(i,m)
     &                  + dmin1(s(i,mw), 0.d0)*wave(i,m,mw)
   90          continue
  100       continue
c           # project momentum component of bmasdq to tangent plane:
            do i=ix1, ixm1
               if (imp.eq.1) then
                  i1 = i-1
               else
                  i1 = i
               endif
               erx = aux1(i1,14)
               ery = aux1(i1,15)
               erz = aux1(i1,16)
               bn = erx*bmasdq(i,2)+ery*bmasdq(i,3)+erz*bmasdq(i,4)
               bmasdq(i,2) = bmasdq(i,2) - bn*erx
               bmasdq(i,3) = bmasdq(i,3) - bn*ery
               bmasdq(i,4) = bmasdq(i,4) - bn*erz

            enddo
c
c

            return
            end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
HANSARD REVISE * NUMERO 82
Le vendredi 27 mars 1998
LE MATCH DE HOCKEY OPPOSANT LES DEPUTES AUX PAGES
Mme Rose-Marie Ur
M. Guy St-Julien
LA COMMISSION SUR L'AVENIR DE L'AEROPORT DE MIRABEL
LE SYSTEME JUDICIAIRE DANS LA REGION DU RESTIGOUCHE
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
L'hon. Ralph E. Goodale
LA COMMISSION SUR L'AVENIR DE L'AEROPORT DE MIRABEL
L'hon. David M. Collenette
L'hon. David M. Collenette
L'hon. Ralph E. Goodale
L'hon. Arthur C. Eggleton
L'hon. Arthur C. Eggleton
L'hon. Ralph E. Goodale
L'hon. Arthur C. Eggleton
L'hon. David M. Collenette
L'hon. David M. Collenette
L'hon. Ralph E. Goodale
LA STATION DE RADIO DE LA GARDE COTIERE
ENERGIE ATOMIQUE DU CANADA LIMITEE
L'hon. Ralph E. Goodale
L'hon. Arthur C. Eggleton
L'hon. David M. Collenette
REPONSE DU GOUVERNEMENT A DES PETITIONS
Justice et droits de la personne
La taxe sur les produits et services
HANSARD REVISE * NUMERO 82
Le vendredi 27 mars 1998
La seance est ouverte a 10 heures.
M. John Harvard (secretaire parlementaire du ministre de l'Agriculture et de l'Agroalimentaire, Lib.):
Cette industrie connait d'ailleurs une croissance impressionnante.
Nos exportations continuent de s'accroitre, tout comme notre marche interieur.
Nous ferons tout pour repondre a l'evolution des besoins de l'industrie.
Ensemble, nous travaillons a l'elaboration de strategies en vue de parvenir a cet objectif.
Je suis donc tres heureux de discuter de ce projet de loi aujourd'hui.
L'industrie du grain divise les cultures en cultures ordinaires ou speciales.
La premiere serie de consultations fut entreprise par un groupe de producteurs.
Un consensus en la matiere s'est degage de ces discussions.
En mars 1995, la commission a publie un rapport refletant ce consensus.
Le gros de la production de cultures speciales vient des trois provinces des Prairies.
En 1984-1985, la production etait d'environ un million de tonnes.
En 1997-98, elle a grimpe a environ trois millions de tonnes.
Les marche des cultures speciales continue de croitre.
Les grands et les petits producteurs contribuent a cette remarquable croissance.
Le secteur des cultures speciales exporte beaucoup.
L'entree de nouveaux intervenants stimulera l'innovation et le developpement de nouveaux marches.
Je souligne que ce regime est facultatif.
Les producteurs pourront choisir d'y adherer ou non.
La contribution initiale equivaudra a un faible pourcentage de la valeur des cereales livrees.
Pour l'industrie, cette mesure permet de reduire la paperasse au strict minimum.
Apres deduction des droits administratifs, la commission enverra l'argent a un assureur.
Je vais revenir un peu plus longuement sur ce comite dans quelques instants.
On s'attend a ce que la plupart des producteurs y adherent.
Une fois la periode de demarrage terminee, le regime d'assurance s'autofinancera.
Il y a environ 800 usines de nettoyage de graines dans les Prairies.
Seuls quelques exploitants d'usines de nettoyage de graines ont une licence.
Comme ils en sont incapables, leur expansion est limitee.
Les regles seront les memes pour les vendeurs grands et petits.
Chacune des provinces des Prairies et chaque secteur de l'industrie seront representes.
La majorite des membres seront des producteurs.
La LSAP a ete adoptee en 1995.
Ces sanctions ont toutefois pour effet de suspendre ou d'arreter les activites commerciales.
La loi met ainsi l'accent sur l'observation, et non sur le chatiment.
Cette mesure legislative devrait etre proclamee sous peu.
On a recu plusieurs messages categoriques par suite du processus de consultation.
Elle prevoit uniquement la reglementation des contrats a terme de grain.
Des discussions preliminaires ont deja eu lieu.
Je recommande aux deputes d'appuyer le projet de loi C-26.
M. Jay Hill (Prince George-Peace River, Ref.):
Au nom des agriculteurs, j'encourage le gouvernement a saisir cette occasion.
Apres tout, n'est-ce pas la le role des deputes?
Les reformistes consultent, ecoutent et consultent de nouveau.
Le temps est venu pour les liberaux de consulter la population.
Je dis pour l'instant, car il est encore tres tot.
Le debat et l'analyse de cette mesure legislative ne font que commencer.
Comme ce sont les membres de ces groupes qui sont les specialistes dans ce domaine.
Cette demarche de consultation doit avoir la chance de produire des resultats.
Ce n'est malheureusement pas quelque chose auquel on peut s'attendre tout simplement.
Ils ont juge qu'ils constituaient une perte de temps.
Le gouvernement n'a pas voulu s'embarrasser de consultations dans son echeancier serre.
Je me preoccupe des agriculteurs.
Je suis sceptique, et c'est comprehensible.
Si la majorite s'y oppose, nous le saurons aussi.
Les interets des agriculteurs doivent etre proteges.
Par exemple, les primes d'assurance doivent etre limitees.
La bureaucratie et l'agriculture ne font pas bon menage.
Tout d'abord, la Commission canadienne des grains administrera le regime d'assurance.
Le regime d'assurance risque de s'enliser dans la bureaucratie gouvernementale.
Pourquoi ne pas confier la fonction d'assurance a une compagnie d'assurance privee?
Il faut etudier cette possibilite.
Il faut reduire l'ingerence gouvernementale dans l'agriculture.
Nous devons soutenir l'accroissement du role du secteur prive.
La question de la contribution suscitera certainement quelques preoccupations.
Le projet de loi a suscite des reactions.
Les deputes sont peu familiers avec le projet de loi.
Nous avons aussi recu des commentaires d'autres producteurs.
La lettre dit ensuite ceci:
L'auteur de la lettre fait reference a la concurrence.
Il devrait y avoir libre concurrence.
Elle concerne egalement le projet de loi.
Ce qui nous inquiete surtout, c'est le regime d'assurance.
Cela se poursuit ainsi.
J'aurais voulu faire part a la Chambre de quelques autres anecdotes pertinentes.
Quand ca arrive, l'agriculteur se retrouve en mauvaise posture.
Malheureusement, c'est arrive a quelques reprises.
Leur souci premier etait le grain entrepose dans les locaux de la societe.
Le president suppleant (M. McClelland):
Je suis desole, mais je dois interrompre le depute de Prince George-Peace River.
Toutefois, il disposera de 16 minutes apres la periode des questions.
LE MATCH DE HOCKEY OPPOSANT LES DEPUTES AUX PAGES
M. David Pratt (Nepean-Carleton, Lib.):
Il s'agissait du match annuel de hockey entre les deputes et les pages.
M. Paul Steckle (Huron-Bruce, Lib.):
A ce jour, les investissements depassent 57 millions de dollars.
Ensemble, nous travaillons a appuyer l'innovation et a promouvoir la croissance economique.
M. Derrek Konrad (Prince Albert, Ref.):
On ne reclame pas de subventions, on reclame seulement la liberte voulue pour etre concurrentiel.
M. John McKay (Scarborough-Est, Lib.):
Nos plus grandes ressources sont les ressources humaines.
La deputee de Lambton-Kent-Middlesex a la parole.
Mme Rose-Marie Ur (Lambton-Kent-Middlesex, Lib.):
M. John Duncan (Ile de Vancouver-Nord, Ref.):
Pourtant, les elements de preuve abondent.
Les liberaux ont fait leur temps.
Mais ils esperent changer leur destin.
En faisant campagne aux partielles,
ils esperent se remettre en selle.
Mais apres lundi, ils auront compris
que leur temps ici est bel et bien fini.
M. Guy St-Julien (Abitibi, Lib.):
M. Guy St-Julien:
L'homme de Sherbrooke a la competence.
Il a des qualites de rassembleur.
Il est courageux, franc et represente un exemple de diplomatie.
Nous sommes fiers de sa decision.
M. Hec Clouthier (Renfrew-Nippissing-Pembroke, Lib.):
Mme Val Meredith (South Surrey-White Rock-Langley, Ref.):
Cela demontre a quel point les liberaux ne connaissent pas la Colombie-Britannique.
Mme Monique Guay (Laurentides, BQ):
Monsieur le President, l'univers du theatre est en fete ces jours-ci.
Par ailleurs, nous celebrons aujourd'hui la Journee mondiale du theatre.
M. John Harvard (Charleswood-Assiniboine, Lib.):
Mme Bev Desjarlais (Churchill, NPD):
M. Mac Harb (Ottawa-Centre, Lib.):
On dirait un extrait du programme du Parti reformiste.
Quand la mouche cessera-t-elle de faire avancer le coche?
Est-ce qu'on ne peut pas faire mieux?
LA COMMISSION SUR L'AVENIR DE L'AEROPORT DE MIRABEL
M. Maurice Dumas (Argenteuil-Papineau, BQ):
LE SYSTEME JUDICIAIRE DANS LA REGION DU RESTIGOUCHE
M. Jean Dube (Madawaska-Restigouche, PC):
Que font les politiciens liberaux?
Ils se lancent la balle d'un niveau de gouvernement a l'autre.
J'implore la ministre de la Justice a demontrer du leadership dans ce dossier.
M. Joseph Volpe (Eglinton-Lawrence, Lib.):
On evalue a 20 000 le nombre des Canadiens qui souffrent d'insuffisance renale.
Des millions d'autres souffrent de maladies du rein.
Pour de nombreux Canadiens, la greffe de rein est la seule solution.
M. Grant Hill (Macleod, Ref.):
M. Joseph Volpe (secretaire parlementaire du ministre de la Sante, Lib.):
Les mesures qui sont actuellement proposees aux Canadiens temoignent justement de cette approche.
M. Grant Hill (Macleod, Ref.):
M. Joseph Volpe (secretaire parlementaire du ministre de la Sante, Lib.):
Monsieur le President, ce qui est propose constitue une offre equitable et empreinte de compassion.
La parole est au depute de Macleod.
M. Grant Hill (Macleod, Ref.):
M. Joseph Volpe (secretaire parlementaire du ministre de la Sante, Lib.):
C'est effectivement le resultat d'une approche fondee sur des negociations.
Nous nous entendions tous pour dire, y compris les groupes concernes, que...
La parole est au depute d'Elk Island.
M. Ken Epp (Elk Island, Ref.):
Mon oncle est decede il y a deux mois.
M. Joseph Volpe (secretaire parlementaire du ministre de la Sante, Lib.):
Monsieur le President, je plains toutes les victimes de cette terrible tragedie.
Si elles ne souhaitent pas...
Le depute de Elk island a la parole.
M. Ken Epp (Elk Island, Ref.):
M. Joseph Volpe (secretaire parlementaire du ministre de la Sante, Lib.):
Les personnes qui preferent proceder de facon differente ont toujours le loisir de le faire.
Mme Suzanne Tremblay (Rimouski-Mitis, BQ):
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
Mme Suzanne Tremblay (Rimouski-Mitis, BQ):
Monsieur le President, c'est facile de se payer la tete des contribuables.
C'est la seule difference entre les deux rapports.
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
Mme Helene Alarie (Louis-Hebert, BQ):
Monsieur le President, ma question s'adresse a la ministre du Patrimoine canadien.
Nous sommes en presence d'un mystere.
Cette lettre est non datee et a moins de trois lignes.
Que se cache-t-il derriere ce nouveau cafouillage?
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
Monsieur le President, j'ai ici une lettre.
Malheureusement, c'est une lettre qui est datee de Quebec, le 28 novembre 1996.
Actuellement, depuis deux ans...
La chicane, c'est de leur cote.
Mme Helene Alarie (Louis-Hebert, BQ):
Monsieur le President, plus ca va, moins c'est clair.
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
Monsieur le President, encore un parti qui ne cherche qu'a faire de la chicane.
Donnez-nous donc la reponse.
Mme Alexa McDonough (Halifax, NPD):
Monsieur le President, ma question s'adresse au premier ministre.
Le tres hon. Jean Chretien (premier ministre, Lib.):
C'etait tres clair.
Mme Alexa McDonough (Halifax, NPD):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Le comite etudie toute cette question.
M. Mark Muise (West Nova, PC):
J'ai interroge le ministre a ce sujet a la periode des questions.
M. Mark Muise (West Nova, PC):
Est-ce que c'est une possibilite que le ministre des Peches envisage serieusement?
Dans l'affirmative, quand peut-on esperer avoir une reponse?
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
M. Mike Scott (Skeena, Ref.):
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
M. Mike Scott (Skeena, Ref.):
Elle parle du suivi au rapport de la commission royale d'enquete.
Rien ne prouve que les conditions dans les reserves canadiennes s'ameliorent.
Tout porte a croire plutot que les conditions se deteriorent.
Qu'est-ce que le gouvernement tente de cacher?
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
LA COMMISSION SUR L'AVENIR DE L'AEROPORT DE MIRABEL
M. Michel Guimond (Beauport-Montmorency-Orleans, BQ):
Ma question s'adresse au ministre des Transports.
L'hon. David M. Collenette (ministre des Transports, Lib.):
M. Michel Guimond (Beauport-Montmorency-Orleans, BQ):
Monsieur le President, on se passe la balle.
Je rappelle au ministre que la responsabilite du gouvernement federal dans ce gachis est evidente.
Ce qu'on demande, ce n'est pas la charite, c'est l'equite.
L'hon. David M. Collenette (ministre des Transports, Lib.):
M. Jay Hill (Prince George-Peace River, Ref.):
Il a dit que c'etait une police d'assurance contre les deversements petroliers.
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
Il devrait le reconnaitre et s'excuser de...
Le depute de Prince-George-Peace River a la parole.
M. Jay Hill (Prince George-Peace River, Ref.):
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
En Ontario, on parle exclusivement de forces d'intervention du secteur prive.
Il n'est aucunement question d'impot.
M. Stephane Bergeron (Vercheres, BQ):
Monsieur le President, ma question s'adresse a la ministre de l'Environnement.
L'hon. Christine Stewart (ministre de l'Environnement, Lib.):
Nous allons nous pencher sur tous les mecanismes a la disposition du gouvernement federal.
M. Stephane Bergeron (Vercheres, BQ):
Le probleme en Ontario est lie a la gestion d'Hydro Ontario.
Il n'a absolument rien a voir avec la technologie du reacteur Candu.
La parole est a la deputee de South Surrey-White Rock-Langley.
Mme Val Meredith (South Surrey-White Rock-Langley, Ref.):
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
Je suis egalement heureux de dire que personne n'a ete blesse.
Une enquete approfondie est en cours.
Mme Val Meredith (South Surrey-White Rock-Langley, Ref.):
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
Monsieur le President, ce n'est simplement pas vrai.
M. Rene Laurin (Joliette, BQ):
L'hon. Lyle Vanclief (ministre de l'Agriculture et de l'Agroalimentaire, Lib.):
Mme Susan Whelan (Essex, Lib.):
La Colombie-Britannique profite tout particulierement de la venue des touristes asiatiques.
M. Walt Lastewka (secretaire parlementaire du ministre de l'Industrie, Lib.):
M. Darrel Stinson (Okanagan-Shuswap, Ref.):
Cette commission a formule 400 recommandations qui ramassent maintenant la poussiere sur les tablettes.
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
Plus d'argent pour le chef.
Le depute crie: Plus d'argent pour le chef.
Le Parti reformiste n'a pas appuye les recommandations.
M. Derrek Konrad (Prince Albert, Ref.):
Ma question s'adresse a la ministre des Affaires indiennes.
Va-t-elle prendre des mesures avant qu'une autre tragedie ne se produise?
L'hon. Sheila Copps (ministre du Patrimoine canadien, Lib.):
M. Yvon Godin (Acadie-Bathurst, NPD):
Ma question s'adresse au premier ministre.
On peut nettement voir les effets de ce rapport dans le budget de 1997.
M. Gordon Earle (Halifax-Ouest, NPD):
Monsieur le President, ma question s'adresse au ministre de la Defense nationale.
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
Monsieur le President, je ne suis pas d'accord avec le preambule de la question.
C'est de cela qu'il s'agit.
La defense nationale a moins de ressources qu'avant.
M. Andre Harvey (Chicoutimi, PC):
Monsieur le President, ma question s'adresse au ministre des Transports.
Je sais que le ministre travaille a l'elaboration d'une nouvelle politique routiere nationale.
L'hon. David M. Collenette (ministre des Transports, Lib.):
Comme l'honorable depute le sait, les autoroutes canadiennes relevent des provinces.
M. Andre Harvey (Chicoutimi, PC):
L'hon. David M. Collenette (ministre des Transports, Lib.):
M. Benoit Serre (Timiskaming-Cochrane, Lib.):
Monsieur le President, ma question s'adresse au ministre des Ressources naturelles.
L'industrie forestiere de la Colombie-Britannique est en difficulte.
Les activites forestieres et les exportations sont en baisse.
M. Eric Lowther (Calgary-Centre, Ref.):
L'hon. Jim Peterson (secretaire d'Etat (Institutions financieres internationales), Lib.):
Le depute de Matapedia-Matane a la parole.
LA STATION DE RADIO DE LA GARDE COTIERE
M. Rene Canuel (Matapedia-Matane, BQ):
Monsieur le President, ma question s'adresse au ministre des Peches et des Oceans.
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
Nous avons egalement analyse en profondeur la technologie en question.
Dans le cas qui nous occupe, il n'y aura pas de perte d'emplois.
ENERGIE ATOMIQUE DU CANADA LIMITEE
M. Pat Martin (Winnipeg-Centre, NPD):
Le gel des mises a pied a ete leve mercredi.
La vente d'EACL a ete reportee.
Le depute de St. John's-Ouest.
M. Charlie Power (St. John's-Ouest, PC):
Le ministre des Peches va-t-il accepter immediatement cette recommandation?
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
Ce sont les memes accords internationaux qui nous permettent...
Le depute de Pontiac-Gatineau-Labelle.
M. Robert Bertrand (Pontiac-Gatineau-Labelle, Lib.):
Monsieur le President, ma question s'adresse au ministre de la Defense nationale.
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
Le gouvernement est a l'ecoute.
M. Lee Morrison (Cypress Hills-Grasslands, Ref.):
L'hon. David M. Collenette (ministre des Transports, Lib.):
Monsieur le President, le depute a une imagination debridee.
Le gouvernement n'a pas de programme rigide.
Si le depute a des suggestions a faire, il devrait les soumettre au comite.
Je sais que ce changement interesserait nombre de societes.
Mme Claudette Bradshaw (Moncton, Lib.):
Monsieur le President, j'invoque le Reglement.
Big Bob a toujours ete tres gentil avec nous et le personnel.
Comme c'est vendredi, la Chambre pourrait lui souhaiter une bonne retraite.
REPONSE DU GOUVERNEMENT A DES PETITIONS
JUSTICE ET DROITS DE LA PERSONNE
M. John McKay (Scarborough-Est, Lib.):
M. Bill Gilmour (Nanaimo-Alberni, Ref.):
LA TAXE SUR LES PRODUITS ET SERVICES
M. Bill Gilmour (Nanaimo-Alberni, Ref.):
M. John Cummins (Delta-South Richmond, Ref.):
Monsieur le President, on repondra aujourd'hui a la question no 5.
Monsieur le President, je suggere que les autres questions soient reservees.
M. John Cummins (Delta-South Richmond, Ref.):
Monsieur le President, j'ai pose quelques questions dont la reponse tarde a venir.
La reponse se fait toujours attendre.
Monsieur le President, il y a un delai determine pour repondre a ces questions.
Monsieur le President, j'ai note les numeros 31, 56 et 51.
Je ferai de mon mieux pour savoir ce qu'il advient de ces questions.
Le president suppleant (M. McClelland):
Est-ce que les autres questions peuvent rester au Feuilleton ?
Le president suppleant (M. McClelland):
M. Jay Hill (Prince George-Peace River, Ref.):
Ils ont du mal avec cela.
Permettez-moi d'expliquer le processus des paiements differes.
Il a ensuite 30 jours pour toucher son cheque.
En somme, ils ont 120 jours pour se faire payer.
Le reglement a ete modifie et le delai a ete ramene a 120 jours.
Tout cela nous rappelle le scenario ou l'acheteur doit prendre garde.
Il s'inquietent aussi du fait qu'il s'agit d'une contribution.
Ils doivent attendre la fin de l'annee avant de se retirer du programme.
Un agriculteur devra tenir des registres pour savoir exactement combien il a paye.
Nous ne savons pas si nous allons appuyer le projet de loi maintenant.
Celui-ci renferme de bons elements.
J'espere que le gouvernement ne rendra pas ces efforts inutiles.
Mme Helene Alarie (Louis-Hebert, BQ):
Mais quelle est donc l'essence du projet de loi C-26?
Dans le sommaire de ce projet de loi, on y lit, et je cite:
[...] afin de proteger les producteurs de cultures speciales qui detiennent des accuses de reception
N'est-ce pas la une facon de faire des nominations politiques?
Cela pourrait facilement devenir un conflit d'interets.
Autrement dit, ils ne doivent pas etre negociants.
M. Pat Martin (Winnipeg-Centre, NPD):
Le projet de loi modifie deux lois et en abroge une troisieme.
Il s'agit d'un projet de loi tres detaille et tres complexe.
Pour le moment, elle est reglementee par la Commission canadienne des grains.
Selon nous, les agriculteurs reagirent avec creativite a cette situation.
Les questions qu'il faut se poser sont les suivantes:
Dans l'affirmative, quel en sera le cout a long terme?
La reponse a ces questions est pour nous une source d'inquietudes graves.
Le projet de loi C-26 propose un systeme de licences pour ces petites compagnies.
Le fait d'obtenir une licence leur eviterait de deposer un important cautionnement de garantie.
Les agriculteurs doivent payer cette contribution.
Ils n'ont pas le choix.
Ce genre d'abonnement par defaut est tout simplement inacceptable aux yeux des agriculteurs.
Le fonds propose est extremement impopulaire aupres des agriculteurs.
Ils ne peuvent se permettre ce prelevement.
Cette fois-ci, c'est pour les cultures speciales.
Ce sont quelques-uns des problemes que posent ces mesures.
Cela inquiete grandement notre caucus.
Le discours du budget est reste completement muet sur l'agriculture.
Voila, a mon avis, qui est tres revelateur.
Ils ont completement oublie les Canadiens des regions rurales.
La reponse serait sans doute pas grand-chose.
C'est en fait beaucoup plus que cela.
M. Bob Kilger (Stormont-Dundas, Lib.):
Le president suppleant (M. McClelland):
La Chambre a entendu la motion.
Plait-il a la Chambre d'adopter la motion?
M. Roy Bailey (Souris-Moose Mountain, Ref.):
Le meme phenomene se produit partout, a la grandeur de la Saskatchewan.
Ils font tout leur possible.
Ils cherchent par tous les moyens a en tirer un petit revenu.
Voila que le gouvernement decide de presenter le projet de loi C-26.
Les producteurs de l'ouest sont desillusionnes.
Les producteurs ont le sentiment que le gouvernement les a trahis.
Ils se fichent eperdument du projet de loi C-26 en ce moment.
Ce qui leur importe, c'est de trouver le moyen de payer leurs factures.
C'est une cause desesperee.
Le gouvernement ne donne aucun espoir.
Les quatre pistes de curling des alentours sont deja fermees.
Et le pire est a venir.
Il constatera qu'il y avait des hopitaux dans 52 villes.
J'espere que les habitants de ma province sont a l'ecoute.
Le president suppleant (M. McClelland):
C'est la periode des questions et des observations.
Le depute de Winnipeg-Centre peut faire des observations s'il le veut.
Le fait que le budget ne fasse pas mention de l'agriculture est revelateur.
M. John Harvard (secretaire parlementaire du ministre de l'Agriculture et de l'Agroalimentaire, Lib.):
Le depute de Souris-Moose Mountain...
M. Jay Hill (Prince George-Peace River, Ref.):
Monsieur le President, j'invoque le Reglement.
Le president suppleant (M. McClelland):
Le depute de Prince George-Peace River a parfaitement raison, bien sur.
Il m'a semble parfaitement recevable que le secretaire parlementaire s'y reporte.
Le secretaire parlementaire peut avoir 30 secondes pour resumer.
Le depute de Winnipeg-Centre aura ensuite une minute pour repondre.
La question de l'economie rurale nous preoccupe tous.
Il a besoin d'espace et d'innovations.
Voila justement les avantages que le projet de loi C-26 tente de lui procurer.
Il faut mettre les choses au clair.
M. Gerald Keddy (South Shore, PC):
Le projet de loi comporte trois parties.
D'abord, il abroge la Loi sur les marches de grain a terme.
C'est un changement positif pour le secteur agricole.
Le Parti conservateur va appuyer cet aspect du projet de loi.
L'assureur serait la Societe pour l'expansion des exportations.
Je reviendrai tout a l'heure sur le regime d'assurance.
Je vais parler de certaines idees exprimees par des organisations agricoles de l'ouest.
La motion est ainsi libellee:
Attendu que le programme jouit d'un appui douteux de la part des agriculteurs;
Attendu que l'industrie des cultures speciales a su prosperer sans un tel programme;
Cette resolution decrit bien ce que le projet de loi C-26 ne fait pas.
Ce serait tout a fait sense.
J'espere que le gouvernement envisagera de donner ce choix aux producteurs.
En principe, le Parti conservateur appuie cette mesure legislative.
M. Jay Hill (Prince George-Peace River, Ref.):
Monsieur le President, j'ai quelques observations a faire.
C'est la meme chose dans le cas du projet de loi C-26.
J'ai une question a poser au depute conservateur.
Monsieur le President, le depute de Prince George-Peace River pose de tres bonnes questions.
Le projet de loi comporte certains elements sur lesquels nous accrochons.
Notre soutien reste conditionnel.
Le projet de loi a certaines qualites.
M. John Harvard (secretaire parlementaire du ministre de l'Agriculture et de l'Agroalimentaire, Lib.):
Monsieur le President, je voudrais faire deux breves observations.
Nous avons pris un bon depart aujourd'hui, un excellent depart meme.
Nous avons l'appui possible de l'opposition officielle.
Nous avons l'appui de principe du Parti progressiste conservateur.
C'est plutot un bon debut.
Cela me parait etre le bon moment de demander le renvoi au comite.
Le president suppleant (M. McClelland):
Le president suppleant (M. McClelland):
Ce que nous faisons ainsi, c'est permettre a la Chambre de s'ajourner normalement.
Avons-nous le consentement unanime pour dire qu'il est 14 h 30?
Le president suppleant (M. McClelland):
(La seance est levee a 13 h 02.)


text:
! { dg-do run }
! This verifies that the compiler will correctly accpet the name="", write out
! an empty string for the binding label to the module file, and then read it
! back in.  Also, during gfc_verify_binding_labels, the name="" will prevent
! any verification (since there is no label to verify).
module one
contains
  subroutine foo() bind(c)
  end subroutine foo
end module one

module two
contains
  ! This procedure is only used accessed in C
  ! as procedural pointer
  subroutine foo() bind(c, name="")
  end subroutine foo
end module two

use one, only: foo_one => foo
use two, only: foo_two => foo
end

! { dg-final { cleanup-modules "one two" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!--------------------------------------------------------------------
!     Created by Mahdi Esmaily Moghadam
!     contact memt63@gmail.com for reporting the bugs.
!--------------------------------------------------------------------
!
!     UC Copyright Notice
!     This software is Copyright ©2012 The Regents of the University of
!     California. All Rights Reserved.
!
!     Permission to copy and modify this software and its documentation
!     for educational, research and non-profit purposes, without fee,
!     and without a written agreement is hereby granted, provided that
!     the above copyright notice, this paragraph and the following three
!     paragraphs appear in all copies.
!
!     Permission to make commercial use of this software may be obtained
!     by contacting:
!     Technology Transfer Office
!     9500 Gilman Drive, Mail Code 0910
!     University of California
!     La Jolla, CA 92093-0910
!     (858) 534-5815
!     invent@ucsd.edu
!
!     This software program and documentation are copyrighted by The
!     Regents of the University of California. The software program and
!     documentation are supplied "as is", without any accompanying
!     services from The Regents. The Regents does not warrant that the
!     operation of the program will be uninterrupted or error-free. The
!     end-user understands that the program was developed for research
!     purposes and is advised not to rely exclusively on the program for
!     any reason.
!
!     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY
!     PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
!     DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS
!     SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
!     CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY
!     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
!     OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
!     SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE
!     UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
!     MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
!
!--------------------------------------------------------------------
!     This is to construct the structure for LS.
!--------------------------------------------------------------------

      SUBROUTINE FSILS_LS_CREATE(ls, LS_type, relTol, absTol, maxItr,   &
     &   dimKry, relTolIn, absTolIn, maxItrIn)
      INCLUDE "FSILS_STD.h"
      TYPE(FSILS_lsType), INTENT(INOUT) :: ls
      INTEGER(KIND=LSIP), INTENT(IN) :: LS_type
      REAL(KIND=LSRP), INTENT(IN), OPTIONAL :: relTol, absTol,          &
     &   relTolIn(2), absTolIn(2)
      INTEGER(KIND=LSIP), INTENT(IN), OPTIONAL :: maxItr, dimKry,       &
     &   maxItrIn(2)

      IF (ls%foC) THEN
         PRINT *, "FSILS: LS is not free You may use FSILS_LS_FREE to", &
     &      " free this structure"
         STOP "FSILS: FATAL ERROR"
      END IF

      ls%foC     = .TRUE.
      ls%LS_type = LS_type

      SELECT CASE (LS_type)
         CASE (LS_TYPE_NS)
            ls%RI%relTol = 0.4_LSRP
            ls%GM%relTol = 1.E-2_LSRP
            ls%CG%relTol = 0.2_LSRP
            ls%RI%mItr = 10
            ls%GM%mItr = 2
            ls%CG%mItr = 500
            ls%GM%sD   = 100
            ls%RI%sD   = 100
         CASE (LS_TYPE_GMRES)
            ls%RI%relTol = 0.1_LSRP
            ls%RI%mItr   = 4
            ls%RI%sD     = 250
         CASE (LS_TYPE_CG)
            ls%RI%reltol = 1.E-2_LSRP
            ls%RI%mItr   = 1000
         CASE (LS_TYPE_BICGS)
            ls%RI%reltol = 1.E-2_LSRP
            ls%RI%mItr   = 500
         CASE DEFAULT
            PRINT *, 'FSILS: LS_TYPE is not defined'
            STOP "FSILS: FATAL ERROR"
      END SELECT
      ls%RI%absTol = 1.E-10_LSRP
      ls%GM%absTol = 1.E-10_LSRP
      ls%CG%absTol = 1.E-10_LSRP

      IF (PRESENT(relTol)) ls%RI%relTol = relTol
      IF (PRESENT(absTol)) ls%RI%absTol = absTol
      IF (PRESENT(maxItr)) ls%RI%mItr   = maxItr

      IF (PRESENT(dimKry)) THEN
         ls%RI%sD = dimKry
         ls%GM%sD = dimKry
      END IF
      IF (PRESENT(relTolIn)) THEN
         ls%GM%relTol = relTolIn(1)
         ls%CG%relTol = relTolIn(2)
      END IF
      IF (PRESENT(absTolIn)) THEN
         ls%GM%absTol = absTolIn(1)
         ls%CG%absTol = absTolIn(2)
      END IF
      IF (PRESENT(maxItrIn)) THEN
         ls%GM%mItr = maxItrIn(1)
         ls%CG%mItr = maxItrIn(2)
      END IF

      RETURN
      END SUBROUTINE FSILS_LS_CREATE
!####################################################################
      SUBROUTINE FSILS_LS_FREE (ls)
      INCLUDE "FSILS_STD.h"
      TYPE(FSILS_lsType), INTENT(INOUT) :: ls

      IF (.NOT.ls%foC) THEN
         PRINT *, 'FSILS: LS is not created yet to be freed'
         STOP "FSILS: FATAL ERROR"
      END IF

      ls%foC  = .FALSE.

      RETURN
      END SUBROUTINE FSILS_LS_FREE
!####################################################################

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


text:
      subroutine nextrg(ibound,x,y,net,nodes,nb,xv,yv)
C...01/09/78
      common /gftrg1/job,xloc,yloc,step,dop,itrgl(3)
      dimension ibound(2,1),x(1),y(1),net(1),nodes(1)
      if(job.eq.0) goto 105
      i=iabs(job)
      if(nodes(i)) 105,105,106
  105 k=iabs(ibound(1,1))
      r=(xv-x(k))**2+(yv-y(k))**2
      do 76 i=1,nb
      l=iabs(ibound(1,i))
      ir=ibound(2,i)
      vx=y(l)-y(ir)
      vy=x(ir)-x(l)
      r1=cosin1(yv-y(ir),x(ir)-xv,vx,vy)
      r2=cosin1(-vy,vx,x(l)-xv,y(l)-yv)
      if(r1*r2.lt.0) goto 73
      r1=cosin1(x(l)-xv,y(l)-yv,vx,vy)**2/(vx**2+vy**2)
      if(r1.ge.r) goto 76
      r=r1
      k=ir
      go to 76
   73 do 75 j=1,2
      k1=iabs(ibound(j,i))
      r1=(xv-x(k1))**2+(yv-y(k1))**2
      if(r1.ge.r) goto75
      r=r1
      k=k1
   75 continue
   76 continue
      job=0
      goto 107
  106 k=iabs(job)
  107 l2=k
      ir1=k
      k1=nodes(k)
      k2=iabs(nodes(k+1))-1
      m=0
      i=k1
   10 l=net(i)
      if(l) 11,7,7
    7 i=i+1
      if(i.gt.k2) i=k1
      go to 10
   11 l1=iabs(l)
   83 if(i.ne.k2) goto 13
      i=k1
      go to 84
   13 i=i+1
   84 if(net(i).eq.0) goto 83
      ir2=net(i)
      r1=cosin1(y(l1)-y(l2),x(l2)-x(l1),x(l1)-xv,y(l1)-yv)
      if(r1.lt.0) goto 15
      r1=cosin1(y(ir1)-y(ir2),x(ir2)-x(ir1),x(ir1)-xv,y(ir1)-yv)
      if(r1) 15,17,17
   15 if(m.ne.0) goto 17
      m=1
      go to 7
   17 l=l1
      ir=l2
      call trg(l,ir,k,kl,kr,net,nodes,ibound,nb,x,y,xv,yv,ind)
      if(ind.ne.0) goto 20
      l=ir1
      ir=ir2
      if(job.le.0) job=0
      i=0
      call trg(l,ir,k,kl,kr,net,nodes,ibound,nb,x,y,xv,yv,ind)
      if(ind.ne.0) goto 20
      k1=nodes(l1)
      k2=iabs(nodes(l1+1))-1
      if(l1.eq.ir2) goto 105
      if(l1.eq.ir1) goto 105
      i=k1
   23 m1=iabs(net(i))
      if(m1.eq.l2) goto 24
      i=i+1
      go to 23
   24 l2=l1
   93 if(i.ne.k1) goto 90
      i=k2
      go to 92
   90 i=i-1
   92 if(net(i).eq.0) goto 93
      l1=iabs(net(i))
      k1=nodes(ir2)
      k2=iabs(nodes(ir2+1))-1
      i=k1
   29 m1=iabs(net(i))
      if(m1.eq.ir1) goto 30
      i=i+1
      go to 29
   30 ir1=ir2
   99 if(i.ne.k2) goto 95
      i=k1
      go to 97
   95 i=i+1
   97 if(net(i).eq.0) goto 99
      ir2=net(i)
      job=0
      go to 17
   20 itrgl(1)=l
      itrgl(2)=ir
      itrgl(3)=k
      l=iabs(l)
      ir=iabs(ir)
      k=iabs(k)
       job=-k
      if(kl.ge.0) goto 102
      if(kr.ge.0) goto 103
      job=k
      go to 103
  102 job=-ir
  103 call ilimtn(l,ir,ibound,nb)
      if(kl.le.0) goto 35
      call ilimtn(k,l,ibound,nb)
      go to 36
   35 j=iabs(kl)
      if(net(j).le.0) goto 36
      nb=nb+1
      ibound(1,nb)=-l
      ibound(2,nb)=k
   36 if(kr.le.0) goto 37
      call ilimtn(ir,k,ibound,nb)
       go to 39
   37 j=iabs(kr)
      if(net(j).lt.0) goto 39
      nb=nb+1
      ibound(1,nb)=-k
      ibound(2,nb)=ir
   39 m=0
   58 k1=nodes(l)
      k2=iabs(nodes(l+1))-1
      i=k1
   42 m1=net(i)
      if(iabs(m1).eq.ir) goto 41
      i=i+1
      go to 42
   41 if(i.ne.k1) goto 43
      m1=net(k2)
      go to 45
   43 m1=net(i-1)
   45 if(m1.gt.0) goto 47
      net(i)=0
      go to 50
   47 net(i)=-net(i)
   50 if(i.ne.k2) goto 48
      i=k1
      go to 51
   48 i=i+1
   51 if(net(i).le.0) net(i)=0
      nodes(l)=-nodes(l)
      do 54 i=k1,k2
      if(net(i).eq.0) goto 54
      nodes(l)=-nodes(l)
      go to 57
   54 continue
   57 m1=l
      l=ir
      ir=k
      k=m1
      m=m+1
      if(m.lt.3) goto 58
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
*DECK DERF
      DOUBLE PRECISION FUNCTION DERF (X)
C***BEGIN PROLOGUE  DERF
C***PURPOSE  Compute the error function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C8A, L5A1E
C***TYPE      DOUBLE PRECISION (ERF-S, DERF-D)
C***KEYWORDS  ERF, ERROR FUNCTION, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DERF(X) calculates the double precision error function for double
C precision argument X.
C
C Series for ERF        on the interval  0.          to  1.00000E+00
C                                        with weighted error   1.28E-32
C                                         log weighted error  31.89
C                               significant figures required  31.05
C                                    decimal places required  32.55
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, DERFC, INITDS
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900727  Added EXTERNAL statement.  (WRB)
C   920618  Removed space from variable name.  (RWC, WRB)
C***END PROLOGUE  DERF
      DOUBLE PRECISION X, ERFCS(21), SQEPS, SQRTPI, XBIG, Y, D1MACH,
     1  DCSEVL, DERFC
      LOGICAL FIRST
      EXTERNAL DERFC
      SAVE ERFCS, SQRTPI, NTERF, XBIG, SQEPS, FIRST
      DATA ERFCS(  1) / -.4904612123 4691808039 9845440333 76 D-1     /
      DATA ERFCS(  2) / -.1422612051 0371364237 8247418996 31 D+0     /
      DATA ERFCS(  3) / +.1003558218 7599795575 7546767129 33 D-1     /
      DATA ERFCS(  4) / -.5768764699 7674847650 8270255091 67 D-3     /
      DATA ERFCS(  5) / +.2741993125 2196061034 4221607914 71 D-4     /
      DATA ERFCS(  6) / -.1104317550 7344507604 1353812959 05 D-5     /
      DATA ERFCS(  7) / +.3848875542 0345036949 9613114981 74 D-7     /
      DATA ERFCS(  8) / -.1180858253 3875466969 6317518015 81 D-8     /
      DATA ERFCS(  9) / +.3233421582 6050909646 4029309533 54 D-10    /
      DATA ERFCS( 10) / -.7991015947 0045487581 6073747085 95 D-12    /
      DATA ERFCS( 11) / +.1799072511 3961455611 9672454866 34 D-13    /
      DATA ERFCS( 12) / -.3718635487 8186926382 3168282094 93 D-15    /
      DATA ERFCS( 13) / +.7103599003 7142529711 6899083946 66 D-17    /
      DATA ERFCS( 14) / -.1261245511 9155225832 4954248533 33 D-18    /
      DATA ERFCS( 15) / +.2091640694 1769294369 1705002666 66 D-20    /
      DATA ERFCS( 16) / -.3253973102 9314072982 3641600000 00 D-22    /
      DATA ERFCS( 17) / +.4766867209 7976748332 3733333333 33 D-24    /
      DATA ERFCS( 18) / -.6598012078 2851343155 1999999999 99 D-26    /
      DATA ERFCS( 19) / +.8655011469 9637626197 3333333333 33 D-28    /
      DATA ERFCS( 20) / -.1078892517 7498064213 3333333333 33 D-29    /
      DATA ERFCS( 21) / +.1281188399 3017002666 6666666666 66 D-31    /
      DATA SQRTPI / 1.772453850 9055160272 9816748334 115D0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DERF
      IF (FIRST) THEN
         NTERF = INITDS (ERFCS, 21, 0.1*REAL(D1MACH(3)))
         XBIG = SQRT(-LOG(SQRTPI*D1MACH(3)))
         SQEPS = SQRT(2.0D0*D1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS(X)
      IF (Y.GT.1.D0) GO TO 20
C
C ERF(X) = 1.0 - ERFC(X)  FOR  -1.0 .LE. X .LE. 1.0
C
      IF (Y.LE.SQEPS) DERF = 2.0D0*X/SQRTPI
      IF (Y.GT.SQEPS) DERF = X*(1.0D0 + DCSEVL (2.D0*X*X-1.D0,
     1  ERFCS, NTERF))
      RETURN
C
C ERF(X) = 1.0 - ERFC(X) FOR ABS(X) .GT. 1.0
C
 20   IF (Y.LE.XBIG) DERF = SIGN (1.0D0-DERFC(Y), X)
      IF (Y.GT.XBIG) DERF = SIGN (1.0D0, X)
C
      RETURN
      END

max_stars_repo_name: nilqed/spadlib

max_issues_repo_name: nilqed/spadlib

max_forks_repo_name: nilqed/spadlib


text:
C
C  Copyright (c) 2013-2020 The University of Tennessee and The University
C                          of Tennessee Research Foundation.  All rights
C                          reserved.
C  $COPYRIGHT$
C
C  Additional copyrights may follow
C
C  $HEADER$
C

C     Simulate the BLANK mode in FTMPI. Only peers communicating
C     directly with a dead process will notice the failures, and no
C     corrective actions will be taken globally (in other terms
C     communicator revocation is not an option). While this approach can
C     minimize the overhead, it makes difficult to recover the failure
C     of the master.

      subroutine error_handler(communicator, error_code)

C     Only the master will react to the failure. Every other process
C     will ignore all MPI_ERR_PROC_FAILED except if the dead process is
C     the master.

      implicit none
      include "mpif.h"
      include "fsolvergen.inc"

C input Args
      integer communicator, error_code !, stdargs(*)

C Local Variables
      integer rc, myrank, num_fails
C     Failed group handle and members array
      integer f_group, f_group_rank(maxworkers)
C     Original group handle and members array
      integer w_group, w_group_rank(maxworkers)

c     Who I was on the original communicator (for debugging purposes)
      call MPI_Comm_rank(communicator, myrank, rc)

c     A failed process: verify the error code
      if (error_code .eq. MPI_ERR_PROC_FAILED) then
         write(*,*) "MPI_ERR_PROC_FAILED triggered @ rank", myrank
         call MPIX_Comm_failure_ack(communicator, rc)

         if (myrank .eq. 0) then
c     Get failed procs group
            call MPIX_Comm_failure_get_acked(communicator, f_group, rc)
c     Get no. of failed procs
            call MPI_Group_size(f_group, num_fails, rc)
c     Get group from current communicator
            call MPI_Comm_group(communicator, w_group, rc)
c     Get ranks of failed procs
        do rc = 1, num_fails
            f_group_rank(rc) = rc-1
        enddo
        call MPI_Group_translate_ranks(f_group, num_fails,
     $           f_group_rank, w_group, w_group_rank, rc)
c     And mark them dead
            do rc = 1, num_fails
               write(*,*) "[", myrank, "] Peer",
     $              w_group_rank(rc), "dead"
               call mark_dead(w_group_rank(rc))
            enddo
            rc = 0
         endif
c     Strange things happen
      else
         write(*,*)"Something weird is happening.",
     $        "Error code", error_code, "on rank", myrank
      endif

      return
      end