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

LAPACK calls the unblocked versions for our small matrices 

so a work that holds exactly the data seems to be the best choice
This commit is contained in:
Martin Diehl 2020-04-10 23:54:38 +02:00
parent ff5da8fb60
commit a6d1e02b32
1 changed files with 9 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
src/math.f90
View File

@ -486,15 +486,14 @@ function math_invSym3333(A)
real(pReal),dimension(3,3,3,3),intent(in) :: A
integer, dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66
real(pReal), dimension(6*(64+2)) :: work
integer :: ierr_i, ierr_f
integer, dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66
real(pReal), dimension(6*6) :: work
integer :: ierr_i, ierr_f
temp66 = math_sym3333to66(A)
call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
if (ierr_i /= 0 .or. ierr_f /= 0) then
call IO_error(400, ext_msg = 'math_invSym3333')
else
@ -513,9 +512,9 @@ subroutine math_invert(InvA, error, A)
real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
logical, intent(out) :: error
integer, dimension(size(A,1)) :: ipiv
real(pReal), dimension(size(A,1)*(64+2)) :: work
integer :: ierr
integer, dimension(size(A,1)) :: ipiv
real(pReal), dimension(size(A,1)**2) :: work
integer :: ierr
invA = A
call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
@ -879,7 +878,7 @@ subroutine math_eigh(m,w,v,error)
logical, intent(out) :: error
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(size(m,1)**2) :: work
v = m ! copy matrix to input (doubles as output) array
call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr)
@ -1034,7 +1033,7 @@ function math_eigvalsh(m)
real(pReal), dimension(size(m,1),size(m,1)) :: m_
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(size(m,1)**2) :: work
m_= m ! copy matrix to input (will be destroyed)
call dsyev('N','U',size(m,1),m_,size(m,1),math_eigvalsh,work,size(work,1),ierr)