submodule(forlab_linalg) forlab_linalg_cat

    implicit none

contains

    module procedure horzcat_r_1_sp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_r_1_sp

    module procedure horzcat_r_2_sp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_r_2_sp

    module procedure horzcat_r_21_sp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_r_21_sp

    module procedure horzcat_r_12_sp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_r_12_sp
    module procedure horzcat_r_1_dp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_r_1_dp

    module procedure horzcat_r_2_dp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_r_2_dp

    module procedure horzcat_r_21_dp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_r_21_dp

    module procedure horzcat_r_12_dp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_r_12_dp
    module procedure horzcat_c_1_sp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_c_1_sp

    module procedure horzcat_c_2_sp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_c_2_sp

    module procedure horzcat_c_21_sp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_c_21_sp

    module procedure horzcat_c_12_sp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_c_12_sp
    module procedure horzcat_c_1_dp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_c_1_dp

    module procedure horzcat_c_2_dp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_c_2_dp

    module procedure horzcat_c_21_dp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_c_21_dp

    module procedure horzcat_c_12_dp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_c_12_dp
    module procedure horzcat_i_1_int8
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_i_1_int8

    module procedure horzcat_i_2_int8
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_i_2_int8

    module procedure horzcat_i_21_int8
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_i_21_int8

    module procedure horzcat_i_12_int8
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_i_12_int8
    module procedure horzcat_i_1_int16
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_i_1_int16

    module procedure horzcat_i_2_int16
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_i_2_int16

    module procedure horzcat_i_21_int16
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_i_21_int16

    module procedure horzcat_i_12_int16
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_i_12_int16
    module procedure horzcat_i_1_int32
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_i_1_int32

    module procedure horzcat_i_2_int32
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_i_2_int32

    module procedure horzcat_i_21_int32
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_i_21_int32

    module procedure horzcat_i_12_int32
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_i_12_int32
    module procedure horzcat_i_1_int64
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(max(m1, m2), 2)
    result(1:m1, 1) = x1
    result(1:m2, 2) = x2

    end procedure horzcat_i_1_int64

    module procedure horzcat_i_2_int64
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(max(m1, m2), n1 + n2)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1:) = A2

    end procedure horzcat_i_2_int64

    module procedure horzcat_i_21_int64
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(max(m1, m2), n1 + 1)
    result(1:m1, 1:n1) = A1
    result(1:m2, n1 + 1) = x2
    return
    end procedure horzcat_i_21_int64

    module procedure horzcat_i_12_int64
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(max(m1, m2), n2 + 1)
    result(1:m1, 1) = x1
    result(1:m2, 2:) = A2

    end procedure horzcat_i_12_int64
    module procedure vertcat_r_1_sp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_r_1_sp

    module procedure vertcat_r_2_sp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_r_2_sp

    module procedure vertcat_r_21_sp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_r_21_sp

    module procedure vertcat_r_12_sp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_r_12_sp
    module procedure vertcat_r_1_dp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_r_1_dp

    module procedure vertcat_r_2_dp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_r_2_dp

    module procedure vertcat_r_21_dp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_r_21_dp

    module procedure vertcat_r_12_dp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_r_12_dp
    module procedure vertcat_c_1_sp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_c_1_sp

    module procedure vertcat_c_2_sp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_c_2_sp

    module procedure vertcat_c_21_sp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_c_21_sp

    module procedure vertcat_c_12_sp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_c_12_sp
    module procedure vertcat_c_1_dp
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_c_1_dp

    module procedure vertcat_c_2_dp
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_c_2_dp

    module procedure vertcat_c_21_dp
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_c_21_dp

    module procedure vertcat_c_12_dp
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_c_12_dp
    module procedure vertcat_i_1_int8
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_i_1_int8

    module procedure vertcat_i_2_int8
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_2_int8

    module procedure vertcat_i_21_int8
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_i_21_int8

    module procedure vertcat_i_12_int8
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_12_int8
    module procedure vertcat_i_1_int16
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_i_1_int16

    module procedure vertcat_i_2_int16
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_2_int16

    module procedure vertcat_i_21_int16
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_i_21_int16

    module procedure vertcat_i_12_int16
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_12_int16
    module procedure vertcat_i_1_int32
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_i_1_int32

    module procedure vertcat_i_2_int32
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_2_int32

    module procedure vertcat_i_21_int32
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_i_21_int32

    module procedure vertcat_i_12_int32
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_12_int32
    module procedure vertcat_i_1_int64
    integer :: m1, m2

    m1 = size(x1)
    m2 = size(x2)
    result = zeros(m1 + m2, 1)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1) = x2

    end procedure vertcat_i_1_int64

    module procedure vertcat_i_2_int64
    integer :: m1, n1, m2, n2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(A2, 1)
    n2 = size(A2, 2)

    result = zeros(m1 + m2, max(n1, n2))
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_2_int64

    module procedure vertcat_i_21_int64
    integer :: m1, n1, m2

    m1 = size(A1, 1)
    n1 = size(A1, 2)
    m2 = size(x2)
    result = zeros(m1 + m2, n1)
    result(1:m1, 1:n1) = A1
    result(m1 + 1:m1 + m2, 1) = x2
    return
    end procedure vertcat_i_21_int64

    module procedure vertcat_i_12_int64
    integer :: m1, m2, n2

    m1 = size(x1)
    m2 = size(A2, 1)
    n2 = size(A2, 2)
    result = zeros(m1 + m2, n2)
    result(1:m1, 1) = x1
    result(m1 + 1:m1 + m2, 1:n2) = A2

    end procedure vertcat_i_12_int64

end submodule forlab_linalg_cat


text:
 Biography
As a first year undergraduate student at UC Davis, Mitty first hit the http://media.www.californiaaggie.com/media/storage/paper981/news/2006/08/07/CampusNews/Yearbook.Support.Group.Takes.Shape.On.Facebook2143670.shtml?sourcedomainwww.californiaaggie.com&MIIHostmedia.collegepublisher.com news stream of Davis with the facebook group Aggies for a Yearbook. This facebook group was created to start an effort on creating a new yearbook organization, as ASUCD had eliminated the Yearbook during the Spring of 2006. After seeing many supporters of the creation of a new yearbook organization, Mitty started the club, The Mustang, as the yearbook club. He later dissolved the club after realizing that he would have to devote much of his life towards yearbook, rather than his true passion: service.

Mitty also was actively involved with with ASUCD student government, serving as a member of the Internal Affairs Commission in 2007, and as the Service Representative of the ASUCD Outreach Assembly in 2007 and 2008.

Mitty is now the President of the Campus Rotaract Club of Davis.

Whups  I uploaded a logo just as you were uploading one.  :)  Feel free to use yours instead.  I didnt mean to duplicate your efforts. wiki:gnome:Users/JabberWokky Evan JabberWokky Edwards

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
! other Shroud Project Developers.
! See the top-level COPYRIGHT file for details.
!
! SPDX-License-Identifier: (BSD-3-Clause)
! #######################################################################
!
! Test Fortran API generated from pointers.yaml.
! Used with pointers-c and pointers-cxx.
!
program tester
  use fruit
  use iso_c_binding
  use pointers_mod
  implicit none
  real(C_DOUBLE), parameter :: pi = 3.1415926_C_DOUBLE
  logical ok

  call init_fruit

  call test_functions
  call test_functions2
  call test_swig
  call test_char_arrays
  call test_out_ptrs
  call test_out_alloc
  call test_nested_ptrs
  call test_dimension

  call fruit_summary
  call fruit_finalize

  call is_all_successful(ok)
  if (.not. ok) then
     call exit(1)
  endif

contains

  subroutine test_functions
    integer(c_int) iargin, iarginout, iargout
    real(c_double) :: in_double(5), out_double(5)
    integer(c_int) :: out_int(4)
    integer(c_int) :: nvalues, values1(3), values2(3)

    call set_case_name("test_functions")

    ! set global_int.
    call intargs_in(5)

    iargout = 0
    call intargs_out(iargout)
    call assert_equals(5, iargout) ! get global_int

    iarginout = 6
    call intargs_inout(iarginout)  ! set global_int
    call intargs_out(iargout)      ! get global_int
    call assert_equals(6, iargout)
    call assert_equals(7, iarginout) ! incremented iarginout

    iargin    = 1
    iarginout = 2
    iargout   = -1
    call intargs(iargin, iarginout, iargout)
    call assert_true(iarginout == 1)
    call assert_true(iargout   == 2)

    in_double = [0.0*pi, 0.5*pi, pi, 1.5*pi, 2.0*pi]
    out_double = 0.0
    call cos_doubles(in_double, out_double)
    call assert_true(all(abs(out_double - cos(in_double)) < 1.e-08 ))

    out_int = 0
    call truncate_to_int([1.2d0, 2.3d0, 3.4d0, 4.5d0], out_int)
    call assert_true(all(out_int == [1, 2, 3, 4]))

    values1 = 0
    call get_values(nvalues, values1)
    call assert_equals(3, nvalues)
    call assert_true(all(values1(1:3) == [1, 2, 3]))

    values1 = 0
    values2 = 0
    call get_values2(values1, values2)
    call assert_true(all(values1(1:3) == [1, 2, 3]))
    call assert_true(all(values2(1:3) == [11, 12, 13]))

    values1 = 0
    call iota_dimension(nvalues, values1)
    call assert_equals(3, nvalues)
    call assert_true(all(values1(1:3) == [1, 2, 3]))

  end subroutine test_functions

  subroutine test_functions2
    integer(c_int) rv_int, out(3), values(5)

    call set_case_name("test_functions2")

    call sum([1,2,3,4,5], rv_int)
    call assert_true(rv_int .eq. 15, "sum")

    out = 0
    call fill_int_array(out)
    call assert_true(all(out(1:3) == [1, 2, 3]), "fillIntArray")

    values = [1, 2, 3, 4, 5]
    call increment_int_array(values)
    call assert_true(all(values(1:5) == [2, 3, 4, 5, 6]), "incrementIntArray")

  end subroutine test_functions2

  subroutine test_swig

    real(C_DOUBLE) zero(10)
    integer(C_INT) sum, count(5)

    call set_case_name("test_swig")

    call fill_with_zeros(zero)

    count = [1, 2, 3, 4, 5]
    sum = accumulate(count)
    call assert_equals(15, sum)

  end subroutine test_swig

  subroutine test_char_arrays
    integer nchar
    character(10) :: in(3) = [ &
         "dog       ", &
         "cat       ", &
         "monkey    "  &
         ]
    character(10), target :: word1, word2, word3
    type(C_PTR)  cin(4)

    call set_case_name("test_char_arrays")

    ! Call the bufferify function.
    ! It will copy strings to create char ** variable.
    nchar = accept_char_array_in(in)
    call assert_equals(len_trim(in(1)), nchar, "acceptCharArrayIn")

    ! Build up a native char ** variable and pass to C.
    ! Caller is responsibile for explicilty NULL terminating.
    word1 = "word1" // C_NULL_CHAR
    word2 = "word2+" // C_NULL_CHAR
    word3 = "word3long" // C_NULL_CHAR
    cin(1) = c_loc(word1)
    cin(2) = c_loc(word2)
    cin(3) = c_loc(word3)
    cin(4) = C_NULL_PTR
    nchar = c_accept_char_array_in(cin)
    call assert_equals(5, nchar, "acceptCharArrayIn") ! 5 = len(word1) - trailing NULL

  end subroutine test_char_arrays

  subroutine test_out_ptrs
    integer(C_INT) :: ivalue, narray
    integer(C_INT), target :: ivalue1, ivalue2
    integer(C_INT), pointer :: iscalar, irvscalar
    integer(C_INT), pointer :: iarray(:)
    integer(C_INT), pointer :: irvarray(:)
    type(C_PTR) :: cptr_scalar, cptr_array
    type(C_PTR) :: void
    type(C_PTR) :: cptr_arrays(2)

    call set_case_name("test_out_ptrs")

    call set_global_int(0)

    nullify(iscalar)
    call get_ptr_to_scalar(iscalar)
    call assert_true(associated(iscalar))
    call assert_equals(0, iscalar)

    ! iscalar points to global_int in pointers.c.
    call set_global_int(5)
    call assert_equals(5, iscalar)

    nullify(iarray)
    call get_ptr_to_fixed_array(iarray)
    call assert_true(associated(iarray))
    call assert_equals(10, size(iarray))
    call assert_true(lbound(iarray,1) == 1, &
         "getPtrToFixedArray - lbound")
    call assert_true(ubound(iarray,1) == 10, &
         "getPtrToFixedArray - ubound")
    iarray = 0
    call assert_equals(0, sum_fixed_array(), &
         "sumFixedArray - initial")
    ! Make sure we're assigning to global_array.
    iarray(1) = 1
    iarray(10) = 2
    call assert_equals(3, sum_fixed_array(), &
         "sumFixedArray - changes")

    ! Returns global_array in pointers.c.
    nullify(iarray)
    call get_ptr_to_dynamic_array(iarray)
    call assert_true(associated(iarray))
    call assert_true(size(iarray) == 10)

    ! Call C version directly.
    cptr_array = C_NULL_PTR
    call c_get_ptr_to_dynamic_array(cptr_array, narray)
    call assert_true(c_associated(cptr_array))
    call assert_true(narray == 10)

    ! Returns global_array in pointers.c.
    ! iarray is used later for deref(raw) tests. Do not reset.
    nullify(iarray)
    call get_ptr_to_func_array(iarray)
    call assert_true(associated(iarray))
    call assert_true(size(iarray) == 10)
    call assert_true(lbound(iarray,1) == 1)
    call assert_true(ubound(iarray,1) == 10)

    call get_raw_ptr_to_scalar(cptr_scalar)
    call assert_true(c_associated(cptr_scalar), &
         "getRawPtrToScalar - c_associated")
    ! associated with global_int in pointers.c
    call assert_true(c_associated(cptr_scalar, c_loc(iscalar)), &
         "getRawPtrToScalar - c_associated(iscalar)")

    call get_raw_ptr_to_scalar_force(cptr_scalar)
    call assert_true(c_associated(cptr_scalar), &
         "getRawPtrToScalarForce - c_associated")
    ! associated with global_int in pointers.c
    call assert_true(c_associated(cptr_scalar, c_loc(iscalar)), &
         "getRawPtrToScalarForce - c_associated(iscalar)")

    cptr_array = C_NULL_PTR
    call get_raw_ptr_to_fixed_array(cptr_array)
    call assert_true(c_associated(cptr_array), &
         "getRawPtrToFixedArray - c_associated")
    ! associated with global_fixed_array in pointers.c
    call assert_true(c_associated(cptr_array, c_loc(iarray(1))), &
         "getRawPtrToFixedArray - c_associated(iarray)")

    cptr_array = C_NULL_PTR
    call get_raw_ptr_to_fixed_array_force(cptr_array)
    call assert_true(c_associated(cptr_array), &
         "getRawPtrToFixedArrayForce - c_associated")
    ! associated with global_fixed_array in pointers.c
    call assert_true(c_associated(cptr_array, c_loc(iarray(1))), &
         "getRawPtrToFixedArrayForce - c_associated(iarray)")

    ! Return pointer to global_int as a type(C_PTR).
    ! via interface
    void = C_NULL_PTR
    void = return_address1(1)
    call assert_true(c_associated(void, cptr_scalar))
    ! via wrapper
    void = C_NULL_PTR
    void = return_address2(1)
    call assert_true(c_associated(void, cptr_scalar))
    ! via argument
    void = C_NULL_PTR
    call fetch_void_ptr(void)
    call assert_true(c_associated(void, cptr_scalar))

    ! Pass array of pointers  (void **)
    ivalue1 = 10
    ivalue2 = 4
    cptr_arrays(1) = c_loc(ivalue1)
    cptr_arrays(2) = c_loc(ivalue2)
    call assert_equals(14, void_ptr_array(cptr_arrays))

    ! ***** Non-const results
    ! Return pointer to global_int as a fortran pointer.
    nullify(irvscalar)
    irvscalar => return_int_ptr_to_scalar()
    call assert_true(associated(irvscalar, iscalar))
    call set_global_int(7)
    call assert_equals(7, irvscalar)

    ! ivalue is not a pointer.
    call set_global_int(8)
    ivalue = return_int_scalar()
    call assert_equals(8, ivalue)

    ! Return pointer to global_fixed_int as a fortran pointer.
    nullify(irvarray)
    irvarray => return_int_ptr_to_fixed_array()
    call assert_true(associated(irvarray))
    call assert_true(size(irvarray) == 10)
    call assert_true(associated(irvscalar, iscalar))

    ! ***** const results
    ! Return pointer to global_int as a fortran pointer.
    nullify(irvscalar)
    irvscalar => return_int_ptr_to_const_scalar()
    call assert_true(associated(irvscalar, iscalar))

    ! Return pointer to global_fixed_int as a fortran pointer.
    nullify(irvarray)
    irvarray => return_int_ptr_to_fixed_const_array()
    call assert_true(associated(irvarray))
    call assert_true(size(irvarray) == 10)
    call assert_true(associated(irvscalar, iscalar))

    ! +deref(scalar)
    ivalue = return_int_scalar()

    ! Return pointer to global_int.
    void = return_int_raw()
    call assert_true(c_associated(void, c_loc(irvscalar)))
    call assert_true(c_associated(void, c_loc(iscalar)))

    void = return_int_raw_with_args("with args")
    call assert_true(c_associated(void, c_loc(irvscalar)))
    call assert_true(c_associated(void, c_loc(iscalar)))

  end subroutine test_out_ptrs

  subroutine test_out_alloc
    !-- XXX - copied from test_out_ptrs, need to fill in.
!-    integer(C_INT) :: ivalue, narray
!-    integer(C_INT), target :: ivalue1, ivalue2
!-    integer(C_INT), pointer :: iscalar, irvscalar
    integer(C_INT), allocatable :: iarray(:)
    integer(C_INT), allocatable :: irvarray(:)
!-    type(C_PTR) :: cptr_scalar, cptr_array
!-    type(C_PTR) :: void
!-    type(C_PTR) :: cptr_arrays(2)

    call set_case_name("test_out_alloc")

!-    call set_global_int(0)
!-
!-    nullify(iscalar)
!-    call get_ptr_to_scalar(iscalar)
!-    call assert_true(associated(iscalar))
!-    call assert_equals(0, iscalar)
!-
!-    ! iscalar points to global_int in pointers.c.
!-    call set_global_int(5)
!-    call assert_equals(5, iscalar)
!-
!-    nullify(iarray)
    call get_alloc_to_fixed_array(iarray)
    call assert_true(allocated(iarray))
    call assert_equals(10, size(iarray))
    call assert_true(lbound(iarray,1) == 1, &
         "getPtrToFixedArray - lbound")
    call assert_true(ubound(iarray,1) == 10, &
         "getPtrToFixedArray - ubound")
    iarray = 0
!-    call assert_equals(0, sum_fixed_array(), &
!-         "sumFixedArray - initial")
!-    ! Make sure we're assigning to global_array.
!-    iarray(1) = 1
!-    iarray(10) = 2
!-    call assert_equals(3, sum_fixed_array(), &
!-         "sumFixedArray - changes")
!-
!-    ! Returns global_array in pointers.c.
!-    nullify(iarray)
!-    call get_ptr_to_dynamic_array(iarray)
!-    call assert_true(associated(iarray))
!-    call assert_true(size(iarray) == 10)
!-
!-    ! Call C version directly.
!-    cptr_array = C_NULL_PTR
!-    call c_get_ptr_to_dynamic_array(cptr_array, narray)
!-    call assert_true(c_associated(cptr_array))
!-    call assert_true(narray == 10)
!-
!-    ! Returns global_array in pointers.c.
!-    ! iarray is used later for deref(raw) tests. Do not reset.
!-    nullify(iarray)
!-    call get_ptr_to_func_array(iarray)
!-    call assert_true(associated(iarray))
!-    call assert_true(size(iarray) == 10)
!-    call assert_true(lbound(iarray,1) == 1)
!-    call assert_true(ubound(iarray,1) == 10)
!-
!-    call get_raw_ptr_to_scalar(cptr_scalar)
!-    call assert_true(c_associated(cptr_scalar), &
!-         "getRawPtrToScalar - c_associated")
!-    ! associated with global_int in pointers.c
!-    call assert_true(c_associated(cptr_scalar, c_loc(iscalar)), &
!-         "getRawPtrToScalar - c_associated(iscalar)")
!-
!-    call get_raw_ptr_to_scalar_force(cptr_scalar)
!-    call assert_true(c_associated(cptr_scalar), &
!-         "getRawPtrToScalarForce - c_associated")
!-    ! associated with global_int in pointers.c
!-    call assert_true(c_associated(cptr_scalar, c_loc(iscalar)), &
!-         "getRawPtrToScalarForce - c_associated(iscalar)")
!-
!-    cptr_array = C_NULL_PTR
!-    call get_raw_ptr_to_fixed_array(cptr_array)
!-    call assert_true(c_associated(cptr_array), &
!-         "getRawPtrToFixedArray - c_associated")
!-    ! associated with global_fixed_array in pointers.c
!-    call assert_true(c_associated(cptr_array, c_loc(iarray(1))), &
!-         "getRawPtrToFixedArray - c_associated(iarray)")
!-
!-    cptr_array = C_NULL_PTR
!-    call get_raw_ptr_to_fixed_array_force(cptr_array)
!-    call assert_true(c_associated(cptr_array), &
!-         "getRawPtrToFixedArrayForce - c_associated")
!-    ! associated with global_fixed_array in pointers.c
!-    call assert_true(c_associated(cptr_array, c_loc(iarray(1))), &
!-         "getRawPtrToFixedArrayForce - c_associated(iarray)")
!-
!-    ! Return pointer to global_int as a type(C_PTR).
!-    ! via interface
!-    void = C_NULL_PTR
!-    void = return_address1(1)
!-    call assert_true(c_associated(void, cptr_scalar))
!-    ! via wrapper
!-    void = C_NULL_PTR
!-    void = return_address2(1)
!-    call assert_true(c_associated(void, cptr_scalar))
!-    ! via argument
!-    void = C_NULL_PTR
!-    call fetch_void_ptr(void)
!-    call assert_true(c_associated(void, cptr_scalar))
!-
!-    ! Pass array of pointers  (void **)
!-    ivalue1 = 10
!-    ivalue2 = 4
!-    cptr_arrays(1) = c_loc(ivalue1)
!-    cptr_arrays(2) = c_loc(ivalue2)
!-    call assert_equals(14, void_ptr_array(cptr_arrays))
!-
!-    ! ***** Non-const results
!-    ! Return pointer to global_int as a fortran pointer.
!-    nullify(irvscalar)
!-    irvscalar => return_int_ptr_to_scalar()
!-    call assert_true(associated(irvscalar, iscalar))
!-    call set_global_int(7)
!-    call assert_equals(7, irvscalar)
!-
!-    ! ivalue is not a pointer.
!-    call set_global_int(8)
!-    ivalue = return_int_scalar()
!-    call assert_equals(8, ivalue)

    ! Return copy of global_fixed_int as a fortran allocatable.
    call assert_false(allocated(irvarray))
    irvarray = return_int_alloc_to_fixed_array()
    call assert_true(allocated(irvarray))
    call assert_true(size(irvarray) == 10)
!-  XXX - check values of irvarray
    deallocate(irvarray)

!-    ! ***** const results
!-    ! Return pointer to global_int as a fortran pointer.
!-    nullify(irvscalar)
!-    irvscalar => return_int_ptr_to_const_scalar()
!-    call assert_true(associated(irvscalar, iscalar))
!-
!-    ! Return pointer to global_fixed_int as a fortran pointer.
!-    nullify(irvarray)
!-    irvarray => return_int_ptr_to_fixed_const_array()
!-    call assert_true(associated(irvarray))
!-    call assert_true(size(irvarray) == 10)
!-    call assert_true(associated(irvscalar, iscalar))
!-
!-    ! +deref(scalar)
!-    ivalue = return_int_scalar()
!-
!-    ! Return pointer to global_int.
!-    void = return_int_raw()
!-    call assert_true(c_associated(void, c_loc(irvscalar)))
!-    call assert_true(c_associated(void, c_loc(iscalar)))
!-
!-    void = return_int_raw_with_args("with args")
!-    call assert_true(c_associated(void, c_loc(irvscalar)))
!-    call assert_true(c_associated(void, c_loc(iscalar)))

  end subroutine test_out_alloc

  subroutine test_nested_ptrs
    type(C_PTR) addr, rvaddr
    type(C_PTR), pointer :: array2d(:)
    integer(C_INT), pointer :: row1(:), row2(:)
    integer total

    call set_case_name("test_nested_ptrs")

    addr = C_NULL_PTR
    call get_raw_ptr_to_int2d(addr)
    call assert_equals(15, check_int2d(addr), "getRawPtrToInt2d")

    call c_f_pointer(addr, array2d, [2])
    call c_f_pointer(array2d(1), row1, [3])
    call c_f_pointer(array2d(2), row2, [2])

    total = row1(1) + row1(2) + row1(3) + row2(1) + row2(2)
    call assert_equals(15, total)

    ! function result
    rvaddr = return_raw_ptr_to_int2d()
    call assert_true(c_associated(rvaddr, addr), "returnRawPtrToInt2d")

  end subroutine test_nested_ptrs

  subroutine test_dimension
    ! Test +dimension(10,20) +intent(in)  together.
!    integer(C_INT) arg(2,3)
    integer(C_INT) arg2(20,30)

    call set_case_name("test_dimension")

    ! gcc Warning: Actual argument contains too few elements for dummy argument 'arg'
    ! intel error #7983: The storage extent of the dummy argument exceeds that of the actual argument.   [ARG]

!    call dimension_in(arg)

    ! compilers seem ok with too much space
    call dimension_in(arg2)

  end subroutine test_dimension

end program tester

max_stars_repo_name: ExternalRepositories/shroud

max_issues_repo_name: ExternalRepositories/shroud

max_forks_repo_name: ExternalRepositories/shroud


text:
MODULE m_slabdim
  USE m_juDFT
CONTAINS
  SUBROUTINE slab_dim(atoms,nsld)
    !***********************************************************************
    !     This subroutine calculates  the number of layers in the slab
    !
    !                                   Yury Koroteev 2003-09-30
    !***********************************************************************
    !                     ABBREVIATIONS
    !
    ! natd                  : in, the number of atoms in the film
    ! pos(3,natd)         : in, the coordinates of atoms in the film
    ! ntypd,ntype           : in, the number of mt-sphere types
    ! neq(ntypd)            : in, the number of mt-spheres of the same type
    !-----------------------------------------------------------------------
    ! nsld                  : out, the number of layers in the film
    !-----------------------------------------------------------------------
    ! znz(nsl)              : work, the z-ordinate of mt-spheres in
    !                               the nsl-layer
    !-----------------------------------------------------------------------
    !
    USE m_types_setup
    IMPLICIT NONE

    TYPE(t_atoms),INTENT(IN)   :: atoms
    !	..
    !       ..Scalar Argument
    INTEGER, INTENT (OUT) :: nsld
    !       ..
    !       ..Array Arguments
    !       ..
    !       ..Local Scalars
    INTEGER  iz,i,j,na,nz
    REAL    zs
    !       ..
    !       ..Local Arrays
    REAL    znz(atoms%nat)
    !       ..
    !    ----------------------------------------------
    REAL,PARAMETER:: epsz=1.e-3
    !    ----------------------------------------------
    !
    ! --->  Calculate the number of the film layers (nsld)
    !
    znz(1) = atoms%pos(3,1)
    nz = 1
    na = 0
    DO i=1,atoms%ntype
       equivAtomsLoop: DO j=1,atoms%neq(i)
          na = na + 1
          zs = atoms%pos(3,na)
          DO iz=1,nz
             IF(ABS(zs-znz(iz)).LT.epsz) CYCLE equivAtomsLoop
          END DO
          nz = nz+1
          znz(nz) = zs
       END DO equivAtomsLoop
    END DO
    nsld = nz
    IF(nsld>atoms%nat)   CALL juDFT_error("nsld.GT.atoms%nat ",calledby="slab_dim")
    !
  END SUBROUTINE slab_dim
END MODULE m_slabdim

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
Ive been a resident of the city of Davis since early in 2004.  I figured out how nice a town this is when, of all the crazy things Ive done, I completed my tour of Borders Books in the greater Sacramento area with a stop off at Borders Davis.  Im a liberal/progressive/whatever.  I think Im probably way too intense for my own good sometimes but hey, gotta be something, ya know.  I have come to have a deep affinity for our town and the university adjoining it.  I spend time both on campus and in town.  When it comes to all things that float and are powered, either by propeller or sail, you can classify me as obsessed (even though I dont own a boatsomebody want company on the water.....well here I am).  Im a firm believer in the Constitution especially the Bill of Rights.  I was born in the 70s under the sign of the bull.  My desire for the city of Davis is that it become an island of beauty, descency, and kindness in a sea of productive fields.

Vitals:

63

200+ lbs

Blue Eyes

Brown Hair(long haired and proud of it)

Right Handed


text:
C *******************************************************************
C COPYRIGHT (c) 1994 Council for the Central Laboratory
*                    of the Research Councils
C All rights reserved.
C
C None of the comments in this Copyright notice between the lines
C of asterisks shall be removed or altered in any way.
C
C This Package is intended for compilation without modification,
C so most of the embedded comments have been removed.
C
C ALL USE IS SUBJECT TO LICENCE. For full details of the ACADEMIC
C SOFTWARE LICENCE, see http://hsl.rl.ac.uk/hsl2007/cou/academic.html
C
C Please note that for an ACADEMIC Licence:
C
C 1. The Packages may only be used for academic research or teaching
C    purposes by the Licensee, and must not be copied by the Licensee for
C    use by any other persons. Use of the Packages in any commercial
C    application shall be subject to prior written agreement between
C    Hyprotech UK Limited and the Licensee on suitable terms and
C    conditions, which will include financial conditions.
C 2. All information on the Package is provided to the Licensee on the
C    understanding that the details thereof are confidential.
C 3. All publications issued by the Licensee that include results obtained
C    with the help of one or more of the Packages shall acknowledge the
C    use of the Packages. The Licensee will notify the Numerical Analysis
C    Group at Rutherford Appleton Laboratory (STFC) of any such publication.
C 4. The Packages may be modified by or on behalf of the Licensee
C    for such use in research applications but at no time shall such
C    Packages or modifications thereof become the property of the
C    Licensee. The Licensee shall make available free of charge to the
C    copyright holder for any purpose all information relating to
C    any modification.
C 5. Neither STFC nor Hyprotech UK Limited shall be liable for any
C    direct or consequential loss or damage whatsoever arising out of
C    the use of Packages by the Licensee.
C *******************************************************************
C
C Original date 30 November 1995

C 12th July 2004 Version 1.0.0. Version numbering added.

      SUBROUTINE MA51AD(M,N,LA,IRN,KEEP,RANK,ROWS,COLS,W)
      INTEGER M,N,LA,IRN(LA),KEEP(*),RANK,ROWS(M),COLS(N),W(*)
      INTEGER I,IQB(1),IPTRD,IPTRL,IPTRO,IPTRU,J,JB,J1,K1,K2,
     +        KBLOCK,MBLOCK,NB,NBLOCK,NC,NE,NP,NR,RANKB
      EXTERNAL MA51ZD
      IQB(1) = 0
      IPTRL = M + N
      IPTRU = IPTRL + N
      IPTRD = IPTRU + N
      IPTRO = IPTRD + N + 1
      NBLOCK = IPTRO + N - 1
      MBLOCK = NBLOCK + 1
      KBLOCK = MBLOCK + 1
      NB = KEEP(KBLOCK+3)
      NE = KEEP(IPTRO+N+1)-1
      K2 = 0
      RANK = 0
      DO 20 JB = 1,NB
         NC = KEEP(NBLOCK+3*JB)
         NR = NC
         IF(NB.EQ.1) NR = M
         K1 = K2 + 1
         K2 = K2 + NC
         NP = KEEP(MBLOCK+3*JB)
         IF (NP.LT.0) THEN
            DO 10 J = K1,K2
               ROWS(J) = 1
               COLS(J) = 1
   10       CONTINUE
            RANK = RANK + NC
         ELSE
            J1 = 1
            IF (JB.GT.1) J1 = KEEP(KBLOCK+3*JB)
            CALL MA51ZD(NR,NC,IQB,NP,LA+1-NE-J1,IRN(NE+J1),
     +                 KEEP(IPTRL+K1),KEEP(IPTRU+K1),RANKB,
     +                 ROWS(K1),COLS(K1),W)
            RANK = RANK + RANKB
         END IF
   20 CONTINUE
      DO 30 I = 1,M
         W(I) = ROWS(KEEP(I))
   30 CONTINUE
      K1 = 1
      K2 = M
      DO 40 I = 1,M
         IF(W(I).EQ.1)THEN
            ROWS(K1) = I
            K1 = K1 + 1
         ELSE
            ROWS(K2) = I
            K2 = K2 - 1
         END IF
   40 CONTINUE
      DO 50 I = 1,N
         W(KEEP(M+I)) = COLS(I)
   50 CONTINUE
      K1 = 1
      K2 = N
      DO 60 I = 1,N
         IF(W(I).EQ.1)THEN
            COLS(K1) = I
            K1 = K1 + 1
         ELSE
            COLS(K2) = I
            K2 = K2 - 1
         END IF
   60 CONTINUE
      END
      SUBROUTINE MA51BD(M,N,IQ,NP,LFACT,IRNF,IPTRL,
     +                  IPTRU,RANK,ROWS,COLS,W)
      INTEGER M,N,IQ(*),NP,LFACT,IRNF(LFACT),IPTRL(N),IPTRU(N)
      INTEGER RANK,ROWS(M),COLS(N),W(*)
      INTEGER I,K1,K2
      EXTERNAL MA51ZD
      CALL MA51ZD(M,N,IQ,NP,LFACT,IRNF,IPTRL,
     +                  IPTRU,RANK,ROWS,COLS,W)
      DO 10 I = 1,M
         W(I) = ROWS(I)
   10 CONTINUE
      K1 = 1
      K2 = M
      DO 20 I = 1,M
         IF(W(I).EQ.1)THEN
            ROWS(K1) = I
            K1 = K1 + 1
         ELSE
            ROWS(K2) = I
            K2 = K2 - 1
         END IF
   20 CONTINUE
      DO 30 I = 1,N
         W(I) = COLS(I)
   30 CONTINUE
      K1 = 1
      K2 = N
      DO 40 I = 1,N
         IF(W(I).EQ.1)THEN
            COLS(K1) = I
            K1 = K1 + 1
         ELSE
            COLS(K2) = I
            K2 = K2 - 1
         END IF
   40 CONTINUE
      END
      SUBROUTINE MA51CD(M,N,LA,A,IRN,KEEP,SGNDET,LOGDET,W)
      INTEGER M,N,LA
      DOUBLE PRECISION A(LA)
      INTEGER IRN(LA),KEEP(*),SGNDET
      DOUBLE PRECISION LOGDET
      INTEGER W(N)
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D+0)
      INTEGER I,IQB(1),IPTRD,IPTRL,IPTRO,IPTRU,J,JB,J1,K,K1,K2,
     +        KBLOCK,L,MBLOCK,NB,NBLOCK,NC,NE,NP,NR,SGNDT
      DOUBLE PRECISION PIV,LOGDT
      EXTERNAL MA51DD
      IQB(1) = 0
      IPTRL = M + N
      IPTRU = IPTRL + N
      IPTRD = IPTRU + N
      IPTRO = IPTRD + N + 1
      NBLOCK = IPTRO + N - 1
      MBLOCK = NBLOCK + 1
      KBLOCK = MBLOCK + 1
      NB = KEEP(KBLOCK+3)
      NE = KEEP(IPTRO+N+1)-1
      LOGDET = ZERO
      SGNDET = 1
      IF (M.NE.N) THEN
          SGNDET = 0
          RETURN
      END IF
      DO 10 I = 1,N
         W(I) = 1
   10 CONTINUE
      DO 40 I = 1,N
         K = KEEP(I)
         IF (W(K).EQ.1) THEN
            DO 30 J = 1,N
               L = KEEP(K)
               W(K) = 0
               IF(K.EQ.I)GO TO 40
               SGNDET = -SGNDET
               K = L
   30       CONTINUE
         END IF
   40 CONTINUE
      DO 60 I = 1,N
         K = KEEP(N+I)
         IF (W(K).EQ.0) THEN
            DO 50 J = 1,N
               L = KEEP(N+K)
               W(K) = 1
               IF(K.EQ.I)GO TO 60
               SGNDET = -SGNDET
               K = L
   50       CONTINUE
         END IF
   60 CONTINUE
      K2 = 0
      DO 80 JB = 1,NB
         NC = KEEP(NBLOCK+3*JB)
         NR = NC
         IF(NB.EQ.1) NR = M
         K1 = K2 + 1
         K2 = K2 + NC
         NP = KEEP(MBLOCK+3*JB)
         IF (NP.LT.0) THEN
            DO 70 J = K1,K2
               PIV = A(KEEP(IPTRD+J))
               IF(PIV.LT.ZERO)THEN
                  SGNDET  = -SGNDET
                  LOGDET  = LOGDET + LOG(-PIV)
               ELSE
                  LOGDET  = LOGDET + LOG(PIV)
               END IF
   70       CONTINUE
         ELSE
            J1 = 1
            IF (JB.GT.1) J1 = KEEP(KBLOCK+3*JB)
            CALL MA51DD(NR,NC,IQB,NP,LA+1-J1-NE,A(NE+J1),IRN(NE+J1),
     +                  KEEP(IPTRL+K1),KEEP(IPTRU+K1),SGNDT,LOGDT,W)
            SGNDET = SGNDET*SGNDT
            LOGDET = LOGDET+LOGDT
            IF (SGNDET.EQ.0) THEN
               LOGDET = 0
               RETURN
            END IF
         END IF
   80 CONTINUE
      END
      SUBROUTINE MA51DD(M,N,IQ,NP,LFACT,FACT,IRNF,IPTRL,
     +                  IPTRU,SGNDET,LOGDET,W)
      INTEGER K,L,M,N,IQ(*),NP,LFACT
      DOUBLE PRECISION FACT(LFACT)
      INTEGER IRNF(LFACT),IPTRL(N),IPTRU(N),SGNDET
      DOUBLE PRECISION LOGDET
      INTEGER W(N)
      DOUBLE PRECISION A
      INTEGER I,IA1,IF1,J,MF,NF
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D+0)
      IF1 = IPTRL(N) + 1
      MF = IRNF(2)
      NF = N - NP
      LOGDET = 0
      SGNDET = 1
      IF (MF.GT.0) THEN
         IF (M.NE.N .OR. MF.NE.NF .OR. IRNF(IF1+MF+NF-1).LT.0) THEN
            SGNDET = 0
            RETURN
         END IF
         CALL MA51XD(MF,FACT(IF1),MF,IRNF(IF1+MF),SGNDET,LOGDET)
      END IF
      DO 10 I = 1,NP
         IA1 = IPTRU(I) + 1
         A = FACT(IA1)
         IF(A.LT.ZERO)THEN
            SGNDET  = -SGNDET
            LOGDET  = LOGDET - LOG(-A)
         ELSE
            LOGDET  = LOGDET - LOG(A)
         END IF
         W(I) = IRNF(IA1)
   10 CONTINUE
      DO 20 I = 1,MF
         W(NP+I) = IRNF(IF1+I-1)
   20 CONTINUE
      DO 40 I = 1,N
         K = W(I)
         IF (K.GT.0) THEN
            DO 30 J = 1,N
               L = W(K)
               W(K) = 0
               IF(K.EQ.I)GO TO 40
               SGNDET = -SGNDET
               K = L
   30       CONTINUE
         END IF
   40 CONTINUE
      IF(IQ(1).LE.0)RETURN
      DO 60 I = 1,N
         K = IQ(I)
         IF (W(K).EQ.0) THEN
            DO 50 J = 1,N
               L = IQ(K)
               W(K) = 1
               IF(K.EQ.I)GO TO 60
               SGNDET = -SGNDET
               K = L
   50       CONTINUE
         END IF
   60 CONTINUE
      END
      SUBROUTINE MA51XD(N,A,LDA,IPIV,SGNDET,LOGDET)
      IMPLICIT NONE
      INTEGER LDA,N
      DOUBLE PRECISION A(LDA,N)
      INTEGER IPIV(N),SGNDET
      DOUBLE PRECISION LOGDET
      INTEGER J
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D+0)
      INTRINSIC LOG
      SGNDET  = 1
      LOGDET  = ZERO
      DO 30 J = 1,N
         IF (IPIV(J).NE.J) SGNDET  = -SGNDET
         IF (A(J,J).LT.ZERO) THEN
            SGNDET  = -SGNDET
            LOGDET  = LOGDET + LOG(-A(J,J))
         ELSE
            LOGDET  = LOGDET + LOG(A(J,J))
         END IF
   30 CONTINUE
      END
      SUBROUTINE MA51YD(M,N,IPIV,RANK,ROWS,COLS)
      INTEGER M, N, IPIV(N), RANK, ROWS(M), COLS(N)
      INTEGER I,J,K
      DO 10 K = 1,N
         IF (IPIV(K).LT.0) GO TO 20
         ROWS(K) = 1
         COLS(K) = 1
   10 CONTINUE
   20 RANK = K - 1
      DO 30 K = RANK+1, M
          ROWS(K) = 0
   30 CONTINUE
      DO 40 K = RANK+1, N
         COLS(K) = 0
   40 CONTINUE
      DO 50 I = RANK + 1,N
         K = -IPIV(I)
         J = COLS(I)
         COLS(I) = COLS(K)
         COLS(K) = J
   50 CONTINUE
      DO 60 I = RANK,1,-1
         K = IPIV(I)
         J = ROWS(I)
         ROWS(I) = ROWS(K)
         ROWS(K) = J
   60 CONTINUE
      END
      SUBROUTINE MA51ZD(M,N,IQ,NP,LFACT,IRNF,IPTRL,
     +                  IPTRU,RANK,ROWS,COLS,W)
      INTEGER M,N,IQ(*),NP,LFACT,IRNF(LFACT),IPTRL(N),IPTRU(N)
      INTEGER RANK,ROWS(M),COLS(N),W(*)
      INTEGER I,IA1,IF1,J,MF,NF
      EXTERNAL MA51YD
      IF1 = IPTRL(N) + 1
      MF = IRNF(2)
      NF = N - NP
      IF (MF.GT.0 .AND. NF.GT.0) THEN
         CALL MA51YD(MF,NF,IRNF(IF1+MF),RANK,W,COLS(NP+1))
      ELSE
         RANK = 0
         DO 10 I = 1,MF
            W(I) = 0
   10    CONTINUE
         DO 20 I = 1,NF
            COLS(NP+I) = 0
   20    CONTINUE
      END IF
      DO 30 I = 1,M
         ROWS(I) = 1
   30 CONTINUE
      DO 40 I = MF,1,-1
         J = IRNF(IF1+I-1)
         ROWS(J) = W(I)
   40 CONTINUE
      RANK = RANK + M - MF
      DO 220 J = NP,1,-1
         IA1 = IPTRU(J)
         IF (IA1.GE.IPTRL(J)) THEN
            COLS(J) = 0
         ELSE
            COLS(J) = 1
         END IF
  220 CONTINUE
      IF (IQ(1).GT.0) THEN
         DO 230 I = 1,N
            W(I) = COLS(I)
  230    CONTINUE
         DO 240 I = 1,N
            COLS(IQ(I)) = W(I)
  240    CONTINUE
      END IF
      END

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
The YellowBilled Magpie (Pica nuttalli) live only in the Central Valley of California. Because of its rarity some birders travel from around the country just to see our special magpie. The highly social Birds and Bird Watching birds prefer to build nests high up in trees, normally oaks. Magpies are not picky about their diets, eating earthworms, insects, eggs and young from other birds nests, and roadkill. Around Davis and Central California magpies take full advantage of the fresh fruits, grains and berries grown by farmers. In the Weather Fall, acorns are also an important food source.

In the Arboretum and around Davis you can easily spot magpie nests high up in the trees.  They are huge and domeshaped, constructed from loosely thatched sticks, mud, and grass.  In late February the birds lay six to eight eggs in these nests.

There are a lot of these hanging around Olive Tree Lane/on Russel by Arlington. The walnut trees, random roadkill, and freshly tilled land makes for an excellent meal I bet.

Magpies are a close relative of Crow Migration crows and WesternScrub Jays jays and share their gregarious behavior with them.

A recent Sacramento Bee article (http://www.sacbee.com/101/story/840247.html Magpies threatened by West Nile virus  SacBee, April 6, 2008) revealed that the magpies are particularly vulnerable to the West Nile Virus. One of the interviewees speculates that if the birds continue to die from the virus at current rates, they could be extinct in 34 years.

If you see a dead bird, please report it here: http://townsend.ucdavis.edu/reportadeadbird/

Pictures


text:
module solvers
use types, only: dp
use lapack, only: dlamch, dsygvx, dggev, dsysv, dsygvd
use utils, only: stop_error
use sorting, only: sortpairs
implicit none
private
public solve_eig_range, spurious, filter_spurious_states, count_nodes, &
    solve_eig_range2, solve_sym, solve_eig_range3

CONTAINS

subroutine solve_eig_range(Am, Bm, Emin, Emax, lam, c)
! solves system equations
! simple dense solver for now, sparse solver better suited for low-order FE/SE ...
real(dp), intent(in) :: Am(:,:)   ! LHS matrix: Am c = lam Bm c
real(dp), intent(in) :: Bm(:,:)   ! RHS matrix: Am c = lam Bm c
real(dp), intent(in) :: Emin, Emax
real(dp), allocatable, intent(out) :: lam(:)   ! eigenvalues: Am c = lam Bm c
real(dp), allocatable, intent(out) :: c(:,:)   ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
integer i, n
! lapack variables
integer lwork, info, m
integer, allocatable:: iwork(:), ifail(:)
real(dp) abstol
real(dp), allocatable:: Amt(:,:), Bmt(:,:), lamt(:), work(:)
real(dp), allocatable:: z(:,:)

! solve
n=size(Am,1)
lwork=8*n
allocate(Amt(n,n),Bmt(n,n),lamt(n),work(lwork),iwork(5*n),ifail(n),z(n,n))
Amt=Am; Bmt=Bm  ! Amt,Bmt temporaries overwritten by dsygvx
abstol=2*dlamch('S')
call dsygvx(1,'V','V','L',n,Amt,n,Bmt,n,Emin,Emax,0,0,abstol,m, &
    lamt,z,n,work,lwork,iwork,ifail,info)
if (info/=0) then
    print *, "dsygvx returned info =", info
    if (info <= n) then
        print *, info, "eigenvectors failed to converge. Their indices are:"
        print *, ifail
    else
        print *, "The leading minor of order ", info-n, &
            "of B is not positive definite. The factorization of B could ", &
            "not be completed and no eigenvalues or eigenvectors were computed."
    end if
    call stop_error('DSYGVX ERROR')
end if
! Copy the found eigenvalues/vectors into lam/c
allocate(lam(m), c(n, m))
lam(:) = lamt(:m)
c(:, :) = z(:, :m)
do i = 1, m ! normalize such that initial values positive
   if (c(2,i)<0) c(:,i)=-c(:,i)
end do
!write(*,*) "Eigenvector residuals ||Am c - lam Bm c||: "
!allocate(r(n))
!do i = 1, m
!   r=matmul(Am-lam(i)*Bm,c(:,i))
!   write(*,'(1x,i5,a,es18.11)') i, ": ", sqrt(dot_product(r,r))
!   print *, dot_product(c(:, i), matmul(Bm, c(:, i)))
!end do
end subroutine

subroutine solve_eig_range2(Am, Bm, Emin, Emax, lam, c)
! solves generalized nonsymmetric eigenvalue problem
real(dp), intent(in) :: Am(:,:)   ! LHS matrix: Am c = lam Bm c
real(dp), intent(in) :: Bm(:,:)   ! RHS matrix: Am c = lam Bm c
real(dp), intent(in) :: Emin, Emax
real(dp), allocatable, intent(out) :: lam(:)   ! eigenvalues: Am c = lam Bm c
real(dp), allocatable, intent(out) :: c(:,:)   ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
integer i, j, n
! lapack variables
integer lwork, info
real(dp), allocatable:: Amt(:,:), Bmt(:,:), work(:)
real(dp), allocatable:: r(:)
real(dp), allocatable :: vl(:, :)
real(dp), allocatable :: vr(:, :), alphar(:), alphai(:), beta(:)
real(dp), allocatable :: lambda(:)
integer, allocatable :: sel(:)
!complex(dp), parameter :: J = (0, 1)

! solve
n=size(Am,1)
lwork=8*n
allocate(Amt(n,n), Bmt(n,n), work(lwork), vr(n, n), alphar(n), alphai(n), &
    beta(n), vl(n, n))
Amt=Am; Bmt=Bm  ! Amt,Bmt temporaries overwritten by dsygvx
call dggev('N', 'V', n, Amt, n, Bmt, n, alphar, alphai, beta, vl, n, &
    vr, n, work, lwork, info)
if (info/=0) then
    print *, "dsygvx returned info =", info
    if (info <= n) then
        print *, "The", info, " iteration failed."
    else
        print *, "Other error"
    end if
    call stop_error('DGGEV ERROR')
end if
allocate(lambda(n), sel(n))
where (abs(beta) > 1e-16_dp) lambda = alphar / beta
j = 0
do i = 1, n
    sel(i) = 0
    if (abs(alphai(i)) < tiny(1._dp)) then
        if (abs(beta(i)) < tiny(1._dp)) then
            ! Maybe we can just ignore the beta==0 cases
            call stop_error("Got beta==0")
        end if
        if (lambda(i) > Emin .and. lambda(i) < Emax) then
            j = j + 1
            sel(i) = j
        end if
    end if
end do
allocate(lam(maxval(sel)), c(n, maxval(sel)))
do i = 1, size(lambda)
    if (sel(i) /= 0) then
        lam(sel(i)) = lambda(i)
        c(:, sel(i)) = vr(:, i)
    end if
end do
call sortpairs(lam, c)
do i = 1, size(lam) ! normalize such that initial values positive
   if (c(2,i)<0) c(:,i)=-c(:,i)
end do
! Normalize the eigenvectors:
do i = 1, size(lam)
   c(:, i) = c(:, i) / sqrt(dot_product(c(:, i), matmul(Bm, c(:, i))))
end do
write(*,*) "Eigenvector residuals ||Am c - lam Bm c||: "
allocate(r(n))
do i = 1, size(lam)
   r=matmul(Am-lam(i)*Bm,c(:,i))
   write(*,'(1x,i5,a,es18.11)') i, ": ", sqrt(dot_product(r,r))
   !print *, dot_product(c(:, i), matmul(Bm, c(:, i)))
end do
end subroutine

subroutine solve_eig_range3(Am, Bm, Emin, Emax, lam, c)
! solves system equations
! simple dense solver for now, sparse solver better suited for low-order FE/SE ...
real(dp), intent(in) :: Am(:,:)   ! LHS matrix: Am c = lam Bm c
real(dp), intent(in) :: Bm(:,:)   ! RHS matrix: Am c = lam Bm c
real(dp), intent(in) :: Emin, Emax
real(dp), allocatable, intent(out) :: lam(:)   ! eigenvalues: Am c = lam Bm c
real(dp), allocatable, intent(out) :: c(:,:)   ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
integer i, n
! lapack variables
integer lwork, liwork, info, m, m1, m2
integer, allocatable:: iwork(:)
real(dp), allocatable:: Amt(:,:), Bmt(:,:), lamt(:), work(:)
real(dp), allocatable:: r(:)

! solve
n=size(Am,1)
lwork=1+6*n+2*n**2
liwork=3+5*n
allocate(Amt(n,n),Bmt(n,n),lamt(n),work(lwork),iwork(liwork))
Amt=Am; Bmt=Bm  ! Amt,Bmt temporaries overwritten by dsygvx
call dsygvd(1,'V','L',n,Amt,n,Bmt,n,lamt,work,lwork,iwork,liwork,info)
if (info/=0) then
    print *, "dsygvd returned INFO =", info
    if (info < 0) then
        print *, "the", -info, "-th argument had an illegal value"
    else if (info <= n) then
        print *, "the algorithm failed to compute an eigenvalue while working"
        print *, "on the submatrix lying in rows and columns", 1.0_dp*info/(n+1)
        print *, "through", mod(info, n+1)
    else
        print *, "The leading minor of order ", info-n, &
            "of B is not positive definite. The factorization of B could ", &
            "not be completed and no eigenvalues or eigenvectors were computed."
    end if
    call stop_error('DSYGVD ERROR')
end if
! TODO: make sure things work with regards m1==n and m2==n:
m1 = 1
do while (lamt(m1) < Emin .and. m1 < n)
    m1 = m1 + 1
end do
! m1 is the first eigenvalue in the range
m2 = m1
do while (lamt(m2) <= Emax .and. m2 < n)
    m2 = m2 + 1
end do
m2 = m2 - 1
! m2 is the last eigenvalue in the range
m = m2-m1+1 ! number of eigenvalues in the range
! Copy the found eigenvalues/vectors into lam/c
allocate(lam(m), c(n, m))
lam(:) = lamt(m1:m2)
c(:, :) = Amt(:, m1:m2)
do i = 1, m ! normalize such that initial values positive
   if (c(2,i)<0) c(:,i)=-c(:,i)
end do
write(*,*) "Eigenvector residuals ||Am c - lam Bm c||: "
allocate(r(n))
do i = 1, m
   r=matmul(Am-lam(i)*Bm,c(:,i))
   write(*,'(1x,i5,a,es18.11)') i, ": ", sqrt(dot_product(r,r))
   print *, dot_product(c(:, i), matmul(Bm, c(:, i)))
end do
end subroutine

integer function count_nodes(u) result(nn)
! Returns the number of nodes of the function
real(dp), intent(in) :: u(:) ! State "u"

integer :: i
logical :: s, last_s, ok
real(dp) :: thres
thres = maxval(abs(u)) * 1e-5_dp
nn = 0
ok = .false.
last_s = .false. ! To fix compiler uninitialized variable warning
do i = 1, size(u)
    ! skip small "tails" of the wavefunction whose oscillations we don't want to
    ! count
    if (abs(u(i)) < thres) cycle
    s = u(i) > 0
    if (ok) then
        if (s .neqv. last_s) nn = nn + 1
    end if
    last_s = s
    ok = .true.
end do
end function

logical function spurious_old2(u)
! Returns .true. if the state "u" is spurious, otherwise .false.
! The state is declared spurious if it has (on average) less nodes per
! oscillation than allowed (see the nn_per_oscillation below), determined by
! counting the number of oscillations and comparing to the length of the vector.
real(dp), intent(in) :: u(:) ! State "u"
! The minimum allowed nodes per oscillation
integer, parameter :: nn_per_oscillation = 4

real(dp) :: osc(size(u)-1)
integer :: i, nn
logical :: s, last_s
real(dp) :: thres
integer :: last_oscillation
osc(:) = u(2:) - u(:size(u)-1) ! Calculate differences
thres = (maxval(u) - minval(u)) * 1e-5_dp ! Determine treshold
!print *, "threshold:", thres
nn = 0
last_s = osc(1) > 0
!print *, "u:"
!print *, u
!print *, "osc:"
!print *, osc
last_oscillation = 1
do i = 1, size(osc)
    ! skip small "tails" of the wavefunction whose oscillations we don't want to
    ! count
    if (abs(osc(i)) < thres) cycle
    s = osc(i) > 0
    if (s .neqv. last_s) then
        nn = nn + 1
        !if (i==99) then
        !    print *, u(i-3:)
        !end if
!        print *, "found oscillation at:", i, size(u), &
!            abs(u(i) - u(last_oscillation))
        if (i - last_oscillation < nn_per_oscillation) then
            spurious_old2 = .true.
            return
        end if
        last_oscillation = i
    end if
    last_s = s
end do
if (nn > size(u) / nn_per_oscillation) then
    spurious_old2 = .true.
else
    spurious_old2 = .false.
end if
!print *, "oscillations:", nn, size(u) / nn_per_oscillation
!print *, maxval(abs(u)), u(:10)
!print *, thres, osc(:10)
!spurious = .false.
end function

logical function spurious(u, ib)
! Returns .true. if the state "u" is spurious, otherwise .false.
! The state is declared spurious if it has (on average) less nodes per
! oscillation than allowed (see the nn_per_oscillation below), determined by
! counting the number of oscillations and comparing to the length of the vector.
real(dp), intent(in) :: u(:) ! State "u"
integer, intent(in) :: ib(:, :)   ! basis connectivity: ib(i,j) = index of basis
   ! function associated with local basis function i of element j. 0 = no
   ! associated basis fn.
! The minimum allowed nodes per oscillation
!integer, parameter :: nn_per_oscillation = 4

real(dp) :: osc(size(u) - 1)
integer :: i, ie, ip, nn
logical :: s, last_s
real(dp) :: thres
integer :: last_oscillation
osc(:) = u(2:) - u(:size(u)-1) ! Calculate differences
thres = (maxval(u) - minval(u)) * 1e-5_dp ! Determine treshold
!print *, "threshold:", thres
nn = 0
last_s = osc(1) > 0
!print *, "u:"
!print *, u
!print *, "osc:"
!print *, osc
if (size(ib, 1)-1 == 1) call stop_error("spurious check doesn't work for p=1")
last_oscillation = 1
do ie = 1, size(ib, 2)
    do ip = 1, size(ib, 1)
        i = ib(ip, ie)
        if (i == 0) cycle
        if (i < 0) cycle
        if (i > size(osc)) exit
        s = osc(i) > 0
        if (s .neqv. last_s) then
            if (abs(u(i)-u(last_oscillation)) < thres) then
                last_oscillation = i
                last_s = s
                cycle
            end if
            nn = nn + 1
            last_oscillation = i
        end if
        last_s = s
    end do
    !print *, "nn = ", nn, size(ib, 1), (size(ib, 1)-1)/3 + 1
    if (nn > (size(ib, 1)-1)/3 + 1) then
        print *, "spurious, nodes:", nn, "element:", ie
        spurious = .true.
        return
    end if
    nn = 0
end do
spurious = .false.
end function

subroutine filter_spurious_states(ib, lam, u, lam2, u2, skipped)
integer, intent(in) :: ib(:,:,:)   ! basis connectivity: ib(i,j) = index of basis
   ! function associated with local basis function i of element j. 0 = no
   ! associated basis fn.
real(dp), intent(in) :: u(:, :), lam(:)
real(dp), intent(out), allocatable :: u2(:, :), lam2(:)
logical, intent(out) :: skipped
integer :: m(size(lam))
integer :: i, j, c, c_begin, c_end, n_eig
logical :: sp
j = 0
print *, "Filtering high oscillatory states:"
do i = 1, size(lam)
    ! Loop over the vector components:
    sp = .false.
    c_begin = 1
    do c = 1, size(ib, 3)
        c_end = maxval(ib(:, :, c))
        if (spurious(u(c_begin:c_end, i), ib(:, :, c))) then
            sp = .true.
            exit
        end if
        c_begin = c_end+1
    end do
    if (sp) then
        m(i) = 0
    else
        j = j + 1
        m(i) = j
    end if
    print "(i4,i4,f12.5)", i, m(i), lam(i)
end do
n_eig = maxval(m)
allocate(u2(size(u, 1), n_eig), lam2(n_eig))
do i = 1, size(lam)
    if (m(i) /= 0) then
        u2(:, m(i)) = u(:, i)
        lam2(m(i)) = lam(i)
    end if
end do
! If we skipped any state (we filtered some in between), then these two numbers
! will mismatch:
!skipped = m(n_eig) /= n_eig
! But we actually want to see, whether we filtered anything at all:
if (size(m) == 0) then
    skipped = .false.
else
    skipped = m(size(m)) /= size(m)
end if
end subroutine

function solve_sym(Am, bv) result(c)
! solves symmetric dense system of equations
real(dp), intent(in) :: Am(:,:)   ! system matrix: Am c = bv
real(dp), intent(in) :: bv(:)     ! source vector: Am c = bv
real(dp) :: c(size(bv))     ! solution vector: Am c = bv
!real(dp) :: r(size(bv))
integer :: n
! lapack variables
integer :: lwork, info
integer, allocatable :: ipiv(:)
real(dp), allocatable :: Amt(:,:),bm(:,:),work(:)

n = size(c)
lwork = n
allocate(Amt(n,n), bm(n,1), ipiv(n), work(lwork))
Amt=Am; bm(:,1)=bv   ! temporaries for dsysv
call dsysv('L', n, 1, Amt, n, ipiv, bm, n, work, lwork, info)
if (info < 0) then
    print *, "The", -info, "-th argument had illegal value"
    call stop_error('DSYSV ERROR.')
end if
if (info > 0) then
    print *, "D(", info, ",", info, ") is exactly zero."
    print *, "The factorization has been completed, but the block diagonal"
    print *, "matrix D is exactly singular, so the solution could not be"
    print *, "computed."
    call stop_error('DSYSV ERROR.')
end if
c = bm(:, 1)
! error
!r=matmul(Am, c)-bv
!write(*,'(1x,a,es18.11)') "Solution vector residual ||Am c - bv||/||bv||: ", &
!   sqrt(dot_product(r,r)/dot_product(bv,bv))
end function

end module

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
C    %W% %G%
      subroutine int3fd(a)
C
C  THIS SUBROUTINE CALCULATES INITIAL STATE VECTORS AND PAST VALUES
C  FOR EXCITER TYPE FD
C
      include 'tspinc/param.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/params.inc'
      include 'tspinc/gentbla.inc'
      include 'tspinc/pointr.inc'
      include 'tspinc/gentblb.inc'
      include 'tspinc/int3.inc'
C     -  Local variables
      dimension array(40),a(40)
       equivalence        (array(1), vrmax),(array(2), vrmin),
     1 (array(3),  vref), (array(4),efdmax),(array(5),   cki),
     2 (array(6),   hbr), (array(7),  cckp),(array(8),   cka),
     3 (array(9),    dr), (array(10),  ckc),(array(11), hbf),
     4 (array(12),   df),(array(13),dk),(array(14),    de),
     5 (array(15),  cke), (array(24),   ckpi),
     6 (array(16),   vao),(array(17),   vbo),(array(18),   vco),
     7 (array(19), vino),(array(20),    xl),(array(21),  rcex),
     8 (array(22),   xcex),(array(23),  ckpr),(array(25),    da),
     9 (array(26),   hba),(array(27),   hbe),(array(28), theta),
     * (array(29),   efd),(array(30),   a1),(array(31),     a2),
     * (array(32),    a3),(array(33),   vio),(array(34),  bnf)
      logical debug                                                     !dem
C     -     Begin     Begin     Begin     Begin     Begin     Begin
      debug = .false.                                                   !dem
      l6 = 6
C
C TRANSFER CITER DATA TO ARRAY
C
      do 100 itr=1,40
100   array(itr) = a(itr)

C
CCC CK FOR LIMIT VIOLATION
C
          if(efdo .gt. efdmax) then
          write(errbuf,2000) name,base,id,efdo
 2000     format(1h0,'  INITIAL FIELD VOLTS VIOLATES LIMITS FOR ',
     1 a8,2x,f5.1,2x,a1,5x,' INITIAL VOLTS = ',f6.3)
          call prterr('E',1)
          iabort = 1
          go to 2300
          endif
          efd = efdo
C
CCCC CALCULATE VOLTAGE TRANSDUCER CKT
C
 2031     if(xt .gt. 0.) go to 2033
          vtr = vhd
          vti = vhq
          go to 2034
 2033     vtr = vhd + rt*agend - xt*agenq
          vti = vhq + rt*agenq + xt*agend
 2034     crt = agend
          cit = agenq
          zcr = rcex*crt - xcex*cit
          zci = rcex*cit + xcex*crt
          vcr = zcr + vtr
          vci = zci + vti
          vc = sqrt(vcr**2 + vci**2)
          vco = vc
C
C R BLOCK
C
          dr = 2.0*dr/dt + 1.0
          hbr = (dr - 2.0)*vc + vc
C
CCC CALCULATE THEVENIN VOLTS
C
          oir = agend
          oii = agenq
          ckiir = cki*oir
          ckiii = cki*oii
          ckvr = cckp*vtr
          ckvi = cckp*vti
          verr = ckvr-ckiii
          veri = ckvi + ckiir
          ve = sqrt(verr**2 + veri**2)
C
C INITIAL FIELD CURRENT,CFD, EQUALS INITIAL FIELD VOLTS,EFDO.
C
          cfd = efdo
          cin = ckc*cfd/ve
          if(cin .le. 0.51) fex = 1.0 - 0.58*cin
          if(0.51 .lt. cin .and. cin .lt. 0.715)
     1    fex = -0.865*(cin + 0.00826)**2 + 0.93233
          if(cin .ge. 0.715)fex = 1.68 - 1.714*cin
          if(fex .lt. 0.0) fex = 0.0
          vb = fex*ve
          vbo = vb
C
C E BLOCK
C
          de = 2.*de/dt
          hbe = efdo*de
C
C A BLOCK
C
          vr = cke*efdo-vb
          x3 = vr
          da = 2.0*da/dt + 1.0
          va = vr/cka
          vao = va
          hba = vr*(da - 2.0) + cka*va
          if(vr .lt. vrmax .and. vr .gt. vrmin) go to 2038
          write (errbuf, 2036) name, base, id, vr
          call prterr ('E',1)
 2036     format(1h0, ' INITIAL VR VIOLATES LIMITS FOR ', a8,
     1           2x, f5.1, 2x, a1, 5x, 'INIT VR = ', f6.3)
          iabort = 1
          go to 2300
 2038     vref = va + vc
C
C F BLOCK
C
      vf = 0.0
      dk = 2.*dk/dt
      df = 2.*df/dt + 1.
      hbf = -dk*efdo
      vino = vr + vb - cke*efdo
c     -
      if (debug) then                                                   !dem
        call dbgeko ('INT3FD - at end')                                 !dem
        call dbgwrf ('  VC /Vterm xducer/ = ',vc)                       !dem
        call dbgwrf ('  VE /lin drop comp/ = ',ve)                      !dem
        call dbgwrf ('  EFDO /field volts/ = ',efdo)                    !dem
        call dbgwrf ('  CFD /field curr/ = ',cfd)                       !dem
        call dbgwrf ('  FEX /V_e multplr/ = ',fex)                      !dem
        call dbgwrf ('  FEX /V_e multplr/ = ',fex)                      !dem
        call dbgwrf ('  VB /source blk out/ = ',vb)                     !dem
        call dbgwrf ('  VF /regultr feedbk/ = ',vf)                     !dem
        call dbgwrf ('  VA /regultr input/ = ',va)                      !dem
        call dbgwrf ('  VR /regultr output/ = ',vr)                     !dem
        rvef = cke * efdo                                               !dem
        call dbgwrf ('  RVEF /exciter feedbk/ = ',rvef)                 !dem
        call dbgwrf ('  VIN /exc integrtr in/ = ',vino)                 !dem
      endif                                                             !dem
C
C TRANSFER DATA TO CITER AND RETURN TO INITAL3
C
 2300     do 2400 itr = 1,40
 2400     a(itr) = array(itr)
          return
          end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
      subroutine bspp2d ( t, bcoef, n, k, m, scrtch, break, coef, l )
c  from  * a practical guide to splines *  by c. de boor
calls  bsplvb
c  this is an extended version of  bsplpp  for use with tensor products
c
converts the b-representation  t, bcoef(.,j), n, k  of some spline into
c  its pp-representation  break, coef(j,.,.), l, k , j=1, ..., m  .
c
c******  i n p u t  ******
c  t.....knot sequence, of length  n+k
c  bcoef(.,j) b-spline coefficient sequence, of length	n ,j=1,...,m
c  n.....length of  bcoef  and	dimension of spline space  spline(k,t)
c  k.....order of the spline
c
c  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)	is impo-
c	 sed by the arbitrary dimension statement for  biatx  below, but
c	 is  n o w h e r e   c h e c k e d   for.
c
c  m	 number of data sets
c
c******  w o r k   a r e a  ******
c  scrtch   of size  (k,k,m), needed to contain bcoeffs of a piece of
c	 the spline and its  k-1  derivatives	for each of the m sets
c
c******  o u t p u t  ******
c  break.....breakpoint sequence, of length  l+1, contains (in increas-

c	 ing order) the distinct points in the sequence  t(k),...,t(n+1)
c  coef(mm,.,.)  array of size (k,n), with  coef(mm,i,j) = (i-1)st der-
c	 ivative of  mm-th  spline at break(j) from the right, mm=1,.,m
c  l.....number of polynomial pieces which make up the spline in the in-
c	 terval  (t(k), t(n+1))
c
c******  m e t h o d  ******
c     for each breakpoint interval, the  k  relevant b-coeffs of the
c  spline are found and then differenced repeatedly to get the b-coeffs

c  of all the derivatives of the spline on that interval. the spline and
c  its first  k-1  derivatives are then evaluated at the left end point

c  of that interval, using  bsplvb  repeatedly to obtain the values of
c  all b-splines of the appropriate order at that point.
c
c     parameter kmax = 20
      integer k,l,m,n,	 i,j,jp1,kmj,left,lsofar,mm
      real bcoef(n,m),break(1),coef(m,k,1),t(1),   scrtch(k,k,m)
     *					      ,biatx(20),diff,fkmj,sum
c
c    *					      ,biatx(kmax),diff,fkmj,sum
c     dimension break(l+1),coef(k,l),t(n+k)
current fortran standard makes it impossible to specify the length of
c  break , coef  and  t  precisely without the introduction of otherwise
c  superfluous additional arguments.
      lsofar = 0
      break(1) = t(k)
      do 50 left=k,n
c				 find the next nontrivial knot interval.
	 if (t(left+1) .eq. t(left))	go to 50
	 lsofar = lsofar + 1
	 break(lsofar+1) = t(left+1)
	 if (k .gt. 1)			go to 9
	 do 5 mm=1,m
    5	    coef(mm,1,lsofar) = bcoef(left,mm)
					go to 50
c	 store the k b-spline coeff.s relevant to current knot interval

c			      in  scrtch(.,1) .
    9	 do 10 i=1,k
	    do 10 mm=1,m
   10	       scrtch(i,1,mm) = bcoef(left-k+i,mm)
c
c	 for j=1,...,k-1, compute the  k-j  b-spline coeff.s relevant to
c	 current knot interval for the j-th derivative by differencing
c	 those for the (j-1)st derivative, and store in scrtch(.,j+1) .

	 do 20 jp1=2,k
	    j = jp1 - 1
	    kmj = k - j
	    fkmj = float(kmj)
	    do 20 i=1,kmj
	       diff = (t(left+i) - t(left+i - kmj))/fkmj
	       if (diff .le. 0.)	 go to 20
	       do 15 mm=1,m
   15		  scrtch(i,jp1,mm) =
     *		  (scrtch(i+1,j,mm) - scrtch(i,j,mm))/diff
   20	       continue
c
c	 for  j = 0, ..., k-1, find the values at  t(left)  of the  j+1

c	 b-splines of order  j+1  whose support contains the current
c	 knot interval from those of order  j  (in  biatx ), then comb-

c	 ine with the b-spline coeff.s (in scrtch(.,k-j) ) found earlier
c	 to compute the (k-j-1)st derivative at  t(left)  of the given
c	 spline.
c	    note. if the repeated calls to  bsplvb  are thought to gene-
c	 rate too much overhead, then replace the first call by
c	    biatx(1) = 1.
c	 and the subsequent call by the statement
c	    j = jp1 - 1
c	 followed by a direct copy of the lines
c	    deltar(j) = t(left+j) - x
c		   ......
c	    biatx(j+1) = saved
c	 from  bsplvb . deltal(kmax)  and  deltar(kmax)  would have to
c	 appear in a dimension statement, of course.
c
	 call bsplvb ( t, 1, 1, t(left), left, biatx )
	 do 25 mm=1,m
   25	    coef(mm,k,lsofar) = scrtch(1,k,mm)
	 do 30 jp1=2,k
	    call bsplvb ( t, jp1, 2, t(left), left, biatx )
	    kmj = k+1 - jp1
	    do 30 mm=1,m
	       sum = 0.
	       do 28 i=1,jp1
   28		  sum = biatx(i)*scrtch(i,kmj,mm) + sum
   30	       coef(mm,kmj,lsofar) = sum
   50	 continue
      l = lsofar
					return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!    Date:       22 Oct 1992
!    Time:       22:01
!    Program:    LP.FOR
!    Version:    6.00.00
!    Programmer: Nicolaas M de Rooij
!
      subroutine lp(mon)
!           4-69
!           sets up a linear programming problem and calls subroutine
!        simple to find a feasible solution.
!           lp uses subprograms
!              bar, leave, simple
!           lp uses double precision version of
!              alog, amin1, float
      include  'char1.inc'
      equivalence (cc(1),xmf(1)),(xx(1),x2(1)),(p(1),v1(1)),
     2   (jh(1),v2(1))
      dimension xx(1), kout(7), cc(1), p(1), jh(1)
      mon = 0
      if (xstart.le.0.0) xstart = 1.e-6
!        clear v4 storage for new column of matrix
      do 10 i=1,m
         p(i) = b(i)
         v4(i) = 0.0
   10 continue
      do 20 j=1,naij
         i = irow(j)
         v4(i) = v4(i) + aij(j)
   20 continue
!        test for space for new column
      if (naij+m.le.maxaij) go to 30
!        no room for new column
      write (not,99999) naij
      call leave(ilp, not)
      return
!        move summation column into matrix area
   30 j = ntot + 1
      do 40 i=1,m
         ii = naij + i
         irow(ii) = i
         jcol(ii) = j
         aij(ii) = v4(i)
   40 continue
      naij = naij + m
      do 50 j=1,ntot
         cc(j) = 0.0
   50 continue
      cc(n+1) = -1.0
!        zero-th simplex is to determine feasibility
      call simple(0, m, n+1, naij, aij, irow, jcol, p, cc, kout, xx,
     1 pie, jh, v3, v4, x3, r)
      zt = xx(n+1)
      if (pf.ge.0) write (not,99998) kout(2), zt, kout(1)
      zzt = dmin1(zt/2.0,xstart)
      k = naij
   60 if (jcol(k).ne.n+1) go to 70
      i = irow(k)
      p(i) = p(i) - zzt*aij(k)
      k = k - 1
      go to 60
!        restore element count to original value
   70 naij = naij - m
      do 80 j=1,ntot
         x(j) = xx(j)
         xmf(j) = 1.0
         x1(j) = 0.0
   80 continue
      if (zt.le.0. .or. kout(1).ne.0) go to 160
!        simplex loop
      fr2 = 1.e+20
      do 130 nn=1,ncomp
         do 90 j=1,ntot
            cc(j) = c(j) + xmf(j) - 1.0
   90    continue
         fn = dble(nn) - 1.0
!     write(not,1234)nn
!1234  format(' loop in lp, simplex: ',i4)
         call simple(1, m, n, naij, aij, irow, jcol, p, cc, kout, xx,
     1    pie, jh, v3, v4, x3, r)
!     write(not,1235)kout(1)
!1235  format(' loop in lp, kout(1): ',i4)
         if (kout(1).ne.0) go to 240
         do 110 j=1,ntot
            x(j) = xx(j)
            x(j) = (fn*x1(j)+x(j))/(fn+1.0)
            x1(j) = x(j)
  110    continue
         call barx(x, xbar)
         k = 1
         fr = 0.0
         do 120 j=1,n
            if (j.ge.kl(k+1)) k = k + 1
            if (j.eq.kl(k) .and. xbar(k).gt.0.0) fr = fr -
     1       xbar(k)*dlog(xbar(k))
            if (x(j).gt.0.0) fr = fr + x(j)*(dlog(x(j))+c(j))
            xmf(j) = 0.
            if (xbar(k).ne.0.) xmf(j) = x(j)/xbar(k)
  120    continue
         if (pf.ge.0) write (not,99997) nn, kout(2), fr
         if (fr.ge.fr2) go to 140
         fr2 = fr
  130 continue
  140 do 150 j=1,n
         x(j) = x(j) + zzt
  150 continue
      return
  160 if (kout(1).gt.1) go to 240
      write (not,99996)
      do 170 i=1,m
         if (pie(i).ne.0.) write (not,99995) pie(i), nr(i,1), nr(i,2)
  170 continue
      write (not,99994)
      ii = 1
      do 210 k=1,ncomp
         mta = kl(k)
         mtb = kl(k+1) - 1
         do 200 j=mta,mtb
            d = 0.
  180       if (jcol(ii).ne.j) go to 190
            i = irow(ii)
            d = pie(i)*aij(ii) + d
            ii = ii + 1
            go to 180
  190       if (d.ne.0.) write (not,99993) d, kn(j), nam(k,1), nam(k,2)
  200    continue
  210 continue
      d = 0.
      do 220 i=1,m
         d = pie(i)*b(i) + d
  220 continue
      write (not,99992) d
  230 mon = 1
      return
  240 if (kout(1).ne.2) go to 250
!           simplex declared an unbounded solution.  species number is
!        in kout(6).  lp will continue with current solution.  no
!        message output.
!     go to 100
!     write(not,99990)
!99990 format(' lp: the microsoft compiler could not handle this, apologies')
      mon = -kout(6)
      write(not,1001)kn(-mon)
1001  format(' simplex declared an unbounded solution for species: ',a6)
      return
  250 write (not,99991)
      go to 230
99999 format (' matrix storage has been exceeded in lp subroutine, wh',
     1 'en new species created for first simplex.  total entries =',i6/
     2 1h )
99998 format (' simplex  0,', i4, ' iterations, max min element=',
     1 1pe15.8, ', condition ', i3)
99997 format (' simplex', i3, 1h,, i4, ' iterations ', ' fr eng=',
     1 1pe15.8)
99996 format (' this problem is infeasible.  the following linear com',
     1 'bination of components'/1x)
99995 format (10x, 3h+ (, f15.8, 5h ) * , 2a6)
99994 format ('   leads to the following infeasible equation,  '/1x)
99993 format (10x, 3h+ (, f15.8, 5h ) * , a6, ' in ', 2a6)
99992 format (' ', 15x, 7h+ 0.0 =, f15.8)
99991 format (' simplex routine has failed due to excessive round-off',
     1 ' error')
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
MODULE PublicPar
!
! Stress histogram parameters
   DOUBLE PRECISION, ALLOCATABLE :: Spect(:,:,:) ! Stress histograms
   INTEGER                       :: nccmax       ! size of 2nd column of stress histograms (number of rows in stress range input files)
   DOUBLE PRECISION              :: R            ! Stress ratio
   INTEGER                       :: nmaxCsr      ! Maximum number of stress range uncertainty values
!
! Inspection parameters
   INTEGER                       :: ninsp        ! number of inspections
   INTEGER                       :: ttd          ! 0 = surface flaws can be detected
					    		                 ! 1 = only through thickness flaws can be detected
   INTEGER, ALLOCATABLE          :: nd(:)        ! number of detected cracks during each inspection
!
! Number of detail locations
   INTEGER                       :: nloc         ! number of potential crack locations
!
! Material parameters
   DOUBLE PRECISION              :: fynom        ! nominal yield stress
   INTEGER                       :: law          ! Crack growth equation
                                                 ! 1 = Forman Mettu
												 ! 2 = BS 7910:2013 in air
!
! Geometry parameters
   INTEGER                       :: LoadType     ! 1 = Membrane loading
                                                 ! 2 = Bending loading
												 ! 3 = Both membrane and bending
   INTEGER                       :: FlawType     ! 1 = Through-thickness flaws in plates
                                                 ! 2 = Surface flaws in plates
												 ! 3 = Through-width flaws in plates
									             ! 4 = semi-circular surface flaw in round bar
												 ! 5 = Corner flaws in plates
												 ! 6 = Semi-circular surface flaw in bolts, solution 2
												 ! 7 = Edge flaws in plates
												 ! 8 = Through-thickness deck plate crack (orthotropic bridge deck)
												 ! 9 = Cope longitudinal girger (raveling) ProRail
   INTEGER                       :: WeldGeometry ! 0 = No weld
                                                 ! 1 = According to fig. M.23
                                                 ! 2 = According to fig. M.24
   DOUBLE PRECISION              :: Bnom         ! nominal plate thickness
   DOUBLE PRECISION              :: p            ! exponent for the dependency of stress on plate thickness (1 for membrane, 2 for bending)
   DOUBLE PRECISION              :: cmax         ! Semi crack length at which end-of-life is assumed if not already fractured
   DOUBLE PRECISION              :: alpha        ! Constraint constant
!   INTEGER, PARAMETER :: ndof = 6
!
END MODULE PublicPar

max_stars_repo_name: monicanicoreac/ProRail-ravelingen

max_issues_repo_name: monicanicoreac/ProRail-ravelingen

max_forks_repo_name: monicanicoreac/ProRail-ravelingen


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ###############################################################
c     ##                                                           ##
c     ##  subroutine estrbnd3  --  stretch-bend energy & analysis  ##
c     ##                                                           ##
c     ###############################################################
c
c
c     "estrbnd3" calculates the stretch-bend potential energy;
c     also partitions the energy among the atoms
c
c
      subroutine estrbnd3
      implicit none
      include 'sizes.i'
      include 'action.i'
      include 'analyz.i'
      include 'angle.i'
      include 'angpot.i'
      include 'atmtyp.i'
      include 'atoms.i'
      include 'bond.i'
      include 'bound.i'
      include 'energi.i'
      include 'group.i'
      include 'inform.i'
      include 'iounit.i'
      include 'math.i'
      include 'strbnd.i'
      include 'usage.i'
      integer i,j,k,istrbnd
      integer ia,ib,ic
      real*8 e,dr1,dr2,dt
      real*8 fgrp,angle
      real*8 force1,force2
      real*8 dot,cosine
      real*8 xia,yia,zia
      real*8 xib,yib,zib
      real*8 xic,yic,zic
      real*8 xab,yab,zab
      real*8 xcb,ycb,zcb
      real*8 rab,rab2
      real*8 rcb,rcb2
      logical proceed
      logical header,huge
c
c
c     zero out the energy component and partitioning terms
c
      neba = 0
      eba = 0.0d0
      do i = 1, n
         aeba(i) = 0.0d0
      end do
      header = .true.
c
c     calculate the stretch-bend energy term
c
      do istrbnd = 1, nstrbnd
         i = isb(1,istrbnd)
         ia = iang(1,i)
         ib = iang(2,i)
         ic = iang(3,i)
         force1 = sbk(1,istrbnd)
         force2 = sbk(2,istrbnd)
c
c     decide whether to compute the current interaction
c
         proceed = .true.
         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,0,0,0)
         if (proceed)  proceed = (use(ia) .or. use(ib) .or. use(ic))
c
c     get the coordinates of the atoms in the angle
c
         if (proceed) then
            xia = x(ia)
            yia = y(ia)
            zia = z(ia)
            xib = x(ib)
            yib = y(ib)
            zib = z(ib)
            xic = x(ic)
            yic = y(ic)
            zic = z(ic)
c
c     compute the value of the bond angle
c
            xab = xia - xib
            yab = yia - yib
            zab = zia - zib
            xcb = xic - xib
            ycb = yic - yib
            zcb = zic - zib
            if (use_polymer) then
               call image (xab,yab,zab)
               call image (xcb,ycb,zcb)
            end if
            rab2 = xab*xab + yab*yab + zab*zab
            rcb2 = xcb*xcb + ycb*ycb + zcb*zcb
            if (rab2.ne.0.0d0 .and. rcb2.ne.0.0d0) then
               rab = sqrt(rab2)
               rcb = sqrt(rcb2)
               dot = xab*xcb + yab*ycb + zab*zcb
               cosine = dot / (rab*rcb)
               cosine = min(1.0d0,max(-1.0d0,cosine))
               angle = radian * acos(cosine)
               dt = angle - anat(i)
c
c     get the stretch-bend interaction energy
c
               j = isb(2,istrbnd)
               k = isb(3,istrbnd)
               dr1 = rab - bl(j)
               dr2 = rcb - bl(k)
               e = stbnunit * (force1*dr1+force2*dr2) * dt
c
c     scale the interaction based on its group membership
c
               if (use_group)  e = e * fgrp
c
c     increment the total stretch-bend energy
c
               neba = neba + 1
               eba = eba + e
               aeba(ib) = aeba(ib) + e
c
c     print a message if the energy of this interaction is large
c
               huge = (abs(e) .gt. 2.0d0)
               if (debug .or. (verbose.and.huge)) then
                  if (header) then
                     header = .false.
                     write (iout,10)
   10                format (/,' Individual Stretch-Bend',
     &                          ' Interactions :',
     &                       //,' Type',18x,'Atom Names',18x,'dSB 1',
     &                          5x,'dSB 2',6x,'Energy',/)
                  end if
                  write (iout,20)  ia,name(ia),ib,name(ib),
     &                             ic,name(ic),dr1*dt,dr2*dt,e
   20             format (' StrBend',3x,3(i7,'-',a3),2x,2f10.4,f12.4)
               end if
            end if
         end if
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib