|
|
|
@ -44,11 +44,11 @@
|
|
|
|
|
!
|
|
|
|
|
! This routine computes
|
|
|
|
|
!
|
|
|
|
|
! Y = beta*Y + alpha*op(M^(-1))*X,
|
|
|
|
|
! Y = beta*Y + alpha*op(ML^(-1))*X,
|
|
|
|
|
! where
|
|
|
|
|
! - M is a multilevel domain decomposition (Schwarz) preconditioner associated
|
|
|
|
|
! - ML is a multilevel preconditioner associated
|
|
|
|
|
! to a certain matrix A and stored in p,
|
|
|
|
|
! - op(M^(-1)) is M^(-1) or its transpose, according to the value of trans,
|
|
|
|
|
! - op(ML^(-1)) is ML^(-1) or its transpose, according to the value of trans,
|
|
|
|
|
! - X and Y are vectors,
|
|
|
|
|
! - alpha and beta are scalars.
|
|
|
|
|
!
|
|
|
|
@ -59,10 +59,14 @@
|
|
|
|
|
! A multilevel preconditioner is regarded as an array of 'one-level' data structures,
|
|
|
|
|
! each containing the part of the preconditioner associated to a certain level
|
|
|
|
|
! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90).
|
|
|
|
|
! For each level ilev, the 'base preconditioner' K(ilev) is stored in
|
|
|
|
|
! p%precv(ilev)%prec
|
|
|
|
|
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
|
|
|
|
|
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
|
|
|
|
|
! For each level lev, there is a smoother stored in
|
|
|
|
|
! p%precv(lev)%sm
|
|
|
|
|
! which in turn contains a solver
|
|
|
|
|
! p$precv(lev)%sm%sv
|
|
|
|
|
! Typically the solver acts only locally, and the smoother applies any required
|
|
|
|
|
! parallel communication/action.
|
|
|
|
|
! Each level has a matrix A(lev), obtained by 'tranferring' the original
|
|
|
|
|
! matrix A (i.e. the matrix to be preconditioned) to the level lev, through smoothed
|
|
|
|
|
! aggregation.
|
|
|
|
|
!
|
|
|
|
|
! The levels are numbered in increasing order starting from the finest one, i.e.
|
|
|
|
@ -85,35 +89,27 @@
|
|
|
|
|
! The multilevel preconditioner data structure containing the
|
|
|
|
|
! local part of the preconditioner to be applied.
|
|
|
|
|
! Note that nlev = size(p%precv) = number of levels.
|
|
|
|
|
! p%precv(ilev)%prec - type(psb_cbaseprec_type)
|
|
|
|
|
! The 'base preconditioner' for the current level
|
|
|
|
|
! p%precv(ilev)%ac - type(psb_cspmat_type)
|
|
|
|
|
! The local part of the matrix A(ilev).
|
|
|
|
|
! p%precv(ilev)%desc_ac - type(psb_desc_type).
|
|
|
|
|
! p%precv(lev)%sm - type(psb_cbaseprec_type)
|
|
|
|
|
! The 'smoother' for the current level
|
|
|
|
|
! p%precv(lev)%ac - type(psb_cspmat_type)
|
|
|
|
|
! The local part of the matrix A(lev).
|
|
|
|
|
! p%precv(lev)%parms - type(psb_sml_parms)
|
|
|
|
|
! Parameters controllin the multilevel prec.
|
|
|
|
|
! p%precv(lev)%desc_ac - type(psb_desc_type).
|
|
|
|
|
! The communication descriptor associated to the sparse
|
|
|
|
|
! matrix A(ilev)
|
|
|
|
|
! p%precv(ilev)%map - type(psb_inter_desc_type)
|
|
|
|
|
! Stores the linear operators mapping level (ilev-1)
|
|
|
|
|
! to (ilev) and vice versa. These are the restriction
|
|
|
|
|
! matrix A(lev)
|
|
|
|
|
! p%precv(lev)%map - type(psb_inter_desc_type)
|
|
|
|
|
! Stores the linear operators mapping level (lev-1)
|
|
|
|
|
! to (lev) and vice versa. These are the restriction
|
|
|
|
|
! and prolongation operators described in the sequel.
|
|
|
|
|
! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable.
|
|
|
|
|
! The integer parameters defining the multilevel
|
|
|
|
|
! strategy
|
|
|
|
|
! p%precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable.
|
|
|
|
|
! The real parameters defining the multilevel strategy
|
|
|
|
|
! p%precv(ilev)%mlia - integer, dimension(:), allocatable.
|
|
|
|
|
! The aggregation map (ilev-1) --> (ilev).
|
|
|
|
|
! p%precv(ilev)%nlaggr - integer, dimension(:), allocatable.
|
|
|
|
|
! The number of aggregates (rows of A(ilev)) on the
|
|
|
|
|
! various processes.
|
|
|
|
|
! p%precv(ilev)%base_a - type(psb_cspmat_type), pointer.
|
|
|
|
|
! p%precv(lev)%base_a - type(psb_cspmat_type), pointer.
|
|
|
|
|
! Pointer (really a pointer!) to the base matrix of
|
|
|
|
|
! the current level, i.e. the local part of A(ilev);
|
|
|
|
|
! the current level, i.e. the local part of A(lev);
|
|
|
|
|
! so we have a unified treatment of residuals. We
|
|
|
|
|
! need this to avoid passing explicitly the matrix
|
|
|
|
|
! A(ilev) to the routine which applies the
|
|
|
|
|
! A(lev) to the routine which applies the
|
|
|
|
|
! preconditioner.
|
|
|
|
|
! p%precv(ilev)%base_desc - type(psb_desc_type), pointer.
|
|
|
|
|
! p%precv(lev)%base_desc - type(psb_desc_type), pointer.
|
|
|
|
|
! Pointer to the communication descriptor associated
|
|
|
|
|
! to the sparse matrix pointed by base_a.
|
|
|
|
|
!
|
|
|
|
@ -134,11 +130,9 @@
|
|
|
|
|
! info - integer, output.
|
|
|
|
|
! Error code.
|
|
|
|
|
!
|
|
|
|
|
! Note that when the LU factorization of the matrix A(ilev) is computed instead of
|
|
|
|
|
! Note that when the LU factorization of the matrix A(lev) is computed instead of
|
|
|
|
|
! the ILU one, by using UMFPACK or SuperLU, the corresponding L and U factors
|
|
|
|
|
! are stored in data structures provided by UMFPACK or SuperLU and pointed by
|
|
|
|
|
! p%precv(ilev)%prec%iprcparm(mld_umf_ptr) or p%precv(ilev)%prec%iprcparm(mld_slu_ptr),
|
|
|
|
|
! respectively.
|
|
|
|
|
! are stored in data structures provided by UMFPACK or SuperLU.
|
|
|
|
|
!
|
|
|
|
|
! This routine is formulated in a recursive way, so it is very compact.
|
|
|
|
|
! In the original code the recursive formulation was explicitly unrolled.
|
|
|
|
@ -147,7 +141,7 @@
|
|
|
|
|
! recursive formulation.
|
|
|
|
|
!
|
|
|
|
|
! This routine computes
|
|
|
|
|
! Y = beta*Y + alpha*op(M^(-1))*X,
|
|
|
|
|
! Y = beta*Y + alpha*op(ML^(-1))*X,
|
|
|
|
|
! where
|
|
|
|
|
! - M is a multilevel domain decomposition (Schwarz) preconditioner
|
|
|
|
|
! associated to a certain matrix A and stored in p,
|
|
|
|
@ -162,138 +156,74 @@
|
|
|
|
|
! The multilevel preconditioner is regarded as an array of 'one-level' data structures,
|
|
|
|
|
! each containing the part of the preconditioner associated to a certain level
|
|
|
|
|
! (for more details see the description of mld_Tonelev_type in mld_prec_type.f90).
|
|
|
|
|
! For each level ilev, the 'base preconditioner' K(ilev) is stored in
|
|
|
|
|
! p%precv(ilev)%prec
|
|
|
|
|
! and is associated to a matrix A(ilev), obtained by 'tranferring' the original
|
|
|
|
|
! matrix A (i.e. the matrix to be preconditioned) to the level ilev, through smoothed
|
|
|
|
|
! For each level ilev, the 'base preconditioner' K(lev) is stored in
|
|
|
|
|
! p%precv(lev)%prec
|
|
|
|
|
! and is associated to a matrix A(lev), obtained by 'tranferring' the original
|
|
|
|
|
! matrix A (i.e. the matrix to be preconditioned) to the level lev, through smoothed
|
|
|
|
|
! aggregation.
|
|
|
|
|
! The levels are numbered in increasing order starting from the finest one, i.e.
|
|
|
|
|
! level 1 is the finest level and A(1) is the matrix A.
|
|
|
|
|
!
|
|
|
|
|
! This routine applies one of the following multilevel strategies:
|
|
|
|
|
!
|
|
|
|
|
! Additive multilevel
|
|
|
|
|
! This is additive both within the levels and among levels.
|
|
|
|
|
!
|
|
|
|
|
! For details on the additive multilevel Schwarz preconditioner, see
|
|
|
|
|
! Algorithm 3.1.1 in the book:
|
|
|
|
|
! B.F. Smith, P.E. Bjorstad & W.D. Gropp,
|
|
|
|
|
! Domain decomposition: parallel multilevel methods for elliptic partial
|
|
|
|
|
! differential equations, Cambridge University Press, 1996.
|
|
|
|
|
!
|
|
|
|
|
! (P(ilev) denotes the smoothed prolongator from level ilev to level
|
|
|
|
|
! ilev-1, while PT(ilev) denotes its transpose, i.e. the corresponding
|
|
|
|
|
! restriction operator from level ilev-1 to level ilev).
|
|
|
|
|
!
|
|
|
|
|
! 1. Transfer the outer vector Xest to x(1) (inner X at level 1)
|
|
|
|
|
!
|
|
|
|
|
!
|
|
|
|
|
! 2. Apply the base preconditioner at the current level:
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev)
|
|
|
|
|
! y(ilev) = (K(ilev)^(-1))*x(ilev)
|
|
|
|
|
! - Additive multilevel
|
|
|
|
|
! - V-cycle
|
|
|
|
|
! - W-cycle
|
|
|
|
|
! - K-cycle
|
|
|
|
|
!
|
|
|
|
|
! 3. If ilev < nlevel
|
|
|
|
|
! a. Transfer x(ilev) to the next level:
|
|
|
|
|
! x(ilev+1) = PT(ilev+1)*x(ilev)
|
|
|
|
|
! b. Call recursively itself
|
|
|
|
|
! c. Transfer y(ilev+1) to the current level:
|
|
|
|
|
! y(ilev) = y(ilev) + P(ilev+1)*y(ilev+1)
|
|
|
|
|
! For details of the algorithms, see:
|
|
|
|
|
!
|
|
|
|
|
! 4. if ilev == 1 Transfer the inner y to the external:
|
|
|
|
|
! Yext = beta*Yext + alpha*y(1)
|
|
|
|
|
!
|
|
|
|
|
!
|
|
|
|
|
!
|
|
|
|
|
! Hybrid multiplicative---pre-smoothing
|
|
|
|
|
!
|
|
|
|
|
! The preconditioner M is hybrid in the sense that it is multiplicative through the
|
|
|
|
|
! levels and additive inside a level.
|
|
|
|
|
!
|
|
|
|
|
! For details on the pre-smoothed hybrid multiplicative multilevel Schwarz
|
|
|
|
|
! preconditioner, see Algorithm 3.2.1 in the book:
|
|
|
|
|
! B.F. Smith, P.E. Bjorstad & W.D. Gropp,
|
|
|
|
|
! - B.F. Smith, P.E. Bjorstad & W.D. Gropp,
|
|
|
|
|
! Domain decomposition: parallel multilevel methods for elliptic partial
|
|
|
|
|
! differential equations, Cambridge University Press, 1996.
|
|
|
|
|
!
|
|
|
|
|
! - W. L. Briggs, V. E. Henson, S. F. McCormick,
|
|
|
|
|
! A Multigrid Tutorial, Second Edition
|
|
|
|
|
! SIAM, 2000.
|
|
|
|
|
!
|
|
|
|
|
! 1 Transfer the outer vector Xest to x(1) (inner X at level 1)
|
|
|
|
|
! - Y. Notay, P. S. Vassilevski,
|
|
|
|
|
! Recursive Krylov-based multigrid cycles
|
|
|
|
|
! Numerical Linear Algebra with Applications, 15 (5), 2008, 473--487.
|
|
|
|
|
!
|
|
|
|
|
! 2. Apply the base preconditioner at the current level:
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev).
|
|
|
|
|
! y(ilev) = (K(ilev)^(-1))*x(ilev)
|
|
|
|
|
!
|
|
|
|
|
! 3. If ilev < nlevel
|
|
|
|
|
! a. Compute the residual:
|
|
|
|
|
! r(ilev) = x(ilev) - A(ilev)*y(ilev)
|
|
|
|
|
! b. Transfer r(ilev) to the next level:
|
|
|
|
|
! x(ilev+1) = PT(ilev+1)*r(ilev)
|
|
|
|
|
! c. Call recursively
|
|
|
|
|
! d. Transfer y(ilev+1) to the current level:
|
|
|
|
|
! y(ilev) = y(ilev) + P(ilev+1)*y(ilev+1)
|
|
|
|
|
! The V-cycle can be described as follows, where
|
|
|
|
|
! P(lev) denotes the smoothed prolongator from level lev to level
|
|
|
|
|
! lev-1, while R(lev) denotes the corresponding restriction operator
|
|
|
|
|
! (normally its transpose) from level lev-1 to level lev.
|
|
|
|
|
! M(lev) is the smoother at the current level.
|
|
|
|
|
!
|
|
|
|
|
! 4. if ilev == 1 Transfer the inner y to the external:
|
|
|
|
|
! Yext = beta*Yext + alpha*y(1)
|
|
|
|
|
! In the code below, the recursive procedure is inner_ml_aply, which
|
|
|
|
|
! in turn makes use of mld_inner_mult (for V-cycle) or similar for
|
|
|
|
|
! the other cycles.
|
|
|
|
|
!
|
|
|
|
|
! 1. Transfer the outer vector Xest to u(1) (inner X at level 1)
|
|
|
|
|
!
|
|
|
|
|
! 2. Invoke V-=cycle(1,M,P,R,A,b,u)
|
|
|
|
|
!
|
|
|
|
|
! Hybrid multiplicative, post-smoothing variant
|
|
|
|
|
! procedure V-cycle(lev,M,P,R,A,b,u)
|
|
|
|
|
! if (lev < nlev) then
|
|
|
|
|
!
|
|
|
|
|
! 1. Transfer the outer vector Xest to x(1) (inner X at level 1)
|
|
|
|
|
! u(lev) = u(lev) + M(lev)*(b(lev)-A(lev)*u(lev)
|
|
|
|
|
!
|
|
|
|
|
! 2. If ilev < nlev
|
|
|
|
|
! a. Transfer x(ilev) to the next level:
|
|
|
|
|
! x(ilev+1) = PT(ilev+1)*x(ilev)
|
|
|
|
|
! b. Call recursively
|
|
|
|
|
! c. Transfer y(ilev+1) to the current level:
|
|
|
|
|
! y(ilev) = P(ilev+1)*y(ilev+1)
|
|
|
|
|
! d. Compute the residual:
|
|
|
|
|
! x(ilev) = x(ilev) - A(ilev)*y(ilev)
|
|
|
|
|
! e. Apply the base preconditioner to the residual at the current level:
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev)
|
|
|
|
|
! y(ilev) = y(ilev) + (K(ilev)^(-1))*x(ilev)
|
|
|
|
|
! b(lev+1) = R(lev+1)*(b(lev)-A(lev)*u(lev)
|
|
|
|
|
!
|
|
|
|
|
! 3. If ilev == nlev apply y(ilev) = (K(ilev)^(-1))*x(ilev)
|
|
|
|
|
! u(lev+1) = V-cycle(lev+1,M,P,R,A,b,u)
|
|
|
|
|
!
|
|
|
|
|
! 4. if ilev == 1 Transfer the inner Y to the external:
|
|
|
|
|
! Yext = beta*Yext + alpha*Y(1)
|
|
|
|
|
! u(lev) = u(lev) + P(lev+1) * u(lev+1)
|
|
|
|
|
!
|
|
|
|
|
! u(lev) = u(lev) + M(lev)*(b(lev)-A(lev)*u(lev)
|
|
|
|
|
!
|
|
|
|
|
! else
|
|
|
|
|
!
|
|
|
|
|
! Hybrid multiplicative, pre- and post-smoothing (two-side) variant
|
|
|
|
|
! solve A(lev)*u(lev) = b(lev)
|
|
|
|
|
!
|
|
|
|
|
! For details on the symmetrized hybrid multiplicative multilevel Schwarz
|
|
|
|
|
! preconditioner, see Algorithm 3.2.2 in the book:
|
|
|
|
|
! B.F. Smith, P.E. Bjorstad & W.D. Gropp,
|
|
|
|
|
! Domain decomposition: parallel multilevel methods for elliptic partial
|
|
|
|
|
! differential equations, Cambridge University Press, 1996.
|
|
|
|
|
! end if
|
|
|
|
|
!
|
|
|
|
|
! return u(lev)
|
|
|
|
|
! end
|
|
|
|
|
!
|
|
|
|
|
! 1. Transfer the outer vector Xest to x(1) (inner X at level 1)
|
|
|
|
|
!
|
|
|
|
|
! 2. Apply the base preconditioner at the current level:
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev)
|
|
|
|
|
! y(ilev) = (K(ilev)^(-1))*x(ilev)
|
|
|
|
|
!
|
|
|
|
|
! 3. If ilev < nlevel
|
|
|
|
|
! a. Compute the residual:
|
|
|
|
|
! r(ilev) = x(ilev) - A(ilev)*y(ilev)
|
|
|
|
|
! b. Transfer r(ilev) to the next level:
|
|
|
|
|
! x(ilev+1) = PT(ilev+1)*r(ilev)
|
|
|
|
|
! c. Call recursively
|
|
|
|
|
! d. Transfer y(ilev+1) to the current level:
|
|
|
|
|
! y(ilev) = y(ilev) + P(ilev+1)*y(ilev+1)
|
|
|
|
|
! d. Compute the residual:
|
|
|
|
|
! r(ilev) = x(ilev) - A(ilev)*y(ilev)
|
|
|
|
|
! e. Apply the base preconditioner at the current level to the residual:
|
|
|
|
|
! ! The sum over the subdomains is carried out in the
|
|
|
|
|
! ! application of K(ilev)
|
|
|
|
|
! y(ilev) = y(ilev) + (K(ilev)^(-1))*r(ilev)
|
|
|
|
|
!
|
|
|
|
|
! 4. if ilev == 1 Transfer the inner Y to the external:
|
|
|
|
|
! Yext = beta*Yext + alpha*Y(1)
|
|
|
|
|
! 3. Transfer u(1) to the external:
|
|
|
|
|
! Yext = beta*Yext + alpha*u(1)
|
|
|
|
|
!
|
|
|
|
|
!
|
|
|
|
|
subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
|
|
|
|
|