C IPLNU3 SOURCE CHAT 05/01/13 00:43:40 5004 SUBROUTINE IPLNU3(MA) c c fonction c Cette routine sert à resoudre un systeme lineaire non symmetrqiue c par la methode LU ( cf Numerical recipes in fortran 2nd edition) c c variables (E:entrée / S:sortie) c a (E/S) est un tableau nXn dont les npXnp 1er membre contienne la c matrice. A la sortie a est sous la forme LU c indx (S) est un tableau de dimension n qui contiendra une tableau c d'indice de permuation c c d (S) est la signature de la permutation stockee dans indx (-/+1) c c date 18/07/94 c c langage fortran pur c IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (TINY=1.0D-20) SEGMENT MA * ce segment contient la matrice du syteme lineaire * a contient les termes et indx contient les indice de permutation * une fois la decomposition effectué REAL*8 A(NN+1,NN+1),D INTEGER INDX(NN+1) ENDSEGMENT SEGMENT MV REAL*8 VV(NN+1) ENDSEGMENT SEGACT,MA*MOD N = A(/1) NN = N-1 SEGINI MV D=1.D0 DO 12 I=1,N AAMAX=0.D0 DO 11 J=1,N IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J)) 11 CONTINUE IF (AAMAX.EQ.0.) THEN c la matrice contient une colonne de termes nuls CALL ERREUR(-292) AAMAX = TINY ENDIF VV(I)=1.D0/AAMAX 12 CONTINUE DO 19 J=1,N IF (J.GT.1) THEN DO 14 I=1,J-1 SUM=A(I,J) IF (I.GT.1)THEN DO 13 K=1,I-1 SUM=SUM-A(I,K)*A(K,J) 13 CONTINUE A(I,J)=SUM ENDIF 14 CONTINUE ENDIF AAMAX=0.D0 DO 16 I=J,N SUM=A(I,J) IF (J.GT.1)THEN DO 15 K=1,J-1 SUM=SUM-A(I,K)*A(K,J) 15 CONTINUE A(I,J)=SUM ENDIF DUM=VV(I)*ABS(SUM) IF (DUM.GE.AAMAX) THEN IMAX=I AAMAX=DUM ENDIF 16 CONTINUE IF (J.NE.IMAX)THEN DO 17 K=1,N DUM=A(IMAX,K) A(IMAX,K)=A(J,K) A(J,K)=DUM 17 CONTINUE D=-D VV(IMAX)=VV(J) ENDIF INDX(J)=IMAX IF(J.NE.N)THEN IF(A(J,J).EQ. 0.D0) THEN * & 'Nuage pathologique pour la méthode d interpolation' * matrice singulière dans l'algorithme d'interpolation A(J,J)=TINY CALL ERREUR(-292) ENDIF DUM=1.D0/A(J,J) DO 18 I=J+1,N A(I,J)=A(I,J)*DUM 18 CONTINUE ENDIF 19 CONTINUE IF(A(N,N).EQ. 0.D0) THEN * & 'Nuage pathologique pour la méthode d interpolation' * matrice singulière dans l'algorithme d'interpolation CALL ERREUR(-292) A(N,N)=TINY ENDIF SEGSUP MV RETURN END