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

substituted hand written matrix inversion by LAPACK version with precision selection. 

also introduced check for inversion into DAMASK_spectral_Utilities.f90 for the stress BC calculation. This part is further improved by using 5% of the reference stiffness to avoid trouble in the fully plastic regime (where the stiffness is underestimated)

Test for Marc 2010 is updated because the new inversion give slightly different results near 0 (order of e-13)
This commit is contained in:
Martin Diehl 2012-08-28 16:59:45 +00:00
parent 73349d02f5
commit 0959ff3299
5 changed files with 114 additions and 242 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
code/DAMASK_spectral.f90
View File

@ -111,14 +111,6 @@ program DAMASK_spectral
materialpoint_results materialpoint_results
implicit none implicit none
#ifdef PETSC
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscsnes.h>
#include <finclude/petscvec.h90>
#include <finclude/petscsnes.h90>
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
real(pReal), dimension(9) :: & real(pReal), dimension(9) :: &
@ -530,6 +522,7 @@ program DAMASK_spectral
close (777) close (777)
coordinates = 0.0 ! change it later!!! coordinates = 0.0 ! change it later!!!
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
if (debugRestart) write(6,'(a)') 'Data read in'
endif endif
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
@ -544,6 +537,7 @@ program DAMASK_spectral
0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF) 0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF)
C = C + dPdF C = C + dPdF
enddo; enddo; enddo enddo; enddo; enddo
if (debugGeneral) write(6,'(a)') 'First call to CPFEM finished'
C = C * wgt C = C * wgt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -596,6 +590,7 @@ program DAMASK_spectral
if (appendToOutFile) then if (appendToOutFile) then
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',& open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',&
form='UNFORMATTED', position='APPEND', status='OLD') form='UNFORMATTED', position='APPEND', status='OLD')
if (debugRestart) write(6,'(a)') 'Result File opened for appending'
else else
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',& open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',&
form='UNFORMATTED',status='REPLACE') form='UNFORMATTED',status='REPLACE')
@ -713,7 +708,7 @@ program DAMASK_spectral
j = j + 1_pInt j = j + 1_pInt
c_reduced(k,j) = temp99_Real(n,m) c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness call math_invert(size_reduced,c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=400_pInt) if(errmatinv) call IO_error(error_ID=400_pInt)
temp99_Real = 0.0_pReal ! build full compliance temp99_Real = 0.0_pReal ! build full compliance
k = 0_pInt k = 0_pInt

39
code/DAMASK_spectral_Utilities.f90
View File

@ -209,7 +209,6 @@ subroutine Utilities_updateGamma(C)
memory_efficient memory_efficient
implicit none implicit none
real(pReal), dimension(3,3,3,3), intent(in) :: C real(pReal), dimension(3,3,3,3), intent(in) :: C
real(pReal), dimension(3,3) :: temp33_Real, xiDyad real(pReal), dimension(3,3) :: temp33_Real, xiDyad
real(pReal) :: filter real(pReal) :: filter
@ -242,6 +241,8 @@ end subroutine Utilities_updateGamma
subroutine Utilities_forwardFFT(row,column) subroutine Utilities_forwardFFT(row,column)
use mesh, only : & use mesh, only : &
virt_dim virt_dim
use math, only: &
math_divergenceFFT
implicit none implicit none
integer(pInt), intent(in), optional :: row, column integer(pInt), intent(in), optional :: row, column
@ -298,7 +299,6 @@ end subroutine Utilities_forwardFFT
subroutine Utilities_backwardFFT(row,column) subroutine Utilities_backwardFFT(row,column)
implicit none implicit none
integer(pInt), intent(in), optional :: row, column integer(pInt), intent(in), optional :: row, column
integer(pInt) :: i, j, k, m, n integer(pInt) :: i, j, k, m, n
@ -347,7 +347,6 @@ subroutine Utilities_fourierConvolution(fieldAim)
memory_efficient memory_efficient
implicit none implicit none
real(pReal), dimension(3,3), intent(in) :: fieldAim real(pReal), dimension(3,3), intent(in) :: fieldAim
real(pReal), dimension(3,3) :: xiDyad, temp33_Real real(pReal), dimension(3,3) :: xiDyad, temp33_Real
real(pReal) :: filter real(pReal) :: filter
@ -393,7 +392,7 @@ end subroutine Utilities_fourierConvolution
!> @brief calculate root mean square of divergence of field_fourier !> @brief calculate root mean square of divergence of field_fourier
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) function Utilities_divergenceRMS() real(pReal) function Utilities_divergenceRMS()
implicit none
integer(pInt) :: i, j, k integer(pInt) :: i, j, k
real(pReal) :: err_div_RMS, err_real_div_RMS, err_post_div_RMS,& real(pReal) :: err_div_RMS, err_real_div_RMS, err_post_div_RMS,&
err_div_max, err_real_div_max err_div_max, err_real_div_max
@ -459,23 +458,26 @@ end function Utilities_divergenceRMS
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function Utilities_maskedCompliance(rot_BC,mask_stressVector,C) function Utilities_maskedCompliance(rot_BC,mask_stressVector,C)
implicit none
real(pReal), dimension(3,3,3,3) :: Utilities_maskedCompliance real(pReal), dimension(3,3,3,3) :: Utilities_maskedCompliance
real(pReal), dimension(3,3,3,3), intent(in) :: C real(pReal), dimension(3,3,3,3), intent(in) :: C
integer(pInt) :: i, j, k, m, n integer(pInt) :: j, k, m, n
real(pReal), dimension(3,3), intent(in) :: rot_BC real(pReal), dimension(3,3), intent(in) :: rot_BC
logical, dimension(9), intent(in) :: mask_stressVector logical, dimension(9), intent(in) :: mask_stressVector
real(pReal), dimension(3,3,3,3) :: C_lastInc real(pReal), dimension(3,3,3,3) :: C_lastInc
real(pReal), dimension(9,9) :: temp99_Real real(pReal), dimension(9,9) :: temp99_Real
integer(pInt) :: size_reduced = 0_pInt integer(pInt) :: size_reduced = 0_pInt
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC) real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced, sTimesC ! reduced compliance and stiffness (only for stress BC)
logical :: errmatinv logical :: errmatinv
character(len=1024):: formatString
size_reduced = count(mask_stressVector) size_reduced = count(mask_stressVector)
if(size_reduced > 0_pInt )then if(size_reduced > 0_pInt )then
allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal)
C_lastInc = math_rotate_forward3333(C,rot_BC) ! calculate stiffness from former inc C_lastInc = math_rotate_forward3333(C*0.95_pReal+0.5_pReal*C_ref,rot_BC) ! calculate stiffness from former inc
temp99_Real = math_Plain3333to99(C_lastInc) temp99_Real = math_Plain3333to99(C_lastInc)
k = 0_pInt ! build reduced stiffness k = 0_pInt ! build reduced stiffness
do n = 1_pInt,9_pInt do n = 1_pInt,9_pInt
@ -487,7 +489,7 @@ function Utilities_maskedCompliance(rot_BC,mask_stressVector,C)
j = j + 1_pInt j = j + 1_pInt
c_reduced(k,j) = temp99_Real(n,m) c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=400_pInt) if(errmatinv) call IO_error(error_ID=400_pInt)
temp99_Real = 0.0_pReal ! build full compliance temp99_Real = 0.0_pReal ! build full compliance
k = 0_pInt k = 0_pInt
@ -500,12 +502,26 @@ function Utilities_maskedCompliance(rot_BC,mask_stressVector,C)
j = j + 1_pInt j = j + 1_pInt
temp99_Real(n,m) = s_reduced(k,j) temp99_Real(n,m) = s_reduced(k,j)
endif; enddo; endif; enddo endif; enddo; endif; enddo
sTimesC = matmul(c_reduced,s_reduced)
do m=1_pInt, size_reduced
do n=1_pInt, size_reduced
if(m==n .and. abs(sTimesC(m,n)) > (1.0_pReal + 10.0e-12_pReal)) errmatinv = .true.
if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true.
enddo
enddo
if(debugGeneral .or. errmatinv) then
write(formatString, '(I16.16)') size_reduced
formatString = '(a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
write(6,trim(formatString),advance='no') 'C * S', transpose(matmul(c_reduced,s_reduced))
write(6,trim(formatString),advance='no') 'S', transpose(s_reduced)
endif
if(errmatinv) call IO_error(error_ID=400_pInt)
deallocate(c_reduced) deallocate(c_reduced)
deallocate(s_reduced) deallocate(s_reduced)
deallocate(sTimesC)
else else
temp99_real = 0.0_pReal temp99_real = 0.0_pReal
endif endif
Utilities_maskedCompliance = math_Plain99to3333(temp99_Real) Utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
end function Utilities_maskedCompliance end function Utilities_maskedCompliance
@ -520,7 +536,6 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
use FEsolving, only: restartWrite use FEsolving, only: restartWrite
implicit none implicit none
real(pReal), dimension(res(1),res(2),res(3)) :: temperature real(pReal), dimension(res(1),res(2),res(3)) :: temperature
real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates
@ -578,9 +593,8 @@ end subroutine Utilities_constitutiveResponse
subroutine Utilities_forwardField(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field) subroutine Utilities_forwardField(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field)
implicit none
real(pReal), intent(in), dimension(3,3) :: delta_aim real(pReal), intent(in), dimension(3,3) :: delta_aim
real(pReal), intent(in) :: timeinc, timeinc_old, guessmode real(pReal), intent(in) :: timeinc, timeinc_old, guessmode
real(pReal), intent(inout), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field real(pReal), intent(inout), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field
@ -620,7 +634,6 @@ end function Utilities_getFilter
subroutine Utilities_destroy() subroutine Utilities_destroy()
implicit none implicit none
if (debugDivergence) call fftw_destroy_plan(plan_divergence) if (debugDivergence) call fftw_destroy_plan(plan_divergence)
if (debugFFTW) then if (debugFFTW) then

75
code/crystallite.f90
View File

@ -666,7 +666,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
if (crystallite_todo(g,i,e) & if (crystallite_todo(g,i,e) &
.and. iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & .and. iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,f12.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& write(6,'(a,f12.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e) crystallite_subStep(g,i,e)
write(6,*) write(6,*)
@ -2771,49 +2771,50 @@ real(pReal), optional, intent(in) :: timeFraction ! fraction of
logical crystallite_integrateStress ! flag indicating if integration suceeded logical crystallite_integrateStress ! flag indicating if integration suceeded
!*** local variables ***! !*** local variables ***!
real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep
Fp_current, & ! plastic deformation gradient at start of timestep Fp_current, & ! plastic deformation gradient at start of timestep
Fp_new, & ! plastic deformation gradient at end of timestep Fp_new, & ! plastic deformation gradient at end of timestep
Fe_new, & ! elastic deformation gradient at end of timestep Fe_new, & ! elastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new invFp_new, & ! inverse of Fp_new
invFp_current, & ! inverse of Fp_current invFp_current, & ! inverse of Fp_current
Lpguess, & ! current guess for plastic velocity gradient Lpguess, & ! current guess for plastic velocity gradient
Lpguess_old, & ! known last good guess for plastic velocity gradient Lpguess_old, & ! known last good guess for plastic velocity gradient
Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
residuum, & ! current residuum of plastic velocity gradient residuum, & ! current residuum of plastic velocity gradient
residuum_old, & ! last residuum of plastic velocity gradient residuum_old, & ! last residuum of plastic velocity gradient
deltaLp, & ! direction of next guess deltaLp, & ! direction of next guess
gradientR, & ! derivative of the residuum norm gradientR, & ! derivative of the residuum norm
Tstar,& ! 2nd Piola-Kirchhoff Stress Tstar,& ! 2nd Piola-Kirchhoff Stress
A,& A,&
B, & B, &
Fe ! elastic deformation gradient Fe ! elastic deformation gradient
real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real(pReal), dimension(9,9):: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
dFe_dLp, & ! partial derivative of elastic deformation gradient real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law
dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law
inv_dR_dLp ! inverse of dRdLp dFe_dLp, & ! partial derivative of elastic deformation gradient
real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
dFe_dLp3333 ! partial derivative of elastic deformation gradient inv_dR_dLp ! inverse of dRdLp
real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
det, & ! determinant dFe_dLp3333 ! partial derivative of elastic deformation gradient
real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress
det, & ! determinant
expectedImprovement, & expectedImprovement, &
steplength0, & steplength0, &
steplength, & steplength, &
steplength_max, & steplength_max, &
dt, & ! time increment dt, & ! time increment
aTol aTol
logical error ! flag indicating an error logical error ! flag indicating an error
integer(pInt) NiterationStress, & ! number of stress integrations integer(pInt) NiterationStress, & ! number of stress integrations
dummy, &
k, & k, &
l, & l, &
m, & m, &
n, & n, &
o, & o, &
p, & p, &
jacoCounter ! counter to check for Jacobian update jacoCounter ! counter to check for Jacobian update
integer(pLongInt) tick, & integer(pLongInt) tick, &
tock, & tock, &
tickrate, & tickrate, &
@ -3025,15 +3026,23 @@ LpLoop: do
if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
dFe_dLp3333 = 0.0_pReal dFe_dLp3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dFe_dLp3333(p,o,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(j,l) dFe_dLp3333(p,o,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(j,l)
enddo; enddo enddo; enddo
dFe_dLp3333 = -dt * dFe_dLp3333 dFe_dLp3333 = -dt * dFe_dLp3333
dFe_dLp = math_Plain3333to99(dFe_dLp3333) dFe_dLp = math_Plain3333to99(dFe_dLp3333)
dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333) dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333)
dR_dLp = math_identity2nd(9_pInt) - & dR_dLp = math_identity2nd(9_pInt) - &
math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp)) math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp))
inv_dR_dLp = 0.0_pReal inv_dR_dLp = dR_dLp ! will be changed in first call to LAPACK
call math_invert(9_pInt,dR_dLp,inv_dR_dLp,dummy,error) ! invert dR/dLp --> dLp/dR #if(FLOAT==8)
call dgetrf(9,9,inv_dR_dLp,9,ipiv,error) ! invert dR/dLp --> dLp/dR
call dgetri(9,inv_dR_dLp,9,ipiv,work,9,error)
#elif(FLOAT==4)
call sgetrf(9,9,inv_dR_dLp,9,ipiv,error) ! invert dR/dLp --> dLp/dR
call sgetri(9,inv_dR_dLp,9,ipiv,work,9,error)
#else
NO SUITABLE PRECISION SELECTED, COMPILATION ABORTED
#endif
if (error) then if (error) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then

4
code/homogenization_RGC.f90
View File

@ -375,7 +375,7 @@ function homogenization_RGC_updateState(&
integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID
integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc
integer(pInt), dimension (2) :: residLoc integer(pInt), dimension (2) :: residLoc
integer(pInt) homID,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain,nGrain integer(pInt) homID,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain
real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD
real(pReal), dimension (3,homogenization_maxNgrains) :: NN,pNN real(pReal), dimension (3,homogenization_maxNgrains) :: NN,pNN
real(pReal), dimension (3) :: normP,normN,mornP,mornN real(pReal), dimension (3) :: normP,normN,mornP,mornN
@ -740,7 +740,7 @@ function homogenization_RGC_updateState(&
!* ------------------------------------------------------------------------------------------------------------- !* -------------------------------------------------------------------------------------------------------------
!*** Computing the update of the state variable (relaxation vectors) using the Jacobian matrix !*** Computing the update of the state variable (relaxation vectors) using the Jacobian matrix
allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot)); jnverse = 0.0_pReal allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot)); jnverse = 0.0_pReal
call math_invert(3_pInt*nIntFaceTot,jmatrix,jnverse,ival,error) ! Compute the inverse of the overall Jacobian matrix call math_invert(size(jmatrix,1),jmatrix,jnverse,error) ! Compute the inverse of the overall Jacobian matrix
!* Debugging the inverse Jacobian matrix !* Debugging the inverse Jacobian matrix
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then

215
code/math.f90
View File

@ -184,8 +184,7 @@ real(pReal), dimension(4,36), parameter, private :: &
#endif #endif
private :: math_partition, & private :: math_partition, &
math_delta, & math_delta
Gauss
contains contains
@ -287,7 +286,6 @@ subroutine math_init
end subroutine math_init end subroutine math_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Quicksort algorithm for two-dimensional integer arrays !> @brief Quicksort algorithm for two-dimensional integer arrays
! Sorting is done with respect to array(1,:) ! Sorting is done with respect to array(1,:)
@ -919,190 +917,44 @@ function math_invSym3333(A)
end function math_invSym3333 end function math_invSym3333
!--------------------------------------------------------------------------------------------------
!> @brief Gauss elimination to invert matrix of arbitrary dimension
!--------------------------------------------------------------------------------------------------
pure subroutine math_invert(dimen,A, InvA, AnzNegEW, error)
! Invertieren einer dimen x dimen - Matrix !--------------------------------------------------------------------------------------------------
! A = Matrix A !> @brief invert matrix of arbitrary dimension
! InvA = Inverse of A !--------------------------------------------------------------------------------------------------
! AnzNegEW = Number of negative Eigenvalues of A subroutine math_invert(myDim,A, InvA, error)
! error = false: Inversion done.
! = true: Inversion stopped in SymGauss because of dimishing
! Pivotelement
implicit none implicit none
integer(pInt), intent(in) :: dimen integer(pInt), intent(in) :: myDim
real(pReal), dimension(dimen,dimen), intent(in) :: A real(pReal), dimension(myDim,myDim), intent(in) :: A
real(pReal), dimension(dimen,dimen), intent(out) :: InvA
integer(pInt), intent(out) :: AnzNegEW
logical, intent(out) :: error
real(pReal) :: LogAbsDetA
real(pReal), dimension(dimen,dimen) :: B
InvA = math_identity2nd(dimen)
B = A integer(pInt) :: ierr
CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,error) integer(pInt), dimension(myDim) :: ipiv
real(pReal), dimension(myDim) :: work
real(pReal), dimension(myDim,myDim), intent(out) :: invA
logical, intent(out) :: error
invA = A
#if(FLOAT==8)
call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#elif(FLOAT==4)
call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr)
call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr)
#else
NO SUITABLE PRECISION SELECTED, COMPILATION ABORTED
#endif
if (ierr == 0_pInt) then
error = .false.
else
error = .true.
endif
end subroutine math_invert end subroutine math_invert
!--------------------------------------------------------------------------------------------------
!> @brief Solves a linear EQS A * X = B with the GAUSS-Algorithm
! For numerical stabilization using a pivot search in rows and columns
!--------------------------------------------------------------------------------------------------
pure subroutine Gauss (dimen,A,B,LogAbsDetA,NegHDK,error)
! input parameters
! A(dimen,dimen) = matrix A
! B(dimen,dimen) = right side B
!
! output parameters
! B(dimen,dimen) = Matrix containing unknown vectors X
! LogAbsDetA = 10-Logarithm of absolute value of determinatns of A
! NegHDK = Number of negative Maindiagonal coefficients resulting
! Vorwaertszerlegung
! error = false: EQS is solved
! = true : Matrix A is singular.
!
! A and B will be changed!
implicit none
logical, intent(out) :: error
integer(pInt), intent(in) :: dimen
integer(pInt), intent(out) :: NegHDK
real(pReal), intent(out) :: LogAbsDetA
real(pReal), intent(inout), dimension(dimen,dimen) :: A, B
logical :: SortX
integer(pInt) :: PivotZeile, PivotSpalte, StoreI, I, IP1, J, K, L
integer(pInt), dimension(dimen) :: XNr
real(pReal) :: AbsA, PivotWert, EpsAbs, Quote
real(pReal), dimension(dimen) :: StoreA, StoreB
error = .true.; NegHDK = 1_pInt; SortX = .false.
! Unbekanntennumerierung
DO I = 1_pInt, dimen
XNr(I) = I
ENDDO
! Genauigkeitsschranke und Bestimmung des groessten Pivotelementes
PivotWert = ABS(A(1,1))
PivotZeile = 1_pInt
PivotSpalte = 1_pInt
do I = 1_pInt, dimen; do J = 1_pInt, dimen
AbsA = ABS(A(I,J))
IF (AbsA .GT. PivotWert) THEN
PivotWert = AbsA
PivotZeile = I
PivotSpalte = J
ENDIF
enddo; enddo
IF (PivotWert .LT. 0.0000001_pReal) RETURN ! Pivotelement = 0?
EpsAbs = PivotWert * 0.1_pReal ** PRECISION(1.0_pReal)
! V O R W A E R T S T R I A N G U L A T I O N
DO I = 1_pInt, dimen - 1_pInt
! Zeilentausch?
IF (PivotZeile .NE. I) THEN
StoreA(I:dimen) = A(I,I:dimen)
A(I,I:dimen) = A(PivotZeile,I:dimen)
A(PivotZeile,I:dimen) = StoreA(I:dimen)
StoreB(1:dimen) = B(I,1:dimen)
B(I,1:dimen) = B(PivotZeile,1:dimen)
B(PivotZeile,1:dimen) = StoreB(1:dimen)
SortX = .TRUE.
ENDIF
! Spaltentausch?
IF (PivotSpalte .NE. I) THEN
StoreA(1:dimen) = A(1:dimen,I)
A(1:dimen,I) = A(1:dimen,PivotSpalte)
A(1:dimen,PivotSpalte) = StoreA(1:dimen)
StoreI = XNr(I)
XNr(I) = XNr(PivotSpalte)
XNr(PivotSpalte) = StoreI
SortX = .TRUE.
ENDIF
! Triangulation
DO J = I + 1_pInt, dimen
Quote = A(J,I) / A(I,I)
DO K = I + 1_pInt, dimen
A(J,K) = A(J,K) - Quote * A(I,K)
ENDDO
DO K = 1_pInt, dimen
B(J,K) = B(J,K) - Quote * B(I,K)
ENDDO
ENDDO
! Bestimmung des groessten Pivotelementes
IP1 = I + 1_pInt
PivotWert = ABS(A(IP1,IP1))
PivotZeile = IP1
PivotSpalte = IP1
DO J = IP1, dimen
DO K = IP1, dimen
AbsA = ABS(A(J,K))
IF (AbsA .GT. PivotWert) THEN
PivotWert = AbsA
PivotZeile = J
PivotSpalte = K
ENDIF
ENDDO
ENDDO
IF (PivotWert .LT. EpsAbs) RETURN ! Pivotelement = 0?
ENDDO
! R U E C K W A E R T S A U F L O E S U N G
DO I = dimen, 1_pInt, -1_pInt
DO L = 1_pInt, dimen
DO J = I + 1_pInt, dimen
B(I,L) = B(I,L) - A(I,J) * B(J,L)
ENDDO
B(I,L) = B(I,L) / A(I,I)
ENDDO
ENDDO
! Sortieren der Unbekanntenvektoren?
IF (SortX) THEN
DO L = 1_pInt, dimen
StoreA(1:dimen) = B(1:dimen,L)
DO I = 1_pInt, dimen
J = XNr(I)
B(J,L) = StoreA(I)
ENDDO
ENDDO
ENDIF
! Determinante
LogAbsDetA = 0.0_pReal
NegHDK = 0_pInt
DO I = 1_pInt, dimen
IF (A(I,I) .LT. 0.0_pReal) NegHDK = NegHDK + 1_pInt
AbsA = ABS(A(I,I))
LogAbsDetA = LogAbsDetA + LOG10(AbsA)
ENDDO
error = .false.
end subroutine Gauss
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief symmetrize a 33 matrix !> @brief symmetrize a 33 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1272,7 +1124,6 @@ pure function math_Plain9to33(v9)
end function math_Plain9to33 end function math_Plain9to33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 33 matrix into Mandel vector 6 !> @brief convert symmetric 33 matrix into Mandel vector 6
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2110,7 +1961,11 @@ subroutine math_spectralDecompositionSym33(M,values,vectors,error)
real(pReal), dimension((64+2)*3) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension((64+2)*3) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
vectors = M ! copy matrix to input (doubles as output) array vectors = M ! copy matrix to input (doubles as output) array
call DSYEV('V','U',3,vectors,3,values,work,(64+2)*3,info) #if(FLOAT==8)
call dsyev('V','U',3,vectors,3,values,work,(64+2)*3,info)
#elif(FLOAT==4)
call ssyev('V','U',3,vectors,3,values,work,(64+2)*3,info)
#endif
error = (info == 0_pInt) error = (info == 0_pInt)
end subroutine end subroutine