OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
invert.f
Go to the documentation of this file.
1 c data set invert (entered from hardcopy, 10 Jan 93, JWhiting)
2  subroutine invert(a,b,n,l,c,ier)
3 
4 c invert -- inverts a matrix in place and solves a set of
5 c simultaneous linear equations.
6 
7 c input:
8 
9 c a(l,l) matrix to be inverted (inverse is returned in this
10 c same array).
11 c b(n) vector of right hand sides of linear equations
12 c n number of equations to be solved (1-20)
13 c l dimension of a, as declared in calling prog (l.ge.n)
14 
15 c output:
16 
17 c a(l,l) inverse of original matrix
18 c c(n) vector of solution to linear equations
19 c ier return code
20 c =0, normal return
21 c =1, singular matrix (matrix a may be altered)
22 c =2, invalid input (n.lt.1, or n.gt.20, or l.lt.n)
23 c (matrix a unaltered)
24 
25 c external references: none
26 
27 c M. Shear, Computer Sciences Corp. 3-30-75
28 c Minor revision to subroutine LINEAR. Name changed and error
29 c return ier=2 inserted. January 1976
30 
31 c Method
32 c Gauss-Jordan elimination with normalization and optimal pivoting.
33 
34 c Reference
35 c Shan S. Kuo, Numerical Methods and Computers,
36 c Addison-Wesley, 1965. pp. 168-169.
37 
38 c Notes
39 c 1. If only matrix inversion is desired, the vector of right hand
40 c sides may be set to zeroes.
41 c 2. Execution time varies approx. as n**3. For n=20, execution
42 c time is about .1 seconds (IBM S/360-95)
43 
44  implicit real*8 (a-h,o-z)
45  dimension a(l,l), b(1), c(1)
46 
47 c Internal declarations
48  dimension ipvot(20), index(20,2)
49 
50 c Initialize
51  if (n.lt.1 .or. n.gt.20 .or. l.lt.n) goto 902
52  do 10 i=1,n
53  c(i)=b(i)
54  10 ipvot(i)=0
55 
56 c Loop on eliminations
57  do 500 npass=1,n
58 
59 c Determine pivot element
60  amax=0.d0
61  do 30 j=1,n
62  if (ipvot(j).ne.0) goto 30
63  do 20 k=1,n
64  if (ipvot(k).ne.0) goto 20
65  if (dabs(amax).ge.dabs(a(j,k))) goto 20
66  irow=j
67  icol=k
68  amax=a(j,k)
69  20 continue
70  30 continue
71  if (amax.eq.0.d0) goto 900
72  ipvot(icol)=1
73 
74 c Normalization and exchange of rows
75  a(irow,icol)=1.d0
76  do 40 i=1,n
77  t=a(irow,i)
78  a(irow,i)=a(icol,i)
79  40 a(icol,i)=t/amax
80  t=c(irow)
81  c(irow)=c(icol)
82  c(icol)=t/amax
83  index(npass,1)=irow
84  index(npass,2)=icol
85 
86 c Perform elimination
87  do 100 i=1,n
88  if (i.eq.icol) goto 100
89  t=a(i,icol)
90  a(i,icol)=0.d0
91  c(i)=c(i)-c(icol)*t
92  do 90 j=1,n
93  90 a(i,j)=a(i,j)-a(icol,j)*t
94  100 continue
95  500 continue
96 
97 c Reorder the inverted matrix by exchanging columns
98  do 600 i=1,n
99  k=n-i+1
100  irow=index(k,1)
101  icol=index(k,2)
102  do 550 j=1,n
103  t=a(j,irow)
104  a(j,irow)=a(j,icol)
105  550 a(j,icol)=t
106  600 continue
107  ier=0
108  goto 999
109 
110  900 ier=1
111  goto 999
112 
113  902 ier=2
114 
115  999 return
116  end
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
subroutine invert(a, b, n, l, c, ier)
Definition: invert.f:3