-
Notifications
You must be signed in to change notification settings - Fork 26
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Added a vector data type to facilitate modular vector arithmetic in a…
…ll representation.
- Loading branch information
Showing
4 changed files
with
249 additions
and
10 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,90 @@ | ||
module vector_allreps_module | ||
|
||
use precision, only: i64, r64 | ||
use misc, only: mux_vector, demux_vector | ||
|
||
implicit none | ||
|
||
private | ||
public :: vector_allreps, & | ||
vector_allreps_add, vector_allreps_sub, & | ||
vector_allreps_print | ||
|
||
type vector_allreps | ||
!! A container for a (3-)vector and relevant arithmetic operations. | ||
!! It is convenient to use under certain circumstances, for example, | ||
!! when low level vector arithmetic is repetitive and error prone. | ||
!! However, the low level method will be more performant. | ||
|
||
integer(i64) :: muxed_index = -1 | ||
integer(i64) :: int(3) = 0 | ||
real(r64) :: frac(3) = 0.0 | ||
real(r64) :: cart(3) = 0.0 | ||
end type vector_allreps | ||
|
||
interface vector_allreps | ||
module procedure create | ||
end interface vector_allreps | ||
|
||
contains | ||
|
||
function create(ind, grid, primitive_vecs) result(vector_obj) | ||
integer(i64), intent(in) :: ind, grid(3) | ||
real(r64), intent(in) :: primitive_vecs(3, 3) | ||
type(vector_allreps) :: vector_obj | ||
|
||
!vector_obj%grid = grid | ||
|
||
!vector_obj%primitive_vecs = primitive_vecs | ||
|
||
vector_obj%muxed_index = ind | ||
|
||
call demux_vector(ind, vector_obj%int, & | ||
grid, 0_i64) | ||
|
||
vector_obj%frac = dble(vector_obj%int)/grid | ||
|
||
vector_obj%cart = matmul(primitive_vecs, vector_obj%frac) | ||
end function create | ||
|
||
subroutine vector_allreps_print(v) | ||
!! Printer | ||
|
||
type(vector_allreps), intent(in) :: v | ||
|
||
print*, v%muxed_index | ||
print*, v%int | ||
print*, v%frac | ||
print*, v%cart | ||
end subroutine vector_allreps_print | ||
|
||
pure function vector_allreps_add(v1, v2, grid, primitive_vecs) result(v3) | ||
!! Adder | ||
|
||
type(vector_allreps), intent(in) :: v1 | ||
type(vector_allreps), intent(in) :: v2 | ||
integer(i64), intent(in) :: grid(3) | ||
real(r64), intent(in) :: primitive_vecs(3, 3) | ||
type(vector_allreps) :: v3 | ||
|
||
v3%int = modulo(v1%int + v2%int, grid) | ||
v3%muxed_index = mux_vector(v3%int, grid, 0_i64) | ||
v3%frac = dble(v3%int)/grid | ||
v3%cart = matmul(primitive_vecs, v3%frac) | ||
end function vector_allreps_add | ||
|
||
pure function vector_allreps_sub(v1, v2, grid, primitive_vecs) result(v3) | ||
!! Subtracter | ||
|
||
type(vector_allreps), intent(in) :: v1 | ||
type(vector_allreps), intent(in) :: v2 | ||
integer(i64), intent(in) :: grid(3) | ||
real(r64), intent(in) :: primitive_vecs(3, 3) | ||
type(vector_allreps) :: v3 | ||
|
||
v3%int = modulo(v1%int - v2%int, grid) | ||
v3%muxed_index = mux_vector(v3%int, grid, 0_i64) | ||
v3%frac = dble(v3%int)/grid | ||
v3%cart = matmul(primitive_vecs, v3%frac) | ||
end function vector_allreps_sub | ||
end module vector_allreps_module |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,145 @@ | ||
program test_vector_allreps | ||
|
||
use precision, only : r64, i64 | ||
use testify_m, only : testify | ||
use vector_allreps_module, only : & | ||
vector_allreps, vector_allreps_add, vector_allreps_sub, & | ||
vector_allreps_print | ||
|
||
implicit none | ||
|
||
integer :: itest | ||
integer, parameter :: num_tests = 17 | ||
type(testify) :: test_array(num_tests), tests_all | ||
|
||
integer(i64) :: imuxed, grid(3) | ||
real(r64) :: primitive_vecs(3, 3) | ||
!integer(i64), parameter :: N = 3 | ||
type(vector_allreps) :: v0, v1, v2, v3 | ||
|
||
print*, '<<module vector_allreps unit tests>>' | ||
|
||
primitive_vecs(:, 1) = [-0.5_r64, 0.0_r64, 0.5_r64] | ||
primitive_vecs(:, 2) = [ 0.0_r64, 0.5_r64, 0.5_r64] | ||
primitive_vecs(:, 3) = [-0.5_r64, 0.5_r64, 0.0_r64] | ||
|
||
grid = [4, 4, 4] | ||
|
||
!Null vector | ||
v0 = vector_allreps(1_i64, grid, primitive_vecs) | ||
|
||
itest = 1 | ||
test_array(itest) = testify("null vector, integer rep") | ||
call test_array(itest)%assert(& | ||
v0%int, & | ||
[0, 0, 0]*1_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("null vector, fractional rep") | ||
call test_array(itest)%assert(& | ||
v0%frac, & | ||
[0.0, 0.0, 0.0]*1.0_r64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("null vector, Cartesian rep") | ||
call test_array(itest)%assert(& | ||
v0%cart, & | ||
[0.0, 0.0, 0.0]*1.0_r64) | ||
|
||
!A non-null vector | ||
v1 = vector_allreps(2_i64, grid, primitive_vecs) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("a non-null vector, integer rep") | ||
call test_array(itest)%assert(& | ||
v1%int, & | ||
[1, 0, 0]*1_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("a non-null vector, fractional rep") | ||
call test_array(itest)%assert(& | ||
v1%frac, & | ||
[0.25, 0.0, 0.0]*1.0_r64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("a non-null vector, Cartesian rep") | ||
call test_array(itest)%assert(& | ||
v1%cart, & | ||
[-0.125, 0.0, 0.125]*1.0_r64) | ||
|
||
!The last vector | ||
v2 = vector_allreps(product(grid), grid, primitive_vecs) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("another non-null vector, integer rep") | ||
call test_array(itest)%assert(& | ||
v2%int, & | ||
[3, 3, 3]*1_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("another non-null vector, fractional rep") | ||
call test_array(itest)%assert(& | ||
v2%frac, & | ||
[0.75, 0.75, 0.75]*1.0_r64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("another non-null vector, Cartesian rep") | ||
call test_array(itest)%assert(& | ||
v2%cart, & | ||
[-0.75, 0.75, 0.75]*1.0_r64) | ||
|
||
!Addition with Umklapp | ||
v3 = vector_allreps_add(v1, v2, grid, primitive_vecs) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("add vectors, muxed index") | ||
call test_array(itest)%assert(v3%muxed_index, 61_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("add vectors, integer rep") | ||
call test_array(itest)%assert(& | ||
v3%int, & | ||
[0, 3, 3]*1_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("add vectors, fractional rep") | ||
call test_array(itest)%assert(& | ||
v3%frac, & | ||
[0.0, 3.0, 3.0]/4.0_r64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("add vectors, Cartesian rep") | ||
call test_array(itest)%assert(& | ||
v3%cart, & | ||
[-3.0, 6.0, 3.0]/8.0_r64) | ||
|
||
!Subtraction with Umklapp | ||
v3 = vector_allreps_sub(v1, v2, grid, primitive_vecs) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("subtract vectors, muxed index") | ||
call test_array(itest)%assert(v3%muxed_index, 23_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("subtract vectors, integer rep") | ||
call test_array(itest)%assert(& | ||
v3%int, & | ||
[2, 1, 1]*1_i64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("subtract vectors, fractional rep") | ||
call test_array(itest)%assert(& | ||
v3%frac, & | ||
[0.5, 0.25, 0.25]*1.0_r64) | ||
|
||
itest = itest + 1 | ||
test_array(itest) = testify("subtract vectors, Cartesian rep") | ||
call test_array(itest)%assert(& | ||
v3%cart, & | ||
[-3.0, 2.0, 3.0]/8.0_r64) | ||
|
||
tests_all = testify(test_array) | ||
call tests_all%report | ||
|
||
if(tests_all%get_status() .eqv. .false.) error stop -1 | ||
end program test_vector_allreps |