interface checking for LAPACK

This commit is contained in:
Martin Diehl 2020-04-10 12:52:27 +02:00
parent a86f00d827
commit a2e70612ff
5 changed files with 61 additions and 14 deletions

59
src/LAPACK_interface.f90 Normal file
View File

@ -0,0 +1,59 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Fortran interfaces for LAPACK routines
!> @details https://www.netlib.org/lapack/
!--------------------------------------------------------------------------------------------------
module LAPACK_interface
interface
subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
use prec
character, intent(in) :: jobvl,jobvr
integer, intent(in) :: n,lda,ldvl,ldvr,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
real(pReal), intent(out), dimension(n) :: wr,wi
real(pReal), intent(out), dimension(ldvl,n) :: vl
real(pReal), intent(out), dimension(ldvr,n) :: vr
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dgeev
subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
use prec
integer, intent(in) :: n,nrhs,lda,ldb
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(n) :: ipiv
real(pReal), intent(out), dimension(ldb,nrhs) :: b
integer, intent(out) :: info
end subroutine dgesv
subroutine dgetrf(m,n,a,lda,ipiv,info)
use prec
integer, intent(in) :: m,n,lda
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(min(m,n)) :: ipiv
integer, intent(out) :: info
end subroutine dgetrf
subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
use prec
integer, intent(in) :: n,lda,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(n) :: ipiv
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dgetri
subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
use prec
character, intent(in) :: jobz,uplo
integer, intent(in) :: n,lda,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
real(pReal), intent(out), dimension(n) :: w
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dsyev
end interface
end module LAPACK_interface

View File

@ -9,6 +9,7 @@
#include "list.f90"
#include "future.f90"
#include "config.f90"
#include "LAPACK_interface.f90"
#include "math.f90"
#include "quaternions.f90"
#include "Lambert.f90"

View File

@ -835,8 +835,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update
logical :: error
external :: &
dgesv
integrateStress = .false.

View File

@ -9,6 +9,7 @@ module math
use prec
use IO
use numerics
use LAPACK_interface
implicit none
public
@ -489,9 +490,6 @@ function math_invSym3333(A)
real(pReal), dimension(6,6) :: temp66
real(pReal), dimension(6*(64+2)) :: work
integer :: ierr_i, ierr_f
external :: &
dgetrf, &
dgetri
temp66 = math_sym3333to66(A)
call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
@ -518,9 +516,6 @@ subroutine math_invert(InvA, error, A)
integer, dimension(size(A,1)) :: ipiv
real(pReal), dimension(size(A,1)*(64+2)) :: work
integer :: ierr
external :: &
dgetrf, &
dgetri
invA = A
call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
@ -885,8 +880,6 @@ 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
external :: &
dsyev
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)
@ -1042,8 +1035,6 @@ 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
external :: &
dsyev
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)

View File

@ -596,8 +596,6 @@ function om2ax(om) result(ax)
real(pReal), dimension(3,3) :: VR, devNull, om_
integer :: ierr, i
external :: dgeev
om_ = om
! first get the rotation angle