- Email: [email protected]
FORTRAN IV subroutines for the inversion of MAGSAT data using an algorithm of one-dimensional arrays
Computers & Geosciences Vol. l 1, No. 1, pp. 79-83, 1985
0098-3004185 $3.00 + .00 © 1985 Pergamon Press Ltd.
Printed in the U.S.A.
F O R T R A N IV S U B R O U T I N E S F O R THE I N V E R S I O N OF M A G S A T D A T A U S I N G A N A L G O R I T H M OF ONE-DIMENSIONAL ARRAYS K. N. N. RAo, N. K. THAKUR, a n d P. K. AGRAWAL National Geophysical Research Institute, Hyderabad 500 007 India (Received 13 November 1983; revised 18 May 1984)
INTRODUCTION
v = ~
The magnetic field satellite (MAGSAT) was launched by NASA specifically to investigate large- and intermediate-scale lithospheric anomalies. Quantitative interpretation and modeling of such data are yet in its infancy. Attempts have been made (Mayhew, Johnson, and I.angel, 1980) to invert the anomalies to arrive at causative magnetization distributions within the lithospheric region. Efforts in this direction generally are foiled as the inversion procedure requires large computer memory. In this study, a suitable algorithm has been developed to store the leastsquares symmetric matrix as a I-D array, which considerably economizes the computer memory. These least-squares equations are solved using the Gauss-Seidel iteration procedure to arrive at the magnetization of the causative body.
[m,~7a - 7') - mo,rb + m~,rc]/P.
(2)
The magnetic induction B at P is B = -grad v
_( Or ' 07'
av 7nO'
av ) 7si~Oax
'
(3)
where the component 0v/by can be expressed (Meyer and others, 1983) as
07 = ~
m,
7~ + 37 - - F - - / -
37'
METHODOLOGY An equal area grid is laid out on the surface of the Earth to cover the region under investigation. A set of dipoles w~re assumed to be situated at the center of each square mesh. The magnetic moment of these dipoles is changed to give the least-squares best fit to the observed field in the space above the grid. In practice, the magnetic moment is divided by the volume (determined by the dipole spacing and an arbitrary depth (40 kin) below the Earth's surface) to compute apparent magnetization values. The contribution of any c n m a l dipole with magnetic moment m(7', 0', ),') to the potential at a point P(y, O, X) is given by v(7, 0 , ),) = m m l cos (rn, 1) 4~r Is '
where 1 = ('r: + 7 '2 - 2"rT'a) m a = cos 0 cos 0 ' + sin 0 sin 0 ' c o s (>, - X') b = cos O sin O' - sin O cos O' cos O, - ),') c = sin O sin (~ - ~') m~,= - - m s i n l me, = - m cos I cos D rnx, = rn cos I sin D. D and I are the declination and inclination of the magnetic field, respectively. The anomaly field due to a set of k dipoles is compared with the measured anomaly at the n satellite points. The computed field at the jth field point can be written as
(1)
where l is the vector distance b e t w ~ n the positions of rn and P and v.o is the magnetic permeability of vacuum. The direction cosines can be expressed in terms of the known coordinates of m and P. The potential is expressed as
k
as= ~ P~F~j j - - 1 , 2 , 3 . . . . . n i=!
79
(5)
K. N. N. RAO, N. K. THAKUR,and P. K. AGRAWAL
80
where Pi = rnivi = magnetic moment of the ith pole, Fij = the dipole source function and is given by eq. (4). The problem is to solve for a set of parameters re,which minimizes the squared sum of residuals between the observed (Sj) and computed anomaly over all data points and the solution is given by n
o
k
{j~[~ ~P,F,.j] 2)
o,
i
k.(6)
1,2,
i= t
ALGORITHM In practice, generally for n observed data values and k dipoles we need an n × k array for the generation of the source geometry. Using this n × k array, the least-squares solution for the unknowns requires another k × k matrix. Because the leastsquares matrix is symmetrical, only k × (k + 1)/2 elements of the coefficient matrix need to be stored. Knowing the geometry of the source function from eq. (4), the values of the function also are generated
so °
\
ss"
96 °
rJ \
7'2 °
as a one-dimensional array which forms the input for the least-squares solution. The least-squares equations are solved using the Gauss=Seidel iterative procedure storing the coefficient matrix as a I-D array. The two subroutines COEFF and ITERN are given in the Appendix, These subroutines were developed in FORTRAN IV for a microprocessor-based desktop system with a maximum memory of 64 kB. These programs can be used on other computer systems with a m i n i m u m of modification. The subroutine COEFF generates the least=squares solution matrix as a I-D array. The subroutine ITERN solves the system of least-squares equations. Using the two subroutines the MAGSAT data inversion problem has been solved for an array of 42 dipoles with 63 observed data points. An example. MAGSAT data after filtering core field at 400 km altitude over the Indian subcontinent were analyzed using the subroutines and the results are given in Figures 1 and 2. The program reproduces the observed field values (Fig. 1) with high accuracy (error < 5%). The magnetization distribution (Fig. 2) for the observed anomalies correlates well with known
r" 80 °
8~ °
96 ~
Figure 1. Residual vertical intensity anomaly map averaged at 4° square grid and at 400 km elevation. Contour interval, is 1 nT.
FORTRAN IV subroutines for the inversion of MAGSAT data 7~
6 4 ='
/
80 °
o_
/,"
//
88 =
/
/
II
~
II
,"
.,.,""~.~ ~"°?
/ / / ~ \ \ \
k.,.-"~ 'xlll I I / I / ~
."
96 °
""--.'~-'~
:-~ ...
I
81
:llllllll/l
I !
..-L~
.... ~<~.,_.-.
\
} ] IIIII ,._.,~,x~v//
l
II/
1 ,:
Ii
1
J[
I
II I
le01
72°
80=
BB°
Figure 2. Apparent magnetization distribution obtained through inversion of residual vertical intensity data. Contour interval is 0. l A/re.
geological features of the Indian landmass and its environment.
Acknowledgmants--The authors are thankful to Dr. J. G. Negi, for his encouragement and critically reading the manuscript and to Mr. Ramesh K.hanna, Mr. J. P. Sastry> and Mr. P. V. Swamy for their help in the preparation of the manuscript. The permi~on accorded by the Director, NGRI, Hyderahad to publish the work is gratefully acknowledged.
REFERENCES Mayhew, M. A., Johnson, B. D., and Langel, R. A., 1980, An equivalent source model of the satellite-altitude magnetic anomaly field over Australia: Earth PI. ScL Lett., v. 51, no. 1, p. 189-198. Meyer, J., Hufen, J. H., Sibert, M., and Hahn, A., 1983, Investigation of the internal geomagnetic field by means of a global model of the Earth's crust: Jour. Geophys., v. 52, no. 2, p. 71-84.
82
K.N.N.R.AO, N. K. THAKUR, and P. K. AGRAWAL APPENDIX
PROGRAM LISTINGS 01 C 02 C 03 C 04 C 05 C 06 C 07 C 08 C 09 C 10 C 11 C 12 C 13 14 15 16 17 C 18 C 19 C 20 C 21 C 22 23 C 24 C 25 26 27 28 C 29 30 31 32 33 C 34 35 36 30 37 38 20 39 C 40 41 42 43 44 45 46 40 47 48 10 49 C 50 51 52 53
SUBROUTINE COEFF . . . . GENERATES COEFFICIENT MATRIX FOR THE GENERALISED INVERSE PROBLEM BY LEAST SQUARES METHOD INPUTS - - -ARRAYS X ( I ) AND F ( I ) THE ARRAY X ( I ) CORRESPONDS TO THE SOURCE DIPOLE CONFIGURATION. IN FACT I T IS A TWO DIMENSIONAL ARRAY OF N ROWS AND K COLUMNS WHICH HAS BEEN STORED AS A ONE DIMENSIONAL ARRAY. .... F ( I ) THE OBSERVED FIELD DATA VALUES. OUTPUT . . . . A ( I ) THE ELEMENTS OF THE UPPER TRIANGULAR MATRIX (SYMMETRIC) AND Y ( K ) THE RHS CONSTANT VECTOR, N .... NUMBER OF FIELD DATA POINTS. K .... NUMBER OF DIPOLES. SUBROUTINE COEFF (N, K)
DIMENSION A(42) COMMON /FIELD/ F(63) COMMON /INPUT/ X(2646), Y(42) IN THE NEXT STATEMENT CALL OPEN FACILITATES THE CREATION OF MEMORY IMAGE FILE OF THE MATRIX A(1) AND THE VECTOR Y(1) ON THE SECONDARY MEMORY DEVICE. THESE TWO VECTORS ARE READ INTO THE MEMORY LOCATIONS OCCUPIED BY X(1) AND Y(1) BEFORE ENTERING INTO THE SUBROUTINE ITERN. CALL OPEN (6, 'COEFFYYYDAT', I) GENERATION OF THE COEFFICIENTS OF THE LEAST SQUARES SYMMETRIC MATRIX, ONLY UPPER TRIANGULAR ELEMENTS NEED BE GENERATED. DO 10 J = 1, K DO 20 JJ = J, K SUM = O. 0 INITIALISES (II) AND (KK) TO PICK UP THE REQUIRED ELEMENT II = J KK = JJ DO 30 I = I, N S U M = SUM+X(II)tX(KK) IN THE NEXT TWO STATEMENTS (If) AND (KK) ARE INCREMENTED II = II + K KK = KK + K CONTINUE A(JJ) = SUM CONTINUE
WRITES THE GENERATED ELEMENTS ONTO THE SECONDARY MEMORY DEVICE WRITE(6) (A(M), M=J, K) SUMY = 0.0 II = J DO 40 I = I, N SUMY = SUMY + F(!)*X(II) II = II + K CONTINUE Y ( J ) = SUMY CONTINUE WRITES THE RHS CONSTANT VECTOR ON THE SECONDARY MEMORY DEVICE WRITE(6) (Y(M), M=I, K) ENDFILE 6 RETURN END
FORTRAN IV subroulines for the inversion of MAGSAT data
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 50 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
C C C C C C C C C C C
C 757
C C C 556
202 101 C
21 20
121 122
SUBROUTINE ITERN . . . . SOLVES SYSTEM OF LINEAR EQUATIONS BY GAUSS-SEIDEL ITERATION METHOD. INPUT . . . . THE ARRAY A ( K ) THE TWO DIMENSIONAL SYMMETRIC COEFFICIENT MATRIX I S STORED AS A ONE DIMENSIONAL ARRAY B(I) .... RHS CONSTANT VECTOR. X(I) .... THE SOLUTION VECTOR. INC(I) .... A LOCAL ARRAY. E(I) .... ERROR ESTIMATES. APPROX(I) .... I N I T I A L APPROXIMATION TO THE SOLUTION VECTOR. N .... NUMBER OF ITERATIONS AND ACRY: THE ACCURACY FACTOR. K .... ORDER OF THE MATRIX. SUBROUTINE ITERN (K)
DIMENSION E(42), INC(42) COMMON /INPUT/ A(2646), B(42) COMMON /OUTPUT/ X(42) COMMON /SOLN/ N, ACRY, APPROX(42) INITIALISES THE SOLUTION VECTOR WITH STARTING APPROXIMATION. DO 757 I = I, K X(1) = APPROX(1) NN = K - I DO 555 IT = I, N ERROR = 0.0 THE FOLLOWING ELEVEN STATEMENTS HAVE BEEN DESIGNED TO GENERATE THE SUBSCRIPTS OF THE TWO DIMENSIONAL SYMMETRIC MATRIX WHICH HAS BEEN STORED AS A ONE DIMENSIONAL ARRAY. DO 556 1 = I, K INC(1) = I DO 10 1 = I, K II = I - I IF(I.EQ.I .OR. I.EQ.K) GO TO 101 L = NN DO 202 J = I, II INC(J) = L L=L-1 CONTINUE KX = I APPLICATION OF THE GAUSS-SEIDEL METHOD. SUM = 0.0 DO 20 J = I, K IF(I.EQ.J) GO TO 21 SUM = SUM + A(KX)*X(J) KX = K X + INC(J) CONTINUE KY = I I F ( I . E Q . I ) GOTO 122 DO 121 I X = 1, II KY = K Y + INC(IX) CONTINUE TEMP = (B(1)-SUM)/A(KY) DIFFER = X ( 1 ) - TEMP E(1) = DIFFER DIFFER
=
ABS(DIFFER)
IF(DIFFER.GT.ERROR)
10 555 666 C C C 667 668
83
ERROR = DIFFER
X(1) = TEMP CONTINUE IF(ERROR.LT.ACRY) GOTO666 CONTINUE CONTINUE THE FOLLOWING FOURSTATEMENTS ARE OPTIONAL. IF ONE WANTS THEY MAY BE SKIPPED. IN THESE STATEMENTS ERROR ESTIMATES FOR THE SOLUTION VECTOR AFTER N ITERATIONS ARE PRINTED, WRITE(2,668)
WRITE(2,667) (E(1), I = I , K) FORMAT(8EIO.2) FORMAT ( / l ERRORVECTOR ' / ) RETURN END
0098-3004185 $3.00 + .00 © 1985 Pergamon Press Ltd.
Printed in the U.S.A.
F O R T R A N IV S U B R O U T I N E S F O R THE I N V E R S I O N OF M A G S A T D A T A U S I N G A N A L G O R I T H M OF ONE-DIMENSIONAL ARRAYS K. N. N. RAo, N. K. THAKUR, a n d P. K. AGRAWAL National Geophysical Research Institute, Hyderabad 500 007 India (Received 13 November 1983; revised 18 May 1984)
INTRODUCTION
v = ~
The magnetic field satellite (MAGSAT) was launched by NASA specifically to investigate large- and intermediate-scale lithospheric anomalies. Quantitative interpretation and modeling of such data are yet in its infancy. Attempts have been made (Mayhew, Johnson, and I.angel, 1980) to invert the anomalies to arrive at causative magnetization distributions within the lithospheric region. Efforts in this direction generally are foiled as the inversion procedure requires large computer memory. In this study, a suitable algorithm has been developed to store the leastsquares symmetric matrix as a I-D array, which considerably economizes the computer memory. These least-squares equations are solved using the Gauss-Seidel iteration procedure to arrive at the magnetization of the causative body.
[m,~7a - 7') - mo,rb + m~,rc]/P.
(2)
The magnetic induction B at P is B = -grad v
_( Or ' 07'
av 7nO'
av ) 7si~Oax
'
(3)
where the component 0v/by can be expressed (Meyer and others, 1983) as
07 = ~
m,
7~ + 37 - - F - - / -
37'
METHODOLOGY An equal area grid is laid out on the surface of the Earth to cover the region under investigation. A set of dipoles w~re assumed to be situated at the center of each square mesh. The magnetic moment of these dipoles is changed to give the least-squares best fit to the observed field in the space above the grid. In practice, the magnetic moment is divided by the volume (determined by the dipole spacing and an arbitrary depth (40 kin) below the Earth's surface) to compute apparent magnetization values. The contribution of any c n m a l dipole with magnetic moment m(7', 0', ),') to the potential at a point P(y, O, X) is given by v(7, 0 , ),) = m m l cos (rn, 1) 4~r Is '
where 1 = ('r: + 7 '2 - 2"rT'a) m a = cos 0 cos 0 ' + sin 0 sin 0 ' c o s (>, - X') b = cos O sin O' - sin O cos O' cos O, - ),') c = sin O sin (~ - ~') m~,= - - m s i n l me, = - m cos I cos D rnx, = rn cos I sin D. D and I are the declination and inclination of the magnetic field, respectively. The anomaly field due to a set of k dipoles is compared with the measured anomaly at the n satellite points. The computed field at the jth field point can be written as
(1)
where l is the vector distance b e t w ~ n the positions of rn and P and v.o is the magnetic permeability of vacuum. The direction cosines can be expressed in terms of the known coordinates of m and P. The potential is expressed as
k
as= ~ P~F~j j - - 1 , 2 , 3 . . . . . n i=!
79
(5)
K. N. N. RAO, N. K. THAKUR,and P. K. AGRAWAL
80
where Pi = rnivi = magnetic moment of the ith pole, Fij = the dipole source function and is given by eq. (4). The problem is to solve for a set of parameters re,which minimizes the squared sum of residuals between the observed (Sj) and computed anomaly over all data points and the solution is given by n
o
k
{j~[~ ~P,F,.j] 2)
o,
i
k.(6)
1,2,
i= t
ALGORITHM In practice, generally for n observed data values and k dipoles we need an n × k array for the generation of the source geometry. Using this n × k array, the least-squares solution for the unknowns requires another k × k matrix. Because the leastsquares matrix is symmetrical, only k × (k + 1)/2 elements of the coefficient matrix need to be stored. Knowing the geometry of the source function from eq. (4), the values of the function also are generated
so °
\
ss"
96 °
rJ \
7'2 °
as a one-dimensional array which forms the input for the least-squares solution. The least-squares equations are solved using the Gauss=Seidel iterative procedure storing the coefficient matrix as a I-D array. The two subroutines COEFF and ITERN are given in the Appendix, These subroutines were developed in FORTRAN IV for a microprocessor-based desktop system with a maximum memory of 64 kB. These programs can be used on other computer systems with a m i n i m u m of modification. The subroutine COEFF generates the least=squares solution matrix as a I-D array. The subroutine ITERN solves the system of least-squares equations. Using the two subroutines the MAGSAT data inversion problem has been solved for an array of 42 dipoles with 63 observed data points. An example. MAGSAT data after filtering core field at 400 km altitude over the Indian subcontinent were analyzed using the subroutines and the results are given in Figures 1 and 2. The program reproduces the observed field values (Fig. 1) with high accuracy (error < 5%). The magnetization distribution (Fig. 2) for the observed anomalies correlates well with known
r" 80 °
8~ °
96 ~
Figure 1. Residual vertical intensity anomaly map averaged at 4° square grid and at 400 km elevation. Contour interval, is 1 nT.
FORTRAN IV subroutines for the inversion of MAGSAT data 7~
6 4 ='
/
80 °
o_
/,"
//
88 =
/
/
II
~
II
,"
.,.,""~.~ ~"°?
/ / / ~ \ \ \
k.,.-"~ 'xlll I I / I / ~
."
96 °
""--.'~-'~
:-~ ...
I
81
:llllllll/l
I !
..-L~
.... ~<~.,_.-.
\
} ] IIIII ,._.,~,x~v//
l
II/
1 ,:
Ii
1
J[
I
II I
le01
72°
80=
BB°
Figure 2. Apparent magnetization distribution obtained through inversion of residual vertical intensity data. Contour interval is 0. l A/re.
geological features of the Indian landmass and its environment.
Acknowledgmants--The authors are thankful to Dr. J. G. Negi, for his encouragement and critically reading the manuscript and to Mr. Ramesh K.hanna, Mr. J. P. Sastry> and Mr. P. V. Swamy for their help in the preparation of the manuscript. The permi~on accorded by the Director, NGRI, Hyderahad to publish the work is gratefully acknowledged.
REFERENCES Mayhew, M. A., Johnson, B. D., and Langel, R. A., 1980, An equivalent source model of the satellite-altitude magnetic anomaly field over Australia: Earth PI. ScL Lett., v. 51, no. 1, p. 189-198. Meyer, J., Hufen, J. H., Sibert, M., and Hahn, A., 1983, Investigation of the internal geomagnetic field by means of a global model of the Earth's crust: Jour. Geophys., v. 52, no. 2, p. 71-84.
82
K.N.N.R.AO, N. K. THAKUR, and P. K. AGRAWAL APPENDIX
PROGRAM LISTINGS 01 C 02 C 03 C 04 C 05 C 06 C 07 C 08 C 09 C 10 C 11 C 12 C 13 14 15 16 17 C 18 C 19 C 20 C 21 C 22 23 C 24 C 25 26 27 28 C 29 30 31 32 33 C 34 35 36 30 37 38 20 39 C 40 41 42 43 44 45 46 40 47 48 10 49 C 50 51 52 53
SUBROUTINE COEFF . . . . GENERATES COEFFICIENT MATRIX FOR THE GENERALISED INVERSE PROBLEM BY LEAST SQUARES METHOD INPUTS - - -ARRAYS X ( I ) AND F ( I ) THE ARRAY X ( I ) CORRESPONDS TO THE SOURCE DIPOLE CONFIGURATION. IN FACT I T IS A TWO DIMENSIONAL ARRAY OF N ROWS AND K COLUMNS WHICH HAS BEEN STORED AS A ONE DIMENSIONAL ARRAY. .... F ( I ) THE OBSERVED FIELD DATA VALUES. OUTPUT . . . . A ( I ) THE ELEMENTS OF THE UPPER TRIANGULAR MATRIX (SYMMETRIC) AND Y ( K ) THE RHS CONSTANT VECTOR, N .... NUMBER OF FIELD DATA POINTS. K .... NUMBER OF DIPOLES. SUBROUTINE COEFF (N, K)
DIMENSION A(42) COMMON /FIELD/ F(63) COMMON /INPUT/ X(2646), Y(42) IN THE NEXT STATEMENT CALL OPEN FACILITATES THE CREATION OF MEMORY IMAGE FILE OF THE MATRIX A(1) AND THE VECTOR Y(1) ON THE SECONDARY MEMORY DEVICE. THESE TWO VECTORS ARE READ INTO THE MEMORY LOCATIONS OCCUPIED BY X(1) AND Y(1) BEFORE ENTERING INTO THE SUBROUTINE ITERN. CALL OPEN (6, 'COEFFYYYDAT', I) GENERATION OF THE COEFFICIENTS OF THE LEAST SQUARES SYMMETRIC MATRIX, ONLY UPPER TRIANGULAR ELEMENTS NEED BE GENERATED. DO 10 J = 1, K DO 20 JJ = J, K SUM = O. 0 INITIALISES (II) AND (KK) TO PICK UP THE REQUIRED ELEMENT II = J KK = JJ DO 30 I = I, N S U M = SUM+X(II)tX(KK) IN THE NEXT TWO STATEMENTS (If) AND (KK) ARE INCREMENTED II = II + K KK = KK + K CONTINUE A(JJ) = SUM CONTINUE
WRITES THE GENERATED ELEMENTS ONTO THE SECONDARY MEMORY DEVICE WRITE(6) (A(M), M=J, K) SUMY = 0.0 II = J DO 40 I = I, N SUMY = SUMY + F(!)*X(II) II = II + K CONTINUE Y ( J ) = SUMY CONTINUE WRITES THE RHS CONSTANT VECTOR ON THE SECONDARY MEMORY DEVICE WRITE(6) (Y(M), M=I, K) ENDFILE 6 RETURN END
FORTRAN IV subroulines for the inversion of MAGSAT data
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 50 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
C C C C C C C C C C C
C 757
C C C 556
202 101 C
21 20
121 122
SUBROUTINE ITERN . . . . SOLVES SYSTEM OF LINEAR EQUATIONS BY GAUSS-SEIDEL ITERATION METHOD. INPUT . . . . THE ARRAY A ( K ) THE TWO DIMENSIONAL SYMMETRIC COEFFICIENT MATRIX I S STORED AS A ONE DIMENSIONAL ARRAY B(I) .... RHS CONSTANT VECTOR. X(I) .... THE SOLUTION VECTOR. INC(I) .... A LOCAL ARRAY. E(I) .... ERROR ESTIMATES. APPROX(I) .... I N I T I A L APPROXIMATION TO THE SOLUTION VECTOR. N .... NUMBER OF ITERATIONS AND ACRY: THE ACCURACY FACTOR. K .... ORDER OF THE MATRIX. SUBROUTINE ITERN (K)
DIMENSION E(42), INC(42) COMMON /INPUT/ A(2646), B(42) COMMON /OUTPUT/ X(42) COMMON /SOLN/ N, ACRY, APPROX(42) INITIALISES THE SOLUTION VECTOR WITH STARTING APPROXIMATION. DO 757 I = I, K X(1) = APPROX(1) NN = K - I DO 555 IT = I, N ERROR = 0.0 THE FOLLOWING ELEVEN STATEMENTS HAVE BEEN DESIGNED TO GENERATE THE SUBSCRIPTS OF THE TWO DIMENSIONAL SYMMETRIC MATRIX WHICH HAS BEEN STORED AS A ONE DIMENSIONAL ARRAY. DO 556 1 = I, K INC(1) = I DO 10 1 = I, K II = I - I IF(I.EQ.I .OR. I.EQ.K) GO TO 101 L = NN DO 202 J = I, II INC(J) = L L=L-1 CONTINUE KX = I APPLICATION OF THE GAUSS-SEIDEL METHOD. SUM = 0.0 DO 20 J = I, K IF(I.EQ.J) GO TO 21 SUM = SUM + A(KX)*X(J) KX = K X + INC(J) CONTINUE KY = I I F ( I . E Q . I ) GOTO 122 DO 121 I X = 1, II KY = K Y + INC(IX) CONTINUE TEMP = (B(1)-SUM)/A(KY) DIFFER = X ( 1 ) - TEMP E(1) = DIFFER DIFFER
=
ABS(DIFFER)
IF(DIFFER.GT.ERROR)
10 555 666 C C C 667 668
83
ERROR = DIFFER
X(1) = TEMP CONTINUE IF(ERROR.LT.ACRY) GOTO666 CONTINUE CONTINUE THE FOLLOWING FOURSTATEMENTS ARE OPTIONAL. IF ONE WANTS THEY MAY BE SKIPPED. IN THESE STATEMENTS ERROR ESTIMATES FOR THE SOLUTION VECTOR AFTER N ITERATIONS ARE PRINTED, WRITE(2,668)
WRITE(2,667) (E(1), I = I , K) FORMAT(8EIO.2) FORMAT ( / l ERRORVECTOR ' / ) RETURN END