Télécharger la plaquette Cast3M


fr en

Télécharger miltp3.eso

Retour à la liste

Numérotation des lignes : 

miltp3

  1. C MILTP3    SOURCE    GOUNAND   05/02/16    21:16:19     5029
  2.       SUBROUTINE MILTP3(NTT,NNZA,A,JA,IA,
  3.      $     XLFIL,XDTOL,XSPIV,NNZMLU,
  4.      $     ALU,JLU,JU,NNZLU,
  5.      $     IWORK,JWORK,RWORK,LWORK,
  6.      $     KWORK,QWORK,
  7.      $     IPERM,JPERM,
  8.      $     IMPR,IRET)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10.       IMPLICIT INTEGER (I-N)
  11. C***********************************************************************
  12. C NOM         : MILTP3
  13. C DESCRIPTION :
  14. C     Calcul du préconditionneur ILUTP d'une matrice Morse
  15. C     ILUT : Incomplete LU factorization with dual truncation strategy
  16. C     La matrice Morse est au format CSR (Column Sparse Row)
  17. C     Le préconditionneur est au format MSR (Modified Sparse Row,
  18. C     stockage de l'inverse de la diagonale) (ALU, JLU).
  19. C     JU est un tableau supplémentaire pour le repérage des diagonales.
  20. C     (voir la doc. de Sparskit version 2)
  21. C
  22. C LANGAGE     : FORTRAN 77 (sauf E/S)
  23. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  24. C               mél : gounand@semt2.smts.cea.fr
  25. C REFERENCE (bibtex-like) :
  26. C  @BOOK{templates,
  27. C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  28. C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  29. C                H. Van der Vorst},
  30. C        TITLE={Templates for the Solution of Linear Systems :
  31. C               Building Blocks for Iterative Methods},
  32. C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  33. C  -> URL : http://www.netlib.org/templates/Templates.html
  34. C  Sparskit : a basic tool kit for sparse matrix computations
  35. C  Version 2 (Youcef Saad)
  36. C  -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  37. C REMARQUES :
  38. C Quelques parties pourraient être améliorés :
  39. C - l'algorithme de quicksplit (cf.qsplit.eso)
  40. C - la partie où l'on doit effectuer l'élimination dans le bon ordre
  41. C   pourrait utiliser un algorithme de tri type quicksort
  42. C
  43. C
  44. C On garde absolument tous les termes qui seraient dans ILU(0).
  45. C
  46. C
  47. C***********************************************************************
  48. C APPELES          : QSPLIT
  49. C APPELE PAR       : PRILUT
  50. C***********************************************************************
  51. C ENTREES            : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
  52. C ENTREES/SORTIES    : IWORK,JWORK,RWORK (tableaux de travail)
  53. C SORTIES            : ALU,JLU,JU,NNZLU
  54. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  55. C***********************************************************************
  56. C VERSION    : v1, 26/04/04, version initiale
  57. C HISTORIQUE : v1, 26/04/04, création
  58. C HISTORIQUE :
  59. C HISTORIQUE :
  60. C***********************************************************************
  61. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  62. C en cas de modification de ce sous-programme afin de faciliter
  63. C la maintenance !
  64. C***********************************************************************
  65.  
  66. -INC PPARAM
  67. -INC CCOPTIO
  68. -INC CCREEL
  69. *     .. Scalar Arguments ..
  70.       INTEGER  NTT,NNZA,NNZMLU,NNZLU
  71.       INTEGER  IMPR,IRET
  72.       REAL*8 XDTOL,XLFIL
  73. *     ..
  74. *     .. Array Arguments ..
  75.       REAL*8  A(NNZA),ALU(NNZMLU+1)
  76.       INTEGER JA(NNZA),IA(NTT+1),JU(NTT),JLU(NNZMLU+1)
  77.       INTEGER IWORK(NTT)
  78.       INTEGER JWORK(NTT)
  79.       REAL*8 RWORK(NTT)
  80.       INTEGER KWORK(NTT)
  81.       REAL*8 QWORK(NTT)
  82.       LOGICAL LWORK(NTT)
  83.       INTEGER IPERM(NTT)
  84.       INTEGER JPERM(NTT)
  85. *
  86. *     .. Variables locales
  87. *     .. Parameters
  88.       REAL*8     ZERO      ,ONE
  89.       PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  90. *     ..
  91. C     Nombre de pivots petits
  92.       INTEGER NBPIVP
  93. C     Nombre de pivots inférieurs à 0
  94.       INTEGER NBPIVI
  95. *
  96.       INTEGER ILFIL,ILFILL,ILFILU
  97.       INTEGER COL,COLLU
  98.       INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
  99.       INTEGER JJL,JPOS,JROW,KJL,KLU
  100.       INTEGER LEN,LENL,LENU
  101.       INTEGER JU0,ITMP
  102.       REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
  103.       LOGICAL LTMP,LFACT,LPIV,LFOUND
  104.       REAL*8 XNTLIG
  105. *      REAL*8 XPOND1,XPOND2
  106. *
  107. * Executable statements
  108. *
  109.       IF (XLFIL.LT.1.D0) THEN
  110.          WRITE(IOIMP,*) 'XLFIL.LT.1.D0 non acceptable'
  111.          GOTO 9999
  112.       ENDIF
  113.       PETIT=XZPREC**(0.75D0)
  114.       NBPIVP=0
  115.       NBPIVI=0
  116. c-----------------------------------------------------------------------
  117. c     initialize ju0 (points to next element to be added to alu,jlu)
  118. c     and pointer array.
  119. c-----------------------------------------------------------------------
  120.       JU0=NTT+2
  121.       JLU(1)=JU0
  122. *
  123. * Initialize permutation arrays and lwork
  124. *
  125.       DO ITT=1,NTT
  126.          JWORK(ITT)=0
  127.          IPERM(ITT)=ITT
  128.          JPERM(ITT)=ITT
  129. *!!!!!!!!!!!!!!
  130.          LWORK(ITT)=.FALSE.
  131.       ENDDO
  132. *
  133. * On boucle sur les ddl
  134. *
  135.       DO 1 ITT=1,NTT
  136. * Calcul de la norme de la ittème ligne
  137.          ISTRT=IA(ITT)
  138.          ISTOP=IA(ITT+1)-1
  139.          TNORM=ZERO
  140.          DO 11 IJA=ISTRT,ISTOP
  141.             TNORM=TNORM+ABS(A(IJA))
  142.  11      CONTINUE
  143.          TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  144.          IF (TNORM.LE.XPETIT) THEN
  145.             WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  146.             GOTO 9999
  147.          ENDIF
  148. * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
  149. * parties triangulaires inférieures et supérieurs.
  150. *
  151. * On peut le programmer autrement compte tenu du fait que les
  152. * colonnes sont ordonnées
  153. *
  154.          LENU=1
  155.          LENL=0
  156.          IWORK(ITT)=ITT
  157.          RWORK(ITT)=ZERO
  158.          JWORK(ITT)=ITT
  159.          LWORK(ITT)=.TRUE.
  160.  
  161.          DO 12 IJA=ISTRT,ISTOP
  162. *!            COL=JA(IJA)
  163.             COL=JPERM(JA(IJA))
  164.             VAL=A(IJA)
  165.             IF (COL.LT.ITT) THEN
  166.                LENL=LENL+1
  167.                IWORK(LENL)=COL
  168.                RWORK(LENL)=VAL
  169. *!!!!!!!!!!!!!!
  170.                LWORK(LENL)=.TRUE.
  171.                JWORK(COL) =LENL
  172.             ELSEIF (COL.EQ.ITT) THEN
  173.                RWORK(ITT)=VAL
  174.             ELSE
  175.                LENU=LENU+1
  176.                JPOS=ITT+LENU-1
  177.                IWORK(JPOS)=COL
  178.                RWORK(JPOS)=VAL
  179. *!!!!!!!!!!!!!!
  180.                LWORK(JPOS)=.TRUE.
  181.                JWORK(COL) =JPOS
  182.             ENDIF
  183.  12      CONTINUE
  184. * On a enlevé les stratégies bizarres
  185.          XNTLIG=DBLE(NNZA)/DBLE(NTT)
  186. *!!         ILFIL=INT(XNTLIG*XLFIL/2.D0)
  187.          ILFIL=INT(XNTLIG*(XLFIL-1.D0)/2.D0)
  188.          ILFILL=ILFIL
  189.          ILFILU=ILFIL
  190. * Elimination avec les lignes précédentes
  191.          JJL=0
  192.          LEN=0
  193.  13      CONTINUE
  194.          JJL=JJL+1
  195.          IF (JJL.LE.LENL) THEN
  196. c-----------------------------------------------------------------------
  197. c     in order to do the elimination in the correct order we must select
  198. c     the smallest column index among jw(k), k=jj+1, ..., lenl.
  199. c-----------------------------------------------------------------------
  200.             JROW=IWORK(JJL)
  201.             KJL=JJL
  202. c     determine smallest column index
  203.             DO 132 ILENL=JJL+1,LENL
  204.                IF (IWORK(ILENL).LT.JROW) THEN
  205.                   JROW=IWORK(ILENL)
  206.                   KJL=ILENL
  207.                ENDIF
  208.  132        CONTINUE
  209.             IF (KJL.NE.JJL) THEN
  210.                ITMP=IWORK(JJL)
  211.                IWORK(JJL)=IWORK(KJL)
  212.                IWORK(KJL)=ITMP
  213.                JWORK(JROW)=JJL
  214.                JWORK(ITMP)=KJL
  215.                RTMP=RWORK(JJL)
  216.                RWORK(JJL)=RWORK(KJL)
  217.                RWORK(KJL)=RTMP
  218. *!!!!!!!!!!!!!!
  219.                LTMP=LWORK(JJL)
  220.                LWORK(JJL)=LWORK(KJL)
  221.                LWORK(KJL)=LTMP
  222.             ENDIF
  223. c     zero out element in row by setting jw(n+jrow) to zero.
  224.             JWORK(JROW)=0
  225. c     get the multiplier for row to be eliminated (jrow).
  226.             FACT=RWORK(JJL)*ALU(JROW)
  227.             LFACT=LWORK(JJL)
  228. *!!!!!!!!!!!!!!
  229.             IF (ABS(FACT).LE.XDTOL.AND.(.NOT.LFACT)) THEN
  230. *            IF (ABS(FACT).LE.XDTOL) THEN
  231.                GOTO 13
  232.             ENDIF
  233. c
  234. c     combine current row and row jrow
  235. c
  236.             DO 134 KLU=JU(JROW),JLU(JROW+1)-1
  237.                TERMLU=FACT*ALU(KLU)
  238. *!               COLLU =JLU(KLU)
  239.                COLLU =JPERM(JLU(KLU))
  240.                JPOS  =JWORK(COLLU)
  241.                IF (COLLU.GE.ITT) THEN
  242. c     dealing with upper part (including diagonal)
  243.                   IF (JPOS.EQ.0) THEN
  244. c     this is a fill-in element
  245.                      LENU=LENU+1
  246.                      IF (LENU.GT.NTT-ITT+1) THEN
  247.                         WRITE(IOIMP,*) 'pb. matrice U'
  248.                         GOTO 9999
  249.                      ENDIF
  250.                      IWORK(ITT+LENU-1)=COLLU
  251.                      JWORK(COLLU)=ITT+LENU-1
  252.                      RWORK(ITT+LENU-1)=-TERMLU
  253.                      LWORK(ITT+LENU-1)=.FALSE.
  254.                   ELSE
  255. c     this is not a fill-in element
  256.                      RWORK(JPOS)=RWORK(JPOS)-TERMLU
  257.                   ENDIF
  258.                ELSE
  259. c     dealing with lower part.
  260.                   IF (JPOS.EQ.0) THEN
  261. c     this is a fill-in element
  262.                      LENL=LENL+1
  263.                      IF (LENL.GT.ITT-1) THEN
  264.                         WRITE(IOIMP,*) 'pb. matrice L'
  265.                         WRITE(IOIMP,*) 'LENL=',LENL
  266.                         WRITE(IOIMP,*) 'ITT=',ITT
  267.                         GOTO 9999
  268.                      ENDIF
  269.                      IWORK(LENL) =COLLU
  270.                      JWORK(COLLU)=LENL
  271.                      RWORK(LENL) =-TERMLU
  272.                      LWORK(LENL) =.FALSE.
  273.                   ELSE
  274. c     this is not a fill-in element
  275.                      RWORK(JPOS)=RWORK(JPOS)-TERMLU
  276.                   ENDIF
  277.                ENDIF
  278.  134        CONTINUE
  279. c     store this pivot element -- (from left to right -- no danger of
  280. c     overlap with the working elements in L (pivots).
  281.             LEN=LEN+1
  282. * ?????????
  283.             LWORK(LEN)=LFACT
  284.             RWORK(LEN)=FACT
  285.             IWORK(LEN)=JROW
  286.             GOTO 13
  287.          ENDIF
  288. c     reset double-pointer to zero (U-part including diagonal)
  289.          DO 14 ILENU=1,LENU
  290.             JWORK(IWORK(ITT+ILENU-1))=0
  291.  14      CONTINUE
  292. c     update L-matrix
  293.          LENL=LEN
  294. *
  295. * First store elements of the ilu(0) factorization
  296. *
  297.          LEN2=0
  298.          DO ILENL=1,LENL
  299.             IF (LWORK(ILENL)) THEN
  300.                LEN2=LEN2+1
  301.                KWORK(LEN2)=IWORK(ILENL)
  302.                QWORK(LEN2)=RWORK(ILENL)
  303.             ENDIF
  304.          ENDDO
  305.          LENIL0=LEN2
  306. * then the fill-in elements
  307.          DO ILENL=1,LENL
  308.             IF (.NOT.LWORK(ILENL)) THEN
  309.                LEN2=LEN2+1
  310.                KWORK(LEN2)=IWORK(ILENL)
  311.                QWORK(LEN2)=RWORK(ILENL)
  312.             ENDIF
  313.          ENDDO
  314.          IF (LEN2.NE.LENL) THEN
  315.             WRITE(IOIMP,*) 'Erreur de programmation 1'
  316.             GOTO 9999
  317.          ENDIF
  318.          LENFIL=LENL-LENIL0
  319. * sort the fill-in elements by quick-split
  320.          LEN=MIN(LENFIL,ILFILL)
  321. c     sort by quick-split
  322.          CALL QSPLIT(QWORK(LENIL0+1),KWORK(LENIL0+1),LENFIL,LEN,
  323.      $        IMPR,IRET)
  324. c     store L-part -- in original coordinates
  325.          DO 15 ILEN=1,LENIL0+LEN
  326.             ALU(JU0)=QWORK(ILEN)
  327.             JLU(JU0)=IPERM(KWORK(ILEN))
  328.             JU0=JU0+1
  329.  15      CONTINUE
  330. *
  331. c     save pointer to beginning of row ii of U
  332.          JU(ITT)=JU0
  333. c     update U-matrix -- first apply dropping strategy
  334.          LEN=0
  335.          DO 16 ILEN=1,LENU-1
  336.             IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
  337.                LEN=LEN+1
  338.                RWORK(ITT+LEN)=RWORK(ITT+ILEN)
  339.                IWORK(ITT+LEN)=IWORK(ITT+ILEN)
  340. *!!!!!!!!!!!!!!
  341.                LWORK(ITT+LEN)=LWORK(ITT+ILEN)
  342.             ENDIF
  343.  16      CONTINUE
  344.          LENU=LEN+1
  345. *
  346. * First store elements of the ilu(0) factorization
  347. *
  348.          LEN2=ITT-1
  349.          DO ILENU=ITT,ITT+LENU-1
  350.             IF (LWORK(ILENU)) THEN
  351.                LEN2=LEN2+1
  352.                KWORK(LEN2)=IWORK(ILENU)
  353.                QWORK(LEN2)=RWORK(ILENU)
  354.             ENDIF
  355.          ENDDO
  356.          LENIL0=LEN2-ITT+1
  357. * then the fill-in elements
  358.          DO ILENU=ITT,ITT+LENU-1
  359.             IF (.NOT.LWORK(ILENU)) THEN
  360.                LEN2=LEN2+1
  361.                KWORK(LEN2)=IWORK(ILENU)
  362.                QWORK(LEN2)=RWORK(ILENU)
  363.             ENDIF
  364.          ENDDO
  365.          IF (LEN2-ITT+1.NE.LENU) THEN
  366.             WRITE(IOIMP,*) 'Erreur de programmation 2'
  367.             GOTO 9999
  368.          ENDIF
  369.          LENFIL=LENU-LENIL0
  370. * sort the fill-in elements by quicksplit
  371. *!!!!!!!!
  372.          LEN=MIN(LENFIL,ILFILU)
  373. *!!         LEN=MIN(LENU,ILFILU+1)
  374. *!!!         LEN=MIN(LENU,ILFILU)
  375.          CALL QSPLIT(QWORK(ITT+LENIL0),KWORK(ITT+LENIL0),LENFIL,LEN,
  376.      $        IMPR,IRET)
  377. C            WRITE(IOIMP,*) 'NTT=',NTT
  378. C            WRITE(IOIMP,*) 'ITT=',ITT
  379. C            WRITE(IOIMP,*) 'ILFILU=',ILFILU
  380. C            WRITE(IOIMP,*) 'LEN=',LEN
  381. C            WRITE(IOIMP,*) 'LENU=',LENU
  382. C            DO JTT=ITT,NTT
  383. C               WRITE(IOIMP,*) 'JTT=',JTT
  384. C               WRITE(IOIMP,*) 'IWORK(JTT)=',IWORK(JTT)
  385. C               WRITE(IOIMP,*) 'RWORK(JTT)=',RWORK(JTT)
  386. C               WRITE(IOIMP,*) 'LWORK(JTT)=',LWORK(JTT)
  387. C            ENDDO
  388. C            DO JTT=ITT,NTT
  389. C               WRITE(IOIMP,*) 'JTT=',JTT
  390. C               WRITE(IOIMP,*) 'KWORK(JTT)=',KWORK(JTT)
  391. C               WRITE(IOIMP,*) 'QWORK(JTT)=',QWORK(JTT)
  392. C            ENDDO
  393. *
  394. *    Determine next pivot
  395. *
  396.          IMAX=ITT
  397.          XMAX=ABS(QWORK(IMAX))
  398.          IF (XSPIV.GT.ZERO) THEN
  399.             XMAX=XMAX/XSPIV
  400.             LPIV=.TRUE.
  401. * On pivote même si xspiv.EQ.zero si la diagonale est petite
  402.          ELSEIF (XSPIV.EQ.ZERO.AND.XMAX.LE.PETIT*TNORM) THEN
  403.             LPIV=.TRUE.
  404.          ELSE
  405.             LPIV=.FALSE.
  406.          ENDIF
  407. *
  408. *         WRITE(IOIMP,*) 'LPIV=',LPIV
  409.          IF (LPIV) THEN
  410.             DO ILEN=ITT+1,ITT+LENIL0+LEN-1
  411.                XCOU=ABS(QWORK(ILEN))
  412.                IF (XCOU.GT.XMAX) THEN
  413.                   IMAX=ILEN
  414.                   XMAX=XCOU
  415.                ENDIF
  416.             ENDDO
  417. * On ne cherche le max que dans le ilu0 sauf si on ne trouve rien
  418. C            LFOUND=.FALSE.
  419. C            DO ILEN=ITT+1,ITT+LENIL0-1
  420. C               XCOU=ABS(QWORK(ILEN))
  421. C               IF (XCOU.GT.XMAX) THEN
  422. C                  IMAX=ILEN
  423. C                  XMAX=XCOU
  424. C                  LFOUND=.TRUE.
  425. C               ENDIF
  426. C            ENDDO
  427. C            IF (.NOT.LFOUND) THEN
  428. C               DO ILEN=ITT+LENIL0,ITT+LENIL0+LEN-1
  429. C                  XCOU=ABS(QWORK(ILEN))
  430. C                  IF (XCOU.GT.XMAX) THEN
  431. C                     IMAX=ILEN
  432. C                     XMAX=XCOU
  433. C                  ENDIF
  434. C               ENDDO
  435. C            ENDIF
  436. *
  437. * Echange des valeurs
  438. *
  439.             XTMP=QWORK(ITT)
  440.             QWORK(ITT)=QWORK(IMAX)
  441.             QWORK(IMAX)=XTMP
  442. *
  443. * Update de la numérotation ancienne
  444. *
  445.             JMAX=KWORK(IMAX)
  446.             ITMP=IPERM(ITT)
  447.             IPERM(ITT)=IPERM(JMAX)
  448.             IPERM(JMAX)=ITMP
  449. *
  450. * Update de la numérotation nouvelle
  451. *
  452.             JPERM(IPERM(ITT))=ITT
  453.             JPERM(IPERM(JMAX))=JMAX
  454.          ENDIF
  455. *
  456. c     copy U-Part in original coordinates
  457. *
  458.          DO 17 ILEN=ITT+1,ITT+LENIL0+LEN-1
  459.             JLU(JU0)=IPERM(KWORK(ILEN))
  460.             ALU(JU0)=QWORK(ILEN)
  461.             JU0=JU0+1
  462.  17      CONTINUE
  463. c     store inverse of diagonal element of u
  464.          VALPIV=QWORK(ITT)
  465.          IF (VALPIV.EQ.ZERO) THEN
  466.             WRITE(IOIMP,*) 'Pivot',ITT,' nul : ',
  467.      $           'le préconditionnement ILUTPGO2 est impossible.'
  468.             IRET=-2
  469.             GOTO 9999
  470.          ELSE
  471.             IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  472.                NBPIVP=NBPIVP+1
  473.             ENDIF
  474.             IF (VALPIV.LT.ZERO) THEN
  475.                NBPIVI=NBPIVI+1
  476.             ENDIF
  477.          ENDIF
  478.          ALU(ITT)=ONE/VALPIV
  479. c     update pointer to beginning of next row of U.
  480.          JLU(ITT+1)=JU0
  481.  1    CONTINUE
  482.       NNZLU=JLU(NTT+1)-1
  483. C
  484. C Permutation des colonnes de LU
  485. C
  486.       DO KLU=JLU(1),JLU(NTT+1)-1
  487.          JLU(KLU)=JPERM(JLU(KLU))
  488.       ENDDO
  489. C
  490. C Permutation des colonnes de A
  491. C
  492.       DO IJA=IA(1),IA(NTT+1)-1
  493.          JA(IJA)=JPERM(JA(IJA))
  494.       ENDDO
  495. *
  496. * Warning(s)
  497. *
  498.       IF (IMPR.GT.0) THEN
  499.          IF (NBPIVP.GT.0) THEN
  500.             WRITE(IOIMP,*) 'WARNING !',
  501.      $           ' miltp3 :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  502.          ENDIF
  503. *!         IF (NBPIVI.GT.0) THEN
  504. *!            WRITE(IOIMP,*) 'WARNING !',
  505. *!     $           ' miltp3 :',NBPIVI,' pivot(s) < 0'
  506. *!         ENDIF
  507.       ENDIF
  508. *
  509. * Normal termination
  510. *
  511.       IRET=0
  512.       RETURN
  513. *
  514. * Error handling
  515. *
  516.  9999 CONTINUE
  517.       IRET=1
  518.       WRITE(IOIMP,*) 'An error was detected in miltp3.eso'
  519.       RETURN
  520. *
  521. *     End of MILTP3
  522. *
  523.       END
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales