Commit 16559ce9 authored by Andreas Marek's avatar Andreas Marek

Cleanup of Toeplitz matrix creation

parent fa68dca7
......@@ -221,38 +221,21 @@ program test
#endif
#if defined(__EIGENVALUES) || defined(__SOLVE_TRIDIAGONAL)
! set up simple toeplitz matrix
#ifdef TEST_DOUBLE
diagonalElement = 0.45_rk8
subdiagonalElement = 0.78_rk8
#ifdef TEST_REAL
#ifdef TEST_SINGLE
call prepare_toeplitz_matrix_real_single(na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
#else
diagonalElement = 0.45_rk4
subdiagonalElement = 0.78_rk4
call prepare_toeplitz_matrix_real_double(na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
#endif
#endif
#ifdef TEST_COMPLEX
#ifdef TEST_SINGLE
call prepare_toeplitz_matrix_complex_single(na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
#else
call prepare_toeplitz_matrix_complex_double(na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
#endif
#endif
d(:) = diagonalElement
sd(:) = subdiagonalElement
! set up the diagonal and subdiagonals (for general solver test)
do ii=1, na ! for diagonal elements
if (map_global_array_index_to_local_index(ii, ii, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(rowLocal,colLocal) = diagonalElement
endif
enddo
do ii=1, na-1
if (map_global_array_index_to_local_index(ii, ii+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(rowLocal,colLocal) = subdiagonalElement
endif
enddo
do ii=2, na
if (map_global_array_index_to_local_index(ii, ii-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(rowLocal,colLocal) = subdiagonalElement
endif
enddo
ds = d
sds = sd
as = a
#endif
if (elpa_init(CURRENT_API_VERSION) /= ELPA_OK) then
......
......@@ -55,8 +55,25 @@ module test_prepare_matrix
#endif
end interface
interface prepare_toeplitz_matrix
module procedure prepare_toeplitz_matrix_complex_double
module procedure prepare_toeplitz_matrix_real_double
#ifdef WANT_SINGLE_PRECISION_REAL
module procedure prepare_toeplitz_matrix_real_single
#endif
#ifdef WANT_SINGLE_PRECISION_COMPLEX
module procedure prepare_toeplitz_matrix_complex_single
#endif
end interface
private prows, pcols, map_global_array_index_to_local_index
contains
#include "../../src/general/prow_pcol.X90"
#include "../../src/general/map_global_to_local.X90"
#define COMPLEXCASE 1
#define DOUBLE_PRECISION 1
#include "../../src/general/precision_macros.h"
......
......@@ -47,6 +47,7 @@
&PRECISION&
& (na, myid, sc_desc, a, z, as)
use test_util
implicit none
......@@ -204,4 +205,76 @@ subroutine prepare_matrix_&
& (na, myid, sc_desc, a, z, as)
end subroutine
subroutine prepare_toeplitz_matrix_&
&MATH_DATATYPE&
&_&
&PRECISION&
& (na, d, sd, ds, sds, a, as, z, nblk, np_rows, np_cols, my_prow, my_pcol)
use iso_c_binding
use test_util
implicit none
integer, intent(in) :: na, nblk, np_rows, np_cols, my_prow, my_pcol
integer :: diagonalElement, subdiagonalElement
real(kind=C_DATATYPE_KIND) :: d(:), sd(:), ds(:), sds(:)
#if REALCASE == 1
real(kind=C_DATATYPE_KIND) :: z(:,:), a(:,:), as(:,:)
#endif
#if COMPLEXCASE == 1
complex(kind=C_DATATYPE_KIND) :: z(:,:), a(:,:), as(:,:)
#endif
integer :: ii, rowLocal, colLocal
#if REALCASE == 1
#ifdef DOUBLE_PRECISION_REAL
diagonalElement = 0.45_rk8
subdiagonalElement = 0.78_rk8
#else
diagonalElement = 0.45_rk4
subdiagonalElement = 0.78_rk4
#endif
#endif /* REALCASE */
#if COMPLEXCASE == 1
#ifdef DOUBLE_PRECISION_COMPLEX
diagonalElement = 0.45_rk8
subdiagonalElement = 0.78_rk8
#else
diagonalElement = 0.45_rk4
subdiagonalElement = 0.78_rk4
#endif
#endif /* COMPLEXCASE */
d(:) = diagonalElement
sd(:) = subdiagonalElement
! set up the diagonal and subdiagonals (for general solver test)
do ii=1, na ! for diagonal elements
if (map_global_array_index_to_local_index(ii, ii, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(rowLocal,colLocal) = diagonalElement
endif
enddo
do ii=1, na-1
if (map_global_array_index_to_local_index(ii, ii+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(rowLocal,colLocal) = subdiagonalElement
endif
enddo
do ii=2, na
if (map_global_array_index_to_local_index(ii, ii-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
a(rowLocal,colLocal) = subdiagonalElement
endif
enddo
ds = d
sds = sd
as = a
end subroutine
! vim: syntax=fortran
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment