Skip to content

Commit

Permalink
Merge pull request #625 from ElmerCSC/ParallelAVBlock
Browse files Browse the repository at this point in the history 
Parallel av block
  • Loading branch information
@raback
raback authored Jan 12, 2025
2 parents af96d9b + 79c35d8 commit 1106c21
Show file tree
Hide file tree
Showing 5 changed files with 150 additions and 61 deletions.
137 changes: 93 additions & 44 deletions fem/src/BlockSolve.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -577,9 +577,21 @@ SUBROUTINE BlockInitVar( Solver, BlockMatrix, BlockIndex )

IF(PRESENT(BlockIndex)) THEN
B => BlockMatrix % SubMatrix(i,i) % Mat
IF(ASSOCIATED(B % InvPerm ) ) THEN
Var % Values = Solver % Variable % Values(B % InvPerm)

IF( ParEnv % PEs == 1 ) THEN
IF(.NOT. ASSOCIATED(B % InvPerm) ) CALL Fatal('BlockInitVar','InvPerm not present!')
IF(ASSOCIATED(B % InvPerm ) ) THEN
Var % Values = Solver % Variable % Values(B % InvPerm)
END IF
ELSE
IF(.NOT. ASSOCIATED(B % Perm) ) CALL Fatal('BlockInitVar','Perm not present!')
DO j=1,SIZE(B % Perm)
k = B % Perm(j)
IF(k==0) CYCLE
Var % Values(k) = Solver % Variable % Values(Solver % Matrix % Perm(j))
END DO
END IF

END IF

END DO
Expand Down Expand Up @@ -622,15 +634,31 @@ SUBROUTINE BlockBackCopyVar( Solver, BlockMatrix )
END IF

! Copy the block part to the monolithic solution
DO j=1,n
k = Amat % InvPerm(j)
IF( k < 1 .OR. k > m ) THEN
PRINT *,'ijk:',i,j,k
CYCLE
IF(ASSOCIATED(Amat % Perm)) THEN
IF( ParEnv % PEs == 1 ) THEN
CALL Warn('BlockBackCopyVar','Why do we have Amat % Perm in serial?')
END IF
Solver % Variable % Values(k) = Var % Values(j)
END DO
ELSE
IF( ParEnv % PEs > 1 ) THEN
CALL Warn('BlockBackCopyVar','Why do we not have Amat % Perm in parallel?')
END IF
END IF

! Note that confusingly InvPerm has different defintion in serial and parallel.
! In serial it points just to indexes of the original matrix.
! In parallel it points to all possible indexes that could be present.
IF( ParEnv % PEs > 1 ) THEN
DO j=1,SIZE(Amat % Perm)
k = Amat % Perm(j)
IF(k==0) CYCLE
Solver % Variable % Values(Solver % Matrix % Perm(j)) = Var % Values(k)
END DO
ELSE
WHERE(Amat % InvPerm > 0)
Solver % Variable % Values(Amat % InvPerm) = Var % Values
END WHERE
END IF

END DO

BlockMatrix % TotSize = BlockMatrix % Offset( NoVar + 1 )
Expand Down Expand Up @@ -1052,10 +1080,10 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
INTEGER :: Novar
LOGICAL :: DoAddMatrix

INTEGER :: bcol,brow,bi,bk,i,k,j,n,istat,NoBlock
INTEGER :: bcol,brow,bi,bk,i,k,j,n,m,istat,NoBlock,dofs
TYPE(Matrix_t), POINTER :: A, B, C
INTEGER, ALLOCATABLE :: BlockNumbering(:), rowcount(:), offset(:)
LOGICAL :: SplitComplex, Found
LOGICAL :: SplitComplex, CreatePermPar, Found
REAL(KIND=dp) :: Coeff, Coeff0

CALL Info('BlockPickMatrixPerm','Picking indexed block matrix from monolithic one',Level=10)
Expand All @@ -1077,13 +1105,13 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
END IF
TotMatrix % BlockPerm = 0


DO i=1,n
brow = BlockIndex(i)
rowcount(brow) = rowcount(brow) + 1
BlockNumbering(i) = rowcount(brow)
END DO

CreatePermPar = ASSOCIATED(A % InvPerm) .AND. ( ParEnv % PEs > 1 )

DO i = 1, NoVar
B => TotMatrix % SubMatrix(i,i) % Mat
Expand All @@ -1096,14 +1124,28 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
ALLOCATE(B % Rhs(n))
END IF
B % rhs = 0.0_dp


! Create Perm only in parallel
IF( CreatePermPar ) THEN
m = SIZE( Solver % Variable % Perm ) * Solver % Variable % Dofs
IF(ASSOCIATED(B % Perm)) THEN
IF(SIZE(B % Perm) /= m) DEALLOCATE( B % Perm)
END IF
IF(.NOT. ASSOCIATED(B % Perm)) THEN
ALLOCATE(B % Perm(m))
END IF
B % Perm = 0
END IF

! Always generate InvPerm
IF(ASSOCIATED(B % InvPerm)) THEN
IF(SIZE(B % InvPerm) /= n) DEALLOCATE( B % InvPerm)
END IF
IF(.NOT. ASSOCIATED(B % InvPerm)) THEN
ALLOCATE(B % InvPerm(n))
ALLOCATE(B % InvPerm(n))
END IF
B % InvPerm = 0
B % InvPerm = 0

! Add the (n,n) entry since this helps to create most efficiently the full ListMatrix
! CALL AddToMatrixElement(B,n,n,0.0_dp)

Expand All @@ -1118,8 +1160,18 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
B => TotMatrix % SubMatrix(brow,brow) % Mat
B % Rhs(bi) = A % Rhs(i)

B % InvPerm(bi) = i

IF( CreatePermPar ) THEN
j = A % InvPerm(i)
IF(j<0 .OR. j>SIZE(B % Perm)) THEN
PRINT *,'Too big j:',j,SIZE(B % Perm),i,SIZE(A % Perm),SIZE(A % InvPerm)
STOP
END IF
B % Perm(j) = bi
B % InvPerm(bi) = j
ELSE
B % InvPerm(bi) = i
END IF

TotMatrix % BlockPerm(offset(brow)+bi) = i

DO j=A % Rows(i+1)-1,A % Rows(i),-1
Expand All @@ -1134,8 +1186,6 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
END DO
END DO



NoBlock = NoVar
IF( DoAddMatrix ) THEN
CALL Info('BlockPickMatrixPerm','Adding add matrix to submatrices of block system!',Level=10)
Expand All @@ -1154,17 +1204,6 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
ALLOCATE(B % Rhs(n))
END IF
B % rhs = 0.0_dp

#if 0
! Don't know what to do yet...
IF(ASSOCIATED(B % InvPerm)) THEN
IF(SIZE(B % InvPerm) /= n) DEALLOCATE( B % InvPerm)
END IF
IF(.NOT. ASSOCIATED(B % InvPerm)) THEN
ALLOCATE(B % InvPerm(n))
END IF
B % InvPerm = 0
#endif

DO i=1,C % NumberOfRows
IF(i <= A % NumberOfRows ) THEN
Expand Down Expand Up @@ -1207,6 +1246,11 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar, DoAddMatrix )
CALL List_toCRSMatrix(B)
END IF
END DO

IF( i==j .AND. ParEnv % PEs > 1 ) THEN
CALL ParallelInitMatrix( Solver, B )
END IF

END DO

IF( ASSOCIATED(A % PrecValues) ) THEN
Expand Down Expand Up @@ -4082,26 +4126,31 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )

IF (isParallel) THEN
CALL Info(Caller,'Performing parallel initializations!',Level=18)
DO i=1,NoVar
A => TotMatrix % SubMatrix(i,i) % Mat
! ParallelInitSolve expects full vectors
CALL Info(Caller,'Initializing submatrix '//I2S(11*i),Level=20)
IF (ASSOCIATED(A % ParMatrix % SplittedMatrix % InsideMatrix % PrecValues)) THEN
IF (.NOT. ASSOCIATED(A % PrecValues)) &
NULLIFY(A % ParMatrix % SplittedMatrix % InsideMatrix % PrecValues)
END IF
CALL ParallelInitSolve(A, TotMatrix % Subvector(i) % Var % Values, A % rhs, r )
IF( ASSOCIATED(SolverMatrix)) THEN
x(offset(i)+1:offset(i+1)) = TotMatrix % SubVector(i) % Var % Values
IF(ASSOCIATED(A % rhs)) b(offset(i)+1:offset(i+1)) = A % rhs
END IF
CALL Info(Caller,'Done initializing submatrix '//I2S(11*i),Level=20)
END DO
DO i=1,NoVar
DO j=1,NoVar
A => TotMatrix % SubMatrix(i,j) % Mat
! ParallelInitSolve expects full vectors
IF ( i /= j ) THEN
IF(ASSOCIATED(A % ParMatrix)) CALL ParallelInitSolve(A,r,r,r)
ELSE
IF (ASSOCIATED(A % ParMatrix % SplittedMatrix % InsideMatrix % PrecValues)) THEN
IF (.NOT. ASSOCIATED(A % PrecValues)) &
NULLIFY(A % ParMatrix % SplittedMatrix % InsideMatrix % PrecValues)
END IF
CALL ParallelInitSolve(A, TotMatrix % Subvector(i) % Var % Values, A % rhs, r )
IF( ASSOCIATED(SolverMatrix)) THEN
x(offset(i)+1:offset(i+1)) = TotMatrix % SubVector(i) % Var % Values
IF(ASSOCIATED(A % rhs)) b(offset(i)+1:offset(i+1)) = A % rhs
END IF
END IF
IF(i==j) CYCLE
CALL Info(Caller,'Initializing submatrix '//I2S(10*i+j),Level=20)
IF(ASSOCIATED(A % ParMatrix)) CALL ParallelInitSolve(A,r,r,r)
END DO
END DO
IF(ASSOCIATED(SolverMatrix)) THEN
CALL Info(Caller,'Initializing SolverMatrix',Level=20)
CALL ParallelInitSolve( SolverMatrix, x, b, r )
x = 0.0_dp
b = 0.0_dp
Expand Down
6 changes: 5 additions & 1 deletion fem/src/DefUtils.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3557,7 +3557,11 @@ RECURSIVE SUBROUTINE DefaultStart( USolver )
IF( Solver % NewtonActive ) THEN
IF( ListGetLogical( Params,'Nonlinear System Reset Newton', Found) ) Solver % NewtonActive = .FALSE.
END IF


IF( ListGetLogical( Params,'Nonlinear System Nullify Guess', Found ) ) THEN
Solver % Variable % Values = 0.0_dp
END IF

! If we changed the system last time to harmonic one then revert back the real system
IF( ListGetLogical( Params,'Harmonic Mode',Found ) ) THEN
CALL ChangeToHarmonicSystem( Solver, .TRUE. )
Expand Down
52 changes: 40 additions & 12 deletions fem/src/ParallelUtils.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,8 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )
#ifdef PARALLEL_FOR_REAL
IF ( ParEnv % PEs <= 1 .OR. .NOT. ASSOCIATED(Matrix) ) RETURN

CALL Info('ParallelInitMatrix','Creating communication structures for matrix!',Level=15)

Mesh => Solver % Mesh
DOFs = Solver % Variable % DOFs

Expand All @@ -133,7 +135,21 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )

n = SIZE(Perm)
k = n*DOFs + Matrix % ExtraDOFs
ALLOCATE( Matrix % Perm(k), Matrix % InvPerm(k))

i = 0
IF(ASSOCIATED( Matrix % Perm) ) i=i+1
IF(ASSOCIATED( Matrix % InvPerm) ) i=i+1

IF(i==1) THEN
CALL Fatal('ParallelInitMatrix','Only Perm or InvPerm is associated!')
END IF
IF(i==2) THEN
CALL Info('ParallelInitMatrix','Skipping generation of Perm and InvPerm',Level=20)
GOTO 1
END IF

j = MAXVAL(Perm)*DOFs + Matrix % ExtraDOFs
ALLOCATE( Matrix % Perm(k), Matrix % InvPerm(j))

BLOCK
LOGICAL :: DoConf = .FALSE.
Expand All @@ -148,24 +164,29 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )
END IF

IF ( Perm(i) /= 0 ) THEN
DO j=1,DOFs
Matrix % Perm((i-1)*DOFs+j) = DOFs * (Perm(i)-1) + j
END DO
DO j=1,DOFs
Matrix % Perm((i-1)*DOFs+j) = DOFs * (Perm(i)-1) + j
END DO
END IF
END DO
END BLOCK
END BLOCK

DO i=n*DOFs+1,SIZE(Matrix % Perm)
Matrix % Perm(i) = i
END DO
DO i=n*DOFs+1,SIZE(Matrix % Perm)
Matrix % Perm(i) = i
END DO

Matrix % INVPerm = 0
DO i=1,SIZE(Matrix % Perm)
IF ( Matrix % Perm(i) /= 0 ) THEN
Matrix % InvPerm(Matrix % Perm(i)) = i
END IF
IF ( Matrix % Perm(i) /= 0 ) THEN
Matrix % InvPerm(Matrix % Perm(i)) = i
END IF
END DO

1 CONTINUE

IF(ASSOCIATED(Matrix % ParallelInfo)) THEN
CALL Fatal('ParallelInitMatrix','ParallelInfo already created!')
END IF
ALLOCATE( Matrix % ParallelInfo )

IF ( .NOT. Matrix % DGMatrix ) THEN
Expand Down Expand Up @@ -213,7 +234,6 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )
l_beg = Mesh % NumberOfNodes

n = Mesh % NumberOfEdges

edofs = Mesh % MaxEdgeDOFS
maxedofs = ParallelReduction(edofs,2)

Expand All @@ -230,6 +250,7 @@ SUBROUTINE ParallelInitMatrix( Solver, Matrix, inPerm )
l = DOFs*(l_beg + edofs*(i-1)+j-1)+m
l = Matrix % Perm(l)
IF(l==0) CYCLE

Matrix % ParallelInfo % GlobalDOFs(l) = &
DOFs*(g_beg+maxedofs*(Element % GelementIndex-1)+j-1)+m
Matrix % ParallelInfo % GInterface(l) = Mesh % ParallelInfo % EdgeInterface(i)
Expand Down Expand Up @@ -807,7 +828,14 @@ SUBROUTINE ParallelInitSolve( Matrix, x, b, r, Update )
!-------------------------------------------------------------------------------
LOGICAL :: Upd
#ifdef PARALLEL_FOR_REAL
IF(.NOT. ASSOCIATED(Matrix % ParMatrix)) THEN
CALL Fatal('ParallelInitSolve','ParMatrix not associated!')
END IF
ParEnv => Matrix % ParMatrix % ParEnv
IF(.NOT. ASSOCIATED(ParEnv)) THEN
CALL Fatal('ParallelInitSolve','ParEnv not associated!')
END IF

ParEnv % ActiveComm = Matrix % Comm
Upd = .TRUE.
IF ( PRESENT(Update) ) Upd=Update
Expand Down
8 changes: 4 additions & 4 deletions fem/src/SolveHypre.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -812,7 +812,7 @@ void STDCALLBULL FC_FUNC(solvehypre2,SOLVEHYPRE2)
/* which process number am I? */
MPI_Comm_rank(comm, &myid);

myid = 0;
if(0) myid = 0;

if( myid == 0 )
myverb = verbosity;
Expand Down Expand Up @@ -961,9 +961,9 @@ void STDCALLBULL FC_FUNC(solvehypre2,SOLVEHYPRE2)
for( i=0,k=0; i<local_size; i++ )
if ( owner[i] ) xvec[i] = txvec[k++];

if(myverb > 5) fprintf(stdout,"SolveHypre: required iterations %d (method %d) to norm %lg\n",
if(myverb > 5) fprintf(stdout,"SolveHypre: Required iterations %d (method %d) to norm %lg\n",
num_iterations,Container->hypre_method, final_res_norm);
if (myverb > 4) fprintf( stdout, "SolveHypre: solution time (method %d): %g\n",
if (myverb > 4) fprintf( stdout, "SolveHypre: Solution time (method %d): %g\n",
Container->hypre_method, realtime_()-st );
free( txvec );
free( rcols );
Expand Down Expand Up @@ -995,7 +995,7 @@ void STDCALLBULL FC_FUNC(solvehypre4,SOLVEHYPRE4)
hypre_pre = Container->hypre_method % 100;


if(myverb > 6 ) fprintf(stdout,"SolveHypre: Detroying Hypre solver structures!\n");
if(myverb > 10 ) fprintf(stdout,"SolveHypre: Detroying Hypre solver structures!\n");

/* Destroy Hypre preconditioner */
if ( hypre_pre == 1 ) {
Expand Down
8 changes: 8 additions & 0 deletions fem/src/SolverUtils.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15928,11 +15928,19 @@ END SUBROUTINE BlockSolveExt
CALL Info(Caller,'Eliminating constraints before going into block matrix!')
CALL EliminateLinearRestriction( A, bb, A % ConstraintMatrix, Acoll, Solver, .TRUE. )
CALL List_ToCRSMatrix(Acoll)

Acoll % Comm = A % Comm
Acoll % AddMatrix => A % AddMatrix
CALL ParallelInitMatrix(Solver, Acoll)

CALL BlockSolveExt( Acoll, x, Acoll % rhs, Solver )

CALL Info(Caller,'Freeing collection matrix after solution',Level=10)
NULLIFY( Acoll % AddMatrix )

CALL FreeMatrix(Acoll)
ParEnv => A % ParMatrix % ParEnv

Acoll => NULL()
END BLOCK
ELSE
Expand Down

0 comments on commit 1106c21

Please sign in to comment.