!***************************************************** !* Solving a complex linear system A Z = B, * !* where A is complex matrix of size N*N, * !* B is a complex matrix of size N*M, * !* Z is a complex solution matrix of size N*M. * !* ------------------------------------------------- * !* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp * !* and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * !* [BIBLI 05]. * !* ------------------------------------------------- * !* SAMPLE RUN: * !* * !* Solve complex linear system: * !* * !* z1 + z2 + iz3 = 5 * !* z1 + iz2 + z3 = -2i * !* iz1 + z2 + z3 = 1 * !* * !* Here: N=3, M=1 * !* * !* System solution: * !* * !* z1 = 1.500000 i -0.500000 * !* z2 = 1.000000 i 1.000000 * !* z3 = -0.500000 i -2.500000 * !* * !* DET = 20.0000000000000 * !* * !* F90 version by J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !***************************************************** ! Explanations: ! The original complex linear system to solve is: ! ( a11 a12 ... a1N ) ( z1 ) ( b1 ) ! ( a21 a22 ... a2N ) ( z2 ) ( b2 ) ! ( ... ... ... ... ) (... ) = (... ) ! ( aN1 aN2 ... aNN ) ( zN ) ( bN ) ! a k,l = c k,l + i d k,l ! z k = x k + i y k ! b k = e k + i f k ! where c, d, e, f, x, y are real numbers. ! The system is replaced by the following REAL system of size 2*N: ! ( c11 -d11 c12 -d12 ... c1N -d1N ) ( x1 ) ( e1 ) ! ( d11 c11 d12 c12 ... d1N c1N ) ( y1 ) ( f1 ) ! ( c21 -d21 c22 -d22 ... c2N -d2N ) ( x2 ) ( e2 ) ! ( d21 c21 d22 c22 ... d2N c2N ) ( y2 ) = ( f2 ) ! ( ... .... ... .... ... ... .... ) ( ...) ( ...) ! ( cN1 -dN1 cN2 -dN2 ... cNN -dNN ) ( xN ) ( eN ) ! ( dN1 cN1 dN2 cN2 ... dNN cNN ) ( yN ) ( fN ) ! The real system is then solved by a classic Gauss-Jordan method. ! PROGRAM Zmatsys parameter (NMAX=3, MMAX=1) real*8 Xm(1:NMAX,1:NMAX),Ym(1:NMAX,1:NMAX) real*8 Xv(1:NMAX), Yv(1:NMAX) real*8 M(1:2*NMAX,1:2*NMAX) real*8 V(1:2*NMAX,1:MMAX) integer i,j,li,co real*8 det li=NMAX; co=MMAX !read Xm matrix real part of A matrix Xm(1,1)= 1; Xm(1,2)=1; Xm(1,3)=0 ! 1 1 0 Xm(2,1)= 1; Xm(2,2)=0; Xm(2,3)=1 ! 1 0 1 Xm(3,1)= 0; Xm(3,2)=1; Xm(3,3)=1 ! 0 1 1 !read Ym matrix imaginary part of A matrix Ym(1,1)= 0; Ym(1,2)=0; Ym(1,3)=1 ! 0 0 1 Ym(2,1)= 0; Ym(2,2)=1; Ym(2,3)=0 ! 0 1 0 Ym(3,1)= 1; Ym(3,2)=0; Ym(3,3)=0 ! 1 0 0 !read Xv vector real part of 1st column of B matrix Xv(1)=5; Xv(2)=0; Xv(3)=1 ! 5 0 1 !read Yv vector imaginary part of 1st column of B matrix Yv(1)=0; Yv(2)=-2; Yv(3)=0 ! 0 -2 0 !The complex system AZ=B is replaced by real system MX=V of size 2*NMAX !create M matrix do i=1, NMAX do j=1, NMAX M(2*i-1,2*j-1)=Xm(i,j); M(2*i,2*j)=Xm(i,j) M(2*i,2*j-1)=Ym(i,j); M(2*i-1,2*j)=-Ym(i,j) end do end do !create V matrix do i=1, NMAX V(2*i-1,1)=Xv(i); V(2*i,1)=Yv(i) end do ! solve real linear system MX=V of zize 2*NMAX ! by Gauss-Jordan method - See program sysmat.f90 call MATINV(li*2, co, M, V, det) ! write complex solutions stored in vector V in the form: ! x1,y1,x2,y2,...,xn,yn, with zk = xk + i yk... print *,'' print *,' System solution:' print *,'' do i=1, NMAX write(*,100) i, V(2*i-1,1), V(2*i,1) end do print *,'' print *,' DET = ', det print *,'' print *,'' 100 format(' z',I2,' = ',F10.6,' i ',F10.6) END !******************************************* !* SOLVING A LINEAR MATRIX SYSTEM AX = B * !* with Gauss-Jordan method using full * !* pivoting at each step. During the pro- * !* cess, original A and B matrices are * !* destroyed to spare storage location. * !* --------------------------------------- * !* INPUTS: A MATRIX N*N * !* B MATRIX N*M * !* --------------------------------------- * !* OUTPUS: A INVERSE OF A N*N * !* DET DETERMINANT OF A * !* B SOLUTION MATRIX N*M * !* --------------------------------------- * !* NOTA - If M=0 inversion of A matrix * !* only. * !******************************************* SUBROUTINE MATINV(N,M,AA,BB,DET) PARAMETER(EPSMACH=1.D-20) REAL*8 AA(N,N),BB(N,M) REAL*8 PC(N),PL(N),CS(N) REAL*8 PV,PAV,DET,TT DET=1.D0 DO I=1,N PC(I)=0.D0 PL(I)=0.D0 CS(I)=0.D0 END DO !main loop DO K=1,N !Searching greatest pivot : PV=AA(K,K) IK=K JK=K PAV=DABS(PV) DO I=K,N DO J=K,N IF (DABS(AA(I,J)).GT.PAV) THEN PV=AA(I,J) PAV=DABS(PV) IK=I JK=J ENDIF ENDDO ENDDO !Search terminated, the pivot is in location I=IK, J=JK. !Memorizing pivot location: : PC(K)=JK PL(K)=IK !Determinant DET is actualised !If DET=0, ERROR MESSAGE and STOP !Machine dependant EPSMACH equals here 1.DE-20 IF (IK.NE.K) DET=-DET IF (JK.NE.K) DET=-DET DET=DET*PV IF (DABS(DET).LT.EPSMACH) THEN !Error message and Stop PRINT 10 STOP ENDIF !POSITIONNING PIVOT IN K,K: IF(IK.NE.K) THEN DO I=1,N !EXCHANGE LINES IK and K: TT=AA(IK,I) AA(IK,I)=AA(K,I) AA(K,I)=TT ENDDO ENDIF IF (M.NE.0) THEN DO I=1,M TT=BB(IK,I) BB(IK,I)=BB(K,I) BB(K,I)=TT ENDDO ENDIF !Pivot is at correct line IF(JK.NE.K) THEN DO I=1,N !Exchange columns JK and K of matrix AA TT=AA(I,JK) AA(I,JK)=AA(I,K) AA(I,K)=TT ENDDO ENDIF !Pivot is at correct column and located in K,K !Store column K in vector CS !then set column K to zero DO I=1,N CS(I)=AA(I,K) AA(I,K)=0.D0 ENDDO ! CS(K)=0. AA(K,K)=1. !Modify line K : IF(DABS(PV).LT.EPSMACH) THEN WRITE(*,*) ' PIVOT TOO SMALL - STOP' STOP ENDIF DO I=1,N AA(K,I)=AA(K,I)/PV ENDDO IF (M.NE.0) THEN DO I=1,M BB(K,I)=BB(K,I)/PV ENDDO ENDIF !Modify other lines of matrix AA: DO J=1,N IF (J.EQ.K) CONTINUE DO I=1,N !Modify line J of matrix AA : AA(J,I)=AA(J,I)-CS(J)*AA(K,I) ENDDO IF (M.NE.0) THEN DO I=1,M BB(J,I)=BB(J,I)-CS(J)*BB(K,I) ENDDO ENDIF ENDDO !Line K is ready. ENDDO !End of K loop !The matrix AA is inverted - Rearrange AA !Exchange lines DO I=N,1,-1 IK=PC(I) IF (IK.EQ.I) CONTINUE !EXCHANGE LINES I AND PC(I) OF AA: DO J=1,N TT=AA(I,J) AA(I,J)=AA(IK,J) AA(IK,J)=TT ENDDO IF (M.NE.0) THEN DO J=1,M TT=BB(I,J) BB(I,J)=BB(IK,J) BB(IK,J)=TT ENDDO ENDIF !NO MORE EXCHANGE NEEDED !GO TO NEXT LINE ENDDO !EXCHANGE COLUMNS DO J=N,1,-1 JK=PL(J) IF (JK.EQ.J) CONTINUE !EXCHANGE COLUMNS J AND PL(J) OF AA : DO I=1,N TT=AA(I,J) AA(I,J)=AA(I,JK) AA(I,JK)=TT ENDDO !NO MORE EXCHANGE NEEDED !GO TO NEXT COLUMN ENDDO !REARRANGEMENT TERMINATED. RETURN 10 FORMAT(///' DETERMINANT EQUALS ZERO, NO SOLUTION!') END ! end of file zmatsys.f90