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