{****************************************************** * LU decomposition routines called by test_clu.pas * * * * TPW UNIT version by J-P Moreau, Paris * * (for a general complex linear system). * * (www.jpmoreau.fr) * * --------------------------------------------------- * * Reference: * * * * "Numerical Recipes by W.H. Press, B. P. Flannery, * * S.A. Teukolsky and W.T. Vetterling, Cambridge * * University Press, 1986". * * * ******************************************************} UNIT CLU; INTERFACE Uses WinCrt; CONST ISIZE = 1024; TYPE {complex number} COMPLEX = Record R, I: REAL; {algebraic form} End; {complex vector} pCVEC = ^CVEC; CVEC = Array[0..ISIZE] of COMPLEX; {integer vector} pIVEC = ^IVEC; IVEC = Array[0..ISIZE] of Integer; PROCEDURE LUDCMP(VAR A:pCVEC; N:INTEGER; VAR INDX:pIVEC; VAR D:INTEGER; VAR CODE:INTEGER); PROCEDURE LUBKSB(VAR A:pCVEC; N:INTEGER; INDX: pIVEC; VAR B:pCVEC); Procedure CMUL(Z1,Z2:Complex; Var Z:Complex); IMPLEMENTATION Function CABS(Z:Complex): Real; Begin CABS := sqrt(Z.R*Z.R+Z.I*Z.I) End; Procedure CMUL(Z1,Z2:Complex; Var Z:Complex); Begin Z.R := Z1.R*Z2.R - Z1.I*Z2.I; Z.I := Z1.R*Z2.I + Z1.I*Z2.R End; Procedure CDIV(Z1,Z2:Complex; Var Z:Complex); Var d:Real; C:Complex; Begin d := Z2.R*Z2.R+Z2.I*Z2.I; if d<1E-10 then writeln(' Complex Divide by zero !') else begin C.R:=Z2.R; C.I:=-Z2.I; CMUL(Z1,C,Z); Z.R:=Z.R/d; Z.I:=Z.I/d end End; {**************************************************************** * Given a complex N x N matrix A, this routine replaces it by * * the LU decomposition of a rowwise permutation of itself. A * * and N are input. INDX is an output vector which records the * * row permutation effected by the partial pivoting; D is out- * * put as -1 or 1, depending on whether the number of row inter-* * changes was even or odd, respectively. This routine is used * * in combination with LUBKSB to solve linear equations or to * * invert a matrix. Return code is 1, if matrix is singular. * ****************************************************************} PROCEDURE LUDCMP; {Release with a complex matrix stored in a complex vector} Const NMAX=100; TINY=2.2e-16; VAR I,IMAX,J,K,N1: INTEGER; VV: pCVEC; AMAX, CONE, DUM, SUM, TMP: COMPLEX; BEGIN New(VV); CONE.R:=1.0; CONE.I:=0.0; D:=1; CODE:=0; N1:=N+1; for I:=1 to N do begin AMAX.R:=0.0; AMAX.I:=0.0; for J:=1 to N do IF CABS(A^[I*N1+J]) > CABS(AMAX) then AMAX:=A^[I*N1+J]; IF CABS(AMAX) < TINY THEN begin CODE := 1; exit end; {VV^[I] := 1.0 / AMAX} CDIV(CONE,AMAX,VV^[I]) end; { i loop } for J:=1 to N do begin for I:=1 to J-1 do begin SUM := A^[I*N1+J]; for K:=1 to I-1 do begin {SUM = SUM - A^[I*N1+K]*A^[K*N1+J] } CMUL(A^[I*N1+K],A^[K*N1+J],TMP); SUM.R:=SUM.R - TMP.R; SUM.I:=SUM.I - TMP.I end; A^[I*N1+J] := SUM END; { i loop } AMAX.R:=0.0; AMAX.I:=0.0; for I:=J to N do begin SUM := A^[I*N1+J]; for K:=1 to J-1 do begin {SUM = SUM - A^[I*N1+K]*A^[K*N1+J] } CMUL(A^[I*N1+K],A^[K*N1+J],TMP); SUM.R:=SUM.R - TMP.R; SUM.I:=SUM.I - TMP.I end; A^[I*N1+J] := SUM; {DUM = VV^[I]*ABS(SUM) } DUM.R:=VV^[I].R * CABS(SUM); DUM.I:=VV^[I].I * CABS(SUM); IF CABS(DUM) >= CABS(AMAX) THEN BEGIN IMAX := I; AMAX := DUM END END; { i loop } IF J <> IMAX THEN begin for K:=1 to N do begin DUM := A^[IMAX*N1+K]; A^[IMAX*N1+K] := A^[J*N1+K]; A^[J*N1+K] := DUM end; { k loop } D := -D; VV^[IMAX] := VV^[J] end; INDX^[J] := IMAX; IF CABS(A^[J*N1+J]) < TINY then begin A^[J*N1+J].R := TINY; A^[J*N1+J].I := 0.0 end; IF J <> N THEN begin {DUM = 1.0 / A^[J*N1+J] } CDIV(CONE,A^[J*N1+J],DUM); for I:=J+1 to N do begin {A^[I*N1+J] = A^[I*N1+J]*DUM } CMUL(A^[I*N1+J],DUM,TMP); A^[I*N1+J] := TMP end end END; { j loop } Dispose(VV) END; { subroutine LUDCMP } {****************************************************************** * Solves the set of N complex linear equations A . X := B. Here * * A is input, not as the matrix A but rather as its LU decomposi-* * tion, determined by the routine LUDCMP. INDX is input as the * * permutation vector returned by LUDCMP. B is input as the right-* * hand side vector B, and returns with the solution vector X. A, * * N and INDX are not modified by this routine and can be used for* * successive calls with different right-hand sides. This routine * * is also efficient for plain matrix inversion. * ****************************************************************** Release with a complex matrix stored in a complex vector } PROCEDURE LUBKSB; VAR SUM, TMP: COMPLEX; I,II,J,LL,N1: INTEGER; BEGIN II := 0; N1:=N+1; for I:=1 to N do begin LL := INDX^[I]; SUM := B^[LL]; B^[LL] := B^[I]; IF II <>0 THEN for J:=II to I-1 do begin {SUM = SUM - A^[I*N1+J]*B^[J] } CMUL(A^[I*N1+J],B^[J],TMP); SUM.R:=SUM.R - TMP.R; SUM.I:=SUM.I - TMP.I end ELSE IF CABS(SUM) <> 0.0 THEN II := I; B^[I] := SUM end; { i loop } for I:=N downto 1 do begin SUM := B^[I]; IF I < N THEN for J:=I+1 to N do begin {SUM = SUM - A^[I*N1+J]*B^[J] } CMUL(A^[I*N1+J],B^[J],TMP); SUM.R:=SUM.R - TMP.R; SUM.I:=SUM.I - TMP.I end; {B^[I] = SUM / A^[I*N1+I] } CDIV(SUM, A^[I*N1+I], B^[I]) end; { i loop } END; { subroutine LUBKSB } END. {of unit lu.pas End of file clu.pas}