Contents of file 'zgaussjdrv.f':
1 C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
2 C File: ZGAUSSJDRV.F (Fortran 77 source for the ZGAUSSDRV prog)
3 C Author: Fredrik Jonsson <fj@phys.soton.ac.uk>
4 C Date: December 29, 2005
5 C Last change: December 29, 2006
6 C Description: The ZGAUSSJDRV program is simply a driver for the
7 C ZGAUSSJ routine for solving complex-valued systems by
8 C Gauss-Jordan elimination. See the documentation of
9 C the ZGAUSSJ routine, enclosed in file ZGAUSSJ.F, for
10 C further information.
11 C
12 C Copyright (C) 2006, Fredrik Jonsson
13 C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
14 C PROGRAM MAIN
15 IMPLICIT LOGICAL (A-Z)
16 EXTERNAL zgaussj
17 INTEGER J,K,L
18 COMPLEX*8 A(3,3),AINV(3,3),T(3,3),B(3),X(3),D(3)
19
20 C *****************************************************************
21 C Define the complex-valued test matrix A of the system A*x=b.
22 C *****************************************************************
23 A(1,1)=CMPLX(1.4d0,3.3d0)
24 A(1,2)=CMPLX(-2.3d0,1.7d0)
25 A(1,3)=CMPLX(2.1d0,2.9d0)
26 A(2,1)=CMPLX(-3.4d0,2.3d0)
27 A(2,2)=CMPLX(1.3d0,-1.7d0)
28 A(2,3)=CMPLX(-3.2d0,-1.4d0)
29 A(3,1)=CMPLX(1.3d0,-2.3d0)
30 A(3,2)=CMPLX(2.8d0,1.1d0)
31 A(3,3)=CMPLX(-1.7d0,-1.2d0)
32
33 C *****************************************************************
34 C Define the complex-valued right-hand side b of the system A*x=b.
35 C *****************************************************************
36 B(1)=1.3
37 B(2)=-2.1
38 B(3)=0.7
39
40 C *****************************************************************
41 C Calculate the matrix inverse inv(A) using the zgaussj routine.
42 C *****************************************************************
43 DO 100,J=1,3
44 X(J)=B(J)
45 DO 110,K=1,3
46 AINV(J,K)=A(J,K)
47 110 CONTINUE
48 100 CONTINUE
49 CALL zgaussj(AINV,3,3,X,1,1)
50
51 C *****************************************************************
52 C Display the complex-valued system matrix A.
53 C *****************************************************************
54 WRITE (*,*) 'A:'
55 DO 130,J=1,3
56 WRITE(6,30) A(J,1),A(J,2),A(J,3)
57 130 CONTINUE
58
59 C *****************************************************************
60 C Display the inverted matrix inv(A).
61 C *****************************************************************
62 WRITE (*,*) ' '
63 WRITE (*,*) 'inv(A):'
64 DO 140,J=1,3
65 WRITE(6,30) AINV(J,1),AINV(J,2),AINV(J,3)
66 140 CONTINUE
67
68 C *****************************************************************
69 C As a check on the validity of the output of the zgaussj routine,
70 C verify that the matrix A times its inverse inv(A) yields the
71 C identity matrix.
72 C *****************************************************************
73 DO 150,J=1,3
74 DO 160,K=1,3
75 T(J,K)=CMPLX(0.0,0.0)
76 DO 170,L=1,3
77 T(J,K)=T(J,K)+A(J,L)*AINV(L,K)
78 170 CONTINUE
79 160 CONTINUE
80 150 CONTINUE
81 WRITE (*,*) ' '
82 WRITE (*,*) 'A*inv(A):'
83 DO 180,J=1,3
84 WRITE(6,30) T(J,1),T(J,2),T(J,3)
85 180 CONTINUE
86
87 C *****************************************************************
88 C Also perform the analogous check of inv(a)*A.
89 C *****************************************************************
90 DO 190,J=1,3
91 DO 200,K=1,3
92 T(J,K)=CMPLX(0.0,0.0)
93 DO 210,L=1,3
94 T(J,K)=T(J,K)+AINV(J,L)*A(L,K)
95 210 CONTINUE
96 200 CONTINUE
97 190 CONTINUE
98 WRITE (*,*) ' '
99 WRITE (*,*) 'inv(A)*A:'
100 DO 220,J=1,3
101 WRITE(6,30) T(J,1),T(J,2),T(J,3)
102 220 CONTINUE
103
104 C *****************************************************************
105 C Check the obtained solution x of A*x=b by computing the
106 C difference d=inv(A)*b-x.
107 C *****************************************************************
108 DO 230,J=1,3
109 D(J)=-X(J)
110 DO 240,K=1,3
111 D(J)=D(J)+AINV(J,K)*B(K)
112 240 CONTINUE
113 230 CONTINUE
114 WRITE (*,*) ' '
115 WRITE (*,*) 'inv(A)*b-x:'
116 WRITE(6,30) D(1),D(2),D(3)
117 WRITE (*,*) ' '
118
119 30 FORMAT(2X,'(',F11.8,',',F11.8,') (',F11.8,',',F11.8,') (',
120 + F11.8,',',F11.8,')')
121
122 STOP
123 END
124
Generated by ::viewsrc::