Choral
precond_mod.F90
Go to the documentation of this file.
1 
23 
24 
26 
29  use real_type
30  use algebra_set, only: shellsort_dec
31  use basic_tools
32  use csr_mod
33 
34  implicit none
35  private
36 
37 
38  ! %----------------------------------------%
39  ! | |
40  ! | PUBLIC DATA |
41  ! | |
42  ! %----------------------------------------%
43  public :: precond
44  public :: clear
46 
47  ! %----------------------------------------%
48  ! | |
49  ! | DERIVED TYPE |
50  ! | |
51  ! %----------------------------------------%
56  type precond
57 
59  integer :: type
60 
62  real(RP), dimension(:), allocatable :: invd
63 
65  type(csr) :: llt
66 
68  integer , dimension(:), allocatable :: perm
69 
71  integer , dimension(:), allocatable :: perminv
72 
74  character(20) :: name
75 
76  contains
77 
79  final :: precond_clear
80 
81  end type precond
82 
83  ! %----------------------------------------%
84  ! | |
85  ! | GENERIc SUBROUTINES |
86  ! | |
87  ! %----------------------------------------%
88  interface precond
89  module procedure precond_create
90  end interface precond
91 
93  interface clear
94  module procedure precond_clear
95  end interface clear
96 
97 contains
98 
99  subroutine precond_clear(pc)
100  type(precond), intent(inout) :: pc
101 
102  call clear(pc%LLT)
103  call freemem(pc%invD)
104  call freemem(pc%perm)
105  call freemem(pc%permInv)
106 
107  end subroutine precond_clear
108 
109 
120  function precond_create(K, type) result(pc)
121  type(precond) :: pc
122  type(csr), intent(inout) :: K
123  integer , intent(in) :: type
124 
125  call clear(pc)
126  if(.NOT.valid(k)) call quit("precond_mod: precond_create:&
127  & matrix 'K' not valid")
128 
129  pc%type = type
130 
131  !! Name definition
132  !!
133  select case(type)
134 
135  case(pc_0)
136  pc%name = "void"
137 
138  case(pc_jacobi)
139  pc%name = "Jacobi"
140 
141  case(pc_icc0)
142  pc%name = "iCC0"
143 
144  case default
145  call quit( 'precond_mod: precond_create:&
146  & preconditioning type incorrect')
147 
148  end select
149 
150  if (choral_verb>0) write(*,*) &
151  &'precond_mod : create = ',&
152  & trim(pc%name)
153 
154  !! assembling
155  !!
156  select case(type)
157 
158  case(pc_jacobi)
159  call jacobi(pc%invD, k)
160 
161  case(pc_icc0)
162  call reverse_cuthill_mackey(pc%perm, pc%permInv, k)
163  call reorder(k, pc%perm, pc%permInv)
164  call icc0(pc%LLT, k)
165 
166  end select
167 
168  end function precond_create
169 
172  subroutine jacobi(invD, s)
173  real(RP) , dimension(:), allocatable :: invD
174  type(csr) , intent(inout) :: s
175 
176  integer :: ii, jj
177 
178  if (choral_verb>1) write(*,*) &
179  & 'precond_mod : jacobi'
180  call getdiag(s)
181 
182  !! check diagonal entries
183 #if DBG>0
184  if (minval(s%diag) <= 0 ) call quit( &
185  & 'precond_mod: jacobi: negative diagonal entries')
186 #endif
187 
188  call allocmem(invd, s%nl)
189  do ii=1, s%nl
190  jj = s%diag(ii)
191  invd(ii) = 1._rp / s%a(jj)
192  end do
193 
194  end subroutine jacobi
195 
196 
197 
200  subroutine reverse_cuthill_mackey(perm, permInv, s)
201  type(csr) , intent(in) :: s
202  integer , dimension(:), allocatable :: perm, permInv
203 
204  integer, dimension(s%nl) :: deg, P, LS, LS0, degLS, new_i
205  logical, dimension(s%nl) :: marked
206  integer :: ii, jj, jDeb, jFin, kk
207  integer :: nL, nL0, node
208 
209  if (choral_verb>1) write(*,*) &
210  & 'precond_mod : reverse_cuthill_macKey'
211 
212  call allocmem(perm, s%nl)
213  call allocmem(perminv, s%nl)
214 
215  ! Initialisation
216  marked = .false.
217 
218  ! Degre des sommets
219  do ii = 1,s%nl
220  deg(ii) = s%row(ii+1)-s%row(ii)
221  end do
222 
223  ! Initialisation de l'algo
224  kk = 1
225  node = 1
226  jj = minval(deg)
227  do ii=1, s%nl
228  if (deg(ii)==jj) then
229  node = ii ! Premier noeud, degre minimal
230  exit
231  end if
232  end do
233  marked(node) = .true. ! Marque
234 
235  nl0 = 1 ! taille de LS0
236  ls0(1) = node ! Premiere Level-set, 1 elt
237 
238  p(kk) = node ! Mise à jour de la permutation
239 
240 
241  do while ( kk<s%nl )
242  ! Construction de la level-set suivant LS0
243 
244  nl = 0
245  do ii = 1, nl0
246 
247  jdeb = s%row( ls0(ii) )
248  jfin = s%row( ls0(ii)+1 ) - 1
249 
250  do jj = jdeb,jfin
251 
252  node = s%col(jj)
253 
254  if (.not.marked(node)) then
255 
256  nl = nl+1
257  ls(nl) = node
258  marked(node) = .true. ! Marquage du noeud
259 
260  end if
261 
262  end do
263  end do
264 
265  ! actualisation NL0, LS0
266  nl0 = nl
267  ls0(1:nl) = ls(1:nl)
268 
269  do ii = 1,nl
270  degls(ii) = deg(ls(ii)) ! degLS := degre(LS)
271  end do
272 
273  ! On parcourt LS par ordre croissant de degre
274  call shellsort_dec(degls(1:nl), new_i(1:nl), nl) ! Tri ordre decroissant
275  do ii = nl,1,-1
276  kk = kk + 1 ! Noeud suivant
277  p(kk) = ls(new_i(ii)) ! MAJ permutation
278  end do
279 
280  end do
281 
282  do ii = 1,s%nl
283  perm(ii) = p(s%nl+1-ii)
284  end do
285 
286  do ii = 1, s%nl
287  perminv(perm(ii)) = ii
288  end do
289 
290  end subroutine reverse_cuthill_mackey
291 
292 
293 
296  subroutine icc0(LLT, s)
297  type(csr), intent(out) :: LLT
298  type(csr), intent(inout) :: s
299 
300  integer :: i,k,j,m,mm,band, kDeb,kFin
301  real(RP) :: int1,int2,sum, t1, t2
302 
303  if (choral_verb>1) write(*,'(A35,$)') &
304  &' precond_mod : icc0 ='
305  call cpu_time(t1)
306 
307  if (.not.s%sorted) call sortrows(s)
308  call getdiag(s)
309  call copy(llt, s)
310 
311  ! Calcul de la largeur de bande
312  band = 0
313  do i = 1,s%nl
314  band = max( band, s%col(s%row(i+1)-1)-i )
315  end do
316  if (choral_verb>1) write(*,*) &
317  & int(band,4), ' : band width'
318 
319 
320  ! Calcul de la decomposition de Choleski LL^T incomplete. C'est a
321  ! dire que l'on ne rempli L (et L^T) que la ou S avait un
322  ! coefficient non nul. Dans l'algo ci-dessous on raisonne par
323  ! ligne, c'est donc la partie L^T (triangulaire superieure) qui
324  ! est remplie en premier.
325  do i = 1,s%nl
326  do m = s%diag(i),s%row(i+1)-1 ! Termes sup-diagonaux
327  j = s%col(m) ! j>=i
328  sum = s%a(m)
329  do k = max(i-band-1,1),i-1 ! Lignes precedentes, k = 1 a i-1
330  int1 = 0._rp
331  int2 = 0._rp
332  do mm = s%diag(k),s%row(k+1)-1 ! On ne cherche que parmi les sup-diagonaux
333  if (s%col(mm)==i) int1 = llt%a(mm) ! C'est LT_{ki}, k<i
334  if (s%col(mm)==j) int2 = llt%a(mm) ! C'est LT_{kj}, k<i<=j
335  end do
336  sum = sum - int1*int2
337  end do
338  if (j==i) then
339  llt%a(m) = sqrt(sum)
340  else
341  llt%a(m) = sum / llt%a(s%diag(i))
342  end if
343  end do
344  end do
345 
346  ! On recopie les coefficients dans la partie triangulaire
347  ! inferieure pour construire L.
348  do i = 1,llt%nl-1 ! Parcours des lignes 1 a N-1
349  kdeb = llt%diag(i)+1
350  kfin = llt%row(i+1)-1
351  do k = kdeb,kfin ! Parcours des colonnes j de la ligne i avec i<j
352  j = llt%col(k)
353  ! Il faut recopier le coeff L_{ij} dans L_{ji}
354  do m = llt%row(j), llt%diag(j)-1 ! Parcours de la ligne j colonnes m tq m<j
355  if (llt%col(m)==i) then ! On a trouve l'indice de la colonne _i_
356  llt%a(m) = llt%a(k) ! On fait LLT_{ji} := LLT_{ij}
357  exit ! On sort de la boucle sur _m_
358  end if
359  end do
360  end do
361  end do
362 
363  call cpu_time(t2)
364 
365  end subroutine icc0
366 
367 
368 end module precond_mod
369 
subroutine precond_clear(pc)
subroutine, public icc0(LLT, s)
Incomplete Choleski decomposition with 0 fill-in for s.
M2 = M1.
Definition: csr_mod.F90:161
subroutine, public reverse_cuthill_mackey(perm, permInv, s)
Reverse CutHill Mac Kee algorithm.
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
The type precond defines preconditioning for linear systems.
Definition: precond_mod.F90:56
type(precond) function precond_create(K, type)
Constructor for the type precond
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public jacobi(invD, s)
Jacobi preconditioner.
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
CHORAL CONSTANTS
Derived type for sparse matrices in CSR format.
Definition: csr_mod.F90:62
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
integer, parameter pc_jacobi
Jacobi diagonal preconditioning.
allocate memory for real(RP) arrays
Definition: real_type.F90:158
subroutine, public getdiag(A)
Computes the array Adiag of the indices in Acol for the diagonal coefficients.
Definition: csr_mod.F90:835
reorder rows and cols of a matrix following a permutation of the indexes
Definition: csr_mod.F90:156
integer, parameter pc_icc0
Incomplete Cholesky, order 0.
integer choral_verb
Verbosity level.
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
DERIVED TYPE precond: preconditioning for linear systems
Definition: precond_mod.F90:25
subroutine, public shellsort_dec(t, new_i, n)
Sort integer array : shell sort.
DERIVED TYPE csr for sparse matrices
Definition: csr_mod.F90:13
integer, parameter pc_0
void preconditioning