Skip to content
GitLab
Commit 16433f9e authored by Andreas Marek's avatar Andreas Marek
Browse files

Map old cholesky routines to new interface

parent 3c7e1972
Show whitespace changes
Inline Side-by-side
src/elpa1/elpa_cholesky_template.X90
View file @ 16433f9e
......@@ -43,8 +43,9 @@
! the original distribution, the GNU Lesser General Public License.
#include "../sanity.X90"
use elpa1_compute
use elpa_utilities
use elpa_type
! use elpa1_compute
! use elpa_utilities
use elpa_mpi
#ifdef HAVE_DETAILED_TIMINGS
use timings
......@@ -56,35 +57,39 @@
integer(kind=ik) :: na, lda, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols
#if REALCASE == 1
#ifdef USE_ASSUMED_SIZE
real(kind=REAL_DATATYPE) :: a(lda,*)
#else
!#ifdef USE_ASSUMED_SIZE
! real(kind=REAL_DATATYPE) :: a(lda,*)
!#else
real(kind=REAL_DATATYPE) :: a(lda,matrixCols)
#endif
!#endif
#endif
#if COMPLEXCASE == 1
#ifdef USE_ASSUMED_SIZE
complex(kind=COMPLEX_DATATYPE) :: a(lda,*)
#else
!#ifdef USE_ASSUMED_SIZE
! complex(kind=COMPLEX_DATATYPE) :: a(lda,*)
!#else
complex(kind=COMPLEX_DATATYPE) :: a(lda,matrixCols)
#endif
#endif
integer(kind=ik) :: my_prow, my_pcol, np_rows, np_cols, mpierr
integer(kind=ik) :: l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
integer(kind=ik) :: n, nc, i, info
integer(kind=ik) :: lcs, lce, lrs, lre
integer(kind=ik) :: tile_size, l_rows_tile, l_cols_tile
#if REALCASE == 1
real(kind=REAL_DATATYPE), allocatable :: tmp1(:), tmp2(:,:), tmatr(:,:), tmatc(:,:)
#endif
#if COMPLEXCASE == 1
complex(kind=COMPLEX_DATATYPE), allocatable :: tmp1(:), tmp2(:,:), tmatr(:,:), tmatc(:,:)
#endif
!#endif
#endif
! integer(kind=ik) :: my_prow, my_pcol, np_rows, np_cols, mpierr
integer(kind=ik) :: nev
! integer(kind=ik) :: l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
! integer(kind=ik) :: n, nc, i, info
! integer(kind=ik) :: lcs, lce, lrs, lre
! integer(kind=ik) :: tile_size, l_rows_tile, l_cols_tile
!#if REALCASE == 1
! real(kind=REAL_DATATYPE), allocatable :: tmp1(:), tmp2(:,:), tmatr(:,:), tmatc(:,:)
!#endif
!#if COMPLEXCASE == 1
! complex(kind=COMPLEX_DATATYPE), allocatable :: tmp1(:), tmp2(:,:), tmatr(:,:), tmatc(:,:)
!#endif
logical, intent(in) :: wantDebug
logical :: success
integer(kind=ik) :: istat
character(200) :: errorMessage
integer(kind=ik) :: successInternal
! integer(kind=ik) :: istat
! character(200) :: errorMessage
type(elpa_t) :: elpaAPI
call timer%start("elpa_cholesky_&
&MATH_DATATYPE&
......@@ -92,248 +97,39 @@
&PRECISION &
")
call timer%start("mpi_communication")
call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
call mpi_comm_size(mpi_comm_rows,np_rows,mpierr)
call mpi_comm_rank(mpi_comm_cols,my_pcol,mpierr)
call mpi_comm_size(mpi_comm_cols,np_cols,mpierr)
call timer%stop("mpi_communication")
success = .true.
! Matrix is split into tiles; work is done only for tiles on the diagonal or above
tile_size = nblk*least_common_multiple(np_rows,np_cols) ! minimum global tile size
tile_size = ((128*max(np_rows,np_cols)-1)/tile_size+1)*tile_size ! make local tiles at least 128 wide
l_rows_tile = tile_size/np_rows ! local rows of a tile
l_cols_tile = tile_size/np_cols ! local cols of a tile
l_rows = local_index(na, my_prow, np_rows, nblk, -1) ! Local rows of a
l_cols = local_index(na, my_pcol, np_cols, nblk, -1) ! Local cols of a
allocate(tmp1(nblk*nblk), stat=istat, errmsg=errorMessage)
if (istat .ne. 0) then
print *,"elpa_cholesky_&
&MATH_DATATYPE&: error when allocating tmp1 "//errorMessage
stop 1
endif
allocate(tmp2(nblk,nblk), stat=istat, errmsg=errorMessage)
if (istat .ne. 0) then
print *,"elpa_cholesky_&
&MATH_DATATYPE&
&: error when allocating tmp2 "//errorMessage
stop 1
endif
tmp1 = 0
tmp2 = 0
allocate(tmatr(l_rows,nblk), stat=istat, errmsg=errorMessage)
if (istat .ne. 0) then
print *,"elpa_cholesky_&
&MATH_DATATYPE&
&: error when allocating tmatr "//errorMessage
stop 1
endif
allocate(tmatc(l_cols,nblk), stat=istat, errmsg=errorMessage)
if (istat .ne. 0) then
print *,"elpa_cholesky_&
&MATH_DATATYPE&
&: error when allocating tmatc "//errorMessage
stop 1
if (elpa_init(20170403) /= ELPA_OK) then
success = .false.
error stop "ELPA API version not supported"
endif
tmatr = 0
tmatc = 0
do n = 1, na, nblk
! Calculate first local row and column of the still remaining matrix
! on the local processor
l_row1 = local_index(n, my_prow, np_rows, nblk, +1)
l_col1 = local_index(n, my_pcol, np_cols, nblk, +1)
l_rowx = local_index(n+nblk, my_prow, np_rows, nblk, +1)
l_colx = local_index(n+nblk, my_pcol, np_cols, nblk, +1)
if (n+nblk > na) then
! This is the last step, just do a Cholesky-Factorization
! of the remaining block
if (my_prow==prow(n, nblk, np_rows) .and. my_pcol==pcol(n, nblk, np_cols)) then
call timer%start("blas")
call PRECISION_POTRF('U', na-n+1, a(l_row1,l_col1), lda, info)
call timer%stop("blas")
if (info/=0) then
if (wantDebug) write(error_unit,*) "elpa_cholesky_&
&MATH_DATATYPE&
#if REALCASE == 1
&: Error in dpotrf: ",info
#endif
#if COMPLEXCASE == 1
&: Error in zpotrf: ",info
#endif
nev = 10
elpaAPI = elpa_create(na, nev, lda, matrixCols, nblk, successInternal)
if (successInternal .ne. ELPA_OK) then
print *, "Cannot create elpa object"
success = .false.
stop
return
endif
endif
exit ! Loop
endif
if (my_prow==prow(n, nblk, np_rows)) then
if (my_pcol==pcol(n, nblk, np_cols)) then
! The process owning the upper left remaining block does the
! Cholesky-Factorization of this block
call timer%start("blas")
call PRECISION_POTRF('U', nblk, a(l_row1,l_col1), lda, info)
call timer%stop("blas")
if (info/=0) then
if (wantDebug) write(error_unit,*) "elpa_cholesky_&
&MATH_DATATYPE&
call elpaAPI%set_comm_rows(mpi_comm_rows)
call elpaAPI%set_comm_cols(mpi_comm_cols)
#if REALCASE == 1
&: Error in dpotrf 2: ",info
#endif
#if COMPLEXCASE == 1
&: Error in zpotrf 2: ",info
call elpaAPI%cholesky(a, successInternal)
#endif
if (successInternal .ne. ELPA_OK) then
print *, "Cannot run cholesky"
stop
success = .false.
return
else
success =.true.
endif
call elpaAPI%destroy()
nc = 0
do i=1,nblk
tmp1(nc+1:nc+i) = a(l_row1:l_row1+i-1,l_col1+i-1)
nc = nc+i
enddo
endif
#ifdef WITH_MPI
call timer%start("mpi_communication")
call MPI_Bcast(tmp1, nblk*(nblk+1)/2, &
#if REALCASE == 1
MPI_REAL_PRECISION, &
#endif
#if COMPLEXCASE == 1
MPI_COMPLEX_PRECISION, &
#endif
pcol(n, nblk, np_cols), mpi_comm_cols, mpierr)
call timer%stop("mpi_communication")
#endif /* WITH_MPI */
nc = 0
do i=1,nblk
tmp2(1:i,i) = tmp1(nc+1:nc+i)
nc = nc+i
enddo
call timer%start("blas")
if (l_cols-l_colx+1>0) &
#if REALCASE == 1
call PRECISION_TRSM('L', 'U', 'T', 'N', nblk, l_cols-l_colx+1, CONST_1_0, tmp2, ubound(tmp2,dim=1), &
a(l_row1,l_colx), lda)
#endif
#if COMPLEXCASE == 1
call PRECISION_TRSM('L', 'U', 'C', 'N', nblk, l_cols-l_colx+1, CONST_COMPLEX_PAIR_1_0, &
tmp2, ubound(tmp2,dim=1), a(l_row1,l_colx), lda)
#endif
call timer%stop("blas")
endif
do i=1,nblk
#if REALCASE == 1
if (my_prow==prow(n, nblk, np_rows)) tmatc(l_colx:l_cols,i) = a(l_row1+i-1,l_colx:l_cols)
#endif
#if COMPLEXCASE == 1
if (my_prow==prow(n, nblk, np_rows)) tmatc(l_colx:l_cols,i) = conjg(a(l_row1+i-1,l_colx:l_cols))
#endif
#ifdef WITH_MPI
call timer%start("mpi_communication")
if (l_cols-l_colx+1>0) &
call MPI_Bcast(tmatc(l_colx,i), l_cols-l_colx+1, &
#if REALCASE == 1
MPI_REAL_PRECISION, &
#endif
#if COMPLEXCASE == 1
MPI_COMPLEX_PRECISION, &
#endif
prow(n, nblk, np_rows), mpi_comm_rows, mpierr)
call timer%stop("mpi_communication")
#endif /* WITH_MPI */
enddo
! this has to be checked since it was changed substantially when doing type safe
call elpa_transpose_vectors_&
&MATH_DATATYPE&
&_&
&PRECISION &
(tmatc, ubound(tmatc,dim=1), mpi_comm_cols, &
tmatr, ubound(tmatr,dim=1), mpi_comm_rows, &
n, na, nblk, nblk)
do i=0,(na-1)/tile_size
lcs = max(l_colx,i*l_cols_tile+1)
lce = min(l_cols,(i+1)*l_cols_tile)
lrs = l_rowx
lre = min(l_rows,(i+1)*l_rows_tile)
if (lce<lcs .or. lre<lrs) cycle
call timer%start("blas")
#if REALCASE == 1
call PRECISION_GEMM('N', 'T', lre-lrs+1, lce-lcs+1, nblk, -CONST_1_0, &
tmatr(lrs,1), ubound(tmatr,dim=1), tmatc(lcs,1), ubound(tmatc,dim=1), &
CONST_1_0, a(lrs,lcs), lda)
#endif
#if COMPLEXCASE == 1
call PRECISION_GEMM('N', 'C', lre-lrs+1, lce-lcs+1, nblk, -CONST_COMPLEX_PAIR_1_0, &
tmatr(lrs,1), ubound(tmatr,dim=1), tmatc(lcs,1), ubound(tmatc,dim=1), &
CONST_COMPLEX_PAIR_1_0, a(lrs,lcs), lda)
#endif
call timer%stop("blas")
enddo
enddo
deallocate(tmp1, tmp2, tmatr, tmatc, stat=istat, errmsg=errorMessage)
if (istat .ne. 0) then
print *,"elpa_cholesky_&
&MATH_DATATYPE&
&: error when deallocating tmp1 "//errorMessage
stop 1
endif
! Set the lower triangle to 0, it contains garbage (form the above matrix multiplications)
call elpa_uninit()
do i=1,na
if (my_pcol==pcol(i, nblk, np_cols)) then
! column i is on local processor
l_col1 = local_index(i , my_pcol, np_cols, nblk, +1) ! local column number
l_row1 = local_index(i+1, my_prow, np_rows, nblk, +1) ! first row below diagonal
a(l_row1:l_rows,l_col1) = 0
endif
enddo
call timer%stop("elpa_cholesky_&
&MATH_DATATYPE&
&_&
......
src/elpa1/elpa_cholesky_template_new_interface.X90
View file @ 16433f9e
......@@ -56,18 +56,18 @@
integer(kind=ik) :: na, lda, nblk, matrixCols, mpi_comm_rows, mpi_comm_cols
#if REALCASE == 1
#ifdef USE_ASSUMED_SIZE
real(kind=REAL_DATATYPE) :: a(lda,*)
#else
!#ifdef USE_ASSUMED_SIZE
! real(kind=REAL_DATATYPE) :: a(lda,*)
!#else
real(kind=REAL_DATATYPE) :: a(lda,matrixCols)
#endif
!#endif
#endif
#if COMPLEXCASE == 1
#ifdef USE_ASSUMED_SIZE
complex(kind=COMPLEX_DATATYPE) :: a(lda,*)
#else
!#ifdef USE_ASSUMED_SIZE
! complex(kind=COMPLEX_DATATYPE) :: a(lda,*)
!#else
complex(kind=COMPLEX_DATATYPE) :: a(lda,matrixCols)
#endif
!#endif
#endif
integer(kind=ik) :: my_prow, my_pcol, np_rows, np_cols, mpierr
integer(kind=ik) :: l_cols, l_rows, l_col1, l_row1, l_colx, l_rowx
......@@ -89,8 +89,8 @@
call timer%start("elpa_cholesky_&
&MATH_DATATYPE&
&_&
&PRECISION &
")
&PRECISION&
&_new")
call timer%start("mpi_communication")
call mpi_comm_rank(mpi_comm_rows,my_prow,mpierr)
......@@ -338,7 +338,7 @@
&MATH_DATATYPE&
&_&
&PRECISION&
")
&_new")
#undef REALCASE
#undef COMPLEXCASE
......
src/elpa_t.F90
View file @ 16433f9e
......@@ -81,10 +81,10 @@ module elpa_type
elpa_multiply_ah_b_double, &
elpa_multiply_at_b_single, &
elpa_multiply_ah_b_single
generic, public :: cholesky => elpa_cholesky_real_double, &
elpa_cholesky_real_single, &
elpa_cholesky_complex_double, &
elpa_cholesky_complex_single
generic, public :: cholesky => elpa_cholesky_double_real, &
elpa_cholesky_single_real, &
elpa_cholesky_double_complex, &
elpa_cholesky_single_complex
......@@ -104,10 +104,10 @@ module elpa_type
procedure, private :: elpa_multiply_ah_b_double
procedure, private :: elpa_multiply_ah_b_single
procedure, private :: elpa_cholesky_real_double
procedure, private :: elpa_cholesky_real_single
procedure, private :: elpa_cholesky_complex_double
procedure, private :: elpa_cholesky_complex_single
procedure, private :: elpa_cholesky_double_real
procedure, private :: elpa_cholesky_single_real
procedure, private :: elpa_cholesky_double_complex
procedure, private :: elpa_cholesky_single_complex
end type elpa_t
logical :: initDone = .false.
......@@ -739,7 +739,7 @@ module elpa_type
#endif
end subroutine
subroutine elpa_cholesky_real_double (self, a, success)
subroutine elpa_cholesky_double_real (self, a, success)
use iso_c_binding
use elpa1_auxiliary_new
use precision
......@@ -766,6 +766,7 @@ module elpa_type
wantDebugIntern = .false.
endif
success_l = elpa_cholesky_real_double_new (self%na, a, self%local_nrows, self%nblk, &
self%local_ncols, self%mpi_comm_rows, self%mpi_comm_cols, &
wantDebugIntern)
......@@ -780,7 +781,7 @@ module elpa_type
endif
end subroutine
subroutine elpa_cholesky_real_single (self, a, success)
subroutine elpa_cholesky_single_real(self, a, success)
use iso_c_binding
use elpa1_auxiliary_new
use precision
......@@ -823,7 +824,7 @@ module elpa_type
endif
end subroutine
subroutine elpa_cholesky_complex_double (self, a, success)
subroutine elpa_cholesky_double_complex (self, a, success)
use iso_c_binding
use elpa1_auxiliary_new
use precision
......@@ -864,7 +865,7 @@ module elpa_type
endif
end subroutine
subroutine elpa_cholesky_complex_single (self, a, success)
subroutine elpa_cholesky_single_complex (self, a, success)
use iso_c_binding
use elpa1_auxiliary_new
use precision
......
Supports Markdown
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