223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

larger block size seems favorable

This commit is contained in:
Martin Diehl 2019-09-21 21:52:58 -07:00
parent 8b908fb350
commit 25d91c79af
1 changed files with 27 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/math.f90
View File

@ -485,19 +485,22 @@ function math_invSym3333(A)
integer :: ierr integer :: ierr
integer, dimension(6) :: ipiv6 integer, dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6,6) :: temp66
real(pReal), dimension(6) :: work6 real(pReal), dimension(6*(64+2)) :: work
logical :: error
external :: & external :: &
dgetrf, & dgetrf, &
dgetri dgetri
temp66_real = math_sym3333to66(A) temp66 = math_sym3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetrf(6,6,temp66,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) error = (ierr /= 0)
if (ierr == 0) then call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr)
math_invSym3333 = math_66toSym3333(temp66_real) error = error .or. (ierr /= 0)
else if (error) then
call IO_error(400, ext_msg = 'math_invSym3333') call IO_error(400, ext_msg = 'math_invSym3333')
else
math_invSym3333 = math_66toSym3333(temp66)
endif endif
end function math_invSym3333 end function math_invSym3333
@ -513,7 +516,7 @@ subroutine math_invert(InvA, error, A)
logical, intent(out) :: error logical, intent(out) :: error
integer, dimension(size(A,1)) :: ipiv integer, dimension(size(A,1)) :: ipiv
real(pReal), dimension(size(A,1)) :: work real(pReal), dimension(size(A,1)*(64+2)) :: work
integer :: ierr integer :: ierr
external :: & external :: &
dgetrf, & dgetrf, &
@ -521,8 +524,9 @@ subroutine math_invert(InvA, error, A)
invA = A invA = A
call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr) call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
call dgetri(size(A,1),InvA,size(A,1),ipiv,work,size(A,1),ierr) error = (ierr /= 0)
error = merge(.true.,.false., ierr /= 0) call dgetri(size(A,1),InvA,size(A,1),ipiv,work,size(work,1),ierr)
error = error .or. (ierr /= 0)
end subroutine math_invert end subroutine math_invert
@ -933,14 +937,14 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1)), intent(out) :: values
real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors
logical, intent(out) :: error logical, intent(out) :: error
integer :: info integer :: ierr
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: & external :: &
dsyev dsyev
vectors = m ! copy matrix to input (doubles as output) array vectors = m ! copy matrix to input (doubles as output) array
call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,size(work,1),ierr)
error = (info == 0) error = (ierr /= 0)
end subroutine math_eigenValuesVectorsSym end subroutine math_eigenValuesVectorsSym
@ -1184,15 +1188,15 @@ function math_eigenvaluesSym(m)
real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(:,:), intent(in) :: m
real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym
real(pReal), dimension(size(m,1),size(m,1)) :: vectors real(pReal), dimension(size(m,1),size(m,1)) :: m_
integer :: info integer :: ierr
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
external :: & external :: &
dsyev dsyev
vectors = m ! copy matrix to input (doubles as output) array m_= m ! copy matrix to input (will be destroyed)
call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) call dsyev('N','U',size(m,1),m_,size(m,1),math_eigenvaluesSym,work,size(work,1),ierr)
if (info /= 0) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) if (ierr /= 0) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
end function math_eigenvaluesSym end function math_eigenvaluesSym