From 0959ff3299e0f3f2a90abbfb0fd7b3ee0f2d5d84 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 28 Aug 2012 16:59:45 +0000 Subject: [PATCH] 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) --- code/DAMASK_spectral.f90 | 13 +- code/DAMASK_spectral_Utilities.f90 | 43 ++++-- code/crystallite.f90 | 79 ++++++----- code/homogenization_RGC.f90 | 4 +- code/math.f90 | 217 +++++------------------------ 5 files changed, 114 insertions(+), 242 deletions(-) diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index ea7429d56..4bf413311 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -111,14 +111,6 @@ program DAMASK_spectral materialpoint_results implicit none - -#ifdef PETSC -#include -#include -#include -#include -#include -#endif !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file real(pReal), dimension(9) :: & @@ -530,6 +522,7 @@ program DAMASK_spectral close (777) coordinates = 0.0 ! change it later!!! CPFEM_mode = 2_pInt + if (debugRestart) write(6,'(a)') 'Data read in' endif ielem = 0_pInt 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) C = C + dPdF enddo; enddo; enddo + if (debugGeneral) write(6,'(a)') 'First call to CPFEM finished' C = C * wgt !-------------------------------------------------------------------------------------------------- @@ -596,6 +590,7 @@ program DAMASK_spectral if (appendToOutFile) then open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',& form='UNFORMATTED', position='APPEND', status='OLD') + if (debugRestart) write(6,'(a)') 'Result File opened for appending' else open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut',& form='UNFORMATTED',status='REPLACE') @@ -713,7 +708,7 @@ program DAMASK_spectral j = j + 1_pInt c_reduced(k,j) = temp99_Real(n,m) 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) temp99_Real = 0.0_pReal ! build full compliance k = 0_pInt diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index 0d05f0359..28aeadf34 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -23,7 +23,7 @@ module DAMASK_spectral_Utilities IO_error implicit none - + !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW type(C_PTR), private :: plan_forward, plan_backward ! plans for fftw @@ -209,7 +209,6 @@ subroutine Utilities_updateGamma(C) memory_efficient implicit none - real(pReal), dimension(3,3,3,3), intent(in) :: C real(pReal), dimension(3,3) :: temp33_Real, xiDyad real(pReal) :: filter @@ -242,6 +241,8 @@ end subroutine Utilities_updateGamma subroutine Utilities_forwardFFT(row,column) use mesh, only : & virt_dim + use math, only: & + math_divergenceFFT implicit none integer(pInt), intent(in), optional :: row, column @@ -298,7 +299,6 @@ end subroutine Utilities_forwardFFT subroutine Utilities_backwardFFT(row,column) implicit none - integer(pInt), intent(in), optional :: row, column integer(pInt) :: i, j, k, m, n @@ -347,7 +347,6 @@ subroutine Utilities_fourierConvolution(fieldAim) memory_efficient implicit none - real(pReal), dimension(3,3), intent(in) :: fieldAim real(pReal), dimension(3,3) :: xiDyad, temp33_Real real(pReal) :: filter @@ -393,7 +392,7 @@ end subroutine Utilities_fourierConvolution !> @brief calculate root mean square of divergence of field_fourier !-------------------------------------------------------------------------------------------------- real(pReal) function Utilities_divergenceRMS() - + implicit none integer(pInt) :: i, j, k real(pReal) :: err_div_RMS, err_real_div_RMS, err_post_div_RMS,& err_div_max, err_real_div_max @@ -458,24 +457,27 @@ end function Utilities_divergenceRMS !> @brief calculates mask compliance !-------------------------------------------------------------------------------------------------- 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), 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 logical, dimension(9), intent(in) :: mask_stressVector real(pReal), dimension(3,3,3,3) :: C_lastInc real(pReal), dimension(9,9) :: temp99_Real 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 + character(len=1024):: formatString size_reduced = count(mask_stressVector) if(size_reduced > 0_pInt )then allocate (c_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) k = 0_pInt ! build reduced stiffness do n = 1_pInt,9_pInt @@ -487,7 +489,7 @@ function Utilities_maskedCompliance(rot_BC,mask_stressVector,C) j = j + 1_pInt c_reduced(k,j) = temp99_Real(n,m) 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) temp99_Real = 0.0_pReal ! build full compliance k = 0_pInt @@ -500,12 +502,26 @@ function Utilities_maskedCompliance(rot_BC,mask_stressVector,C) j = j + 1_pInt temp99_Real(n,m) = s_reduced(k,j) 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(s_reduced) + deallocate(sTimesC) else temp99_real = 0.0_pReal endif - Utilities_maskedCompliance = math_Plain99to3333(temp99_Real) end function Utilities_maskedCompliance @@ -520,7 +536,6 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti use FEsolving, only: restartWrite implicit none - real(pReal), dimension(res(1),res(2),res(3)) :: temperature 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) - + implicit none real(pReal), intent(in), dimension(3,3) :: delta_aim - real(pReal), intent(in) :: timeinc, timeinc_old, guessmode 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() implicit none - if (debugDivergence) call fftw_destroy_plan(plan_divergence) if (debugFFTW) then diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 56ba8b948..3f147c290 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -666,7 +666,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2 if (crystallite_todo(g,i,e) & .and. iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & .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: ',& crystallite_subStep(g,i,e) write(6,*) @@ -2771,49 +2771,50 @@ real(pReal), optional, intent(in) :: timeFraction ! fraction of logical crystallite_integrateStress ! flag indicating if integration suceeded !*** local variables ***! -real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep - Fp_current, & ! plastic deformation gradient at start of timestep - Fp_new, & ! plastic deformation gradient at end of timestep - Fe_new, & ! elastic deformation gradient at end of timestep - invFp_new, & ! inverse of Fp_new - invFp_current, & ! inverse of Fp_current - Lpguess, & ! current guess for plastic velocity gradient - Lpguess_old, & ! known last good guess for plastic velocity gradient - Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law - residuum, & ! current residuum of plastic velocity gradient - residuum_old, & ! last residuum of plastic velocity gradient - deltaLp, & ! direction of next guess - gradientR, & ! derivative of the residuum norm - Tstar,& ! 2nd Piola-Kirchhoff Stress +real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep + Fp_current, & ! plastic deformation gradient at start of timestep + Fp_new, & ! plastic deformation gradient at end of timestep + Fe_new, & ! elastic deformation gradient at end of timestep + invFp_new, & ! inverse of Fp_new + invFp_current, & ! inverse of Fp_current + Lpguess, & ! current guess for plastic velocity gradient + Lpguess_old, & ! known last good guess for plastic velocity gradient + Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law + residuum, & ! current residuum of plastic velocity gradient + residuum_old, & ! last residuum of plastic velocity gradient + deltaLp, & ! direction of next guess + gradientR, & ! derivative of the residuum norm + Tstar,& ! 2nd Piola-Kirchhoff Stress A,& B, & - Fe ! elastic deformation gradient -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 - dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law - dFe_dLp, & ! partial derivative of elastic deformation gradient - dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) - inv_dR_dLp ! inverse of dRdLp -real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress - dFe_dLp3333 ! partial derivative of elastic deformation gradient -real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress - det, & ! determinant + Fe ! elastic deformation gradient +real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation +real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK +integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK +real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law + dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law + dFe_dLp, & ! partial derivative of elastic deformation gradient + dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) + inv_dR_dLp ! inverse of dRdLp +real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress + dFe_dLp3333 ! partial derivative of elastic deformation gradient +real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress + det, & ! determinant expectedImprovement, & steplength0, & steplength, & steplength_max, & dt, & ! time increment aTol -logical error ! flag indicating an error -integer(pInt) NiterationStress, & ! number of stress integrations - dummy, & +logical error ! flag indicating an error +integer(pInt) NiterationStress, & ! number of stress integrations k, & l, & m, & n, & o, & p, & - jacoCounter ! counter to check for Jacobian update + jacoCounter ! counter to check for Jacobian update integer(pLongInt) tick, & tock, & tickrate, & @@ -2893,7 +2894,7 @@ LpLoop: do #endif return endif - + !* calculate 2nd Piola-Kirchhoff stress tensor @@ -3025,15 +3026,23 @@ LpLoop: do if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then dFe_dLp3333 = 0.0_pReal 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 dFe_dLp3333 = -dt * dFe_dLp3333 dFe_dLp = math_Plain3333to99(dFe_dLp3333) dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333) dR_dLp = math_identity2nd(9_pInt) - & - math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp)) - inv_dR_dLp = 0.0_pReal - call math_invert(9_pInt,dR_dLp,inv_dR_dLp,dummy,error) ! invert dR/dLp --> dLp/dR + math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp)) + inv_dR_dLp = dR_dLp ! will be changed in first call to LAPACK +#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 #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 7d8278800..76fbce7fc 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -375,7 +375,7 @@ function homogenization_RGC_updateState(& integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc 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,homogenization_maxNgrains) :: NN,pNN 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 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 if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then diff --git a/code/math.f90 b/code/math.f90 index 4b320a552..c28bcedf5 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -184,8 +184,7 @@ real(pReal), dimension(4,36), parameter, private :: & #endif private :: math_partition, & - math_delta, & - Gauss + math_delta contains @@ -287,7 +286,6 @@ subroutine math_init end subroutine math_init - !-------------------------------------------------------------------------------------------------- !> @brief Quicksort algorithm for two-dimensional integer arrays ! Sorting is done with respect to array(1,:) @@ -919,190 +917,44 @@ function math_invSym3333(A) 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 -! InvA = Inverse of A -! AnzNegEW = Number of negative Eigenvalues of A -! error = false: Inversion done. -! = true: Inversion stopped in SymGauss because of dimishing -! Pivotelement +!-------------------------------------------------------------------------------------------------- +!> @brief invert matrix of arbitrary dimension +!-------------------------------------------------------------------------------------------------- +subroutine math_invert(myDim,A, InvA, error) implicit none - integer(pInt), intent(in) :: dimen - real(pReal), dimension(dimen,dimen), intent(in) :: A - real(pReal), dimension(dimen,dimen), intent(out) :: InvA - integer(pInt), intent(out) :: AnzNegEW + integer(pInt), intent(in) :: myDim + real(pReal), dimension(myDim,myDim), intent(in) :: A + + + integer(pInt) :: ierr + integer(pInt), dimension(myDim) :: ipiv + real(pReal), dimension(myDim) :: work + + real(pReal), dimension(myDim,myDim), intent(out) :: invA logical, intent(out) :: error - real(pReal) :: LogAbsDetA - real(pReal), dimension(dimen,dimen) :: B - - InvA = math_identity2nd(dimen) - B = A - CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,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 -!-------------------------------------------------------------------------------------------------- -!> @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 !-------------------------------------------------------------------------------------------------- @@ -1272,7 +1124,6 @@ pure function math_Plain9to33(v9) end function math_Plain9to33 - !-------------------------------------------------------------------------------------------------- !> @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 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) end subroutine