text:
MODULE io_tools_module

  implicit none

  PRIVATE
  PUBLIC :: init_io_tools
  PUBLIC :: IOTools_readStringKeyword
  PUBLIC :: IOTools_readIntegerKeyword
  PUBLIC :: IOTools_readReal8Keyword
  PUBLIC :: IOTools_readReal8Keywords
  PUBLIC :: IOTools_readIntegerString
  PUBLIC :: IOTools_findKeyword

  integer,parameter :: max_trial_read = 10000

#ifndef _NOMPI_
  include 'mpif.h'
#endif

  integer :: myrank = 0
  integer :: unit_default = 1
  integer :: unit_input_result = 110
  logical :: flag_init = .false.

  INTERFACE IOTools_readIntegerKeyword
     MODULE PROCEDURE IOTools_readIntegerKeyword_sca &
                     ,IOTools_readIntegerKeyword_vec
  END INTERFACE

CONTAINS

  SUBROUTINE init_io_tools( myrank_in, unit_in )
    implicit none
    integer,intent(IN) :: myrank_in, unit_in
    myrank = myrank_in
    unit_default = unit_in
!    if ( myrank == 0 ) open(unit_input_result,file="input_result")
    flag_init = .true.
  END SUBROUTINE init_io_tools

  SUBROUTINE check_init
    implicit none
    integer :: ierr
#ifdef _NOMPI_
    flag_init = .true.
    return
#else
    include 'mpif.h'
    if ( .not.flag_init ) then
       call MPI_COMM_RANK( MPI_COMM_WORLD, myrank, ierr )
    end if
    flag_init = .true.
#endif
  END SUBROUTINE check_init

  SUBROUTINE IOTools_readStringKeyword( keyword, variable, unit_in )
    implicit none
    character(*),intent(IN) :: keyword
    character(*),intent(INOUT) :: variable
    integer,optional,intent(IN) :: unit_in
    character(10) :: cbuf,ckey
    integer :: i,unit
    call check_init
    unit=unit_default ; if ( present(unit_in) ) unit=unit_in
    if ( myrank == 0 ) then
       rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             backspace(unit)
             read(unit,*) cbuf, variable
             write(*,'(1x,A10," : ",A10)') keyword, variable
             exit
          end if
       end do ! i
999    continue
    end if
#ifndef _NOMPI_
    call MPI_BCAST( variable,len(variable),MPI_CHARACTER,0,MPI_COMM_WORLD,i )
#endif
  END SUBROUTINE IOTools_readStringKeyword

  SUBROUTINE IOTools_readIntegerKeyword_sca( keyword, variable, unit_in )
    implicit none
    character(*),intent(IN) :: keyword
    integer,intent(INOUT) :: variable
    integer,optional,intent(IN) :: unit_in
    character(10) :: cbuf,ckey
    integer :: i,unit
    call check_init
    unit=unit_default ; if ( present(unit_in) ) unit=unit_in
    if ( myrank == 0 ) then
       rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             backspace(unit)
             read(unit,*) cbuf,variable
             write(*,'(1x,A10," : ",3I10)') keyword,variable
             exit
          end if
       end do ! i
999    continue
       !write(unit_input_result,'(a10,3x,i10)') keyword,variable
    end if
#ifndef _NOMPI_
    call MPI_BCAST(variable,1,MPI_INTEGER,0,MPI_COMM_WORLD,i)
#endif
  END SUBROUTINE IOTools_readIntegerKeyword_sca

  SUBROUTINE IOTools_readIntegerKeyword_vec( keyword, variables, unit_in )
    implicit none
    character(*),intent(IN) :: keyword
    integer,intent(INOUT) :: variables(:)
    integer,optional,intent(IN) :: unit_in
    character(10) :: cbuf,ckey
    integer :: i,unit
    call check_init
    unit=unit_default ; if ( present(unit_in) ) unit=unit_in
    if ( myrank == 0 ) then
       rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             backspace(unit)
             read(unit,*) cbuf,variables(:)
             write(*,'(1x,A10," : ",7I10)') keyword,variables(:)
             exit
          end if
       end do ! i
999    continue
    end if
#ifndef _NOMPI_
    call MPI_BCAST(variables,size(variables),MPI_INTEGER,0,MPI_COMM_WORLD,i)
#endif
  END SUBROUTINE IOTools_readIntegerKeyword_vec

  SUBROUTINE IOTools_readReal8Keyword( keyword, variable, unit_in )
    implicit none
    character(*),intent(IN) :: keyword
    real(8),intent(INOUT) :: variable
    integer,optional,intent(IN) :: unit_in
    character(10) :: cbuf,ckey
    integer :: i,unit
    call check_init
    unit=unit_default ; if ( present(unit_in) ) unit=unit_in
    if ( myrank == 0 ) then
       rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             backspace(unit)
             read(unit,*) cbuf,variable
             if ( abs(variable) < 1.d-2 ) then
                write(*,'(1x,A10," : ",ES15.7)') keyword,variable
             else
                write(*,'(1x,A10," : ",F15.10)') keyword,variable
             end if
             exit
          end if
       end do ! i
999    continue
    end if
#ifndef _NOMPI_
    call MPI_BCAST(variable,1,MPI_REAL8,0,MPI_COMM_WORLD,i)
#endif
  END SUBROUTINE IOTools_readReal8Keyword

  SUBROUTINE IOTools_readReal8Keywords( keyword, variables, unit_in )
    implicit none
    character(*),intent(IN) :: keyword
    real(8),intent(INOUT) :: variables(:)
    integer,optional,intent(IN) :: unit_in
    character(10) :: cbuf,ckey
    integer :: i,unit
    call check_init
    unit=unit_default ; if ( present(unit_in) ) unit=unit_in
    if ( myrank == 0 ) then
       rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             backspace(unit)
             read(unit,*) cbuf,variables(:)
             if ( any(abs(variables) < 1.d-2) ) then
                write(*,'(1x,A10," : ",3ES15.7)') keyword,variables(:)
             else
                write(*,'(1x,A10," : ",3F15.10)') keyword,variables(:)
             end if
             exit
          end if
       end do ! i
999    continue
    end if
#ifndef _NOMPI_
    call MPI_BCAST(variables,size(variables),MPI_REAL8,0,MPI_COMM_WORLD,i)
#endif
  END SUBROUTINE IOTools_readReal8Keywords

  SUBROUTINE IOTools_readIntegerString( keyword, variable1, variable2, u, norewind )
    implicit none
    character(*),intent(IN) :: keyword
    integer,intent(INOUT) :: variable1
    character(*),intent(INOUT) :: variable2
    integer,optional,intent(IN) :: u
    logical,optional,intent(IN) :: norewind
    character(10) :: cbuf,ckey
    integer :: i,unit
    call check_init
    unit=unit_default ; if ( present(u) ) unit=u
    if ( myrank == 0 ) then
       if ( .not.present(norewind) ) rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             backspace(unit)
             read(unit,*) cbuf,variable1, variable2
             write(*,'(1x,A10," : ",i4,2x,a20)') keyword,variable1,variable2
             exit
          end if
       end do ! i
999    continue
    end if
#ifndef _NOMPI_
    call MPI_BCAST(variable1,1,MPI_INTEGER,0,MPI_COMM_WORLD,i)
    call MPI_BCAST(variable2,len(variable2),MPI_CHARACTER,0,MPI_COMM_WORLD,i)
#endif
  END SUBROUTINE IOTools_readIntegerString

  SUBROUTINE IOTools_findKeyword( keyword, hasKeyword, unit_out, flag_bcast )
    implicit none
    character(*),intent(IN) :: keyword
    logical,intent(OUT) :: hasKeyword
    integer,optional,intent(OUT) :: unit_out
    logical,optional,intent(IN) :: flag_bcast
    integer :: i,unit
    character(10) :: cbuf,ckey
    call check_init
    unit=unit_default
    if ( present(unit_out) ) unit_out=unit
    hasKeyword = .false.
    if ( myrank == 0 ) then
       rewind unit
       do i=1,max_trial_read
          read(unit,*,END=999) cbuf
          call convertToCapital(cbuf,ckey)
          if ( ckey == keyword ) then
             hasKeyword = .true.
             exit
          end if
       end do ! i
       if ( hasKeyword ) write(*,'(1x,A10)') keyword
999    continue
    end if
    if ( present(flag_bcast) ) then
#ifndef _NOMPI_
       call MPI_BCAST( hasKeyword, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, i )
#endif
    end if
  END SUBROUTINE IOTools_findKeyword

  SUBROUTINE convertToCapital(cbuf,CKEY)
    implicit none
    character(*),intent(IN)  :: cbuf
    character(*),intent(OUT) :: CKEY
    integer :: j,k,n
    n=len_trim(cbuf)
    CKEY=""
    do j=1,n
      k=iachar( cbuf(j:j) )
      if ( k >= 97 ) k=k-32
      CKEY(j:j) = achar(k)
    end do
  END SUBROUTINE convertToCapital

END MODULE io_tools_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
MODULE common_mpi
!=======================================================================
!
! [PURPOSE:] General MPI procedures
!
! [HISTORY:]
!   09/06/2005 Takemasa MIYOSHI  created
!
!=======================================================================
  USE common, only: r_size, r_dble, r_sngl
  IMPLICIT NONE
  PUBLIC
  INCLUDE 'mpif.h'

  INTEGER,SAVE :: nprocs
  INTEGER,SAVE :: myrank
  INTEGER,SAVE :: MPI_r_size

CONTAINS
SUBROUTINE initialize_mpi
  IMPLICIT NONE
  INTEGER :: ierr
  CALL MPI_INIT(ierr)
  CALL MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr)
  CALL MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)
  WRITE(6,'(A,I3.3,A,I3.3)') 'Hello from MYRANK ',myrank,'/',nprocs-1
  IF(r_size == r_dble) THEN
    MPI_r_size = MPI_DOUBLE_PRECISION
  ELSE IF(r_size == r_sngl) THEN
    MPI_r_size = MPI_REAL
  END IF

  RETURN
END SUBROUTINE initialize_mpi

SUBROUTINE finalize_mpi
  IMPLICIT NONE
  INTEGER :: ierr
  CALL MPI_FINALIZE(ierr)

  RETURN
END SUBROUTINE finalize_mpi

END MODULE common_mpi

max_stars_repo_name: skotsuki/speedy-lpf

max_issues_repo_name: skotsuki/speedy-lpf

max_forks_repo_name: skotsuki/speedy-lpf


text:
! sum and other array functions have an optional dim argument
program main
implicit none
integer :: v(2,3)   ! declare 2x3 matrix
v(1,:) = [1,2,3]    ! set 1st row
v(2,:) = 10*v(1,:)  ! set 2nd row to [10,20,30]
print*,sum(v)       ! sum all elements: 66
print*,sum(v,dim=1) ! sum each row: 6 60
print*,sum(v,dim=2) ! sum each col: 11 22 33
end program main

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
C     This File is Automatically generated by ALOHA
C     The process calculated in this file is:
C     P(1,2)*P(2,1) - (13*P(1,3)*P(2,1))/62. - (13*P(1,2)*P(2,3))/62.
C      + (P(1,3)*P(2,3))/62. - (20*P(-1,1)*P(-1,2)*Metric(1,2))/31. +
C      (19*P(-1,2)**2*Metric(1,2))/62. + (2*P(-1,1)*P(-1,3)*Metric(1,2))
C     /31. + (23*P(-1,2)*P(-1,3)*Metric(1,2))/62. -
C      (13*P(-1,3)**2*Metric(1,2))/62.
C
      SUBROUTINE VVS8_0(V1, V2, S3, COUP,VERTEX)
      IMPLICIT NONE
      COMPLEX*16 CI
      PARAMETER (CI=(0D0,1D0))
      COMPLEX*16 V2(*)
      COMPLEX*16 TMP2
      COMPLEX*16 S3(*)
      COMPLEX*16 TMP11
      COMPLEX*16 V1(*)
      COMPLEX*16 TMP12
      COMPLEX*16 P2(0:3)
      COMPLEX*16 P3(0:3)
      COMPLEX*16 TMP16
      COMPLEX*16 TMP15
      COMPLEX*16 TMP14
      COMPLEX*16 P1(0:3)
      COMPLEX*16 VERTEX
      COMPLEX*16 COUP
      COMPLEX*16 TMP9
      COMPLEX*16 TMP17
      COMPLEX*16 TMP13
      COMPLEX*16 TMP8
      P1(0) = V1(1)
      P1(1) = V1(2)
      P1(2) = V1(3)
      P1(3) = V1(4)
      P2(0) = V2(1)
      P2(1) = V2(2)
      P2(2) = V2(3)
      P2(3) = V2(4)
      P3(0) = S3(1)
      P3(1) = S3(2)
      P3(2) = S3(3)
      P3(3) = S3(4)
      TMP9 = (V2(5)*P3(0)-V2(6)*P3(1)-V2(7)*P3(2)-V2(8)*P3(3))
      TMP8 = (V2(5)*P1(0)-V2(6)*P1(1)-V2(7)*P1(2)-V2(8)*P1(3))
      TMP2 = (P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3))
      TMP15 = (P1(0)*P3(0)-P1(1)*P3(1)-P1(2)*P3(2)-P1(3)*P3(3))
      TMP14 = (P2(0)*P2(0)-P2(1)*P2(1)-P2(2)*P2(2)-P2(3)*P2(3))
      TMP17 = (P3(0)*P3(0)-P3(1)*P3(1)-P3(2)*P3(2)-P3(3)*P3(3))
      TMP16 = (P2(0)*P3(0)-P2(1)*P3(1)-P2(2)*P3(2)-P2(3)*P3(3))
      TMP11 = (P2(0)*V1(5)-P2(1)*V1(6)-P2(2)*V1(7)-P2(3)*V1(8))
      TMP13 = (V2(5)*V1(5)-V2(6)*V1(6)-V2(7)*V1(7)-V2(8)*V1(8))
      TMP12 = (P3(0)*V1(5)-P3(1)*V1(6)-P3(2)*V1(7)-P3(3)*V1(8))
      VERTEX = COUP*(-13D0/62D0 * S3(5)*(TMP13*23D0/13D0*(-40D0/23D0 *
     $  CI*(TMP2)+19D0/23D0 * CI*(TMP14)+CI*(TMP16)+4D0/23D0 * CI
     $ *(TMP15)-13D0/23D0 * CI*(TMP17))+(TMP11*(-CI*(TMP9)+62D0/13D0 *
     $  CI*(TMP8))+TMP12*(-CI*(TMP8)+1D0/13D0 * CI*(TMP9)))))
      END

ext: f

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
! This file was automatically generated by SWIG (http://www.swig.org).
! Version 4.0.0
!
! Do not make changes to this file unless you know what you are doing--modify
! the SWIG interface file instead.

! ---------------------------------------------------------------
! Programmer(s): Auto-generated by swig.
! ---------------------------------------------------------------
! SUNDIALS Copyright Start
! Copyright (c) 2002-2021, Lawrence Livermore National Security
! and Southern Methodist University.
! All rights reserved.
!
! See the top-level LICENSE and NOTICE files for details.
!
! SPDX-License-Identifier: BSD-3-Clause
! SUNDIALS Copyright End
! ---------------------------------------------------------------

module fsunmatrix_sparse_mod
 use, intrinsic :: ISO_C_BINDING
 use fsundials_matrix_mod
 use fsundials_types_mod
 use fsundials_nvector_mod
 use fsundials_types_mod
 implicit none
 private

 ! DECLARATION CONSTRUCTS
 integer(C_INT), parameter, public :: CSC_MAT = 0_C_INT
 integer(C_INT), parameter, public :: CSR_MAT = 1_C_INT
 public :: FSUNSparseMatrix
 public :: FSUNSparseFromDenseMatrix
 public :: FSUNSparseFromBandMatrix
 public :: FSUNSparseMatrix_ToCSR
 public :: FSUNSparseMatrix_ToCSC
 public :: FSUNSparseMatrix_Realloc
 public :: FSUNSparseMatrix_Reallocate
 public :: FSUNSparseMatrix_Print
 public :: FSUNSparseMatrix_Rows
 public :: FSUNSparseMatrix_Columns
 public :: FSUNSparseMatrix_NNZ
 public :: FSUNSparseMatrix_NP
 public :: FSUNSparseMatrix_SparseType
 public :: FSUNSparseMatrix_Data
 public :: FSUNSparseMatrix_IndexValues
 public :: FSUNSparseMatrix_IndexPointers
 public :: FSUNMatGetID_Sparse
 public :: FSUNMatClone_Sparse
 public :: FSUNMatDestroy_Sparse
 public :: FSUNMatZero_Sparse
 public :: FSUNMatCopy_Sparse
 public :: FSUNMatScaleAdd_Sparse
 public :: FSUNMatScaleAddI_Sparse
 public :: FSUNMatMatvec_Sparse
 public :: FSUNMatSpace_Sparse

! WRAPPER DECLARATIONS
interface
function swigc_FSUNSparseMatrix(farg1, farg2, farg3, farg4) &
bind(C, name="_wrap_FSUNSparseMatrix") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T), intent(in) :: farg1
integer(C_INT64_T), intent(in) :: farg2
integer(C_INT64_T), intent(in) :: farg3
integer(C_INT), intent(in) :: farg4
type(C_PTR) :: fresult
end function

function swigc_FSUNSparseFromDenseMatrix(farg1, farg2, farg3) &
bind(C, name="_wrap_FSUNSparseFromDenseMatrix") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
real(C_DOUBLE), intent(in) :: farg2
integer(C_INT), intent(in) :: farg3
type(C_PTR) :: fresult
end function

function swigc_FSUNSparseFromBandMatrix(farg1, farg2, farg3) &
bind(C, name="_wrap_FSUNSparseFromBandMatrix") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
real(C_DOUBLE), intent(in) :: farg2
integer(C_INT), intent(in) :: farg3
type(C_PTR) :: fresult
end function

function swigc_FSUNSparseMatrix_ToCSR(farg1, farg2) &
bind(C, name="_wrap_FSUNSparseMatrix_ToCSR") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FSUNSparseMatrix_ToCSC(farg1, farg2) &
bind(C, name="_wrap_FSUNSparseMatrix_ToCSC") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FSUNSparseMatrix_Realloc(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_Realloc") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT) :: fresult
end function

function swigc_FSUNSparseMatrix_Reallocate(farg1, farg2) &
bind(C, name="_wrap_FSUNSparseMatrix_Reallocate") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT64_T), intent(in) :: farg2
integer(C_INT) :: fresult
end function

subroutine swigc_FSUNSparseMatrix_Print(farg1, farg2) &
bind(C, name="_wrap_FSUNSparseMatrix_Print")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
end subroutine

function swigc_FSUNSparseMatrix_Rows(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_Rows") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT64_T) :: fresult
end function

function swigc_FSUNSparseMatrix_Columns(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_Columns") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT64_T) :: fresult
end function

function swigc_FSUNSparseMatrix_NNZ(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_NNZ") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT64_T) :: fresult
end function

function swigc_FSUNSparseMatrix_NP(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_NP") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT64_T) :: fresult
end function

function swigc_FSUNSparseMatrix_SparseType(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_SparseType") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT) :: fresult
end function

function swigc_FSUNSparseMatrix_Data(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_Data") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

function swigc_FSUNSparseMatrix_IndexValues(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_IndexValues") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

function swigc_FSUNSparseMatrix_IndexPointers(farg1) &
bind(C, name="_wrap_FSUNSparseMatrix_IndexPointers") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

function swigc_FSUNMatGetID_Sparse(farg1) &
bind(C, name="_wrap_FSUNMatGetID_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT) :: fresult
end function

function swigc_FSUNMatClone_Sparse(farg1) &
bind(C, name="_wrap_FSUNMatClone_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

subroutine swigc_FSUNMatDestroy_Sparse(farg1) &
bind(C, name="_wrap_FSUNMatDestroy_Sparse")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
end subroutine

function swigc_FSUNMatZero_Sparse(farg1) &
bind(C, name="_wrap_FSUNMatZero_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT) :: fresult
end function

function swigc_FSUNMatCopy_Sparse(farg1, farg2) &
bind(C, name="_wrap_FSUNMatCopy_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FSUNMatScaleAdd_Sparse(farg1, farg2, farg3) &
bind(C, name="_wrap_FSUNMatScaleAdd_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
integer(C_INT) :: fresult
end function

function swigc_FSUNMatScaleAddI_Sparse(farg1, farg2) &
bind(C, name="_wrap_FSUNMatScaleAddI_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FSUNMatMatvec_Sparse(farg1, farg2, farg3) &
bind(C, name="_wrap_FSUNMatMatvec_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
integer(C_INT) :: fresult
end function

function swigc_FSUNMatSpace_Sparse(farg1, farg2, farg3) &
bind(C, name="_wrap_FSUNMatSpace_Sparse") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
integer(C_INT) :: fresult
end function

end interface

contains
 ! MODULE SUBPROGRAMS
function FSUNSparseMatrix(m, n, nnz, sparsetype) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(SUNMatrix), pointer :: swig_result
integer(C_INT64_T), intent(in) :: m
integer(C_INT64_T), intent(in) :: n
integer(C_INT64_T), intent(in) :: nnz
integer(C_INT), intent(in) :: sparsetype
type(C_PTR) :: fresult
integer(C_INT64_T) :: farg1
integer(C_INT64_T) :: farg2
integer(C_INT64_T) :: farg3
integer(C_INT) :: farg4

farg1 = m
farg2 = n
farg3 = nnz
farg4 = sparsetype
fresult = swigc_FSUNSparseMatrix(farg1, farg2, farg3, farg4)
call c_f_pointer(fresult, swig_result)
end function

function FSUNSparseFromDenseMatrix(a, droptol, sparsetype) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(SUNMatrix), pointer :: swig_result
type(SUNMatrix), target, intent(inout) :: a
real(C_DOUBLE), intent(in) :: droptol
integer(C_INT), intent(in) :: sparsetype
type(C_PTR) :: fresult
type(C_PTR) :: farg1
real(C_DOUBLE) :: farg2
integer(C_INT) :: farg3

farg1 = c_loc(a)
farg2 = droptol
farg3 = sparsetype
fresult = swigc_FSUNSparseFromDenseMatrix(farg1, farg2, farg3)
call c_f_pointer(fresult, swig_result)
end function

function FSUNSparseFromBandMatrix(a, droptol, sparsetype) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(SUNMatrix), pointer :: swig_result
type(SUNMatrix), target, intent(inout) :: a
real(C_DOUBLE), intent(in) :: droptol
integer(C_INT), intent(in) :: sparsetype
type(C_PTR) :: fresult
type(C_PTR) :: farg1
real(C_DOUBLE) :: farg2
integer(C_INT) :: farg3

farg1 = c_loc(a)
farg2 = droptol
farg3 = sparsetype
fresult = swigc_FSUNSparseFromBandMatrix(farg1, farg2, farg3)
call c_f_pointer(fresult, swig_result)
end function

function FSUNSparseMatrix_ToCSR(a, bout) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR), target, intent(inout) :: bout
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(a)
farg2 = c_loc(bout)
fresult = swigc_FSUNSparseMatrix_ToCSR(farg1, farg2)
swig_result = fresult
end function

function FSUNSparseMatrix_ToCSC(a, bout) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR), target, intent(inout) :: bout
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(a)
farg2 = c_loc(bout)
fresult = swigc_FSUNSparseMatrix_ToCSC(farg1, farg2)
swig_result = fresult
end function

function FSUNSparseMatrix_Realloc(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_Realloc(farg1)
swig_result = fresult
end function

function FSUNSparseMatrix_Reallocate(a, nnz) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT64_T), intent(in) :: nnz
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT64_T) :: farg2

farg1 = c_loc(a)
farg2 = nnz
fresult = swigc_FSUNSparseMatrix_Reallocate(farg1, farg2)
swig_result = fresult
end function

subroutine FSUNSparseMatrix_Print(a, outfile)
use, intrinsic :: ISO_C_BINDING
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR) :: outfile
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(a)
farg2 = outfile
call swigc_FSUNSparseMatrix_Print(farg1, farg2)
end subroutine

function FSUNSparseMatrix_Rows(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT64_T) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_Rows(farg1)
swig_result = fresult
end function

function FSUNSparseMatrix_Columns(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT64_T) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_Columns(farg1)
swig_result = fresult
end function

function FSUNSparseMatrix_NNZ(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT64_T) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_NNZ(farg1)
swig_result = fresult
end function

function FSUNSparseMatrix_NP(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT64_T) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_NP(farg1)
swig_result = fresult
end function

function FSUNSparseMatrix_SparseType(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_SparseType(farg1)
swig_result = fresult
end function

function FSUNSparseMatrix_Data(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), dimension(:), pointer :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_Data(farg1)
call c_f_pointer(fresult, swig_result, [1])
end function

function FSUNSparseMatrix_IndexValues(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T), dimension(:), pointer :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_IndexValues(farg1)
call c_f_pointer(fresult, swig_result, [1])
end function

function FSUNSparseMatrix_IndexPointers(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T), dimension(:), pointer :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNSparseMatrix_IndexPointers(farg1)
call c_f_pointer(fresult, swig_result, [1])
end function

function FSUNMatGetID_Sparse(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(SUNMatrix_ID) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNMatGetID_Sparse(farg1)
swig_result = fresult
end function

function FSUNMatClone_Sparse(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(SUNMatrix), pointer :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNMatClone_Sparse(farg1)
call c_f_pointer(fresult, swig_result)
end function

subroutine FSUNMatDestroy_Sparse(a)
use, intrinsic :: ISO_C_BINDING
type(SUNMatrix), target, intent(inout) :: a
type(C_PTR) :: farg1

farg1 = c_loc(a)
call swigc_FSUNMatDestroy_Sparse(farg1)
end subroutine

function FSUNMatZero_Sparse(a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(a)
fresult = swigc_FSUNMatZero_Sparse(farg1)
swig_result = fresult
end function

function FSUNMatCopy_Sparse(a, b) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(SUNMatrix), target, intent(inout) :: b
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(a)
farg2 = c_loc(b)
fresult = swigc_FSUNMatCopy_Sparse(farg1, farg2)
swig_result = fresult
end function

function FSUNMatScaleAdd_Sparse(c, a, b) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
real(C_DOUBLE), intent(in) :: c
type(SUNMatrix), target, intent(inout) :: a
type(SUNMatrix), target, intent(inout) :: b
integer(C_INT) :: fresult
real(C_DOUBLE) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c
farg2 = c_loc(a)
farg3 = c_loc(b)
fresult = swigc_FSUNMatScaleAdd_Sparse(farg1, farg2, farg3)
swig_result = fresult
end function

function FSUNMatScaleAddI_Sparse(c, a) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
real(C_DOUBLE), intent(in) :: c
type(SUNMatrix), target, intent(inout) :: a
integer(C_INT) :: fresult
real(C_DOUBLE) :: farg1
type(C_PTR) :: farg2

farg1 = c
farg2 = c_loc(a)
fresult = swigc_FSUNMatScaleAddI_Sparse(farg1, farg2)
swig_result = fresult
end function

function FSUNMatMatvec_Sparse(a, x, y) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: y
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(a)
farg2 = c_loc(x)
farg3 = c_loc(y)
fresult = swigc_FSUNMatMatvec_Sparse(farg1, farg2, farg3)
swig_result = fresult
end function

function FSUNMatSpace_Sparse(a, lenrw, leniw) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(SUNMatrix), target, intent(inout) :: a
integer(C_LONG), dimension(*), target, intent(inout) :: lenrw
integer(C_LONG), dimension(*), target, intent(inout) :: leniw
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(a)
farg2 = c_loc(lenrw(1))
farg3 = c_loc(leniw(1))
fresult = swigc_FSUNMatSpace_Sparse(farg1, farg2, farg3)
swig_result = fresult
end function

end module

max_stars_repo_name: bassenj/sundials

max_issues_repo_name: bassenj/sundials

max_forks_repo_name: bassenj/sundials


text:
      SUBROUTINE DROTG(DA,DB,C,S)

C     CONSTRUCT GIVENS PLANE ROTATION.
C     JACK DONGARRA, LINPACK, 3/11/78.
C                    MODIFIED 9/27/86.

      DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z,ONE,ZERO
      DATA ONE, ZERO /1.0D+00, 0.0D+00/

      ROE = DB
      IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
      SCALE = ABS(DA) + ABS(DB)
      IF( SCALE .NE. ZERO ) GO TO 10
         C = ONE
         S = ZERO
         R = ZERO
         GO TO 20
   10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
      R = SIGN(ONE,ROE)*R
      C = DA/R
      S = DB/R
   20 Z = S
      IF( ABS(C) .GT. ZERO .AND. ABS(C) .LE. S ) Z = ONE/C
      DA = R
      DB = Z
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
Davis Neighbors Night Out is a celebration of community. The event seeks to unite neighborhoods by creating an opportunity for all residents to participate in a block party, barbecue, ice cream social, pot luck or other public event. The goal of Davis Neighbors Night Out is to reduce the potential for misunderstandings or conflicts by promoting familiarity and open communication among neighbors. The event, which is scheduled for Sunday, October 11th, 2009 is timed to enable neighborhoods to include new or returning UC Davis students in the festivities.
Date: October 11th, 2009
Location: Everywhere in Davis
Sponsors:
 Associated Students, University of California, Davis (ASUCD)
 City of Davis
 University of California, Davis
Website: http://www.cityofdavis.org/neighbors/index.cfm
Interactive Find a Party Map: http://maps.ci.davis.ca.us/gis/google/nnout/
Contact: Stacey Winton at 7575661.
Facebook: http://www.facebook.com/event.php?eid131663718481&refsearch&sid536264065.980587515..1

Student Incentives
In an effort to encourage more students to host events and attend there are a number of student incentives.
     Hosting a party puts you in the running for a free Kaplan course ($1,900 value), a CPR course from the Red Cross ($37 value) and a $25 gift certificate to Yolo Berry. All student party hosts will receive a 2009 DNNO tshirt courtesy of Ink Monkey Graphics and a raffle kit for your student guests.
     Student attendees can also enter to win great prizes, including gift certificates from the Yogurt Shack, Woodstock’s, Steve’s Pizza and InNOut. Other prizes include Transportation and Parking Services (TAPS) and Health Education and Promotion giveaway baskets.

There is still time for students to host an event. If you do something like a potluck its easy to organize and city staff will help you.

2009 Planning Committee:
Stacey Winton, City of Davis
Michele Reynolds, City of Davis

Gary Sandy, UC Davis
Michelle Johnston, UC Davis

Steven D. Lee, ASUCD
Chad Roberts, IFC
Molly Bechtel, Panhellenic Council

Articles
http://www.davisenterprise.com/articles/2007/10/14/news/094new1.txt
http://media.www.californiaaggie.com/media/storage/paper981/news/2007/10/12/Features/Hello.Neighbor.A.Citywide.Meet.And.Greet.In.Davis3029145.shtml
http://www.dateline.ucdavis.edu/dl_detail.lasso?id9781

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C
C +-----------------------------------------------------------------+
C |                                                                 |
C |                Copyright (C) 1986 by UCAR                       |
C |        University Corporation for Atmospheric Research          |
C |                    All Rights Reserved                          |
C |                                                                 |
C |                 NCARGRAPHICS  Version 1.00                      |
C |                                                                 |
C +-----------------------------------------------------------------+
C
C
C ---------------------------------------------------------------------
C
      SUBROUTINE AGEXUS (SVAL,ZMIN,ZMAX,ZLOW,ZHGH,
     +                               ZDRA,NVIZ,IIVZ,NEVZ,IIEZ,UMIN,UMAX)
C
      DIMENSION SVAL(2),ZDRA(1)
C
C The routine AGEXUS is used by AGSTUP to determine tentative values of
C the user-window edge coordinates.  Its arguments are as follows:
C
C -- SVAL is the array of special values.
C
C -- ZMIN and ZMAX are user-supplied minimum and maximum values of the
C    data x (or y) coordinates.
C
C -- ZLOW and ZHGH are, respectively, the smallest and largest data
C    values to be considered in choosing the minimum and maximum, if
C    those values, as given by the user, are null.
C
C -- ZDRA, NVIZ, IIVZ, NEVZ, and IIEZ specify the array of x (or y)
C    data coordinates (see AGMAXI or AGMINI for complete description).
C
C -- UMIN and UMAX are returned with tentative minimum and maximum
C    values for use at the appropriate user-window edges (left/right
C    or bottom/top).
C
C The following common block contains AUTOGRAPH variables which are
C not control parameters.  The only one used here is SMRL, which is a
C (machine-dependent) small real which, when added to a number in the
C range (1,10), will round it upward without seriously affecting the
C leading significant digits.  The object of this is to get rid of
C strings of nines.
C
      COMMON /AGORIP/ SMRL , ISLD , MWCL,MWCM,MWCE,MDLA,MWCD,MWDQ ,
     +                INIF
C
C Assume initially that the user has provided actual values to be used.
C
      UMIN=ZMIN
      UMAX=ZMAX
C
C If either of the values is null, replace it by a data-based value.
C
      IF (UMIN.EQ.SVAL(1).OR.UMIN.EQ.SVAL(2))
     +                UMIN=AGMINI(SVAL(1),ZLOW,ZDRA,NVIZ,IIVZ,NEVZ,IIEZ)
      IF (UMAX.EQ.SVAL(1).OR.UMAX.EQ.SVAL(2))
     +                UMAX=AGMAXI(SVAL(1),ZHGH,ZDRA,NVIZ,IIVZ,NEVZ,IIEZ)
C
C Either or both values might still be null (if the user data was null).
C
      IF (UMIN.EQ.SVAL(1)) UMIN=UMAX
      IF (UMAX.EQ.SVAL(1)) UMAX=UMIN
C
C Check the relative values of UMIN and UMAX for problems.
C
      IF (ABS(UMIN-UMAX).LT.50.*SMRL*(ABS(UMIN)+ABS(UMAX))) GO TO 102
      IF (UMAX-UMIN) 101,102,103
  101 IF (ZMIN.NE.SVAL(1).AND.ZMIN.NE.SVAL(2)) UMAX=UMIN
      IF (ZMAX.NE.SVAL(1).AND.ZMAX.NE.SVAL(2)) UMIN=UMAX
C
  102 UMIN=UMIN-.5*ABS(UMIN)
      UMAX=UMAX+.5*ABS(UMAX)
      IF (UMIN.NE.UMAX) GO TO 103
      UMIN=-1.
      UMAX=+1.
C
C If the user wanted these values back-stored, do it.
C
  103 IF (ZMIN.EQ.SVAL(2)) ZMIN=UMIN
      IF (ZMAX.EQ.SVAL(2)) ZMAX=UMAX
C
C Done.
C
      RETURN
C
      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
!!# FUNCTION MODULE: <FUN_xyINTEGRALYY>
MODULE FUN_xyINTEGRALYY

!!## PURPOSE
!! Integrates the function $YY$ over a polygon in 2D.

!*******************************************************************************
!
!! xyINTEGRALYY_Pg integrates the function Y*Y over a polygon in 2D.
!
!  Discussion:
!
!    The polygon is bounded by the points (X(1:N), Y(1:N)).
!
!    INTEGRAL = (1/12) * sum ( 1 <= I <= N )
!      - ( Y(I)**3 + Y(I)**2 * Y(I-1) + Y(I) * Y(I-1)**2 + Y(I-1)**3 )
!      * ( X(I) - X(I-1) )
!
!    where X(0) and Y(0) should be replaced by X(N) and Y(N).
!
!  Modified:
!
!    10 July 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    SF Bockman,
!    Generalizing the Formula for Areas of Polygons to Moments,
!    American Mathematical Society Monthly,
!    1989, pages 131-132.
!
!  Parameters:
!
!    Input, integer N, the number of Pgertices of the polygon.
!    N should be at least 3 for a nonzero INTEGRALYY.
!
!    Input, real ( kind = 8 ) Pg(2,N), the coordinates of the Pgertices
!    of the polygon.  These Pgertices should be giPgen in
!    counter clockwise order.
!
!    Output, real ( kind = 8 ) INTEGRALYY, the Pgalue of the integral.
!

!!## EXTERNAL MODULES
USE FUN_Error                                   !!((04-A-FUN_Error.f90))
USE FUN_xyAVERAGEX                              !!((05-A-FUN_xyAVERAGEX.f90))
USE FUN_xyINTEGRAL1                             !!((06-A-FUN_xyINTEGRAL1.f90))

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## PROCEDURE OVERLOADING
INTERFACE xyINTEGRALYY_Pg
 MODULE PROCEDURE xyINTEGRALYY_Pg_Rsp
 MODULE PROCEDURE xyINTEGRALYY_Pg_Rdp
END INTERFACE

!!## PROCEDURE OVERLOADING
!INTERFACE xyINTEGRALYY_Ls
! MODULE PROCEDURE xyINTEGRALYY_Ls_Rsp
! MODULE PROCEDURE xyINTEGRALYY_Ls_Rdp
!END INTERFACE
!
!!!## PUBLIC ACCESS LIST
!PUBLIC :: xyINTEGRALYY_Ls
!PUBLIC :: xyINTEGRALYY_Ls_Rsp
!PUBLIC :: xyINTEGRALYY_Ls_Rdp

!!## PUBLIC ACCESS LIST
PUBLIC :: xyINTEGRALYY_Pg
PUBLIC :: xyINTEGRALYY_Pg_Rsp
PUBLIC :: xyINTEGRALYY_Pg_Rdp

!!## CONTAINED PROCEDURES
CONTAINS

!
!!!### FUNCTION <xyINTEGRALYY_Ls_Rsp>
!FUNCTION xyINTEGRALYY_Ls_Rsp( Ls ) RESULT(INTEGRALX)
!
!!!#### LOCAL MAPPINGS
!USE PAR_Constants_Rsp                           !!((02-A-PAR_Constants_Rsp.f90))
!USE KND_IntrinsicTypes,ONLY: KIND_R => KIND_Rsp !!((01-A-KND_IntrinsicTypes.f90))
!INCLUDE "06-A-FUN_xyINTEGRALYY_Ls.f90.hdr"
!!!--begin--
!INCLUDE "06-A-FUN_xyINTEGRALYY_Ls.f90.bdy"
!!!--end--
!END FUNCTION
!
!
!
!!!### FUNCTION <xyINTEGRALYY_Ls_Rdp>
!FUNCTION xyINTEGRALYY_Ls_Rdp( Ls ) RESULT(INTEGRALX)
!
!!!#### LOCAL MAPPINGS
!USE PAR_Constants_Rdp                           !!((02-A-PAR_Constants_Rsp.f90))
!USE KND_IntrinsicTypes,ONLY: KIND_R => KIND_Rdp !!((01-A-KND_IntrinsicTypes.f90))
!INCLUDE "06-A-FUN_xyINTEGRALYY_Ls.f90.hdr"
!!!--begin--
!INCLUDE "06-A-FUN_xyINTEGRALYY_Ls.f90.bdy"
!!!--end--
!END FUNCTION

!!### PURE FUNCTION <xyINTEGRALYY_Pg_Rsp>
PURE FUNCTION xyINTEGRALYY_Pg_Rsp( N , Pg ) RESULT(INTEGRALYY)

!!#### LOCAL MAPPINGS
USE PAR_Constants_Rsp                           !!((02-A-PAR_Constants_Rsp.f90))
USE KND_IntrinsicTypes,ONLY: KIND_R => KIND_Rsp !!((01-A-KND_IntrinsicTypes.f90))

INCLUDE "07-A-FUN_xyINTEGRALYY_Pg.f90.hdr"

!!--begin--
INCLUDE "07-A-FUN_xyINTEGRALYY_Pg.f90.bdy"
!!--end--

END FUNCTION

!!### PURE FUNCTION <xyINTEGRALYY_Pg_Rdp>
PURE FUNCTION xyINTEGRALYY_Pg_Rdp( N , Pg ) RESULT(INTEGRALYY)

!!#### LOCAL MAPPINGS
USE PAR_Constants_Rdp                           !!((02-A-PAR_Constants_Rdp.f90))
USE KND_IntrinsicTypes,ONLY: KIND_R => KIND_Rdp !!((01-A-KND_IntrinsicTypes.f90))

INCLUDE "07-A-FUN_xyINTEGRALYY_Pg.f90.hdr"
!!--begin--
INCLUDE "07-A-FUN_xyINTEGRALYY_Pg.f90.bdy"

!!--end--
END FUNCTION

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
According to http://www.sushimasters.com/topspotstoprestaurants.htm SushiMasters.com, Miyagi is one of the top 10 best Sushi places in the Sacramento region. Its located at the corner of Fair Oaks and Munroe, close to the Pavillion Shops of the Arden Arcade area of Sacramento.

Miyagi Bar and Sushi is traditional Japanese cuisine, with a large variety of sushi and sashimi.  Its slightly upscale, but keeps a relaxed friendly atmosphere.  The sushi chefs and owner/manager talk and joke with regulars, making it an enjoyable experience.  Upon being seated, youre given a bowl of edamame while you look over the menu.  Most of the rolls range from 12 to 16 dollars, and the dinner combination entrees are $16.  It has a full bar connected to the main dining area, and patio seating.

The SacBee had a http://www.sacbee.com/154/story/120114.html piece about a Battle of the Bartenders at the Sacramento Hard Rock Cafe, in February of 2007, where Miyagi won the award for the MostCreative Drink, its signature Lycheetini, which was also the crowd favorite.

  See also the wiki:sacramento:Miyagi Bar & Sushi SacWiki entry

20070217 00:03:14 nbsp My girlfriend and I were looking for a great sushi place in Sac, and Miyagi was http://www.yelp.com/search?find_locSacramento%2C+CA&radius10.0&sortbyrating&cfltjapanese top rated on Yelp and elsewhere.  This is an amazing place to take a date!  We went on Valentines Day, and one regular gave a Sushi Chef a pair of Valentines Day Boxers, which was hilarious.  The food was delicious, it was definitely the best sushi weve had in the area.  The fish was fresh, the unagi melted in my mouth, and there were interesting flavor profiles on some rolls.  And just as important, its finally a place with great presentation.  Our rolls were presented on the same platter/plate beautifully, and it was a fun experience.   Overall its pricey, especially with drinks or wine, but very reasonable and well worth it for what you receive.  Their drinks are high quality; the Lycheetini was delicious, and was made with Grey Goose.  Users/EdWins ES


text:
program fibo
  use iso_fortran_env, only : real64
  implicit none

  integer, parameter :: fib = 45

  write(*,*) fibonacci(fib)

contains

  function fibonacci(n) result(res)
    integer, intent(in) :: n
    integer :: res

    res = fibo_naive(n)
    !res = fibo_tail(0,1,n)

  end function fibonacci

  recursive function fibo_naive(n) result(res)
    integer, intent(in) :: n
    integer :: res

    if(n <= 0) then
       res = 0
       return
    end if

    if(n == 1) then
       res = 1
       return
    end if

    res = fibo_naive(n-1) + fibo_naive(n-2)

  end function fibo_naive

  recursive function fibo_tail(prev2,prev1,n) result(res)
    integer,intent(in) :: prev2,prev1,n
    integer :: res

    if(n == 0) then
       res = prev2
       return
    end if

    res = fibo_tail(prev1,prev1+prev2,n-1)

  end function fibo_tail

end program fibo

max_stars_repo_name: afanfa/OptimizedModernFortran

max_issues_repo_name: afanfa/OptimizedModernFortran

max_forks_repo_name: afanfa/OptimizedModernFortran


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     The subroutine opt_qmdff_ser optimizes the noncovalent parameters of a
!     given QMDFF. It sets up the parameters
!
!     The optimized parameters are written to [qmdffprefix].update
!
!
!     part of EVB
!

subroutine opt_qmdff_ser(filegeo,fileenergy,points,fffile1)
use evb_mod
use qmdff
use general
use lm_module
implicit none
logical::has_next
integer::nat,nat3,i,j,k,l,p,m
integer::inc
integer::dg_evb_mode,mat_size
integer::points,mp
integer::coord_mode ! for init_int
! random number generation:
INTEGER SEED,t
INTEGER :: i_seed
INTEGER, DIMENSION(:), ALLOCATABLE :: a_seed
INTEGER, DIMENSION(1:8) :: dt_seed
real(kind=8)::r
character(len=:),allocatable::char_num
character(len=30)::write_format
character(len=60)::filegeo,fileenergy
character(len=60)::fffile1
real(kind=8)::e1,e2,eref,V12,str_e1,str_e2
real(kind=8)::expo,d_p,root  ! short form for exp-part of the gaussians, dot-product
real(kind=8),allocatable::coord(:,:)
real(kind=8),dimension(3*natoms-6)::int_buf
real(kind=8),dimension(3*natoms)::g1,g2,gref
real(kind=8),dimension(3,natoms)::gqmdff1,gqmdff2
real(kind=8),dimension(points)::energies_ref,e_path,ff_e1,ff_e2,e_initial
real(kind=8),dimension(:,:),allocatable::path_xyz
real(kind=8),dimension(:),allocatable::int_path,int_coord  ! for internal coordinates
real(kind=8),dimension(:),allocatable::test_vec  ! TESTTEST
real(kind=8)::deldiscr,delsqrt,ediff,VdV,diffstep
integer::analytic  ! Test for EVB gradient calculation: Do it analytically
logical::bad_run    ! decides, if a levenberg marquardt run should be broken off
! parameters for newton/LM optimizations
real(kind=8)::chi,chi_old,derivative,y1,y2,num_step,deltaE,func
real(kind=8)::s_rest,step
integer::maxstep,o,worse,INFO
real(kind=8),dimension(3*natoms-6)::g_ff1,g_ff2  ! for EVB-QMDFF gradients
real(kind=8),dimension(:,:,:),allocatable::xyz_path  ! xyz structures of the reactionpath
! for multi start local search
real(kind=8)::glob_chi,glob_chi_new
real(kind=8)::coord_test(3,7)
real(kind=8)::dist,ang,oop,dihed
real(kind=8),dimension(:),allocatable::rmsd_path_int
real(kind=8)::minimum,temp
real(kind=8),allocatable::temp3(:)
integer,allocatable::temp2(:)
integer::location,start,ihalf
integer,dimension(:,:),allocatable::coord_tmp,coord_reshape
integer::msls_step,line,readstat
integer,dimension(:),allocatable::coord_types
character(len=40)::coord_line
real(kind=8)::e_one
integer::par_num  ! number of QMDFF modification parameters
real(kind=8),allocatable::par_vec(:)  ! vector of QMDFF modification parameters
real(kind=8),allocatable::par_opt(:)  ! vector with optimized QMDFF modification parameters
real(kind=8)::infinity,db1_prev_var  ! for check if any number is infinity
! for external Levenberg-Maquardt call
real(kind=8),allocatable::fvec(:),wa(:)
real(kind=8)::tol
integer::iflag,lwa
integer,allocatable::iwa(:)
character(len=60)::outfile_name
external::lm_qmdff ! the fitness function to be called in the Levenberg-Marquardt routine
integer::nat_mono  ! number of atoms of the monomer
integer::rank  ! MPI dummy parameter

rank=0
!
!     check, if any variable is infinity (define largest number)
!
infinity=HUGE(db1_prev_var)
!
!    Arrays for all needed input informations: energies (all_ens),
!    gradients (all_grads) and hessians (all_hess)
!    of the reference and the single QMDFFs (f1,f2)
!

!    short forms for variables
nat=natoms
nat3=3*nat
!
!    reduced dimension for internal variables
!
mp=dg_evb_points
allocate(xyz_path(3,nat,points))
allocate(path_xyz(nat*3,points))
allocate(all_ens(dg_evb_points))
allocate(all_xyz(nat*3,dg_evb_points))
allocate(coord(3,natoms))
allocate(fvec(points))
!
!   allocate multi start local search parameters
!    First, read in the newly witten dimer QMDFF
!
call prepare("dimer.qmdff","dummy","dummy",1)
natoms=n_one
nat_mono=natoms/2  ! determine atom number of the monomer
!
!     Open the trainset data and read them into arrays
!
open (unit=126,file=fileenergy,status="old")
open (unit=127,file=filegeo,status='old')

do i=1,points
   read(126,*,iostat=readstat) energies_ref(i)
   if (readstat .ne. 0) then
      write(*,*) "The reference-energies file contains too few lines!"
      call fatal
   end if
   call next_geo(xyz_path(:,:,i),natoms,127,has_next)
   if (.not.has_next) exit
end do
allocate(xyz_ref_lm(3,natoms,points))
xyz_ref_lm=xyz_path/bohr
close(126)
close(127)

!
!    Calculate the initial QMDFF energies for comparison
!
chi=0.d0
par_num=1+7*nat_mono
allocate(mn_par(1+2*(par_num-1)))
mn_par=1.d0
do i=1,points

   call ff_eg(n_one,at,xyz_ref_lm(:,:,i),e_one,g_one)
   call ff_nonb(n_one,at,xyz_ref_lm(:,:,i),q,r0ab,zab,r094_mod,sr42,c6xy,e_one,g_one)
   call ff_hb(n_one,at,xyz_ref_lm(:,:,i),e_one,g_one)
   e_initial(i)=e_one
end do
e_initial=e_initial-(minval(e_initial)-minval(energies_ref))
do i=1,points
   chi=chi+(e_initial(i)-energies_ref(i))*(e_initial(i)-energies_ref(i))
end do

write(*,*) "The intial fitness is:",chi
chi=0.d0
!
!    determine number of parameters to optimize:
!    1 Coulomb (global scaling)
!    3*nat_mono Van der Waals (three parameters per atom)
!    nat_mono Hydrogen bond (one parameter per atom)
!
allocate(par_vec(par_num))
allocate(par_opt(par_num))
!
!   set important parameters
!
Lambda=lm_par
num_step=0.00001d0
glob_chi=0
msls_step=1
write(*,*)  " . ...- -... --.- -- -.. ..-. ..-."
write(*,*) "noncovalent QMDFF parameters are optimized with Multi Start Local Search."
write(*,*) "We will start",maxstart,"single optimization runs."

write(15,*) "--- QMDFF NONCOVALENT PARAMETERIZATION ---"
write(15,*)
write(15,*) "Calculation initiated at: "
call timestamp ( )
write(15,*)
write(15,*) "You have requested a Multi Start local search optimization"
write(15,*) "for noncovalent QMDFF parameters."
write(15,*)

write(15,*)

! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
!   MULTI START LOCAL SEARCH: for all parameters individual values!
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
!   Write out global settings for coupling and optimization
!
write(*,*) "Multi Start Local Search algorithm:"
write(*,*) "---------------------------------------------------"

write(15,*) "Multi Start Local Search algorithm:"
write(15,*)
write(15,*) "SETTINGS:"
write(15,'(a,i5)') " * Number of atoms in the monomer: ", nat_mono
write(15,'(a,i5)') " * Number of optimizable parameters: ", par_num
write(15,'(a,i5)') " * Number of energy points in the trainset: ", points
write(15,'(a,i5)') " * Number of local searches ", maxstart
write(15,'(a,f10.5,f10.5)') " * Bounds for spawning of new params: ", lower_bond, upper_bond
write(15,*)
write(15,*) "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
!
!    define the variables in lm_module to use them for Levenberg-Marquardt
!
allocate(energies_ref_lm(points))
energies_ref_lm=energies_ref
!
!    define other Levenberg-Marquardt parameters
!
lwa=points*par_num+5*par_num+points
allocate(wa(lwa),iwa(par_num))
allocate(f_vec(points))
!
!    The tolerance value for optimizations
!
tol=0.00001D+00

do p=1,maxstart

!   Generate random start number
!
!   Generate random start numbers for all alphas
!   (taken from http://web.ph.surrey.ac.uk/fortweb/glossary/random_seed.html)
! ----- Set up random seed portably -----
   bad_run=.false.
   do i=1,par_num
      CALL RANDOM_SEED(size=i_seed)
      ALLOCATE(a_seed(1:i_seed))
      CALL RANDOM_SEED(get=a_seed)
      CALL DATE_AND_TIME(values=dt_seed)
      a_seed(i_seed)=dt_seed(8); a_seed(1)=dt_seed(8)*dt_seed(7)*dt_seed(6)
      CALL RANDOM_SEED(put=a_seed)
      DEALLOCATE(a_seed)
! ----- Done setting up random seed -----

      CALL RANDOM_NUMBER(r)
!
!    Interval, in which the number should be
!
      r=r*(upper_bond-lower_bond)+lower_bond
      par_vec(i)=r
   end do
!   write(*,*) "hhh alpha manipulated!"
!   alph=40.d0

   write(15,'(a,i4)') " Multi start local search round ",p
!
!    Call the Levenberg-Marquardt routines from MINPACK
!
   call r8vec_print (par_num, par_vec, '  Initial parameter values:' )
   iflag = 1
!   tol = 0.000001D+00
   write(15,*)
   write(15,*) "Do Levenberg-Marquardt optimization..."
!
!     Fitness function is submitted via lm_function
!
   call lmdif1 ( lm_qmdff, points, par_num, par_vec, f_vec, tol, &
         & info, iwa, wa, lwa )
   if (info.eq.0) then
      write(15,*) "FATAL ERROR: The input parameters for Levenberg-Marquardt are improper!"
      write(15,*) "Contact Julien Steffen to solve this problem!"
      call fatal
   else if ((info.eq.1) .or. (info.eq.2) .or. (info .eq.3)) then
      write(15,*) "The Levenberg-Marquardt run has converged!"
   else if (info.eq.4) then
      write(15,*) "ERROR: The alpha-vector is orthogonal to the Jacobian columns to machine"
      write(15,*) "precision in this round!  No convergence could be reached!"
      bad_run=.true.
   else if (info.eq.5) then
      write(15,*) "ERROR: Too many Levenberg-Marquardt runs were needed in this round!"
      write(15,*) "No convergence could be reached!"
      bad_run=.true.
   else if ((info.eq.6) .or. (info.eq.7)) then
      write(15,*) "Levenberg-Marquardt run has converged, but the errorsum was higher than"
      write(15,*) "tolerance! Maybe try a new training set?"
   end if
   call r8vec_print ( par_num, par_vec, '  Optimized parameter values:' )
!
!    Evaluate the results of local optimization
!
!    TEST 25.10.2017: add penalty if wiggles are present..
!
!   if (extrema .gt. 3) then
   chi=sum(f_vec)
!   end if
!
!    catch any exception in optimization runs!
!
   if (chi .ne. chi) then
      write(*,*) "ERROR: The value of the errorsum is NaN! Please check your input!"
      call fatal
   else if (chi .gt. infinity) then
      write(*,*) "ERROR: The value of the errorsum is infinity! Please check your input!"
      call fatal
   end if

   write(15,*) "Error sum for this round:",chi
   write(15,*)
   if (.not. bad_run) then
       if (glob_chi .eq. 0) then
          glob_chi=chi
          par_opt=par_vec
          write(*,*) "After Multi-Start-Local-Search step:",msls_step
          write(*,*) "Local optimization finished, new best fitness:",glob_chi
          write(15,*) ">> Local optimization finished!"
          write(15,*) ">> First best fitness value:",glob_chi
          write(15,*)
       else if (chi .le. glob_chi) then
          glob_chi=chi
          par_opt=par_vec
          write(*,*) "After Multi-Start-Local-Search step:",msls_step
          write(*,*) "Local optimization finished, new best fitness:",glob_chi
          write(15,*) ">> Local optimization finished!"
          write(15,*) ">> New best fitness value:",glob_chi
          write(15,*)
       else
!             write(*,*) "Local optimization finished, the best fitness is still:",glob_chi
          write(15,*) ">> Local optimization finished!"
          write(15,*) ">> No better value found..."
          write(15,*)
       end if
   else
       if (glob_chi .eq. 0) then
!
!      avoid endings of optimizations without optimized fitness value
!
          glob_chi=chi
          par_opt=par_vec
          write(*,*) "After Multi-Start-Local-Search step:",msls_step
          write(*,*) "Local optimization not converged, new best fitness:",glob_chi
          write(15,*) ">> Local optimization not converged..."
          write(15,*) ">> However, this is the first measured fitness value."
          write(15,*) ">> First best fitness value:",glob_chi
          write(15,*)
       else if (chi .le. glob_chi) then
          glob_chi=chi
          par_opt=par_vec
          write(*,*) "After Multi-Start-Local-Search step:",msls_step
          write(*,*) "Local optimization not converged, new best fitness:",glob_chi
          write(15,*) ">> Local optimization not converged..."
          write(15,*) ">> New best fitness value:",glob_chi
          write(15,*)
      end if
   end if
   write(15,*) "-----------------------------------------------------------------"
   msls_step=msls_step+1
end do
write(*,*)  ".. .----. -- -.. --- -. . "
write(*,*) "Multi start local optimization finished!"
write(*,*) "The best fitness-value is:", glob_chi
write(*,*) "Look into evbopt.log for further details."

write(15,*) "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
write(15,*)
write(15,*)  ".. .----. -- -.. --- -. . "
write(15,*)
write(15,*) ">>> Multi start local optimization finished!"
write(15,*) ">>> The best fitness-value is:", glob_chi
write(15,*) ">>> The optimized alpha values are:"
do i=1,par_num
   write ( 15, '(2x,i8,2x,g16.8)' ) i, par_opt(i)
end do
write(15,*)
write(15,*) ">>> File evb_pars.dat was written"
write(15,*) ">>> This file is needed for specification of the coupling term!"
write(15,*)
!
!     At least calculate the path energies for the optimized parameters
!     and plot them with the inital energies for comparison
!
if (maxstart .gt. 1) then
   par_vec=par_opt
end if

par_opt=par_vec
chi=0.d0
open(unit=31,file="en_before_after.dat",status="unknown")
write(31,*) "# Reference energy    QMDFF energy (before)    QMDFF energy (after) "
mn_par(1)=par_vec(1)
mn_par(2:par_num)=par_vec(2:par_num)
mn_par(par_num+1:2*par_num)=par_vec(2:par_num)

do i=1,points

   call ff_eg(n_one,at,xyz_ref_lm(:,:,i),e_one,g_one)
   call ff_nonb(n_one,at,xyz_ref_lm(:,:,i),q,r0ab,zab,r094_mod,sr42,c6xy,e_one,g_one)
   call ff_hb(n_one,at,xyz_ref_lm(:,:,i),e_one,g_one)
   e_path(i)=e_one
end do
e_path=e_path-(minval(e_path)-minval(energies_ref_lm))
do i=1,points
   chi=chi+(e_path(i)-energies_ref_lm(i))*(e_path(i)-energies_ref_lm(i))
   write(31,*) energies_ref(i),e_initial(i),e_path(i)
end do
write(*,*)
write(*,*) "Energies of the reactionpath are written to energies.qmdff,"
write(*,*) "gradients are written to gradients.qmdff"
close(31)
!
!   write the par_opt reference file for further calculations
!
outfile_name=fffile1(1:len(trim(fffile1))-6) // "_mod.dat"
open(unit=67,file=outfile_name,status="unknown")
write(67,*) "** This file contains correction parameters for QMDFF nonbonded interactions **"
write(67,'(f14.8)') par_opt(1)
do i=1,natoms/2
   write(67,'(7f14.8)') par_opt(2+(i-1)*7:1+i*7)
end do
close(67)
write(*,*)
write(*,*) trim(outfile_name)," file written. Use this file for further calculations."
write(*,*) "Exiting normally..."
write(*,*)

close(134)

!
!   deallocate arrays
!

return
end subroutine opt_qmdff_ser

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
!==============================================================================!
  subroutine Probe_1d_Cells_Nodes(grid)
!------------------------------------------------------------------------------!
!   This subroutine finds the coordinate of cell's nodes in non-homogeneous    !
!   direction and write them in file called "name.1d"                          !
!------------------------------------------------------------------------------!
  use Name_Mod, only: problem_name
  use Grid_Mod
  use Sort_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!----------------------------------[Calling]-----------------------------------!
  include "Approx_Real.int"
!-----------------------------------[Locals]-----------------------------------!
  integer           :: n_prob, p, c, n
  real              :: zp(16384)
  character(len=80) :: name_prob
  character(len=80) :: answer
  logical           :: isit
!==============================================================================!

  print *, '#==========================================='
  print *, '# Creating 1d file with the node '
  print *, '# coordinates in non-homogeneous directions '
  print *, '#-------------------------------------------'
  print *, '# Insert non-homogeneous direction '
  print *, '# (x, y, z, rx, ry, rz or skip)'
  print *, '# -------------------------------------------'
  read(*,*) answer
  call To_Upper_Case(answer)
  if(answer .eq. 'SKIP') return

  n_prob = 0
  zp(:)  = 0.0

  !-----------------------------!
  !   Browse through all cells  !
  !-----------------------------!
  do c = -grid % n_bnd_cells, grid % n_cells
    do n = 1, grid % cells_n_nodes(c)

      ! Try to find the cell among the probes
      do p=1, n_prob
        if(answer .eq. 'X') then
          if( Approx_Real(grid % xn(grid % cells_n(n,c)), zp(p)) ) go to 1
        else if(answer .eq. 'Y') then
          if( Approx_Real(grid % yn(grid % cells_n(n,c)), zp(p)) ) go to 1
        else if(answer .eq. 'Z') then
          if( Approx_Real(grid % zn(grid % cells_n(n,c)), zp(p)) ) go to 1
        else if(answer .eq. 'RX') then
          if( Approx_Real( sqrt(grid % zn(grid % cells_n(n,c))**2 +           &
                                grid % yn(grid % cells_n(n,c))**2), zp(p)) )  &
            go to 1
        else if(answer .eq. 'RY') then
          if( Approx_Real( sqrt(grid % xn(grid % cells_n(n,c))**2 +           &
                                grid % zn(grid % cells_n(n,c))**2), zp(p)) )  &
            go to 1
        else if(answer .eq. 'RZ') then
          if( Approx_Real( sqrt(grid % xn(grid % cells_n(n,c))**2 +           &
                                grid % yn(grid % cells_n(n,c))**2), zp(p)) )  &
            go to 1
        end if
      end do

      ! Couldn't find a cell among the probes, add a new one
      n_prob = n_prob + 1
      if(answer .eq. 'X') zp(n_prob) = grid % xn(grid % cells_n(n,c))
      if(answer .eq. 'Y') zp(n_prob) = grid % yn(grid % cells_n(n,c))
      if(answer .eq. 'Z') zp(n_prob) = grid % zn(grid % cells_n(n,c))

      if(answer .eq. 'RX') zp(n_prob) =                                  &
                           sqrt(grid % zn(grid % cells_n(n,c))**2 +      &
                                grid % yn(grid % cells_n(n,c))**2)
      if(answer .eq. 'RY') zp(n_prob) =                                  &
                           sqrt(grid % xn(grid % cells_n(n,c))**2 +      &
                                grid % zn(grid % cells_n(n,c))**2)
      if(answer .eq. 'RZ') zp(n_prob) =                                  &
                           sqrt(grid % xn(grid % cells_n(n,c))**2 +      &
                                grid % yn(grid % cells_n(n,c))**2)

      if(n_prob .eq. 16384) then
        print *, '# Probe 1d: Not a 1d (channel flow) problem.'
        isit = .false.
        return
      end if
    end do
1 end do

  isit = .true.

  !--------------------!
  !   Create 1d file   !
  !--------------------!
  name_prob = problem_name
  name_prob(len_trim(problem_name)+1:len_trim(problem_name)+3) = '.1d'
  print *, '# Creating the file: ', trim(name_prob)
  open(9, file=name_prob)

  ! Write the number of probes
  write(9,'(i8)') n_prob

  call Sort_Mod_Real(zp(1:n_prob))

  ! Write the probe coordinates out
  do p=1, n_prob
    write(9,'(i8,1e17.8)') p, zp(p)
  end do

  close(9)

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
!{\src2tex{textfont=tt}}
!!****f* ABINIT/sg_fftpx
!! NAME
!! sg_fftpx
!!
!! FUNCTION
!! This subroutine is called by the 3-dimensional fft to conduct the
!! "x" transforms for all y and z.
!! Accomodate more optimal treatment of
!! zero padding following the method of fft3dp.
!!
!! COPYRIGHT
!! Copyright by Stefan Goedecker, Ithaca, NY USA, July 14, 1993
!! Copyright (C) 1998-2001 ABINIT group (DCA, XG)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~ABINIT/Infos/copyright
!! or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  fftcache=size of the cache (kB)
!!  mfac = maximum number of factors in 1D FFTs
!!  mg = maximum length of 1D FFTs
!!  mgbound = effective maximum length of 1D FFTs, for dimensioning gbound
!!  nd1=first dimension of (complex) arrays z and zbr (treated as real within
!!   this subroutine)
!!  nd2=second dimension of (complex) arrays z and zbr (treated as real within
!!   this subroutine)
!!  nd3=third dimension of (complex) arrays z and zbr (treated as real within
!!   this subroutine)
!!  n2,n3=actual length of y and z transforms
!!  z(2,nd1,nd2,nd3)=INPUT array; destroyed by transformation
!!  trig, aft, now, bef, ind=provided by previous call to ctrig
!!   Note that in this routine (and in ctrig) the values in array trig are
!!   actually cos and tan, not cos and sin.  Use of tan allows advantageous
!!   use of FMA on the ibm rs6000.
!!  ic=number of (radix) factors of x transform length (from ctrig)
!!  ris=sign of exponential in transform (should be 1 or -1; real)
!!
!! OUTPUT
!!  zbr(2,nd1,nd2,nd3)=OUTPUT transformed array; no scaling applied
!!
!! SIDE EFFECTS
!!
!! NOTES
!! This routine blocks the x transforms
!! so that all transforms under consideration at one step fit within
!! the cache memory, which is crucial for optimal performance.
!! The blocking factor is set by parameter "fftcache" below, which should
!! be adjusted to be somewhat smaller (say 3/4) than the actual cache size
!! of the machine.
!!
!!
!! TODO
!! Use latex for the equation above
!!
!! PARENTS
!!      sg_fftpad,sg_fftrisc
!!
!! SOURCE

 subroutine sg_fftpx(fftcache,mfac,mg,mgbound,nd1,nd2,nd3,n2,n3,&
&    z,zbr,trig,aft,now,bef,ris,ind,ic,gbound)

 use basis_defs
 implicit none

!Arguments ------------------------------------
 integer :: fftcache,ic,mfac,mg,mgbound,nd1,nd2,nd3,n1,n2,n3
!Dimensions of aft, now, bef, ind, and trig should agree with
!those in subroutine ctrig.
 integer :: aft(mfac),bef(mfac),gbound(2*mgbound+4),ind(mg),now(mfac)
 real(dp) :: ris
 real(dp) :: trig(2,mg),z(2,nd1,nd2,nd3),zbr(2,nd1,nd2,nd3)

!Local variables-------------------------------
 integer :: g2,g2min,g2max,gg3,g3,&
& g3max,g3min,i,ia,ib,igb,ihalfy,indx,i3,j,jj,len3,lot,&
& lowlim,ma,mb,ntb,upplim
 real(dp),parameter :: &
& cos2=0.3090169943749474d0,&   !cos(2.d0*pi/5.d0)
& cos4=-0.8090169943749474d0,&  !cos(4.d0*pi/5.d0)
& sin42=0.6180339887498948d0    !sin(4.d0*pi/5.d0)/sin(2.d0*pi/5.d0)
 real(dp) :: bb,cr2,cr2s,cr3,cr3p,cr4,cr5,ct2,ct3,ct4,ct5,&
& factor,r,r1,r2,r25,r3,r34,r4,r5,s,sin2,s1,s2,s25,s3,s34,s4,s5
 character message*500

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

 g3min=gbound(1)
 g3max=gbound(2)
 igb=3
 len3=g3max-g3min+1

!Do x transforms in blocks of size "lot" which is set by how
!many x transform arrays (of size nd1 each) fit into the nominal
!cache size "fftcache".
!Loop over blocks in the loop below.

 factor=0.75d0
 lot=(fftcache*factor*1000d0)/(nd1*8*2)

!Express loop over y, z in terms of separate z and y loops

!$OMP PARALLEL DO DEFAULT(PRIVATE)&
!$OMP&SHARED(aft,bef,gbound,g3max,ic,ind,len3,lot)&
!$OMP&SHARED(n2,n3,nd2,now,ris,trig,z,zbr)
 do gg3=1,len3

  if (gg3<=g3max+1) then
   g3=gg3
  else
!  wrap around for negative gg3
   g3=gg3-len3+n3
  endif

  igb=gg3*2+1
  g2min=gbound(igb)
  g2max=gbound(igb+1)

! Split the y loop into positive and wrapped-around negative parts

  do ihalfy=1,2

!  Start at 1 for ihalfy=1; g2min+1+n2 for ihalfy=2
   lowlim=1+(ihalfy-1)*(g2min+n2)
!  End at g2max+1 for ihalfy=1; n2 for ihalfy=2
   upplim=g2max+1+(ihalfy-1)*(n2-g2max-1)

   do g2=lowlim,upplim,lot

!   Find array starting address ma and ending address mb
!   modified xg 980107
!   ma=g2+(g3-1)*nd2
    ma=g2
!   Perform "lot" transforms at a time (until out of data)
!   mb=min(g2+(lot-1),upplim)+(g3-1)*nd2
    mb=min(g2+(lot-1),upplim)

!   -------------------------------------------------------------------------
!
!   Direct transformation

!   Run over all factors except the last (to ic-1), performing
!   x transform

!   Note: fortran should skip this loop if ic=1; beware "onetrip"
!   compiler option which forces each loop at least once

    do i=1,ic-1
     ntb=now(i)*bef(i)

!    Treat radix 4
     if (now(i)==4) then
      ia=0

!     First step of factor 4
      do ib=1,bef(i)
       do j=ma,mb
        r4=z(1,ia*ntb+3*bef(i)+ib,j,g3)
        s4=z(2,ia*ntb+3*bef(i)+ib,j,g3)
        r3=z(1,ia*ntb+2*bef(i)+ib,j,g3)
        s3=z(2,ia*ntb+2*bef(i)+ib,j,g3)
        r2=z(1,ia*ntb+bef(i)+ib,j,g3)
        s2=z(2,ia*ntb+bef(i)+ib,j,g3)
        r1=z(1,ia*ntb+ib,j,g3)
        s1=z(2,ia*ntb+ib,j,g3)

        r=r1 + r3
        s=r2 + r4
        z(1,ia*ntb+ib,j,g3) = r + s
        z(1,ia*ntb+2*bef(i)+ib,j,g3) = r - s
        r=r1 - r3
        s=s2 - s4
        z(1,ia*ntb+bef(i)+ib,j,g3) = r - s*ris
        z(1,ia*ntb+3*bef(i)+ib,j,g3) = r + s*ris
        r=s1 + s3
        s=s2 + s4
        z(2,ia*ntb+ib,j,g3) = r + s
        z(2,ia*ntb+2*bef(i)+ib,j,g3) = r - s
        r=s1 - s3
        s=r2 - r4
        z(2,ia*ntb+bef(i)+ib,j,g3) = r + s*ris
        z(2,ia*ntb+3*bef(i)+ib,j,g3) = r - s*ris
       enddo
      enddo

!     Second step of factor 4
      do ia=1,aft(i)-1
       indx=ind(ia*4*bef(i)+1)-1
       indx=indx*bef(i)
       cr2=trig(1,indx)
       ct2=trig(2,indx)
       cr3=trig(1,2*indx)
       ct3=trig(2,2*indx)
       cr4=trig(1,3*indx)
       ct4=trig(2,3*indx)
       cr4=cr4/cr2
       cr2s=cr2*ris
       do ib=1,bef(i)
        do j=ma,mb
         r4=z(1,ia*ntb+3*bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+3*bef(i)+ib,j,g3)*ct4
         s4=z(1,ia*ntb+3*bef(i)+ib,j,g3)*ct4 + &
&           z(2,ia*ntb+3*bef(i)+ib,j,g3)
         r3=z(1,ia*ntb+2*bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+2*bef(i)+ib,j,g3)*ct3
         s3=z(1,ia*ntb+2*bef(i)+ib,j,g3)*ct3 + &
&           z(2,ia*ntb+2*bef(i)+ib,j,g3)
         r2=z(1,ia*ntb+bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+bef(i)+ib,j,g3)*ct2
         s2=z(1,ia*ntb+bef(i)+ib,j,g3)*ct2 + &
&           z(2,ia*ntb+bef(i)+ib,j,g3)
         r1=z(1,ia*ntb+ib,j,g3)
         s1=z(2,ia*ntb+ib,j,g3)

         r=r1 + r3*cr3
         s=r2 + r4*cr4
         z(1,ia*ntb+ib,j,g3) = r + s*cr2
         z(1,ia*ntb+2*bef(i)+ib,j,g3) = r - s*cr2
         r=r1 - r3*cr3
         s=s2 - s4*cr4
         z(1,ia*ntb+bef(i)+ib,j,g3) = r - s*cr2s
         z(1,ia*ntb+3*bef(i)+ib,j,g3) = r + s*cr2s
         r=s1 + s3*cr3
         s=s2 + s4*cr4
         z(2,ia*ntb+ib,j,g3) = r + s*cr2
         z(2,ia*ntb+2*bef(i)+ib,j,g3) = r - s*cr2
         r=s1 - s3*cr3
         s=r2 - r4*cr4
         z(2,ia*ntb+bef(i)+ib,j,g3) = r + s*cr2s
         z(2,ia*ntb+3*bef(i)+ib,j,g3) = r - s*cr2s
        enddo
       enddo
      enddo

!    Treat radix 2
     else if (now(i)==2) then
      ia=0

!     First step of factor 2
      do ib=1,bef(i)
       do j=ma,mb
        r1=z(1,ia*ntb+ib,j,g3)
        s1=z(2,ia*ntb+ib,j,g3)
        r2=z(1,ia*ntb+bef(i)+ib,j,g3)
        s2=z(2,ia*ntb+bef(i)+ib,j,g3)
        z(1,ia*ntb+ib,j,g3) =  r2 + r1
        z(2,ia*ntb+ib,j,g3) =  s2 + s1
        z(1,ia*ntb+bef(i)+ib,j,g3) = -r2 + r1
        z(2,ia*ntb+bef(i)+ib,j,g3) = -s2 + s1
       enddo
      enddo

!     Second step of radix 2
      do ia=1,aft(i)-1
       indx=ind(ia*2*bef(i)+1)-1
       indx=indx*bef(i)
       cr2=trig(1,indx)
       ct2=trig(2,indx)
       do ib=1,bef(i)
        do j=ma,mb
         r1=z(1,ia*ntb+ib,j,g3)
         s1=z(2,ia*ntb+ib,j,g3)
         r2=z(1,ia*ntb+bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+bef(i)+ib,j,g3)*ct2
         s2=z(1,ia*ntb+bef(i)+ib,j,g3)*ct2 + &
&           z(2,ia*ntb+bef(i)+ib,j,g3)
         z(1,ia*ntb+ib,j,g3) =  r2*cr2 + r1
         z(2,ia*ntb+ib,j,g3) =  s2*cr2 + s1
         z(1,ia*ntb+bef(i)+ib,j,g3) = -r2*cr2 + r1
         z(2,ia*ntb+bef(i)+ib,j,g3) = -s2*cr2 + s1
        enddo
       enddo
      enddo

!    Treat radix 3
     else if (now(i)==3) then
!     .5d0*sqrt(3.d0)=0.8660254037844387d0
      ia=0
      bb=ris*0.8660254037844387d0

!     First step of radix 3
      do ib=1,bef(i)
       do j=ma,mb
        r1=z(1,ia*ntb+ib,j,g3)
        s1=z(2,ia*ntb+ib,j,g3)
        r2=z(1,ia*ntb+bef(i)+ib,j,g3)
        s2=z(2,ia*ntb+bef(i)+ib,j,g3)
        r3=z(1,ia*ntb+2*bef(i)+ib,j,g3)
        s3=z(2,ia*ntb+2*bef(i)+ib,j,g3)
        r=r2 + r3
        s=s2 + s3
        z(1,ia*ntb+ib,j,g3) = r + r1
        z(2,ia*ntb+ib,j,g3) = s + s1
        r1=r1 - r*.5d0
        s1=s1 - s*.5d0
        r2=r2-r3
        s2=s2-s3
        z(1,ia*ntb+bef(i)+ib,j,g3) = r1 - s2*bb
        z(2,ia*ntb+bef(i)+ib,j,g3) = s1 + r2*bb
        z(1,ia*ntb+2*bef(i)+ib,j,g3) = r1 + s2*bb
        z(2,ia*ntb+2*bef(i)+ib,j,g3) = s1 - r2*bb
       enddo
      enddo

!     Second step of radix 3
      do ia=1,aft(i)-1
       indx=ind(ia*3*bef(i)+1)-1
       indx=indx*bef(i)
       cr2=trig(1,indx)
       ct2=trig(2,indx)
       cr3=trig(1,2*indx)
       ct3=trig(2,2*indx)
       cr2=cr2/cr3
       cr3p=.5d0*cr3
       bb=ris*cr3*0.8660254037844387d0
       do ib=1,bef(i)
        do j=ma,mb
         r1=z(1,ia*ntb+ib,j,g3)
         s1=z(2,ia*ntb+ib,j,g3)
         r2=z(1,ia*ntb+bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+bef(i)+ib,j,g3)*ct2
         s2=z(1,ia*ntb+bef(i)+ib,j,g3)*ct2 + &
&           z(2,ia*ntb+bef(i)+ib,j,g3)
         r3=z(1,ia*ntb+2*bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+2*bef(i)+ib,j,g3)*ct3
         s3=z(1,ia*ntb+2*bef(i)+ib,j,g3)*ct3 + &
&           z(2,ia*ntb+2*bef(i)+ib,j,g3)
         r=cr2*r2 + r3
         s=cr2*s2 + s3
         z(1,ia*ntb+ib,j,g3) = r*cr3 + r1
         z(2,ia*ntb+ib,j,g3) = s*cr3 + s1
         r1=r1 - r*cr3p
         s1=s1 - s*cr3p
         r2=cr2*r2-r3
         s2=cr2*s2-s3
         z(1,ia*ntb+bef(i)+ib,j,g3) = r1 - s2*bb
         z(2,ia*ntb+bef(i)+ib,j,g3) = s1 + r2*bb
         z(1,ia*ntb+2*bef(i)+ib,j,g3) = r1 + s2*bb
         z(2,ia*ntb+2*bef(i)+ib,j,g3) = s1 - r2*bb
        enddo
       enddo
      enddo

!    Treat radix 5
     else if (now(i)==5) then
!     sin(2.d0*pi/5.d0)
      sin2=ris*0.9510565162951536d0
      ia=0

!     First step of radix 5
      do ib=1,bef(i)
       do j=ma,mb
        r1=z(1,ia*ntb+ib,j,g3)
        s1=z(2,ia*ntb+ib,j,g3)
        r2=z(1,ia*ntb+bef(i)+ib,j,g3)
        s2=z(2,ia*ntb+bef(i)+ib,j,g3)
        r3=z(1,ia*ntb+2*bef(i)+ib,j,g3)
        s3=z(2,ia*ntb+2*bef(i)+ib,j,g3)
        r4=z(1,ia*ntb+3*bef(i)+ib,j,g3)
        s4=z(2,ia*ntb+3*bef(i)+ib,j,g3)
        r5=z(1,ia*ntb+4*bef(i)+ib,j,g3)
        s5=z(2,ia*ntb+4*bef(i)+ib,j,g3)
        r25 = r2 + r5
        r34 = r3 + r4
        s25 = s2 - s5
        s34 = s3 - s4
        z(1,ia*ntb+ib,j,g3) = r1 + r25 + r34
        r = r1 + cos2*r25 + cos4*r34
        s = s25 + sin42*s34
        z(1,ia*ntb+bef(i)+ib,j,g3) = r - sin2*s
        z(1,ia*ntb+4*bef(i)+ib,j,g3) = r + sin2*s
        r = r1 + cos4*r25 + cos2*r34
        s = sin42*s25 - s34
        z(1,ia*ntb+2*bef(i)+ib,j,g3) = r - sin2*s
        z(1,ia*ntb+3*bef(i)+ib,j,g3) = r + sin2*s
        r25 = r2 - r5
        r34 = r3 - r4
        s25 = s2 + s5
        s34 = s3 + s4
        z(2,ia*ntb+ib,j,g3) = s1 + s25 + s34
        r = s1 + cos2*s25 + cos4*s34
        s = r25 + sin42*r34
        z(2,ia*ntb+bef(i)+ib,j,g3) = r + sin2*s
        z(2,ia*ntb+4*bef(i)+ib,j,g3) = r - sin2*s
        r = s1 + cos4*s25 + cos2*s34
        s = sin42*r25 - r34
        z(2,ia*ntb+2*bef(i)+ib,j,g3) = r + sin2*s
        z(2,ia*ntb+3*bef(i)+ib,j,g3) = r - sin2*s
       enddo
      enddo

!     Second step of radix 5
      do ia=1,aft(i)-1
       indx=ind(ia*5*bef(i)+1)-1
       indx=indx*bef(i)
       cr2=trig(1,indx)
       ct2=trig(2,indx)
       cr3=trig(1,2*indx)
       ct3=trig(2,2*indx)
       cr4=trig(1,3*indx)
       ct4=trig(2,3*indx)
       cr5=trig(1,4*indx)
       ct5=trig(2,4*indx)
       do ib=1,bef(i)
        do j=ma,mb
         r1=z(1,ia*ntb+ib,j,g3)
         s1=z(2,ia*ntb+ib,j,g3)
         r2=cr2*(z(1,ia*ntb+bef(i)+ib,j,g3) - &
&                z(2,ia*ntb+bef(i)+ib,j,g3)*ct2)
         s2=cr2*(z(1,ia*ntb+bef(i)+ib,j,g3)*ct2 + &
&                z(2,ia*ntb+bef(i)+ib,j,g3))
         r3=cr3*(z(1,ia*ntb+2*bef(i)+ib,j,g3) - &
&                z(2,ia*ntb+2*bef(i)+ib,j,g3)*ct3)
         s3=cr3*(z(1,ia*ntb+2*bef(i)+ib,j,g3)*ct3 + &
&                z(2,ia*ntb+2*bef(i)+ib,j,g3))
         r4=z(1,ia*ntb+3*bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+3*bef(i)+ib,j,g3)*ct4
         s4=z(1,ia*ntb+3*bef(i)+ib,j,g3)*ct4 + &
&           z(2,ia*ntb+3*bef(i)+ib,j,g3)
         r5=z(1,ia*ntb+4*bef(i)+ib,j,g3) - &
&           z(2,ia*ntb+4*bef(i)+ib,j,g3)*ct5
         s5=z(1,ia*ntb+4*bef(i)+ib,j,g3)*ct5 + &
&           z(2,ia*ntb+4*bef(i)+ib,j,g3)
         r25 = r2 + r5*cr5
         r34 = r3 + r4*cr4
         s25 = s2 - s5*cr5
         s34 = s3 - s4*cr4
         z(1,ia*ntb+ib,j,g3) = r1 + r25 + r34
         r = r1 + cos2*r25 + cos4*r34
         s = s25 + sin42*s34
         z(1,ia*ntb+bef(i)+ib,j,g3) = r - sin2*s
         z(1,ia*ntb+4*bef(i)+ib,j,g3) = r + sin2*s
         r = r1 + cos4*r25 + cos2*r34
         s = sin42*s25 - s34
         z(1,ia*ntb+2*bef(i)+ib,j,g3) = r - sin2*s
         z(1,ia*ntb+3*bef(i)+ib,j,g3) = r + sin2*s
         r25 = r2 - r5*cr5
         r34 = r3 - r4*cr4
         s25 = s2 + s5*cr5
         s34 = s3 + s4*cr4
         z(2,ia*ntb+ib,j,g3) = s1 + s25 + s34
         r = s1 + cos2*s25 + cos4*s34
         s = r25 + sin42*r34
         z(2,ia*ntb+bef(i)+ib,j,g3) = r + sin2*s
         z(2,ia*ntb+4*bef(i)+ib,j,g3) = r - sin2*s
         r = s1 + cos4*s25 + cos2*s34
         s = sin42*r25 - r34
         z(2,ia*ntb+2*bef(i)+ib,j,g3) = r + sin2*s
         z(2,ia*ntb+3*bef(i)+ib,j,g3) = r - sin2*s
        enddo
       enddo
      enddo

     else

!     All radices treated
      write(message, '(a,a,a,a)' )ch10,&
&      ' sg_fftpx : BUG -',ch10,&
&      '  called with factors other than 2, 3, and 5'
!      call wrtout(06,message,'PERS')
!      call leave_new('PERS')
     endif

    enddo  ! End of direct transformation (loop over ic)

!   -----------------------------------------------------------------

!   Bitreversal
!   Perform bit reversal on last factor of transformation

!   Treat radix 4
    if (now(ic)==4) then
     ia=0

!    First step of radix 4
     do j=ma,mb
      r4=z(1,ia*4+4,j,g3)
      s4=z(2,ia*4+4,j,g3)
      r3=z(1,ia*4+3,j,g3)
      s3=z(2,ia*4+3,j,g3)
      r2=z(1,ia*4+2,j,g3)
      s2=z(2,ia*4+2,j,g3)
      r1=z(1,ia*4+1,j,g3)
      s1=z(2,ia*4+1,j,g3)

      r=r1 + r3
      s=r2 + r4
      zbr(1,ind(ia*4+1),j,g3) = r + s
      zbr(1,ind(ia*4+3),j,g3) = r - s
      r=r1 - r3
      s=s2 - s4
      zbr(1,ind(ia*4+2),j,g3) = r - s*ris
      zbr(1,ind(ia*4+4),j,g3) = r + s*ris
      r=s1 + s3
      s=s2 + s4
      zbr(2,ind(ia*4+1),j,g3) = r + s
      zbr(2,ind(ia*4+3),j,g3) = r - s
      r=s1 - s3
      s=r2 - r4
      zbr(2,ind(ia*4+2),j,g3) = r + s*ris
      zbr(2,ind(ia*4+4),j,g3) = r - s*ris
     enddo

!    Second step of radix 4
     do ia=1,aft(ic)-1
      indx=ind(ia*4+1)-1
      cr2=trig(1,indx)
      ct2=trig(2,indx)
      cr3=trig(1,2*indx)
      ct3=trig(2,2*indx)
      cr4=trig(1,3*indx)
      ct4=trig(2,3*indx)
      cr4=cr4/cr2
      cr2s=cr2*ris
      do j=ma,mb
       r4=z(1,ia*4+4,j,g3) - z(2,ia*4+4,j,g3)*ct4
       s4=z(1,ia*4+4,j,g3)*ct4 + z(2,ia*4+4,j,g3)
       r3=z(1,ia*4+3,j,g3) - z(2,ia*4+3,j,g3)*ct3
       s3=z(1,ia*4+3,j,g3)*ct3 + z(2,ia*4+3,j,g3)
       r2=z(1,ia*4+2,j,g3) - z(2,ia*4+2,j,g3)*ct2
       s2=z(1,ia*4+2,j,g3)*ct2 + z(2,ia*4+2,j,g3)
       r1=z(1,ia*4+1,j,g3)
       s1=z(2,ia*4+1,j,g3)

       r=r1 + r3*cr3
       s=r2 + r4*cr4
       zbr(1,ind(ia*4+1),j,g3) = r + s*cr2
       zbr(1,ind(ia*4+3),j,g3) = r - s*cr2
       r=r1 - r3*cr3
       s=s2 - s4*cr4
       zbr(1,ind(ia*4+2),j,g3) = r - s*cr2s
       zbr(1,ind(ia*4+4),j,g3) = r + s*cr2s
       r=s1 + s3*cr3
       s=s2 + s4*cr4
       zbr(2,ind(ia*4+1),j,g3) = r + s*cr2
       zbr(2,ind(ia*4+3),j,g3) = r - s*cr2
       r=s1 - s3*cr3
       s=r2 - r4*cr4
       zbr(2,ind(ia*4+2),j,g3) = r + s*cr2s
       zbr(2,ind(ia*4+4),j,g3) = r - s*cr2s
      enddo
     enddo

!   Treat radix 2
    else if (now(ic)==2) then

     ia=0

!    First step of radix 2
     do j=ma,mb
      r1=z(1,ia*2+1,j,g3)
      s1=z(2,ia*2+1,j,g3)
      r2=z(1,ia*2+2,j,g3)
      s2=z(2,ia*2+2,j,g3)
      zbr(1,ind(ia*2+1),j,g3) =  r2 + r1
      zbr(2,ind(ia*2+1),j,g3) =  s2 + s1
      zbr(1,ind(ia*2+2),j,g3) = -r2 + r1
      zbr(2,ind(ia*2+2),j,g3) = -s2 + s1
     enddo

!    Second step of radix 2
     do ia=1,aft(ic)-1
      indx=ind(ia*2+1)-1
      cr2=trig(1,indx)
      ct2=trig(2,indx)
      do j=ma,mb
       r1=z(1,ia*2+1,j,g3)
       s1=z(2,ia*2+1,j,g3)
       r2=z(1,ia*2+2,j,g3) - z(2,ia*2+2,j,g3)*ct2
       s2=z(1,ia*2+2,j,g3)*ct2 + z(2,ia*2+2,j,g3)
       zbr(1,ind(ia*2+1),j,g3) =  r2*cr2 + r1
       zbr(2,ind(ia*2+1),j,g3) =  s2*cr2 + s1
       zbr(1,ind(ia*2+2),j,g3) = -r2*cr2 + r1
       zbr(2,ind(ia*2+2),j,g3) = -s2*cr2 + s1
      enddo
     enddo

!   Treat radix 3
    else if (now(ic)==3) then
!    radix 3
!    .5d0*sqrt(3.d0)=0.8660254037844387d0
     ia=0
     bb=ris*0.8660254037844387d0

!    First step of radix 3
     do j=ma,mb
      r1=z(1,ia*3+1,j,g3)
      s1=z(2,ia*3+1,j,g3)
      r2=z(1,ia*3+2,j,g3)
      s2=z(2,ia*3+2,j,g3)
      r3=z(1,ia*3+3,j,g3)
      s3=z(2,ia*3+3,j,g3)
      r=r2 + r3
      s=s2 + s3
      zbr(1,ind(ia*3+1),j,g3) = r + r1
      zbr(2,ind(ia*3+1),j,g3) = s + s1
      r1=r1 - r*.5d0
      s1=s1 - s*.5d0
      r2=r2-r3
      s2=s2-s3
      zbr(1,ind(ia*3+2),j,g3) = r1 - s2*bb
      zbr(2,ind(ia*3+2),j,g3) = s1 + r2*bb
      zbr(1,ind(ia*3+3),j,g3) = r1 + s2*bb
      zbr(2,ind(ia*3+3),j,g3) = s1 - r2*bb
     enddo

     do ia=1,aft(ic)-1
      indx=ind(ia*3+1)-1
      cr2=trig(1,indx)
      ct2=trig(2,indx)
      cr3=trig(1,2*indx)
      ct3=trig(2,2*indx)
      cr2=cr2/cr3
      cr3p=.5d0*cr3
      bb=ris*cr3*0.8660254037844387d0
      do j=ma,mb
       r1=z(1,ia*3+1,j,g3)
       s1=z(2,ia*3+1,j,g3)
       r2=z(1,ia*3+2,j,g3) - z(2,ia*3+2,j,g3)*ct2
       s2=z(1,ia*3+2,j,g3)*ct2 + z(2,ia*3+2,j,g3)
       r3=z(1,ia*3+3,j,g3) - z(2,ia*3+3,j,g3)*ct3
       s3=z(1,ia*3+3,j,g3)*ct3 + z(2,ia*3+3,j,g3)
       r=cr2*r2 + r3
       s=cr2*s2 + s3
       zbr(1,ind(ia*3+1),j,g3) = r*cr3 + r1
       zbr(2,ind(ia*3+1),j,g3) = s*cr3 + s1
       r1=r1 - r*cr3p
       s1=s1 - s*cr3p
       r2=cr2*r2-r3
       s2=cr2*s2-s3
       zbr(1,ind(ia*3+2),j,g3) = r1 - s2*bb
       zbr(2,ind(ia*3+2),j,g3) = s1 + r2*bb
       zbr(1,ind(ia*3+3),j,g3) = r1 + s2*bb
       zbr(2,ind(ia*3+3),j,g3) = s1 - r2*bb
      enddo
     enddo

!   Treat radix 5
    else if (now(ic)==5) then
!    radix 5
!    sin(2.d0*pi/5.d0)
     sin2=ris*0.9510565162951536d0
     ia=0

!    First step of radix 5
     do j=ma,mb
      r1=z(1,ia*5+1,j,g3)
      s1=z(2,ia*5+1,j,g3)
      r2=z(1,ia*5+2,j,g3)
      s2=z(2,ia*5+2,j,g3)
      r3=z(1,ia*5+3,j,g3)
      s3=z(2,ia*5+3,j,g3)
      r4=z(1,ia*5+4,j,g3)
      s4=z(2,ia*5+4,j,g3)
      r5=z(1,ia*5+5,j,g3)
      s5=z(2,ia*5+5,j,g3)
      r25 = r2 + r5
      r34 = r3 + r4
      s25 = s2 - s5
      s34 = s3 - s4
      zbr(1,ind(ia*5+1),j,g3) = r1 + r25 + r34
      r = r1 + cos2*r25 + cos4*r34
      s = s25 + sin42*s34
      zbr(1,ind(ia*5+2),j,g3) = r - sin2*s
      zbr(1,ind(ia*5+5),j,g3) = r + sin2*s
      r = r1 + cos4*r25 + cos2*r34
      s = sin42*s25 - s34
      zbr(1,ind(ia*5+3),j,g3) = r - sin2*s
      zbr(1,ind(ia*5+4),j,g3) = r + sin2*s
      r25 = r2 - r5
      r34 = r3 - r4
      s25 = s2 + s5
      s34 = s3 + s4
      zbr(2,ind(ia*5+1),j,g3) = s1 + s25 + s34
      r = s1 + cos2*s25 + cos4*s34
      s = r25 + sin42*r34
      zbr(2,ind(ia*5+2),j,g3) = r + sin2*s
      zbr(2,ind(ia*5+5),j,g3) = r - sin2*s
      r = s1 + cos4*s25 + cos2*s34
      s = sin42*r25 - r34
      zbr(2,ind(ia*5+3),j,g3) = r + sin2*s
      zbr(2,ind(ia*5+4),j,g3) = r - sin2*s
     enddo

!    Second step of radix 5
     do ia=1,aft(ic)-1
      indx=ind(ia*5+1)-1
      cr2=trig(1,indx)
      ct2=trig(2,indx)
      cr3=trig(1,2*indx)
      ct3=trig(2,2*indx)
      cr4=trig(1,3*indx)
      ct4=trig(2,3*indx)
      cr5=trig(1,4*indx)
      ct5=trig(2,4*indx)
      do j=ma,mb
       r1=z(1,ia*5+1,j,g3)
       s1=z(2,ia*5+1,j,g3)
       r2=cr2*(z(1,ia*5+2,j,g3) - z(2,ia*5+2,j,g3)*ct2)
       s2=cr2*(z(1,ia*5+2,j,g3)*ct2 + z(2,ia*5+2,j,g3))
       r3=cr3*(z(1,ia*5+3,j,g3) - z(2,ia*5+3,j,g3)*ct3)
       s3=cr3*(z(1,ia*5+3,j,g3)*ct3 + z(2,ia*5+3,j,g3))
       r4=z(1,ia*5+4,j,g3) - z(2,ia*5+4,j,g3)*ct4
       s4=z(1,ia*5+4,j,g3)*ct4 + z(2,ia*5+4,j,g3)
       r5=z(1,ia*5+5,j,g3) - z(2,ia*5+5,j,g3)*ct5
       s5=z(1,ia*5+5,j,g3)*ct5 + z(2,ia*5+5,j,g3)
       r25 = r2 + r5*cr5
       r34 = r3 + r4*cr4
       s25 = s2 - s5*cr5
       s34 = s3 - s4*cr4
       zbr(1,ind(ia*5+1),j,g3) = r1 + r25 + r34
       r = r1 + cos2*r25 + cos4*r34
       s = s25 + sin42*s34
       zbr(1,ind(ia*5+2),j,g3) = r - sin2*s
       zbr(1,ind(ia*5+5),j,g3) = r + sin2*s
       r = r1 + cos4*r25 + cos2*r34
       s = sin42*s25 - s34
       zbr(1,ind(ia*5+3),j,g3) = r - sin2*s
       zbr(1,ind(ia*5+4),j,g3) = r + sin2*s
       r25 = r2 - r5*cr5
       r34 = r3 - r4*cr4
       s25 = s2 + s5*cr5
       s34 = s3 + s4*cr4
       zbr(2,ind(ia*5+1),j,g3) = s1 + s25 + s34
       r = s1 + cos2*s25 + cos4*s34
       s = r25 + sin42*r34
       zbr(2,ind(ia*5+2),j,g3) = r + sin2*s
       zbr(2,ind(ia*5+5),j,g3) = r - sin2*s
       r = s1 + cos4*s25 + cos2*s34
       s = sin42*r25 - r34
       zbr(2,ind(ia*5+3),j,g3) = r + sin2*s
       zbr(2,ind(ia*5+4),j,g3) = r - sin2*s
      enddo
     enddo

    else

!    All radices are treated
     write(message, '(a,a,a,a)' )ch10,&
&     ' sg_fftpx : BUG -',ch10,&
&     '  called with factors other than 2, 3, and 5'
!     call wrtout(06,message,'PERS')
!     call leave_new('PERS')
    endif

!  End of bit reversal

!  -------------------------------------------------------------------
   enddo
  enddo
 enddo
!$OMP END PARALLEL DO

 end subroutine
!!***

max_stars_repo_name: mreineck/benchfft

max_issues_repo_name: mreineck/benchfft

max_forks_repo_name: mreineck/benchfft


text:
C$Procedure      VLCOM3 ( Vector linear combination, 3 dimensions )

      SUBROUTINE VLCOM3 ( A, V1, B, V2, C, V3, SUM )

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

      DOUBLE PRECISION   A
      DOUBLE PRECISION   V1  ( 3 )
      DOUBLE PRECISION   B
      DOUBLE PRECISION   V2  ( 3 )
      DOUBLE PRECISION   C
      DOUBLE PRECISION   V3  ( 3 )
      DOUBLE PRECISION   SUM ( 3 )

C$ Brief_I/O
C
C      VARIABLE  I/O  DESCRIPTION
C      --------  ---  --------------------------------------------------
C      A          I   Coefficient of V1
C      V1         I   Vector in 3-space
C      B          I   Coefficient of V2
C      V2         I   Vector in 3-space
C      C          I   Coefficient of V3
C      V3         I   Vector in 3-space
C      SUM        O   Linear Vector Combination A*V1 + B*V2 + C*V3
C
C$ Detailed_Input
C
C      A     is a double precision number.
C
C      V1    is a double precision 3-dimensional vector.
C
C      B     is a double precision number.
C
C      V2    is a double precision 3-dimensional vector.
C
C      C     is a double precision number.
C
C      V3    is a double precision 3-dimensional vector.
C
C$ Detailed_Output
C
C      SUM   is a double precision 3-dimensional vector which contains
C            the linear combination A*V1 + B*V2 + C*V3
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C      For each index from 1 to 3, this routine implements in FORTRAN
C      code the expression:
C
C      SUM(I) = A*V1(I) + B*V2(I) + C*V3(I)
C
C      No error checking is performed to guard against numeric overflow.
C
C$ Examples
C
C      Often one has the components (A,B,C) of a vector in terms
C      of a basis V1, V2, V3.  The vector represented by (A,B,C) can
C      be obtained immediately from the call
C
C      CALL VLCOM3 ( A, V1, B, V2, C, V3,   VECTOR )
C
C$ Restrictions
C
C      None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber     (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.0, 1-NOV-1990 (WLT)
C
C-&

C$ Index_Entries
C
C     linear combination of three 3-dimensional vectors
C
C-&

      SUM(1) = A*V1(1) + B*V2(1) + C*V3(1)
      SUM(2) = A*V1(2) + B*V2(2) + C*V3(2)
      SUM(3) = A*V1(3) + B*V2(3) + C*V3(3)

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework