Choral
mumps_mod.F90
Go to the documentation of this file.
1 
29 
30 
31 module mumps_mod
32 
35  use real_type
36  use basic_tools
37  use csr_mod
38 
39  implicit none
40  private
41 
42 #ifdef WMUMPS
43  include 'mpif.h'
44  include 'dmumps_struc.h'
45 #endif
46 
47  public :: mumps_init
48  public :: mumps, clear, solve
49 
50  ! %----------------------------------------%
51  ! | |
52  ! | DERIVED TYPE |
53  ! | |
54  ! %----------------------------------------%
58  type mumps
60  logical :: created = .false.
61 
62 #ifdef WMUMPS
63  type(dmumps_struc) :: dmumps
64 #endif
65 
66  contains
67 
69  final :: mumps_clear
70 
71  end type mumps
72 
73  ! %----------------------------------------%
74  ! | |
75  ! | GENERIc SUBROUTINES |
76  ! | |
77  ! %----------------------------------------%
79  interface clear
80  module procedure mumps_clear
81  end interface clear
82 
83  interface mumps
84  module procedure mumps_create
85  end interface mumps
86 
88  interface solve
89  module procedure mumps_solve
90  end interface solve
91 
92 contains
93 
94 
96  subroutine mumps_init()
97 #ifdef WMUMPS
98  logical :: ierr
99 
100  if (choral_verb>1) write(*,*)"mumps_mod : mumps_init"
101 
102  if (rp /= dp) call quit(&
103  & "mumps_mod: mumps_init: double precision only")
104 
105  call mpi_init(ierr)
106 #endif
107  end subroutine mumps_init
108 
109 
111  subroutine mumps_clear(lu)
112  type(mumps), intent(inout) :: lu
113 
114 #ifdef WMUMPS
115 
116  if (lu%created) then
117 
118  if (associated(lu%dmumps%irn)) deallocate(lu%dmumps%irn)
119  lu%dmumps%irn => null()
120  if (associated(lu%dmumps%jcn)) deallocate(lu%dmumps%jcn)
121  lu%dmumps%jcn => null()
122  if (associated(lu%dmumps%rhs)) deallocate(lu%dmumps%rhs)
123  lu%dmumps%rhs => null()
124  if (associated(lu%dmumps%a )) deallocate(lu%dmumps%a)
125  lu%dmumps%a => null()
126 
127  ! Destroy the instance (deallocate internal data structures)
128  lu%dmumps%job = -2
129  call dmumps(lu%dmumps)
130  end if
131 
132 #endif
133 
134  lu%created = .false.
135 
136  end subroutine mumps_clear
137 
138 
139 
140 
152  function mumps_create(g, type, verb) result(lu)
153  type(mumps) :: lu
154  type(csr), intent(in) :: g
155  integer , intent(in) :: type
156  integer , intent(in), optional :: verb
157 
158  integer :: ii, jj, cpt
159  real(DP) :: cpu
160 
161  cpu = clock()
162 
163  if (choral_verb>0) write(*,*)"mumps_mod : create"
164  if(.NOT.valid(g)) call quit( &
165  & "mumps_mod: mumps_create: CSR matrix 'g' not valid" )
166  call clear(lu)
167 
168 #ifdef WMUMPS
169 
170  !!
171  !! MUMPS : INITIALISATION
172  !!
173  if (choral_verb>1) write(*,*)" mumps structure initialisation"
174  lu%dmumps%comm = mpi_comm_world
175  lu%created = .true.
176  lu%dmumps%job = -1
177 
178  select case(type)
179  case(mumps_lu)
180  lu%dmumps%sym = 0
181  if (choral_verb>1) write(*,*)" Non symmetric CSR"
182 
183  case(mumps_ldlt)
184  lu%dmumps%sym = 2
185  if (choral_verb>1) write(*,*)" Symmetric CSR"
186 
187  case(mumps_ldlt_sdp)
188  lu%dmumps%sym = 1
189  if (choral_verb>1) write(*,*)" Symmetric positive definite CSR"
190 
191  case default
192  call quit("mumps_mod: mumps_create: &
193  & unknown parameter 'type' " )
194 
195  end select
196 
197  lu%dmumps%par = 1
198  call dmumps(lu%dmumps)
199  !!
200  !! check error
201  if (lu%dmumps%infog(1)<0) then
202  if (choral_verb>0) then
203  write(*,*) " dmumps error, code =", &
204  & lu%dmumps%infog(1), lu%dmumps%infog(2)
205 
206  end if
207  call quit("mumps_mod: mumps_create: init. error")
208  end if
209  !!
210  !! transfer the CSR matrix entries
211  !!
212  if (lu%dmumps%myid==0) then
213 
214  !! number of lines
215  lu%dmumps%n = g%nl
216 
217  !! symmetric case : count entries a_ij with j <= i only
218  if (lu%dmumps%sym > 0) then
219  cpt = 0
220  do ii=1, g%nl
221  do jj=g%row(ii), g%row(ii+1) - 1
222  if (g%col(jj) <= ii ) cpt = cpt + 1
223  end do
224  end do
225  lu%dmumps%nnz = cpt
226 
227  else !! non symmetric case, all entries are counted
228  lu%dmumps%nnz = g%nnz
229  end if
230 
231  allocate(lu%dmumps%rhs(g%nl ))
232  allocate(lu%dmumps%irn(g%nnz))
233  allocate(lu%dmumps%jcn(g%nnz))
234  allocate(lu%dmumps%a( g%nnz))
235 
236  !! fill in the entries
237  cpt = 0
238  if (lu%dmumps%sym > 0) then !! symmetric case
239 
240  do ii=1, g%nl
241 
242  do jj=g%row(ii), g%row(ii+1) - 1
243 
244  if (g%col(jj) <= ii ) then
245  cpt = cpt + 1
246 
247  lu%dmumps%irn(cpt) = ii
248  lu%dmumps%jcn(cpt) = g%col(jj)
249  lu%dmumps%a(cpt) = g%a(jj)
250 
251  end if
252 
253  end do
254  end do
255 
256  else !! non symmetric case
257 
258  do ii=1, g%nl
259 
260  do jj=g%row(ii), g%row(ii+1) - 1
261  cpt = cpt + 1
262 
263  lu%dmumps%irn(cpt) = ii
264  lu%dmumps%jcn(cpt) = g%col(jj)
265  lu%dmumps%a(cpt) = g%a(jj)
266 
267  end do
268  end do
269 
270  end if
271  end if
272 
273  !!
274  !! MUMPS VERBOSITY
275  !!
276  lu%dmumps%icntl(4) = -1
277  !!
278  if (present(verb)) then
279  lu%dmumps%icntl(4) = verb
280  end if
281 
282  !!
283  !! MUMPS : ANALYSIS PHASE
284  !!
285  if (choral_verb>1) write(*,*)" mumps analysis phase"
286  lu%dmumps%job = 1
287  call dmumps(lu%dmumps)
288  !!
289  !! check error
290  if (lu%dmumps%infog(1)<0) then
291  if (choral_verb>0) then
292  write(*,*) " dmumps error, code =", &
293  & lu%dmumps%infog(1), lu%dmumps%infog(2)
294 
295  end if
296  call quit("mumps_mod: mumps_create: analysis phase error")
297  end if
298 
299  !!
300  !! MUMPS : FACTORISATION PHASE
301  !!
302  if (choral_verb>1) write(*,*)" mumps factorisation phase"
303  lu%dmumps%job = 2
304  call dmumps(lu%dmumps)
305  !!
306  !! check error
307  if (lu%dmumps%infog(1)<0) then
308  if (choral_verb>0) then
309  write(*,*) " dmumps error, code =", &
310  & lu%dmumps%infog(1), lu%dmumps%infog(2)
311 
312  end if
313  call quit("mumps_mod: mumps_create: factorisation phase error")
314  end if
315 
316 #else
317  call quit("mumps_mod: mumps_create: MUMPS not available")
318 
319 #endif
320 
321  cpu = clock() - cpu
322  if (choral_verb>1) write(*,*) &
323  & " CPU =", real(cpu, sp)
324 
325  end function mumps_create
326 
327 
330  subroutine mumps_solve(x, LU, b)
331  real(RP), dimension(:), intent(inout) :: x
332  type(mumps) , intent(inout) :: LU
333  real(RP), dimension(:), intent(in) :: b
334 
335 #ifdef WMUMPS
336 
337  lu%dmumps%rhs = b
338  lu%dmumps%job = 3
339  call dmumps(lu%dmumps)
340  !!
341  !! check error
342  if (lu%dmumps%infog(1)<0) then
343 
344  write(*,*) "mumps_mod : mumps_solve = error ", &
345  & lu%dmumps%infog(1), lu%dmumps%infog(2)
346  call quit("mumps_mod: mumps_solve; error detected")
347  end if
348 
349  !! retrieve solution
350  x = lu%dmumps%rhs
351 
352 #else
353  call quit("mumps_mod: mumps_solve: MUMPS not available")
354 
355 #endif
356 
357  end subroutine mumps_solve
358 
359 end module mumps_mod
is the csr matrix well defined ?
Definition: csr_mod.F90:101
integer, parameter mumps_ldlt
MUMPS LDLT , symm mat.
subroutine mumps_solve(x, LU, b)
SOLVE LU x = b.
Definition: mumps_mod.F90:331
BASIC TOOLS
Definition: basic_tools.F90:8
double precision function, public clock()
Returns the internal clock time
Definition: basic_tools.F90:94
integer, parameter, public rp
real(kind=RP) = real precision in the code REAL_TOL = epsilon to test real equality ...
Definition: real_type.F90:90
subroutine, public quit(message)
Stop program execution, display an error messahe.
type(mumps) function mumps_create(g, type, verb)
Constructor for the type mumps
Definition: mumps_mod.F90:153
integer, parameter, public sp
simple precision for real numbers
Definition: real_type.F90:60
solve LU x = b
Definition: mumps_mod.F90:88
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
subroutine, public mumps_init()
Initialisation.
Definition: mumps_mod.F90:97
integer, parameter mumps_ldlt_sdp
MUMPS LDLT-SDP, sym def pos mat.
integer choral_verb
Verbosity level.
integer, parameter, public dp
double precision for real numbers
Definition: real_type.F90:63
subroutine mumps_clear(lu)
Destructor.
Definition: mumps_mod.F90:112
DERIVED TYPE mumps: to use the library MUMPS
Definition: mumps_mod.F90:31
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
Factorisation of csr matrices with MUMPS.
Definition: mumps_mod.F90:58
DERIVED TYPE csr for sparse matrices
Definition: csr_mod.F90:13
integer, parameter mumps_lu
MUMPS LU.
destructor
Definition: mumps_mod.F90:79