From aa679129bed6857da2dd6ec28067a46e913d1608 Mon Sep 17 00:00:00 2001 From: Peter R Date: Tue, 26 Nov 2024 15:26:47 +0200 Subject: [PATCH] Fix to previous --- fem/src/SolverUtils.F90 | 55 +++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90 index 532a269355..7ee273a3dd 100644 --- a/fem/src/SolverUtils.F90 +++ b/fem/src/SolverUtils.F90 @@ -15917,15 +15917,17 @@ END SUBROUTINE BlockSolveExt IF( BlockMode ) THEN CALL Info(Caller,'Solving linear system with block strategy',Level=10) ! Here activate constraint solve only if constraints are not treated as blocks - IF( RestrictionMode ) THEN + IF( RestrictionMode .AND. & + ListGetLogical( Params, 'Eliminate Linear Constraints', Found) ) THEN BLOCK TYPE(Matrix_t), POINTER :: Acoll Acoll => AllocateMatrix() Acoll % FORMAT = MATRIX_LIST - CALL Warn(Caller,'Eliminating constraints before going into block matrix!') + CALL Info(Caller,'Eliminating constraints before going into block matrix!') CALL EliminateLinearRestriction( A, bb, A % ConstraintMatrix, Acoll, Solver, .TRUE. ) CALL List_ToCRSMatrix(Acoll) CALL BlockSolveExt( Acoll, x, Acoll % rhs, Solver ) + CALL Info(Caller,'Freeing collection matrix after solution',Level=10) CALL FreeMatrix(Acoll) Acoll => NULL() END BLOCK @@ -18527,8 +18529,14 @@ END SUBROUTINE ChangeToHarmonicSystem +!------------------------------------------------------------------------------ +!> Eliminate linear restriction only using the ListMatrix structure. +!> This is placed in a separate routine such that it can be called +!> when solving with and without block matrix being active. +!------------------------------------------------------------------------------ SUBROUTINE EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, & CollectionMatrix, Solver, CopyStiffMatrix ) + IMPLICIT NONE TYPE(Matrix_t) :: StiffMatrix REAL(KIND=dp) :: ForceVector(:) TYPE(Matrix_t), POINTER :: RestMatrix @@ -18869,11 +18877,10 @@ END SUBROUTINE EliminateLinearRestriction !> ConstraintMatrix. !------------------------------------------------------------------------------ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & - Solution, Norm, DOFs, Solver ) + Solution, Norm, DOFs, Solver ) !------------------------------------------------------------------------------ IMPLICIT NONE TYPE(Matrix_t), POINTER :: StiffMatrix !< Linear equation matrix information. - !< The restriction matrix is assumed to be in the EMatrix-field REAL(KIND=dp),TARGET :: ForceVector(:) !< The right hand side of the linear equation REAL(KIND=dp),TARGET :: Solution(:) !< Previous solution as input, new solution as output. REAL(KIND=dp) :: Norm !< The L2 norm of the solution. @@ -18881,37 +18888,28 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & TYPE(Solver_t), TARGET :: Solver !< Linear equation solver options. !------------------------------------------------------------------------------ TYPE(Solver_t), POINTER :: SolverPointer - TYPE(Matrix_t), POINTER :: CollectionMatrix, RestMatrix, AddMatrix, & - RestMatrixTranspose, TMat, XMat - REAL(KIND=dp), POINTER CONTIG :: CollectionVector(:), RestVector(:),& - AddVector(:), Tvals(:), Vals(:) + TYPE(Matrix_t), POINTER :: CollectionMatrix, RestMatrix, AddMatrix, RestMatrixTranspose + REAL(KIND=dp), POINTER CONTIG :: CollectionVector(:), RestVector(:), AddVector(:) REAL(KIND=dp), POINTER :: MultiplierValues(:), pSol(:),DiagScaling(:) - REAL(KIND=dp), ALLOCATABLE, TARGET :: CollectionSolution(:), TotValues(:) + REAL(KIND=dp), ALLOCATABLE, TARGET :: CollectionSolution(:) INTEGER :: NumberOfRows, NumberOfValues, MultiplierDOFs, istat, NoEmptyRows INTEGER :: i, j, k, l, m, n, p,q, ix, Loop, colj, nIter TYPE(Variable_t), POINTER :: MultVar, iterV REAL(KIND=dp) :: scl, rowsum, Relax, val LOGICAL :: Found, ExportMultiplier, NotExplicit, Refactorize, EnforceDirichlet, & - NonEmptyRow, ComplexSystem, ConstraintScaling, UseTranspose, EliminateConstraints, & - SkipConstraints, ResidualMode - SAVE MultiplierValues, SolverPointer - - TYPE(ListMatrix_t), POINTER :: cList - TYPE(ListMatrixEntry_t), POINTER :: cPtr, cPrev, cTmp - - INTEGER, ALLOCATABLE, TARGET :: SlavePerm(:), SlaveIPerm(:), MasterPerm(:), MasterIPerm(:) - INTEGER, POINTER :: UsePerm(:), UseIPerm(:) - REAL(KIND=dp), POINTER :: UseDiag(:), svals(:) - TYPE(ListMatrix_t), POINTER :: Lmat(:) - LOGICAL :: EliminateFromMaster, EliminateSlave, Parallel, UseTreeGauge, & - NeedMassDampValues, DoOwnScaling - REAL(KIND=dp), ALLOCATABLE, TARGET :: SlaveDiag(:), MasterDiag(:), DiagDiag(:) + NonEmptyRow, ComplexSystem, ConstraintScaling, UseTranspose, EliminateConstraints, & + SkipConstraints, ResidualMode + INTEGER, POINTER :: UseIPerm(:) + REAL(KIND=dp), POINTER :: UseDiag(:) + LOGICAL :: Parallel, UseTreeGauge, NeedMassDampValues, DoOwnScaling LOGICAL, ALLOCATABLE :: TrueDof(:) INTEGER, ALLOCATABLE :: Iperm(:) REAL(KIND=dp) :: t0,rt0,st,rst CHARACTER(:), ALLOCATABLE :: str,MultiplierName TYPE(ValueList_t), POINTER :: Params CHARACTER(*), PARAMETER :: Caller = 'SolveWithLinearRestriction' + + SAVE MultiplierValues, SolverPointer !------------------------------------------------------------------------------ CALL Info( Caller, ' ', Level=12 ) @@ -19274,10 +19272,9 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & !------------------------------------------------------------------------------ IF (ASSOCIATED(RestMatrix) .AND. EliminateConstraints) THEN - IF( .NOT. ListGetLogical( Params,'Old Elimination', Found ) ) THEN - CALL EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, CollectionMatrix, Solver) - - ELSE +#if 1 + CALL EliminateLinearRestriction( StiffMatrix, ForceVector, RestMatrix, CollectionMatrix, Solver) +#else CALL Info(Caller,'Eliminating Constraints from CollectionMatrix',Level=12) n = StiffMatrix % NumberOfRows @@ -19555,9 +19552,9 @@ RECURSIVE SUBROUTINE SolveWithLinearRestriction( StiffMatrix, ForceVector, & END IF CALL Info(Caller,'Finished Adding ConstraintMatrix',Level=12) - END IF - +#endif END IF + CALL Info(Caller,'Reverting CollectionMatrix back to CRS matrix',Level=10) IF(CollectionMatrix % FORMAT==MATRIX_LIST) THEN