*IF DEF,C93_1A,OR,DEF,C93_2A GNF0F402.16
C ******************************COPYRIGHT****************************** GTS2F400.5815
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.5816
C GTS2F400.5817
C Use, duplication or disclosure of this code is subject to the GTS2F400.5818
C restrictions as set forth in the contract. GTS2F400.5819
C GTS2F400.5820
C Meteorological Office GTS2F400.5821
C London Road GTS2F400.5822
C BRACKNELL GTS2F400.5823
C Berkshire UK GTS2F400.5824
C RG12 2SZ GTS2F400.5825
C GTS2F400.5826
C If no contract has been raised with this copy of the code, the use, GTS2F400.5827
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.5828
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.5829
C Modelling at the above address. GTS2F400.5830
C ******************************COPYRIGHT****************************** GTS2F400.5831
C GTS2F400.5832
CLL SUBROUTINE MATINV-------------------------------------------- MATINV1A.3
CLL MATINV1A.4
CLL MODEL MODIFICATION HISTORY FROM MODEL VERSION 3.0: MATINV1A.5
CLL VERSION DATE MATINV1A.6
CLL vn4.2 07/11/96:Allow this routine to be used in C93_2A (Farnon) GNF0F402.15
CLL MATINV1A.7
CLL PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4, MATINV1A.8
CLL VERSION 2, DATED 18/01/90 MATINV1A.9
CLL MATINV1A.10
CLL SYSTEM TASK: ???? MATINV1A.11
CLL MATINV1A.12
CLL PURPOSE: INVERTS AN N BY N MATRIX BY ROW REDUCTION MATINV1A.13
CLL MATINV1A.14
CLL DOCUMENTATION: ??? MATINV1A.15
CLL MATINV1A.16
CLLEND------------------------------------------------------------- MATINV1A.17
C MATINV1A.18
C*L ARGUMENTS:----------------------------------------------------- MATINV1A.19
SUBROUTINE MATINV 1MATINV1A.20
& (C,B,N) MATINV1A.21
C MATINV1A.22
IMPLICIT NONE MATINV1A.23
C MATINV1A.24
INTEGER MATINV1A.25
* N ! IN SIZE OF MATRIX TO BE INVERTED MATINV1A.26
C MATINV1A.27
REAL MATINV1A.28
* C(N,N) ! IN MATRIX TO BE INVERTED MATINV1A.29
*, B(N,N) ! OUT INVERTED MATRIX MATINV1A.30
C*----------------------------------------------------------------- MATINV1A.31
C MATINV1A.32
C*L WORKSPACE USAGE------------------------------------------------ MATINV1A.33
REAL MATINV1A.34
* A(N,N) MATINV1A.35
C*----------------------------------------------------------------- MATINV1A.36
C MATINV1A.37
C*L EXTERNAL SUBROUTINES CALLED------------------------------------ MATINV1A.38
C NONE MATINV1A.39
C*----------------------------------------------------------------- MATINV1A.40
C MATINV1A.41
C------------------------------------------------------------------ MATINV1A.42
C DEFINE LOCAL VARIABLES MATINV1A.43
C------------------------------------------------------------------ MATINV1A.44
INTEGER MATINV1A.45
* M1,M2,J,M ! LOOP COUNTERS MATINV1A.46
*, L MATINV1A.47
C MATINV1A.48
REAL MH050593.1
* AMM ! = A(M,M) MATINV1A.50
*, AMJ ! = A(M,J) MATINV1A.51
C MATINV1A.52
C MATINV1A.53
C------------------------------------------------------------------ MATINV1A.54
CL 1. SET UP EXTENDED MATRIX (A|B)=(C|I) TO ROW REDUCE MATINV1A.55
C------------------------------------------------------------------ MATINV1A.56
DO 400 M1=1,N MATINV1A.57
DO 401 M2=1,N MATINV1A.58
IF (M1.EQ.M2) THEN MATINV1A.59
B(M1,M2)=1 MATINV1A.60
ELSE MATINV1A.61
B(M1,M2)=0 MATINV1A.62
ENDIF MATINV1A.63
A(M1,M2)=C(M1,M2) MATINV1A.64
401 CONTINUE MATINV1A.65
400 CONTINUE MATINV1A.66
C MATINV1A.67
C------------------------------------------------------------------ MATINV1A.68
CL 2. FOR EACH COLUMN M FIND THE FIRST NON-ZERO ENTRY. MATINV1A.69
CL ADD ROW WITH NON-ZERO ENTRY TO ROW M SO THAT A(M,M) IS NOT 0. MATINV1A.70
CL IF COLUMN M CONSISTS OF ZEROS THEN MATRIX CANNOT BE INVERTED, MATINV1A.71
CL SO STOP. MATINV1A.72
C------------------------------------------------------------------ MATINV1A.73
DO 2 M=1,N MATINV1A.74
L=M MATINV1A.75
200 CONTINUE MATINV1A.76
IF (A(M,M).EQ.0.) THEN MATINV1A.77
IF (L.EQ.N) STOP 5 MATINV1A.78
L=L+1 MATINV1A.79
DO 402 M1=1,N MATINV1A.80
A(M1,M)=A(M1,M)+A(M1,L) MATINV1A.81
B(M1,M)=B(M1,M)+B(M1,L) MATINV1A.82
402 CONTINUE MATINV1A.83
ENDIF MATINV1A.84
IF (A(M,M).EQ.0.) GOTO 200 MATINV1A.85
C MATINV1A.86
C------------------------------------------------------------------ MATINV1A.87
CL 3. DIVIDE ROW BY A(M,M) IN ORDER TO OBTAIN A(M,M)=1 MATINV1A.88
C------------------------------------------------------------------ MATINV1A.89
AMM=A(M,M) MATINV1A.90
DO 403 M1=1,N MATINV1A.91
A(M1,M)=A(M1,M)/AMM MATINV1A.92
B(M1,M)=B(M1,M)/AMM MATINV1A.93
403 CONTINUE MATINV1A.94
C MATINV1A.95
C------------------------------------------------------------------ MATINV1A.96
CL 4. IN COLUMN M MAKE EVERY VALUE 0 (EXCEPT A(M,M)) BY SUBTRACTING MATINV1A.97
CL ROW M*THE APPROPRIATE VALUE. MATINV1A.98
C------------------------------------------------------------------ MATINV1A.99
DO 3 J=1,N MATINV1A.100
IF (J.NE.M) THEN MATINV1A.101
AMJ=A(M,J) MATINV1A.102
DO 404 M1=1,N MATINV1A.103
A(M1,J)=A(M1,J)-AMJ*A(M1,M) MATINV1A.104
B(M1,J)=B(M1,J)-AMJ*B(M1,M) MATINV1A.105
404 CONTINUE MATINV1A.106
ENDIF MATINV1A.107
3 CONTINUE MATINV1A.108
2 CONTINUE MATINV1A.109
C MATINV1A.110
CL EXTENDED MATRIX (A|B)=(I|C**(-1)) WHERE C**(-1)=INVERSE OF C MATINV1A.111
RETURN MATINV1A.112
END MATINV1A.113
*ENDIF MATINV1A.114