Esope source files

   1 : C IPLNU3    SOURCE    CHAT      05/01/13    00:43:40     5004           
   2 :       SUBROUTINE IPLNU3(MA)
   3 : c
   4 : c fonction
   5 : c     Cette routine sert à resoudre un systeme lineaire non symmetrqiue
   6 : c     par la methode LU ( cf Numerical recipes in fortran 2nd edition)
   7 : c
   8 : c variables (E:entrée / S:sortie)
   9 : c     a (E/S) est un tableau nXn dont les npXnp 1er membre contienne la 
  10 : c             matrice. A la sortie a est sous la forme LU
  11 : c     indx (S) est un tableau de dimension n qui contiendra une tableau
  12 : c              d'indice de permuation
  13 : c
  14 : c     d (S) est la signature de la permutation stockee dans indx (-/+1)
  15 : c
  16 : c date 18/07/94
  17 : c
  18 : c langage fortran pur
  19 : c
  20 :       IMPLICIT INTEGER(I-N)
  21 :       IMPLICIT    REAL*8(A-H,O-Z)     
  22 :       PARAMETER (TINY=1.0D-20)
  23 :       SEGMENT MA 
  24 : *     ce segment contient la matrice du syteme lineaire
  25 : *     a contient les termes et indx contient les indice de permutation
  26 : *     une fois la decomposition effectué  
  27 :          REAL*8 A(NN+1,NN+1),D
  28 :          INTEGER INDX(NN+1)
  29 :       ENDSEGMENT
  30 :       SEGMENT MV
  31 :          REAL*8 VV(NN+1)
  32 :       ENDSEGMENT
  33 :       SEGACT,MA*MOD
  34 :       N = A(/1)
  35 :       NN = N-1
  36 :       SEGINI MV
  37 :       D=1.D0
  38 :       DO 12 I=1,N
  39 :         AAMAX=0.D0
  40 :         DO 11 J=1,N
  41 :           IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
  42 : 11      CONTINUE
  43 :         IF (AAMAX.EQ.0.) THEN
  44 : c          la matrice contient une colonne de termes nuls
  44 :            CALL ERREUR(-292)
  45 :            AAMAX = TINY
  46 :         ENDIF
  47 :         VV(I)=1.D0/AAMAX
  48 : 12    CONTINUE
  49 :       DO 19 J=1,N
  50 :         IF (J.GT.1) THEN
  51 :           DO 14 I=1,J-1
  52 :             SUM=A(I,J)
  53 :             IF (I.GT.1)THEN
  54 :               DO 13 K=1,I-1
  55 :                 SUM=SUM-A(I,K)*A(K,J)
  56 : 13            CONTINUE
  57 :               A(I,J)=SUM
  58 :             ENDIF
  59 : 14        CONTINUE
  60 :         ENDIF
  61 :         AAMAX=0.D0
  62 :         DO 16 I=J,N
  63 :           SUM=A(I,J)
  64 :           IF (J.GT.1)THEN
  65 :             DO 15 K=1,J-1
  66 :               SUM=SUM-A(I,K)*A(K,J)
  67 : 15          CONTINUE
  68 :             A(I,J)=SUM
  69 :           ENDIF
  70 :           DUM=VV(I)*ABS(SUM)
  71 :           IF (DUM.GE.AAMAX) THEN
  72 :             IMAX=I
  73 :             AAMAX=DUM
  74 :           ENDIF
  75 : 16      CONTINUE
  76 :         IF (J.NE.IMAX)THEN
  77 :           DO 17 K=1,N
  78 :             DUM=A(IMAX,K)
  79 :             A(IMAX,K)=A(J,K)
  80 :             A(J,K)=DUM
  81 : 17        CONTINUE
  82 :           D=-D
  83 :           VV(IMAX)=VV(J)
  84 :         ENDIF
  85 :         INDX(J)=IMAX
  86 :         IF(J.NE.N)THEN
  87 :           IF(A(J,J).EQ. 0.D0) THEN
  88 : *     &  'Nuage pathologique pour la méthode d interpolation'
  89 : *        matrice singulière dans l'algorithme d'interpolation
  90 :             A(J,J)=TINY
  90 :             CALL ERREUR(-292)
  91 :           ENDIF
  92 :           DUM=1.D0/A(J,J)
  93 :           DO 18 I=J+1,N
  94 :             A(I,J)=A(I,J)*DUM
  95 : 18        CONTINUE
  96 :         ENDIF
  97 : 19    CONTINUE
  98 :       IF(A(N,N).EQ. 0.D0) THEN
  99 : *     &  'Nuage pathologique pour la méthode d interpolation'
 100 : *        matrice singulière dans l'algorithme d'interpolation
 100 :            CALL ERREUR(-292)
 101 :            A(N,N)=TINY
 102 :       ENDIF
 103 :       SEGSUP MV
 104 :       RETURN
 105 :       END
 106 : 
 107 :  
 108 :