The TELEMAC-MASCARET system  trunk
solve_mumps_par.F
Go to the documentation of this file.
1 ! **************************
2  SUBROUTINE solve_mumps_par
3 ! **************************
4 !
5  &(npoin,nsegb,gloseg,maxseg,da,xa,xinc,rhs,typext,knolg,
6  & npoin_tot)
7 !
8 !***********************************************************************
9 ! BIEF VERSION 7.1
10 !***********************************************************************
11 !
12 !brief PARALLEL DIRECT SYSTEM SOLUTION
13 !+
14 !
15 !history C. DENIS (SINETICS)
16 !+ 26/02/2010
17 !+ V6P0
18 !+ First version
19 !
20 !history J-M HERVOUET (LNHE)
21 !+ 23/11/2015
22 !+ V7P1
23 !+ Correction in the case of non symmetric matrix.
24 !+ Decommenting the deallocation of allocatable arrays.
25 !
26 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
27 !| NPOIN |-->| NUMBER OF POINTS IN THE PROCESSOR.
28 !| ... |-->| ...
29 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
30 !
31 ! PARALLEL SOLVING USING MUMPS
33 !
35  IMPLICIT NONE
36 ! structures MPI et MUMPS
37 #if defined HAVE_MUMPS
38  include 'dmumps_struc.h'
39 #endif
40 !
41 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
42 !
43  INTEGER, INTENT(IN) :: NPOIN,NSEGB,MAXSEG,NPOIN_TOT
44  INTEGER, INTENT(IN) :: GLOSEG(maxseg,2)
45  DOUBLE PRECISION, INTENT(INOUT) :: XA(*),RHS(npoin)
46  DOUBLE PRECISION, INTENT(INOUT) :: XINC(npoin),DA(npoin)
47  CHARACTER(LEN=1), INTENT(IN) :: TYPEXT
48  INTEGER, INTENT(IN) :: KNOLG(*)
49 !
50 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
51 !
52 #if defined HAVE_MUMPS
53 !
54  DOUBLE PRECISION TIME_IN_SECONDS2
55  EXTERNAL time_in_seconds2
56  type(dmumps_struc) mumps_par
57 ! INDIRECT ARRAYS
58  DOUBLE PRECISION ,ALLOCATABLE :: TEMP1(:),TEMP2(:)
59  INTEGER, ALLOCATABLE :: TEMP3(:)
60  INTEGER I,J,K,ERR,NBELEM,IER, IAUX, KTRI, NPOIN2, N2
61 ! variables chantier logiciel ob
62  DOUBLE PRECISION :: RMIN, RMAX, RAUX, EPSBLR
63  INTEGER :: INIV, INPREC, NAUX
64  CHARACTER(LEN=80) :: KVERS, KSTRATEGY, KNOM, KRENUM, KPRE, KMEM,
65  & kpost, kacce
66  LOGICAL :: LAUTOCORREC, LTRI, LPIVOTSTAT, LDEFAULT
67 
68 ! constantes
69  rmin=1.d-100
70  rmax=1.d+100
71 ! temporary parameters
72  iniv=1 ! verbose mode=2, standard mode=1
73  inprec=8 ! 8, 9, 10... or <0
74  epsblr=1.d-9 ! between d-6 and d-15
75 !
76  knom='<solve_mumps_par>'
77  knom=trim(adjustl(knom))
78  kstrategy='user' ! user, vitesse, memoire, robustesse, precision
79 !
80  kacce='FR' ! 'FR', 'FR+', 'LR+'
81  kpre='AUTO' ! 'AUTO' ou 'SANS'
82  krenum='AUTO' ! 'AUTO'/'AMD'/'AMF'/'SCOTCH'/'PTSCOTCH'/'PORD'&
83  & ! 'METIS'/'PARMETIS'/'QAMD'
84  kmem = 'IC' ! 'IC'/'OOC'
85  kpost = 'AUTO' ! 'AUTO' ou 'MINI' ou 'SANS' ou 'FORCE'
86 !
87  lautocorrec=.true. ! .false. less user-friendly
88  lpivotstat =.false. ! .true. less robust but quicker
89  ltri=.true. ! .true. filtering of tiny or huge matrix terms
90 !
91 ! utils/bief PARAMETRIZATION
92  ldefault=.true.
93  IF (ldefault) THEN
94  inprec=8
95  lpivotstat=.false.
96  ltri=.false.
97  kpre='AUTO'
98  krenum='AUTO'
99  kmem='IC'
100  kpost='FORCE'
101  kacce='FR'
102  ENDIF
103 !
104 !MUMPS INITIALISATION
105  mumps_par%COMM = comm
106  mumps_par%JOB = -1
107  IF(typext.EQ.'S') THEN
108  mumps_par%SYM = 2
109  ELSE
110  mumps_par%SYM = 0
111  ENDIF
112  mumps_par%PAR = 1
113  CALL dmumps(mumps_par)
114 
115 !CHECK MUMPS VERSION NUMBER
116  kvers=mumps_par%VERSION_NUMBER
117  kvers=trim(adjustl(kvers))
118  SELECT CASE (kvers)
119  CASE('4.10.0','5.2.1','5.2.1consortium')
120  CASE DEFAULT
121  WRITE(lu,*)'*************************************************'
122  WRITE(lu,*)knom//' * Warning *, MUMPS version not validated'
123  WRITE(lu,*)'version=',kvers
124  WRITE(lu,*)'*************************************************'
125  END SELECT
126 
127 !MUMPS PARAMETRIZATION
128 !-- VERBOSE MODE
129  IF (iniv.EQ.1) THEN
130  mumps_par%ICNTL(1)=lu
131  mumps_par%ICNTL(2)=0
132  mumps_par%ICNTL(3)=0
133  mumps_par%ICNTL(4)=1
134  ELSE
135  mumps_par%ICNTL(1)=lu
136  mumps_par%ICNTL(2)=lu
137  mumps_par%ICNTL(3)=lu
138  mumps_par%ICNTL(4)=2
139  ENDIF
140 !-- FRAMEWORK OF MUMPS COMPUTATION
141  mumps_par%ICNTL(5)=0 ! Input matrix format
142  mumps_par%ICNTL(9)=1 ! AX=b
143  mumps_par%ICNTL(18)=3 ! Distributed MPI mode
144 !-- PRETREATEMENTS
145  IF (kpre(1:4).EQ.'AUTO') THEN
146  mumps_par%ICNTL(6)=7
147  mumps_par%ICNTL(8)=77
148  mumps_par%ICNTL(12)=0
149  ELSE IF (kpre(1:4).EQ.'SANS') THEN
150  mumps_par%ICNTL(6)=0
151  mumps_par%ICNTL(8)=0
152  mumps_par%ICNTL(12)=1
153  ELSE
154  WRITE(lu,*)'***************************************************'
155  WRITE(lu,*)knom//' *** Syntax Error ****, Mauvaise option kpre'
156  WRITE(lu,*)'***************************************************'
157  CALL plante(1)
158  stop
159  ENDIF
160 ! -- BLR COMPRESSIONS
161  SELECT CASE(kacce)
162  CASE('FR')
163  mumps_par%ICNTL(35)=0
164  CASE('FR+')
165  mumps_par%KEEP(370)=1
166  mumps_par%KEEP(371)=1
167  mumps_par%ICNTL(35)=0
168  CASE('LR')
169  mumps_par%ICNTL(35)=2
170  mumps_par%ICNTL(36)=0
171  mumps_par%ICNTL(37)=0
172  mumps_par%CNTL(7)=epsblr
173  CASE('LR+')
174  mumps_par%KEEP(370)=1
175  mumps_par%KEEP(371)=1
176  mumps_par%ICNTL(35)=2
177  mumps_par%ICNTL(36)=0
178  mumps_par%ICNTL(37)=0
179  mumps_par%CNTL(7)=epsblr
180  CASE DEFAULT
181  WRITE(lu,*)'***************************************************'
182  WRITE(lu,*)knom//' ** Syntax Error ***, Mauvaise option kacce'
183  WRITE(lu,*)'***************************************************'
184  CALL plante(1)
185  stop
186  END SELECT
187 ! -- ORDERING
188  SELECT CASE(krenum)
189  CASE('AUTO')
190  mumps_par%ICNTL(7)=7
191  mumps_par%ICNTL(28)=1
192  mumps_par%ICNTL(29)=0
193  CASE('AMD')
194  mumps_par%ICNTL(7)=0
195  mumps_par%ICNTL(28)=1
196  mumps_par%ICNTL(29)=0
197  CASE('AMF')
198  mumps_par%ICNTL(7)=2
199  mumps_par%ICNTL(28)=1
200  mumps_par%ICNTL(29)=0
201  CASE('SCOTCH')
202  mumps_par%ICNTL(7)=3
203  mumps_par%ICNTL(28)=1
204  mumps_par%ICNTL(29)=0
205  CASE('PTSCOTCH')
206  mumps_par%ICNTL(7)=3
207  mumps_par%ICNTL(28)=2
208  mumps_par%ICNTL(29)=1
209  CASE('PORD')
210  mumps_par%ICNTL(7)=4
211  mumps_par%ICNTL(28)=1
212  mumps_par%ICNTL(29)=0
213  CASE('METIS')
214  mumps_par%ICNTL(7)=5
215  mumps_par%ICNTL(28)=1
216  mumps_par%ICNTL(29)=0
217  CASE('PARMETIS')
218  mumps_par%ICNTL(7)=5
219  mumps_par%ICNTL(28)=2
220  mumps_par%ICNTL(29)=2
221  CASE('QAMD')
222  mumps_par%ICNTL(7)=6
223  mumps_par%ICNTL(28)=1
224  mumps_par%ICNTL(29)=0
225  CASE DEFAULT
226  WRITE(lu,*)'***************************************************'
227  WRITE(lu,*)knom//' ** Syntax Error ***, Mauvaise option krenum'
228  WRITE(lu,*)'***************************************************'
229  CALL plante(1)
230  stop
231  END SELECT
232 ! -- NULL PIVOT DETECTION
233  mumps_par%ICNTL(25)=0
234  IF (inprec .GE. 0) THEN
235  mumps_par%ICNTL(13)=1
236  mumps_par%ICNTL(24)=1
237  mumps_par%CNTL(3)=-10.d0**(-inprec)
238  mumps_par%CNTL(5)=1.d+6
239  ELSE
240  mumps_par%ICNTL(13)=0
241  mumps_par%ICNTL(24)=0
242  mumps_par%CNTL(3)=0.d0
243  mumps_par%CNTL(5)=0.d0
244  ENDIF
245 ! -- PIVOT MANAGEMENT
246  IF (lpivotstat) THEN
247  mumps_par%CNTL(1)=-1.d0
248  ENDIF
249 ! -- MEMORY MANAGEMENT
250  mumps_par%ICNTL(14)=50
251  IF (kmem(1:2).EQ.'IC') THEN
252  mumps_par%ICNTL(22)=0
253  mumps_par%ICNTL(23)=0
254  ELSE IF (kmem(1:3).EQ.'OOC') THEN
255  mumps_par%ICNTL(22)=1
256  mumps_par%ICNTL(23)=0
257  mumps_par%OOC_TMPDIR='.'
258  ELSE
259  WRITE(lu,*)'***************************************************'
260  WRITE(lu,*)knom//' *** Syntax Error ***, Mauvaise option kmem'
261  WRITE(lu,*)'***************************************************'
262  CALL plante(1)
263  stop
264  ENDIF
265 ! -- POSTTREATMENT
266  SELECT CASE(kpost)
267  CASE('AUTO')
268  mumps_par%CNTL(2)=1.d-14
269  mumps_par%ICNTL(10)=4
270  mumps_par%ICNTL(11)=2
271  CASE('SANS')
272  mumps_par%CNTL(2)=0.d0
273  mumps_par%ICNTL(10)=0
274  mumps_par%ICNTL(11)=0
275  CASE('MINI')
276  mumps_par%CNTL(2)=0.d0
277  mumps_par%ICNTL(10)=-2
278  mumps_par%ICNTL(11)=0
279  CASE('FORCE')
280  mumps_par%CNTL(2)=1.d-50
281  mumps_par%ICNTL(10)=10
282  mumps_par%ICNTL(11)=1
283  CASE DEFAULT
284  WRITE(lu,*)'***************************************************'
285  WRITE(lu,*)knom//' *** Syntax Error ****, Mauvaise option kpost'
286  WRITE(lu,*)'***************************************************'
287  CALL plante(1)
288  stop
289  END SELECT
290 !
291 ! SIZE OF THE MATRIX
292  mumps_par%N = 2*npoin_tot
293  IF(typext.EQ.'S') THEN
294  mumps_par%NZ_LOC = npoin+nsegb
295  ELSE
296  mumps_par%NZ_LOC = npoin+2*nsegb
297  ENDIF
298 ! IF (.NOT. ASSOCIATED(MUMPS_PAR%IRN_LOC)) THEN
299  ALLOCATE(temp1(mumps_par%N),stat=err)
300  CALL check_allocate(err, "TEMP1")
301  ALLOCATE(temp2(mumps_par%N),stat=err)
302  CALL check_allocate(err, "TEMP2")
303  ALLOCATE(temp3(npoin),stat=err)
304  CALL check_allocate(err, "TEMP3")
305  ALLOCATE(mumps_par%IRN_LOC(mumps_par%NZ_LOC),stat=err)
306  CALL check_allocate(err, "MUMPS_PAR%IRN_LOC")
307  ALLOCATE(mumps_par%JCN_LOC(mumps_par%NZ_LOC),stat=err)
308  CALL check_allocate(err, "MUMPS_PAR%JCN_LOC")
309  ALLOCATE(mumps_par%A_LOC(mumps_par%NZ_LOC),stat=err)
310  CALL check_allocate(err, "MUMPS_PAR%A_LOC")
311  ALLOCATE(mumps_par%RHS(mumps_par%N),stat=err)
312  CALL check_allocate(err, "MUMPS_PAR%RHS")
313 ! END IF
314  temp1(:)=0.d0
315  temp2(:)=0.d0
316  temp3(:)=0
317  mumps_par%IRN_LOC(:)=0
318  mumps_par%JCN_LOC(:)=0
319  mumps_par%A_LOC(:)=0.0
320  mumps_par%RHS(:)=0.0
321 !
322  npoin2=npoin/2
323  n2=mumps_par%N/2
324  IF (ltri) THEN
325  ktri=1
326  DO k = 1,npoin
327  raux = abs(da(k))
328  IF (k.LE.npoin2) THEN
329  iaux = knolg(k)
330  ELSE
331  iaux = knolg(k-npoin2) + n2
332  ENDIF
333  IF ((raux.GT.rmin).AND.(raux.LT.rmax).AND.(iaux.GT.0)) THEN
334  mumps_par%IRN_LOC(ktri) = iaux
335  mumps_par%JCN_LOC(ktri) = iaux
336  mumps_par%A_LOC(ktri) = da(k)
337  ktri=ktri+1
338  END IF
339  temp1(iaux)=rhs(k)
340  temp3(k)=iaux
341  ENDDO
342  ELSE
343  DO k = 1,npoin
344  IF (k .LE. npoin/2) THEN
345  mumps_par%IRN_LOC(k) = knolg(k)
346  mumps_par%JCN_LOC(k) = knolg(k)
347  temp1(knolg(k))=rhs(k)
348  temp3(k)=knolg(k)
349  ELSE
350  mumps_par%IRN_LOC(k) = knolg(k-npoin/2) + mumps_par%N/2
351  mumps_par%JCN_LOC(k) = knolg(k-npoin/2) + mumps_par%N/2
352  temp1(knolg(k-npoin/2)+ mumps_par%N/2)=rhs(k)
353  temp3(k)=knolg(k-npoin/2)+ mumps_par%N/2
354  END IF
355  mumps_par%A_LOC(k) = da(k)
356  ENDDO
357  ENDIF
358 #if defined COMPAD
359  WRITE(lu,*) '(AD) COMPAD :: SOLVE_MUMPS_PAR.F: DIRECT CALL OF ',
360  & 'MPI_ALLREDUCE NOT AD-READY.'
361  WRITE(lu,*) ' PLEASE CONTACT JR @ ADJOINTWARE'
362  CALL plante(1)
363  stop
364 #endif
365 ! GLOBAL REDUCTION OF RHS VECTOR
366  CALL mpi_allreduce(temp1,temp2,mumps_par%N,mpi_double_precision,
367  & mpi_sum,
368  & comm,ier)
369  DO i=1,mumps_par%N
370  mumps_par%RHS(i)=temp2(i)
371  END DO
372 !
373  IF (ltri) THEN
374  nbelem=ktri
375  raux=100.d0*ktri/(npoin*1.0)
376  WRITE(lu,*)'*************************************************'
377  WRITE(lu,*)knom//'MATRIX n: ',mumps_par%N
378  WRITE(lu,*)knom//'MATRIX filtering 1: ',raux,' %'
379  WRITE(lu,*)'*************************************************'
380  ELSE
381  nbelem = npoin
382  ENDIF
383  IF(typext.EQ.'S') THEN
384  IF (ltri) THEN
385  naux=0
386  DO k = 1,nsegb
387  raux = abs(xa(k))
388  i = temp3(gloseg(k,1))
389  j = temp3(gloseg(k,2))
390  IF ((raux.GT.rmin).AND.(raux.LT.rmax).AND.(i.GT.0).
391  & and.(j.GT.0)) THEN
392  nbelem = nbelem + 1
393  IF(i.LT.j) THEN
394  mumps_par%IRN_LOC(nbelem) = i
395  mumps_par%JCN_LOC(nbelem) = j
396  ELSE
397  mumps_par%IRN_LOC(nbelem) = j
398  mumps_par%JCN_LOC(nbelem) = i
399  ENDIF
400  mumps_par%A_LOC(nbelem) = xa(k)
401  ELSE
402  naux=naux+1
403  ENDIF
404  ENDDO
405  raux=100.d0*naux/(1.d0*nsegb)
406  WRITE(lu,*)'*************************************************'
407  WRITE(lu,*)knom//'MATRIX nz_loc: ',mumps_par%NZ_LOC
408  WRITE(lu,*)knom//'MATRIX filtering 2: ',raux,' %'
409  WRITE(lu,*)'*************************************************'
410  ELSE
411  DO k = 1,nsegb
412  i = temp3(gloseg(k,1))
413  j = temp3(gloseg(k,2))
414  nbelem = nbelem + 1
415  IF(i.LT.j) THEN
416  mumps_par%IRN_LOC(nbelem) = i
417  mumps_par%JCN_LOC(nbelem) = j
418  ELSE
419  mumps_par%IRN_LOC(nbelem) = j
420  mumps_par%JCN_LOC(nbelem) = i
421  ENDIF
422  mumps_par%A_LOC(nbelem) = xa(k)
423  ENDDO
424  ENDIF
425  ELSE
426  IF (ltri) THEN
427  naux=0
428  DO k = 1,nsegb
429  raux = abs(xa(k))
430  i = temp3(gloseg(k,1))
431  j = temp3(gloseg(k,2))
432  IF ((raux.GT.rmin).AND.(raux.LT.rmax).AND.(i.GT.0).
433  & and.(j.GT.0)) THEN
434  nbelem = nbelem + 1
435  mumps_par%A_LOC(nbelem) = xa(k)
436  ELSE
437  naux=naux+1
438  ENDIF
439  ENDDO
440  DO k = 1,nsegb
441  raux = abs(xa(k+nsegb))
442  i = temp3(gloseg(k,2))
443  j = temp3(gloseg(k,1))
444  IF ((raux.GT.rmin).AND.(raux.LT.rmax).AND.(i.GT.0).
445  & and.(j.GT.0)) THEN
446  nbelem = nbelem + 1
447  mumps_par%IRN_LOC(nbelem) = i
448  mumps_par%JCN_LOC(nbelem) = j
449  mumps_par%A_LOC(nbelem) = xa(k+nsegb)
450  ELSE
451  naux=naux+1
452  ENDIF
453  ENDDO
454  raux=100.d0*naux/(2.d0*nsegb)
455  WRITE(lu,*)'*************************************************'
456  WRITE(lu,*)knom//'MATRIX nz_loc: ',mumps_par%NZ_LOC
457  WRITE(lu,*)knom//'MATRIX filtering 2: ',raux,' %'
458  WRITE(lu,*)'*************************************************'
459  ELSE
460  DO k = 1,nsegb
461  i = temp3(gloseg(k,1))
462  j = temp3(gloseg(k,2))
463  nbelem = nbelem + 1
464  mumps_par%IRN_LOC(nbelem) = i
465  mumps_par%JCN_LOC(nbelem) = j
466  mumps_par%A_LOC(nbelem) = xa(k)
467  ENDDO
468  DO k = 1,nsegb
469  i = temp3(gloseg(k,2))
470  j = temp3(gloseg(k,1))
471  nbelem = nbelem + 1
472  mumps_par%IRN_LOC(nbelem) = i
473  mumps_par%JCN_LOC(nbelem) = j
474  mumps_par%A_LOC(nbelem) = xa(k+nsegb)
475  ENDDO
476  ENDIF
477  ENDIF
478 !
479 ! -----------
480 ! RESOLUTION
481 ! -----------
482 !
483  mumps_par%JOB = 1
484  CALL dmumps(mumps_par)
485 !
486  mumps_par%JOB = 2
487  CALL dmumps(mumps_par)
488 !
489  mumps_par%JOB = 3
490  CALL dmumps(mumps_par)
491 !
492  temp1(:)=0.0
493  IF(mumps_par%MYID.EQ. 0 ) THEN
494 ! An error occured...
495  IF(mumps_par%INFO(1).LT.0) THEN
496  WRITE(lu,2001) mumps_par%INFO(1)
497  CALL plante(1)
498  stop
499  ENDIF
500  DO k = 1,mumps_par%N
501  temp1(k)=mumps_par%RHS(k)
502  END DO
503  END IF
504  temp2(:)=0.0
505 !
506 !! COMPAD-DCO-MPICHECK BEGIN JR2016
507 #if defined COMPAD
508  WRITE(lu,*) '(AD) COMPAD :: SOLVE_MUMPS_PAR.F: DIRECT CALL OF ',
509  & 'MPI_ALLREDUCE NOT AD-READY.'
510  WRITE(lu,*) ' PLEASE CONTACT JR @ ADJOINTWARE'
511  CALL plante(1)
512  stop
513 #endif
514 !! COMPAD-DCO-MPICHECK END JR2016
515  CALL mpi_bcast(temp1,mumps_par%N,mpi_double_precision,
516  & 0,comm,ier)
517  CALL mpi_barrier(comm,ier)
518  DO k = 1,npoin
519  i=temp3(k)
520  xinc(k)=temp1(i)
521  ENDDO
522 
523 !
524 ! DEALLOCATION
525 !
526  DEALLOCATE(temp1)
527  DEALLOCATE(temp2)
528  DEALLOCATE(temp3)
529  DEALLOCATE(mumps_par%IRN_LOC)
530  DEALLOCATE(mumps_par%JCN_LOC)
531  DEALLOCATE(mumps_par%A_LOC)
532  DEALLOCATE(mumps_par%RHS)
533  CALL mpi_barrier(comm,ier)
534 !
535 ! Killing the current instance
536 !
537  mumps_par%JOB = -2
538  CALL dmumps(mumps_par)
539  RETURN
540 2001 FORMAT(1x,'SOLVE_MUMPS: ERROR DURING SOLVE: '
541  & ,/,1x,'ERROR CODE INFO(1): ',1i6)
542 #else
543  WRITE(lu,2019)
544 2019 FORMAT(1x,'MUMPS_PAR NOT INSTALLED IN THIS SYSTEM',/,1x,
545  & 'CHOOSE OTHER METHOD ',///)
546  CALL plante(1)
547  stop
548 !
549 #endif
550 !
551  END
integer, parameter mpi_double_precision
subroutine solve_mumps_par(NPOIN, NSEGB, GLOSEG, MAXSEG, DA, XA, XINC, RHS, TYPEXT, KNOLG, NPOIN_TOT)