^T 0f >

\ &k

^fAU Of *

CO NBS TECHNICAL NOTE 1126

U.S. DEPARTMENT OF COMMERCE / National Bureau of Standards

Problems Used in Testing the

Efficiency and Accuracy of the

Modified Gram-Schmidt

Least Squares Algorithm
•I

NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards' was established by an act of Congress on March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and technology

and facilitate their effective application for public benefit. To this end, the Bureau conducts

research and provides: (1) a basis for the Nation's physical measurement system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau's technical work is per- formed by the National Measurement Laboratory, the National Engineering Laboratory, and

the Institute for Computer Sciences and Technology.

THE NATIONAL MEASUREMENT LABORATORY provides the national system of

physical and chemical and materials measurement; coordinates the system with measurement

systems of other nations and furnishes essential services leading to accurate and uniform

physical and chemical measurement throughout the Nation's scientific community, industry,

and commerce; conducts materials research leading to improved methods of measurement,

standards, and data on the properties of materials needed by industry, commerce, educational

institutions, and Government; provides advisory and research services to other Government

agencies; develops, produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

Absolute Physical Quantities2 — Radiation Research — Thermodynamics and

Molecular Science — Analytical Chemistry — Materials Science. THE NATIONAL ENGINEERING LABORATORY provides technology and technical ser- vices to the public and private sectors to address national needs and to solve national problems; conducts research in engineering and applied science in support of these efforts; builds and maintains competence in the necessary disciplines required to carry out this research and technical service; develops engineering data and measurement capabilities;

provides engineering measurement traceability services; develops test methods and proposes

engineering standards and code changes; develops and proposes new engineering practices;

and develops and improves mechanisms to transfer results of its research to the ultimate user. The Laboratory consists of the following centers:

Applied Mathematics — Electronics and Electrical Engineering2 — Mechanical

Engineering and Process Technology2 — Building Technology — Fire Research —

Consumer Product Technology — Field Methods.

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts

research and provides scientific and technical services to aid Federal agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and

economy in Government operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing the Federal Information Processing Standards Program, developing Federal ADP standards

guidelines, and managing Federal participation in ADP voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal agencies; and

provides the technical foundation for computer-related policies of the Federal Government.

The Institute consists of the following centers:

Programming Science and Technology — Computer Systems Engineering.

'Headquarters and Laboratories at Gaithersburg, MD, unless otherwise noted; mailing address Washington, DC 20234.

2 Some divisions within the center are located at Boulder, CO 80303.
National Bureau oi Star

Library, F.-01 Admin.

tniciency

and Accuracy of the Modified

Problems Used in Testing the Efficiency

Gram-Schmidt Least Squares Algorithm

/ !

Roy H. Wampler

Center for Applied Mathematics

National Engineering Laboratory

National Bureau of Standards

Washington, D.C. 20234

U.S. DEPARTMENT OF COMMERCE, Philip M. Klutznick, Secretary

Luther H. Hodges, Jr., Deputy Secretary

Jordan J. Baruch, Assistant Secretary for Productivity, Technology and Innovation

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Issued August 1980
National Bureau of Standards Technical Note 1126

Nat. Bur. Stand. (U.S.), Tech. Note 1126, 83 pages (Aug. 1980) CODEN: NBTNAE

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1980

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402

Price $4.00

(Add 25 percent for other than U.S. mailing)
TABLE OF CONTENTS

Page

Abstract 1

1. Introduction 2

2. Testing Routines and Test Results 3

2.1. TEST1: Problems Constructed from Orthonormal Vectors . 4

2.2. TEST2: Problems Having Rank Deficiencies 5 2.3. TEST3: Problems Constructed from Hadamard Matrices . . 6

2.4. TEST4: Problems with Heavy Weights, Unsealed

and Scaled 10

2.5. TEST5: General Least Squares Problems 12

3. Structure Chart (Figure 1) 14

*• How to Omit Pivoting 14

5. Computing a Pseudoinverse 14

6. Miscellaneous Remarks 16

Tables I, II, III, IV, V 18

Acknowledgments 24

References 24

Appendix A. Listing of Program ..... 26

Appendix B. Listing of Data 58

Appendix C. Sample Output from Five Problems Using TEST5

.(Figures 2, 3, 4, 5, 6) 69

in

Problems Used In Testing the Efficiency and Accuracy

of the Modified Gram-Schmidt Least Squares Algorithm

Roy H. Wampler

Center for Applied Matematics

National Engineering Laboratory

National Bureau of Standards

Washington, DC 20234

ABSTRACT

In preparing "Algorithm 544: L2A and L2B, Weighted Least Squares

Solutions by Modified Gram-Schmidt with Iterative Refinement" for publication in ACM Transactions on Mathematical Software (Vol. 5, 1979), the

Fortran computer program was extensively tested. This note describes

the various types of problems which were used to explore the efficiency

and accuracy of this algorithm. The Fortran subprograms which performed

the various tests are listed in an appendix. Also listed are the data

used in executing the testing routines as well as typical output from

several different types of problems. Among the testing routines is one which is suitable for handling general linear least squares problems.

Here, the user has the option of scaling his raw data in order to mitigate the effects of ill-conditioning.

Key words: Algorithms; curve fitting; least squares; modified GramSchmidt; pseudoinverse; regression; statistics; test problems; test results.
1. INTRODUCTION

In preparing "Algorithm 544: L2A and L2B, Weighted Least Squares

Solutions by Modified Gram-Schmidt with Iterative Refinement" [23] for

publication in ACM Transactions on Mathematical Software , the Fortran

computer program was subjected to extensive testing. Section 5 of the

companion paper [22] which appeared in the same issue of ACM Transactions on Mathematical Software gave a very brief account of some of

the test results and also cited previously published relevant findings

of Wampler [21] and other authors. The purpose of the present paper is

to present a more detailed description of the various types of problems

which were used to explore the efficiency and accuracy of the GramSchmidt algorithm.

Attention is given to the execution time required to solve problems

of various sizes, to how the choice of a tolerance parameter can affect

the computed rank of a matrix, to the effect of heavy weights on computed solutions, and to how one can enlarge the class of problems which

this algorithm solves by scaling the raw data. The Fortran subroutines

which performed the various tests are listed in an appendix, together

with a driver program and input data. Several examples of output generated by the computer program are displayed.

The results described in this paper were obtained chiefly on the

Univac 1108 at the National Bureau of Standards under the 1100 Operating

System, Versions 33. R2 and 33R3.A (Exec 8). Some results obtained from

an IBM 360 Model 65 using the Fortran IV G level 21 compiler are also

included. The testing routines, data, and sample output listed in the

appendices may prove to be useful in implementing algorithms L2A and L2B

on other computing systems.

The reader is referred to [22] and [23] for a complete formulation

of the least squares problem, a description of the method of solution,

notes on the calling sequences, and a listing of the introductory com- ments of subroutines L2A and L2B.

The linear least squares problem may be briefly defined [22, p.

458] as follows: Given a real M x N matrix A of rank Nl < min(M,N), and

given a real M-vector b of observations, find the N-vector x of coefficients which minimizes the sum of squares of the residual vector 6 = b - Ax. If M > N and Nl = N, the solution is unique. If M < N or if Nl

< N a unique solution can be obtained by imposing the condition that the

vector x be of minimal Euclidean norm. In some problems, the solution

vector x is subject to Ml linear equality constraints. In others, b is a vector of observations whose components have unequal variances; here

an M-vector W of weights enters the calculations. (It should be noted

that the parameters M, N, Nl and Ml which are here printed in capital

letters were printed in lower case italics with subscripts (m, ja, nj,

HI) in Sections 1, 2 and 3 of [22].)

2
When matrix A is of rank N and M > N, algorithms L2A and L2B furnish the same solutions. When matrix A is of rank Nl < N, the two algorithms follow different paths and furnish different solutions. Properties of both algorithms are explored in the testing routines discussed in this paper.

In algorithm L2A, if the original matrix A is found to be of rank

* Nl < N, the original matrix A is replaced by a smaller matrix A which

consists of Nl columns of A found to be linearly independent. A solution is sought for the smaller system of equations.

In algorithm L2B, if the rank of the M x N matrix A is found to be

Nl < N, a unique solution is obtained by imposing the condition that x,

the N x 1 solution vector, be of minimal Euclidean norm.

The L2A/L2B algorithms are structured so that solutions can be

found for several right-hand sides corresponding to one given matrix A.

It is sometimes convenient to refer to the set of one or more b-vectors

as a matrix B, and to the set of x-vectors as a matrix X.

2. TESTING ROUTINES AND TEST RESULTS

The Fortran program listed in Appendix A consists of a main (driver) program, nine subroutines and one function. Depending upon two

user-furnished input parameters, the main program calls one of five

testing subroutines (named TEST1, TEST2, TEST3, TEST4, and TEST5) each

of which calls either L2A or L2B for computing solutions after generating or reading input data. Subroutines ORTHO and HADMAR and function

GEN are called to construct certain types of input matrices. Subroutine

SCALE is called by TEST4 and TEST5 to scale the raw input data when this

option is requested by the user. Subroutine OUTPUT prints computed so- lutions.

Subroutines L2A, L2B and the seven subroutines which they call are

not listed here. Introductory comments explaining these algorithms

appear in [23]. The complete listing is available in "Collected Algorithms from ACM" [6J. It is also available (as listing, card deck, and

magnetic tape) from ACM Algorithms Distribution Service, c/o International Mathematical and Statistical Libraries, Inc., GNB Building, Sixth

Floor, 7500 Bellaire Boulevard, Houston, Texas 77036.

The driver program reads in two input parameters, KIND and MODE.

The value of KIND (1, 2, 3, 4 or 5) determines which of the five testing

subroutines shall be called — TEST1, TEST2, TEST3, TEST4, or TEST5.

The value of MODE (1 or 2) determines whether subroutine L2A or L2B

shall be called. MODE 1 corresponds to L2A, and MODE 2 corresponds to

L2B.
The first four testing subroutines generate data matrices for

certain specific types of least squares problems. Whereas subroutine

TEST2 requires no input from the user, subroutine TEST1, TEST3 and TEST4

require some input parameters and/or data in order to construct the

matrices A and B which are input to L2A or L2B.

Subroutine TEST5 is a general-purpose testing routine which reads

in matrices A and B and other quantities from card images. Weights are read in when the user wishes to assign unequal weights to different observations. TEST5 is suitable for use in production runs as well as for

testing and development purposes.

A description will now be given of the various types of problems

which were considered and of the results which were obtained from the

five testing routines.

2.1 TEST1: Problems Constructed from Orthonormal Vectors

Subroutines TEST1 and ORTHO construct input matrices A and B for

least squares problems through the use of a set of orthonormal vectors.

The method for generating these vectors is described in Example 2 of

Hastings [9]. These problems were used to obtain data on the time re- quired to solve various-sized problems. The timing data is summarized

in Table I for a set of 60 problems run on the Univac 1108. Besides

giving the total time required for each problem, the table shows the

time spent in subroutines DECOMP, SOLVE, and COVAR. ("DECOMP" here denotes DEC0M1 in the case of L2A and DEC0M2 in the case L2B; "SOLVE" denotes SOLVE 1 in the case of L2A, and S0LVE2 and SOLVE 3 in the case of

L2B. ) Most of the remaining time not otherwise accounted for in the table was consumed by generating the data and printing the results. These

operations required from 36% (large-sized problems) to 77% (small-sized

problems) of the total time.

In all 60 problems the rank was found to equal N. All but two

problems converged to a solution after two iterations; the two cases where M = 150, N = 80, Ml = both required three iterations.

The data of element DATA1A were used as input to TEST1 for the 44

problems which called L2A (MODE 1). Element DATA1B gives the data for

the 16 problems which called L2B (M0DE2). These two sets of data were

processed in two separate runs. In each run it was necessary to modify

the DIMENSION statements and DATA statements in the driver program

(lines 141-145 and 152-158) to handle problems as large as M = 200, N -

80.

In the run which called L2A we used:

INTEGER IPIVOT(80)

REAL A(200,81), B(200,l), 0(1682)

REAL Q(200,80), R(80,80), QR(1,1)

REAL RES(200,1), SDX(80), SF(81), W(200), X(80,l)
LOGICAL FAIL(l)

DATA MM, NN, LL, MMPNN, NNP1, NNPLL, MNL

* / 200, 80, 1, 1, 81, 81, 1682 /

The run which called L2B used the same dimensions and dimensioning

parameters as those used in the L2A run except for

REAL Q(l,l), R(l,l), QR(280,80)

and in the DATA statement MMPNN was set equal to 280 rather than 1.

The timing data given in Table I were obtained from two batch runs

(one for MODE 1, one for MODE 2) made on February 16, 1978, at night un- der N (lowest) priorty on the Univac 1108 at NBS. These two jobs were

repeated on two other occasions in a similar environment, and the details of the timing were quite similar to those given in Table I. Runs

made during the day (prime time) in batch mode and runs made in interactive mode could be expected to require different amounts of time —

probably greater.

The chart below gives the total run time and the central processor

unit (CPU) time required for execution of these two jobs on three different occasions.

MODE 1 (L2A)

44 Problems

MODE 2 (L2B)

16 Problems

Run 1 Date

Total

CPU

time

2-16-78

4 min. 11.268

3 min. 40.461

sec.

sec.

2-16-78

2 min. 54.677 sec.

2 min. 32.982 sec.

Run 2

Date

Total

CPU

time

2-27-80

4 min. 1.003

3 min. 28.232

sec.

sec.

8-28-79

3 min. 0.761 sec.

2 min. 37.009 sec.

Run 3

Date

Total

CPU

time

2-29-80

4 min. 5.124

3 min. 28.907

sec.

sec.

2-27-80

2 min. 53.016 sec.

2 min. 29.238 sec.

2.2 TEST2: Problems Having Rank Deficiencies

Subroutine TEST2 and function GEN generate two sets of test problems as described by Lawson and Hanson [13, pp. 252, 277]. The first

set (18 cases) was constructed from a sequence of integers having a short period so that some of the matrices are mathematically rankdeficient. The second set of 18 cases is the same as the first set ex- cept that "noise," simulating data uncertainty, has been added to all

data values.
These two sets of problems were run using subroutine L2B with two

different values of the tolerance parameter TOL, namely 0.0 and 0.5.

The results obtained from L2B were compared with results from three of

the programs given in [13] which furnish least squares solutions:

PR0G1, PR0G2 and PR0G3. PR0G1 and PR0G2 both use Householder transformation algorithms, whereas PR0G3 uses a singular value decomposition.

The results summarized below were all obtained on a Univac 1108.

For the problems of the first set, whether TOL was set equal to 0.0

or 0.5, the computed rank agreed with that given by PR0G2 and PR0G3. In

six of the 18 cases, Nl < min(M, N). The agreement of computed coefficients obtained by the two programs L2B and PR0G2 ranged from 5 to 8

digits. Agreement between L2B and PR0G3 was also 5 to 8 digits. The

setting of TOL =0.5 in L2B corresponds to TAU =0.5 in PR0G2 and PR0G3.

For the problems of the second set, where noise was added to all

elements of A and B, when TOL was set equal to 0.5, the ranks computed

by L2B again agreed with those computed by PR0G2 and PR0G3. When TOL

was set equal 0.0, the rank in all cases was computed to be min(M, N).

This was in agreement with PR0G1 in all problems for which PR0G1 gave

results. The agreement of coefficients obtained by the two programs L2B

and PR0G2 ranged from 6 to 8 digits. Between L2B and PR0G3 , as well as between L2B and PR0G1, agreement ranged from 1 to 8 digits.

Table II lists the solutions, computed ranks and Euclidean norms of

residuals for the test case M = 7, N = 6, with the two different noise

levels. These results are very similar to those reported in Table C.3,

page 252 of [13]. (It should be noted that 6x6 matrix of coefficients

in Table C.3 of [13] gives incorrect signs in three instances. The

(3,3), (2,5) and (3,5) entries should read, respectively, +0.0833,

-0.0217 and +0.0846.)

No input data is required for subroutine TEST2. One can construct

additional problems of this type by changing one or more of the parameters NOISE, ANOISE, MNl and MN2 in TEST2.

2.3. TEST3: Problems Constructed from Hadamard Matrices

Subroutines TEST3 and HADMAR construct least squares problems from

Hadamard matrices by a method due to Richard J. Hanson [8] which makes

use of the singular value decomposition of a rectangular matrix. In

these examples, the problem size dimensions were varied, the condition

number of the matrix A was varied and the data uncertainty in the righthand sides was varied. Such problems can be devised in the following

manner:

(a) Let A = U S*VT

where A is M x N, U is M x M, S* is M x N, V

T is N x N, with M = 2

l

,

kj - 1,2,3,4, and N = 2 2

, k2 = 1 ,. . .

, k^
(b) The matrices U and V are symmetric Hadamard matrices.

(An Hadamard matrix H of order n is a matrix of +l's and -l's such that

HHT - nlj any two rows of H are orthogonal. Cf. Hall [7, pp. 204-221].

For order 2 one can take

H2 c i

Hadamard matrices of orders 4, 8 and 16 can be constructed from lower

order Hadamard matrices. Let H4 =• H2

x H2 , Hg = H2 x H4 and Hjg = H4 x H4, where x denotes the Kronecker product of two matrices.

(c) Let k3 be a positive integer, and let

n

S*

s

where S is a diagonal N x N matrix whose elements are uniformly

k,

distributed as integers on [1,2 •*]. Set one of the diagonal elements

k.

equal to 2

J and another equal to 1. The largest and smallest singular

k.

values of S are thus 2

J and 1, and hence the condition number of S is

2

J . It is easily shown that the condition number of A is also equal to

J^2 J . Values chosen for k3 depend upon the computer precision available.

(d) Let x = (1,1,... ,1)^ and b = A x. Next, compute a

residual vector

6 = U

_pj > M-N

where p is a vector of integers uniformly distributed on [0, ||b||], and

A A |b I denotes the Euclidean norm of b. With A defined as above, we have

1/2

b - (N,N,..., N) T and | |b| I - NM

(e) Finally, compute

b = b + oc6 for four values of a, namely 0, 0.01, 1, 100

The true solution for all problems constructed in this manner is

x = (1,1,...,1)T.
Sets of Hadamard-type problems were run on the Univac 1108 and the

IBM 360/65. On the Univac 1108, where the machine-dependent parameter

ETA was set equal to 2"~ 26

, we chose k3 = 12, 16, 20, 24, 25. On the IBM

360/65, with ETA set equal to 16~5

, we chose k3 = 12, 16, 20. (ETA is

the smallest positive real number such that 1.0 + ETA > 1.0 in floatingpoint arithmetic.) In all cases we used subroutine L2A with TOL = 0.

The values of p, having been obtained from [15] and [16], were read in

from data elements (DATA3A, DATA3B, DATA3C, DATA3D, DATA3E). DATA3A was used with k3 = 12; DATA3B with k3 = 16; DATA3C with k3 = 20; DATA3D with

k3 = 24; DATA3E with k3 = 25.

Tables III and IV summarize the results for four values of k3 on

the Univac. In Table III, results for ten combinations of M and N are condensed on each row of the table, whereas Table IV gives greater details for the cases where k3 = 24. That the error in the computed residuals was sometimes as great as 8 units in the seventh significant digit

when a = 0.01 is explained by the fact that the elements of b in these cases were not integers and could not, in general, be exactly represented within the computer as octal fractions. In all other examples

covered by Tables III and IV, the elements of A and b were integers

which could be converted exactly to octal fractions.

For the problems where k3 = 25, in seven of the ten (M,N)-combinations, the solution failed to converge for all values of a. For M =

16, N = 16 the computed rank was found to be 15, and a solution was then

given for the problem of reduced size. In the cases where M = 8, N = 8

and M = 16, N = 8, solutions were computed, but here the matrices A as represented in the computer were perturbed versions of the precise problems defined in steps (a) through (e) above. The largest element of A

in both these cases, when computed without rounding error, is 134217729.

By printing out matrix A in octal format, it could be seen that this

number was represented in the Univac as 134217728. The exact

solution for this perturbed problem has all elements of x equal to 8/7.

The computed coefficients (1.1428571) for M = 8, N = 8 and for M = 16,

N = 8, a = 0, 0.01, 1 were thus correct to 8 digits; for M = 16, N = 8, a = 100 the computed coefficients were in error by 5 units in the eighth

significant digit.

Table V summarizes the results for problems run on an IBM 360/65,

with k3 = 12, 16, 20. In all cases where k3 = 20, the matrix A was found to be rank deficient, with Nl = N -1. Solutions were then computed for the problems of reduced size.
For one of the Hadamard-type problems run on the Univac 1108,

further details are now given. Here, M = 8, N = 4; k.3 = 24; the

diagonals of S were sj = 1, S2 5592406, S3 = 11184811, S4 = 16777216;

p = (10, 7, 2, 10) T ; and o - 1.

We btain '"*

' 33554434 -11184810 -22369620 33

-11184810 33554434 -22369620 -1

-22369620 33554434 -11184810 9

A = -22369620 -11184810 33554434 , b = 15

33554434 -11184810 -22369620

9

** -25

-11184810 33554434 -22369620 9

-22369620 33554434 -11184810 -1

-22369620 -11184810 33554434 -7

With ETA = 2" 26 , TOL =0.0, the computed rank was Nl = 4 and a solution was obtained after 8 iterations.

Computed results were:

.99999999

.99999999

.99999999

.99999999

x = diagonals of covariance matrix =

6 =

29.000000

-5.0000000

5.0000000

11.000000

-29.000000

5.0000000

-5.0000000

-11.000000

0071733

,0071733

,0071733

,0071733
Accurately calculated to five significant digits the diagonal terms of

the unsealed covariance matrix are .0078125; the computed values,

.0071733, are correct to only 1.09 digits. This accuracy is close to

the value of DIGITX reported for this problem, 1.33. Recall that DIGITX

[22, p. 463] is an estimate of the accuracy of the initial solution of

x, before iterative refinement. There is no iterative refinement in the

computation of the unsealed covariance matrix.

2.4 TEST4: Problems with Heavy Weights, Unsealed and Scaled

In the problems considered thus far, all components of the

observation vectors were assumed to have equal variances. We now

consider problems where some observations are assumed to have smaller

variances than others, and hence have greater weights. These examples

were constructed by Bjorck [3] from the polynomial data of Wampler [20].

(a) Let A = (a^j) be the 21 x 6 matrix given by

aij = 1 if i = j = 1,

j-l

ajM = (i-1) otherwise.

This corresponds to fitting a fifth-degree polynomial to observations at

the points 0,1,2 20.

(b) Let Aw

- DWA where

Dw

= diag (w,l,l,..., 1, w,l,..., 1, w)

is a diagonal matrix which scales the three rows 1, 11, and 21 of A with

the weight w > 1.

(c) Table 1 of [20] lists a vector of observations called Y2

and a vector of residuals called A. Let y 10^ .Y2. The vector of

observations t^ is chosen as

-1

bw

= Dwy + wDw A.

i-1

(d) Let w - 2 , i = l,2,...,p, where the choice of p depends

upon the computer precision available. We thus have a sequence of

problems Ay^x = t^ with moderate to large residuals (depending on the

magnitute of w), the exact solution of which is

x - (105

, 104, 10 3, 102

, 10, 1) T .

(e) Consider also modifications of the problem defined thus

far in which the columns of Aw are scaled in the following manner.

(i) By G - diag(106

, 105

, 104, 10 3, 102

, 1).

10
Then (Aw G) xq = bw has the solution

-1 -1 -1 -1 -1 nT xG

= (10 ,10 ,10 ,10 ,10 , 1) T

from which one can calculate x = G xq.

(ii) By H = diagdlai-ir1

, I |a2 l I"1 ,... , I

I

6 | I"

1

)

where ||aj|| denotes the Euclidean norm of the jth

column of Aw

. From xjj, the solution to (A^) x^ =

bw , one can obtain x = H Xfj. Under this scaling

all columns of AW

H have unit length. (Cf. Lawson

and Hanson [13, pp. 185-188] and Bj'orck [3].)

Note that for w = 1 the elements of Aw and bw are integers having

up to seven significant digits. As w increases, the number of significant digits, especially in the weighted elements of bw , increases.

Subroutine TEST4 was used to construct the matrices A and bw by

reading in the data of element DATA4 and applying the formulas given in

(a), (b), (c), and (d) above. Subroutine SCALE was used to apply the

scaling described in (e) above.

For the problem defined in steps (a) through (d) above, solutions

from L2A obtained on a Univac 1108 with TOL =0.0 reported the computed

rank Nl to be 6 for w = 1, 2 and 4, but Nl < 5 was reported for w > 8.

With scaling (i) of step (e) introduced, solutions with Nl = 6 were

14

obtained for w < 2 . With scaling (ii), solutions with Nl = 6 were

16

obtained for w < 2

14

Example: w = 2 = 16384, scaling (ii), TOL = 0. Two kinds of perturbations occurred in expressing the data of this problem on the Univac. Firstly, there were two instances of integers having too many significant digits to be exactly expressible in octal represenation. These were the first and last elements of bw , equal to 1,638,400,759 and 103,

219,200,759, which were converted to 1,638,400,752 and 103,219,200,000,

respectively, in octal representation. Secondly, the scaling caused the

elements of AWH to be decimal fractions between zero and one. Such numbers cannot, in general, be converted exactly to octal fractions. The

exact solution to the perturbed problem, obtained with the aid of a double precision program, is listed below, together with the solution to

which L2A converged after four iterations.

11
Exact solution (to 8 digits) Solution computed by L2A

100000.00 99999.998

10035.918 10035.917

948.16813 948.16848

102.26114 102.26109

9.8705418 9.8705446

1.0025746 1.0025745

In discussing weighted problems such as the example considered

here, Bjorck [3] cautioned that the iterative refinement procedure would

fail to converge for extremely large values of w. In this particular

example we note, however, that rank-deficiency was reported for the larger values of w, and reports of failure to converge were not encountered.

2.5 TEST5: General Least Squares Problems

Subroutine TEST5 is a general-purpose routine which reads in matrices A and B as well as vector W if the observations are to be weighted.

Input parameters which are read in by TEST5 are:

ISCALE Denotes whether the raw data are to be scaled (by calling

subroutine SCALE) and if so, what type of scaling is to be

applied.

M Total number of equations.

N Number of unknown coefficients.

Ml Number of linear constraints. (If Ml > 0, the first Ml

equations are to be satisfied exactly.)

L Number of right-hand sides (vectors of observations).

ITYPE Parameter which specifies whether a polynomial type fit or

a multiple regression fit is to be performed.

IWGHT Parameter which specifies whether or not weights are to be

read in.

TOL Parameter determining the type of tolerance to be used in

determining rank.

Further remarks on these input parameters are given in the listing of

TEST5 in Appendix A as well as in Wampler [22; 23].

12
Appendix B lists data elements DATA5A, DATA5B, DATA5C, DATA5D and

DATA5E which were used in conjunction with TEST5.

Through the use of Sande's Fortran execution profiler [17] it was found that on the Univac 1108 each executable statement of L2A and L2B

(and the seven subroutines called by them) was used at least once in

executing the 24 problems listed in DATA5A and DATA5B.

Element DATA5C contains five problems, a subset of DATA5A, for

which the results from L2B are different from those of L2A.

Elements DATA5B, DATA5D and DATA5E illustrate that problems with

zero weights can be easily handled by this algorithm, and that the com- puted rank of the system of equations can be affected by scaling the raw

data.

The first example of DATA5D, taken from Bennett and Franklin [2, pp

379-385] illustrates how orthogonal designs with missing observations

can be treated by assigning weight zero to the missing cells. Whenever

the I-th weight W(I) 0, the I-th rows of A and B are ignored in com- puting the solution vector and the residual sum of squares. Residuals,

however, are computed corresponding to observations given zero weights.

The second example of DATA5D, taken from Hogben, Peavy and Varner

[10, pp. 134-135], is one where 5 of the 12 observations are given zero

weight. It was reported in [10] that the fact that the system of equations is singular was not detected by OMNITAB II in connection with a twoway analysis of variance. OMNITAB 80 [11], however, which uses an adaptation of algorithm L2A for least squares fits and for twoway analysis of variance detects and reports to the user the fact that this system is singular.

The last two examples (Nos. 27 and 28) of DATA5D are, respectively,

sixth and fifth degree polynomial-type problems. On the Univac 1108,

using unsealed data (ISCALE = 0), both problems were reported to be

rank-deficient. The computed ranks were 4 and 3, respectively. If the

user specifies that input matrices A and B should be automatically

scaled (ISCALE = 2), algorithm L2A reports for both problems that matrix

A is of full rank. Computed results for the fifth degree polynomial,

with and without scaling, are shown in Appendix C.

Another example of computed output in Appendix C is problem No. 22

from element DATA5B. This illustrates several features of the L2A/L2B

algorithms:

(1) Linear constraints: The first 3 of the 14 equations are to

be satisfied exactly.

(2) Unequal weights: The observations are assumed to have

different variances, hence different weights are assigned.

13
(3) Zero weights: The last two observations are given zero weights. Hence the last two rows of A and B are ignored in computing

the solution vector of coefficients, but residuals are computed for

these two observations.

3. STRUCTURE CHART

A structure chart showing subprograms and data elements used by the

driver program and each of the five testing subroutines is given in

Figure 1.

4. HOW TO OMIT PIVOTING

In order to increase the accuracy of computations, pivoting is used

in the decomposition subroutines (DEC0M1 and DEC0M2) which are called by

L2A and L2B. In certain applications, however, the user may wish to

omit pivoting. Pivoting was omitted in the adaptation of L2A prepared

for OMNITAB 80 [11 J so that the user can control the order in which

variables enter the regression equation and thus obtain the proper se- quence of certain sums of squares used in an analysis of variance.

Pivoting can be omitted by making the following changes in DEC0M1

and DEC0M2. (DEC0M1 and DEC0M2 are not listed in the present paper, but

complete listings of these subroutines are given in "Collected Algorithms from ACM" [6].)

DEC0M1 :

(1) Delete statements DC1 630 through DC1 990, inclusive.

(2) Replace DC1 1000 by:

130 IK = NS

(3) After DC1 1060 insert:

IF (NS.EQ.l .AND. DS.EQ.0.0) RETURN

DEC0M2 :

(1) Delete statements DC2 670 through DC2 1020, inclusive.

(2) Replace DC2 1030 by:

140 IK = NS

(3) After statement DC2 1160 insert:

IF (NS.EQ.l .AND. DS.EQ.0.0)RETURN

5. COMPUTING A PSEUDOINVERSE

One can compute the pseudoinverse of an M x N matrix A through the use of subroutine L2B. If one chooses MODE =2, L = M, and sets B equal

to the identity matrix of order M, the N x M pseudoinverse of A will be

returned by L2B in array X. That is, the K-th coefficient vector will

be the K-th column of the desired pseudoinverse.

Problem 9 of elements DATA5A and DATA5C is an example given by

Albert [1 ; p. 63] in which A is the 3 by 4 matrix

14
< PQ u Q w m 1/1 m in u-l < < < < <c H H H H H

< < <C < < \zlKJkjkJkJ

15
1 1 1

1 -1

1 1 1

This matrix has rank 2; its pseudoinverse is

3 3

1 15

-15 4

4 -5 -1

3 3

Appendix C gives computed results for this problem from both L2A

(MODE 1) and L2B (MODE 2). Both algorithms find the matrix A to be of

rank 2 and give the same sum of squared residuals. In L2A, columns 3

and 4 of A are ignored in computing the solutions since only columns 1 and 2 were found to be linearly independent. In L2B, however, all four

columns of A are used in obtaining the solutions. Uniqueness is

achieved by requiring that the solution vectors be of minimal norm.

6. MISCELLANEOUS REMARKS

In using the scaling option in connection with subroutine L2B, one should note that if the computed rank Nl is less than N, the solution

vector of coefficients depends upon the choice of scaling. (Recall that

when Nl < N, a unique solution is not obtained until we impose the con- dition that the solution vector be of minimal norm.).

It was mentioned previously that an adaptation of algorithm L2A was prepared for OMNITAB 80 [11]. Special features of this version which

differ from L2A as described in [22] are:

(1) Pivoting is omitted in the decomposition subroutine.

(2) Sums of squares useful in analysis of variance are computed.

(3) Standard deviations of predicted values are computed.

(4) An estimate of the number of correctly computed digits in

each coefficient is obtained by following the user-specified fit by a fit to predicted values and comparing how well the two sets of

coefficients agree. See Wampler [19].

(5) If the desired fit using raw data fails to converge to a solution, or if the system of equations is found to be singular, the

data are automatically scaled and another attempt is made to obtain a solution.

(6) All two-dimensional arrays were treated as one-dimensional

vectors.

16
(7) The feature of allowing Ml linear constraints was omitted.

As was mentioned in Wampler [22], double precision versions of L2A

and L2B (in which the iterative refinement of the solution is omitted)

have been successfully implemented on a Univac 1108 and an IBM 360/65.

Additional results of test problems run on L2A are given in

Wampler ! s "Test Procedures and Test Problems for Least Squares

Algorithms" [24]. This paper considers two types of test problems: (1)

Data due to Longley [14]; (2) Wampler 's modification of a problem due to

Lauchli [12]. Comparative results are given from running these problems

on five different algorithms:

(1) Cholesky decomposition.

(2) Givens transformations.

(3) Modified Gram-Schmidt orthogonalization.

(4) Householder transformations.

(5) Singular value decomposition.

Wampler [24] also discusses and classifies some of the useful least

squares test problems which have appeared in the literature during the

past twenty years.

17
Table I

Time Required to Solve Hastings-type Problems on the Univac 1108

Values of parameters Time required (seconds)

M N Ml Subroutine DECOMP SOLVE COVAR TOTAL

10 2 L2A .004

10 3 L2A .006

10 4 L2A .007

10 6 L2A .009

10 8 L2A .013

10 8 L2B .014

10 8 4 L2A .010

10 8 4 L2B .011

20 2 L2A .005

20 5 L2A .012

20 10 L2A .028

20 15 L2A .056

20 15 L2B .057

20 15 5 L2A .039

20 15 5 L2B .042

30 5 L2A .014

30 10 L2A .038

30 15 L2A .078

30 20 L2A .131

30 20 L2B .134

30 20 10 L2A .086

30 20 10 L2B .091

40 5 L2A .016

40 10 L2A .048

40 20 L2A .169

40 30 L2A .368

40 30 L2B .372

40 30 10 L2A .246

40 30 10 L2B .263

50 10 L2A .058

50 20 L2A .207

50 30 L2A .451

50 40 L2A .791

50 40 L2B .796

50 40 10 L2A .551

50 40 10 L2B .604

100 10 L2A .108

100 20 L2A .396

100 40 L2A 1.524

.008 .003 .349

.010 .004 .251

.011 .005 .268

.015 .006 .300

.018 .009 .337

.037 .008 .495

.018 .008 .333

.036 .008 .376

.012 .004 .317

.019 .005 .365

.033 .012 .472

.048 .024 .635

.104 .025 .734

.047 .024 .617

.102 .025 .709

.024 .005 .410

.043 .011 .558

.064 .024 .752

.087 .047 1.002

.188 .047 1.160

.087 .047 .942

.189 .048 1.109

.030 .005 .470

.055 .012 .642

.107 .047 1.141

.165 .134 1.895

.371 .134 2.209

.161 .133 1.752

.362 .134 2.068

.065 .011 .723

.128 .047 1.284

.195 .134 2.136

.272 .292 3.293

.646 .308 3.850

.262 .292 2.984

.615 .300 3.590

.120 .011 1.144

.231 .048 2.007

.471 .292 4.950

18
Time Requireid to !

Values of parameters

M N Ml

100 60

100 80

100 80

100 80 10

100 80 10

150 10

150 20

150 40

150 60

150 80

150 80

150 80 10

150 80 10

200 10

200 20

Table 1

200 10

200 20

200 40

200 60

200 80

200 80

200 80 10

200 80 10

Time required (seconds)

Subroutine DECOMP SOLVE COVAR TOTAL

L2A 3.389 .746 .915 9.792

L2A 6.156 1.034 2.083 16.912

L2B 6.122 2.404 2.125 18.152

L2A 4.840 1.001 2.115 15.093

L2B 4.986 2.288 2.112 16.619

L2A .159 .174 .011 1.562

L2A .585 .334 .047 2.730

L2A 2.258 .672 .291 6.634

L2A 5.058 1.232 1.044 13.276

L2A 10.255 2.295 2.298 24.691

L2B 8.976 4.426 2.084 24.563

L2A 7.458 1.448 2.190 20.842

L2B 7.473 3.074 2.086 21.580

L2A .213 .234 .012 2.009

L2A .840 .496 .054 3.665

L2A 3.227 .918 .295 8.789

L2A 7.382 1.625 1.049 17.684

L2A 12.879 1.917 2.102 29.626

L2B 11.901 3.963 2.091 29.278

L2A 11.848 1.823 2.299 28.129

L2B 12.419 4.221 2.157 30.669

19
OQCM

hJ 1 m• O

n vO

—i o n

60 o o n ~* •<!•-* co -h c 1—1

<J NMOOON •H II «—< © O O «—i CM

CO II • ••••• 3 H Jo o o o o o o 1 III A <u H vO >

cu

II .J

25 cu

CO

a •H

O

r-~

O • s-**oocom

SS o • ••••• II CI ^ CO fl O O

cu II >J <T vO vO CM

N 2 > m m o o m m

•H

J st «>r

m

m I cu •U o I I CO CO H

«J iH

CJ CU

•u

CO

CU H

<U

m

.C • 00 00 CO CO CM 00

4-1 o o uo o co co -3- o

>*N00« ON

M t-l • II •—< OOO

—< CM

M o o • ••••• UH ij o o o o o o

CU

a

II oH 1 III

ja c iH

cd o CU H S3iH

CO >

cu

O 3 CU • 00 00 CO CO CM 00

T3 CO o m o co co -<r o

•H •H <*N COCOON

CO

O II ~H

OOO -h CM

CU 53 • ••••• Pd .Jo o o o o o o 1 III T3 H C

crjA

.*c2

fj

cu

4-1 3

O.aooA ^h cm co sf in vo

COco X

XXXXX

•H

4-1 3

iHO

CO

\D

CO• OCOmm

CM

-a- m

co

co• oCOm

coco• oCO

in

CO3 M T3 C •H

CO CO U CUM

T3V M-4 4J O a& B o o o S3

20
• o i—< i—I .—t

60 o

•H CO 1—1 U

UM

CO r-1 II

O

O

O • CO 8

vO x: s O

OO

CM 4J -O 1 00 «H i—i i—i i—i

CN CO C CU i—i

M

hh

II •H I-i II8 O

OO < M Mh O O O H O O w u h oooo H 4-» CO •* «*

CU -H i-H oo • 00 o o o o o w • -H * 4-t 4-» 4J —• 4-1 x *o II -4 «H 1 8 O O vO

00 CN ~H CN

O -H

CO H3 >

•H 00

O c H»CO » >> 1 1 1 iH o o o C CU a

i—1 r-1 —I O 4J XXX cfl o <t sf m e a c ~-t <t r-^ 3 -H • • - • oi X iH O o CO II o o o

COh => 8 4-» 4-» 4-> e o. TJ

a) a. •H r^ -* ~*

iH CO CO t—1 ^ ^H

-O cu 1 1 1 o a pej o o o M O

•—1 i-H i—<

P-. -H XXX CO m >a- o\

CU i-t cti cr> cn a o • • • >. CO •

4-> I-l 00 CO i a. •H 4J

O

T3 CO M CU CU •—i o —• o

CO

C ,c -h II

a -h 4-> CJ 8

CO 42 00 «H T3 cj <U i—i

SI

r

in

coef

II8 o o o —1 • o o M U <4-t • o o o M P o II CU 4J 8

CU • <H

O

iH

X

w II o o o

Xi CO «H 8

CO

a TJ

H

COC

<4-l

O -3-

mm

O iH

4-1 o o o • CO 4J 4-t 4J o n Z CU CM CN St

4J •H

CO CN vO

O ^ HHN CN CN CN CN

21
•

si" ("MCM

II ("O ^

CN• \ £> ( NI1CM

II <HW

00 • o CO —1 u

I— •H

00

CJ •H

CO T3 >

•H 00 c& >. .H a CU o U

CO

a

a 3 •H & XO

co

M e a

CO a.

T-i CO Xio c u o

P-, •H

CO

0) •H

p. o

>•. CU •u u 1 a

T3

>-i <u

CO C3 s •H

CO 42

T3 CJ 53• >H

CU .H

43

CO H •60

•H CO

CO rH

CO 4= d

4J •O

00 •H

CO

a CU •H t-l

Vi 14-1 O O UM 4J

<U •H

60 • •H X -a £

CO3

i-lCO Pi

60 CO •H 4-1 CO c

CU 4= •H

4-> o

00 •H

M-l c <4-l •H CUO M CJ Ou (4-4 u O

CU

4-1 • •H X 60 3

•H

COC

<4-l o o •H

4-t • CO O l-l 55 CU

4-i •H

II8

OD 00 00 1 1 1 CO CON III oo

n

nn i i i 1 1 1 O O O f—I t—1 i—4 1 1 1 o o o I—1 —1 -H i i i i o o o o —1 —I 1—1 1-^ o

ii8

»-•XXXO v£> cm m in XXX CO * H

<t r-» cn

XXXX

i—i <—* as «—<

^O

—( CN •—1 o8

on

CM oo m

.—I

II8

i—i o

II8

^Ho o •

II8 * O o o o

II8OO

m

00

00 00 00

CM CM <t

cm <r <r

-4

O

O CO

en

m

CO o o o o o o o

r» \D

r-. oo

oo oo

-H

O •—<

i—i o •—i

00

O CM

CM CO 00 • • • •

-h CM

O

OO

O o o o o o o o

LTnON

scooo

oo oo on

00 00 ON 00 00 ON ON

cn <r oo cm -* oo vo

00 00 00 vO ^O vO v£>

ON

ON

ONON

ONON

ONON

CU

>-l %

CO

4-1 c

cu

•-*-1

<4-l cuo

CJ

cu 4J&o

22
•60

•H CO

CO iH

CO X! 3

4J T3

v£> •H

CO

C CU •H u

1-1 (4-1 o O uu 4J

01 •H

t>0 • •H X T3 s

60 CO •H 4J

CO

a

CU x: •H

4J

O

VO •H

M-l e •4-1 v-i CUo

VJ o o

1-1 <4-l p O

CU

4-» • •H X 00

fl •H a T3 oo

II8

COC

<4-i

O

O "H

4-1 • CO 2 M z cu

CO : CM

O

O

OO

ir> <t I I o o «—• i—

X X

CO CN o o

4J 4J I I o o i—I r-i X X

vO ^

CM ^ O I-** o o o O •H

i—

1 25 o o CO8at0

i-i «* u~t CU o O 4-)

4-1 4-1 3a

CM CM

O

CM ^> o —1 •—I CM CM CM CM

23
ACKNOWLEDGMENTS

The author wishes to thank Prentice-Hall, Inc., Englewood Cliffs,

New Jersey, for permission to reproduce Function GEN and portions of

Program PR0G2 (pages 311, 281-282), from Charles L. Lawson and Richard

J. Hanson, SOLVING LEAST SQUARES PROBLEMS, copyright 1974.

Thanks are also due to Richard J. Hanson for permission to use the

problems constructed from Hadamard matrices, and to Joseph M. Cameron

and James W. Longley for helpful collaborations.

REFERENCES

1. Albert, Arthur. Regression and the Moore-Penrose Pseudoinverse . Academic Press, New York, 1972.

2. Bennett, Carl A., and Franklin, Norman L. Statistical Analysis in

Chemistry and the Chemical Industry . Wiley, New York, 1954.

3. Bjorck, Ake. Comment on the iterative refinement of least-squares

solutions. Journal of the American Statistical Association 73

(1978), 161-166.

4. Bjorck, Ake, and Golub, Gene. ALGOL Programming Contribution No.

22: Iterative refinement of linear least square solutions by

Householder transformation. BIT ^7 (1967), 322-337.

5. Businger, Peter, and Golub, Gene H. Linear least squares by

Householder transformations. Numerische Mathematik 7 (1965),

269-276.

6. Collected Algorithms from ACM. Association for Computing

Machinery, 1133 Avenue of the Americas, New York, N.Y. 10036.

7. Hall, Marshall, Jr. Combinatorial Theory . Blaisdell, Waltham,

Massachusetts, 1967.

8. Hanson, Richard J. Personal communication.

9. Hastings, W. Keith. Test data for statistical algorithms: Least

squares and ANOVA. Journal of the American Statistical Association

67 (1972), 874-879.

10. Hogben, David, Peavy, Sally T. , and Varner, Ruth N. OMNITAB II

User's Reference Manual. National Bureau of Standards Technical

Note 552, U. S. Government Printing Office, Washington, D.C., 1971.

11. Hogben, David, and Peavy, Sally T. OMNITAB 80: A High-Quality

Integrated Programming Language and Statistical Software Computing

System. NBS Handbook, in preparation.

24
12. Lauchli, Peter. Jordan-Elimination und Ausgleichung nach kleinsten

Quadraten. Numerische Mathematik 3 (1961), 226-240.

13. Lawson, Charles L. , and Hanson, Richard J. Solving Least Squares

Problems . Prentice-Hall, Englewood Cliffs, New Jersey, 1974.

14. Longley, James W. An appraisal of least squares programs for the

electronic computer from the point of view of the user. Journal

of the American Statistical Association 62 (1967), 819-841.

15. Moses, L. E. , and Oakford, R. V. Tables of Random Permutations . Stanford University Press, Stanford, California, 1963.

16. The RAND Corporation. A Million Random Digits with 100,000 Normal

Deviates . The Free Press, Glencoe, Illinois, 1955.

17. Sande, G. Program execution profiles. Proceedings of the Computer

Science and Statistics Eighth Annual Symposium on the Interface,

ed. by James W. Frane, University of California, Los Angeles

(1975), 325-326.

18. Seber, G. A. F. Linear Regression Analysis . Wiley, New York,

1977.

19. Wampler, Roy H. An evaluation of linear least squares computer

programs. Journal of Research of the National Bureau of Standards

73B (1969), 59-90.

20. Wampler, Roy H. A report on the accuracy of some widely used least

squares computer programs. Journal of the American Statistical

Association 65 (1970), 549-565.

21. Wampler, Roy H. Some recent developments in linear least-squares

computations. Proceedings of the Computer Science and Statistics

Sixth Annual Symposium on the Interface, ed. by Michael E. Tarter,

University of California, Berkeley (1972), 94-110.

22. Wampler, Roy H. Solutions to weighted least squares problems by

modified Gram-Schmidt with iterative refinement. ACM Transactions

on Mathematical Software , Vol. 5, No. 4 (Dec. 1979), 457-465.

23. Wampler, Roy H. Algorithm 544: L2A and L2B, Weighted least

squares solutions by modified Gram-Schmidt with iterative

refinement. ACM Transactions on Mathematical Software , Vol. 5,

No. 4 (Dec. 1979), 494-499.

24. Wampler, Roy H. Test procedures and test problems for least

squares algorithms. Journal of Econometrics 12 (1980), 3-22.

25
APPENDIX A. LISTING OF PROGRAM

Main program (DRIVER)

Subroutine TEST1

Subroutine TEST2

Subroutine TEST3

Subroutine TEST4

Subroutine TEST5

Subroutine OUTPUT

Subroutine ORTHO

Function GEN

Subroutine HADMAR

Subroutine SCALE

26
THIS MAIN PROGRAM CALLS FIVE TESTING SUBROUTINES FOR SOLVING LINEAR

LEAST SQUARES TEST PROBLEMS, THE TESTING SUBROUTINES. IN TURN, CALL

SUBROUTINES WHICH USE A MODIFIED GRAM-SCHMIDT ALGORITHM WITH

ITERATIVE REFINEMENT IN ORDER TO OBTAIN THE SOLUTIONS.

VERSION OF FEBRUARY 25. 1980.

WRITTEN BY ROY H. WAMPLER» STATISTICAL ENGINEERING DIVISIONS

NATIONAL BUREAU OF STANDARDS. WASHINGTON. D. C. 20234.

REFERENCES

1. ROY H. WAMPLER. SOLUTIONS TO WEIGHTED LEAST SQUARES PROBLEMS BY

MODIFIED GRAM-SCHMIDT WITH ITERATIVE REFINEMENT.

ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE. VOL. 5. 1979.

PP» 457-465.

2. ROY H. WAMPLER. ALGORITHM 544. L2A AND L2B, WEIGHTED LEAST

SQUARES SOLUTIONS BY MODIFIED GRAM-SCHMIDT WITH ITERATIVE

REFINEMENT (F4). ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE.

VOL. 5e 1979. PP, 494-499.

3. ROY H. WAMPLER. PROBLEMS USED IN TESTING THE EFFICIENCY AND

ACCURACY OF THE MODIFIED GRAM-SCHMIDT LEAST SQUARES ALGORITHM.

N.B.S. TECHNICAL NOTE. NATIONAL BUREAU OF STANDARDS.

WASHINGTON. D. C. 20234. (THIS PUBLICATION)

THIS MAIN PROGRAM READS IN TWO INPUT PARAMETERS — KIND AND MODE.

DEPENDING ON THE VALUE OF KIND (1, 2. 3. 4 OR 5). ONE OF FIVE TESTING

SUBROUTINES IS CALLED TO EITHER GENERATE DATA OR READ IN DATA FOR

LEAST SQUARES PROBLEMS. THESE FIVE SUBROUTINES ARE NAMED TEST1.

TEST2. TEST3. TEST4 AND TEST5. THE FIRST FOUR GENERATE SPECIFIC TYPE

OF TEST PROBLEMS. WHEREAS TEST5 READS IN DATA FOR GENERAL LEAST

SQUARES PROBLEMS.

EACH OF THESE FIVE SUBROUTINES WILL CALL EITHER SUBROUTINE L2A OR L2B

TO COMPUTE THE LEAST SQUARES SOLUTIONS. SUBSEQUENTLY. EACH OF THE

FIVE TESTING SUBROUTINES CALLS SUBROUTINE OUTPUT FOR PRINTING

RESULTS.

THE VALUE OF THE PARAMETER MODE (1 OR 2 ) DETERMINES WHETHER SUBROUTIN

L2A OR L2B SHALL BE CALLED. THESE TWO SUBROUTINES WILL FURNISH THE

SAME SOLUTIONS WHENEVER THE COMPUTED RANK OF THE SYSTEM OF M EQUATION

IN N UNKNOWNS EQUALS N AND M IS GREATER THAN OR EQUAL TO N. IN CASES

WHERE THE COMPUTED RANK Nl IS LESS THAN N. THE USER SPECIFIES THE TYP

OF SOLUTION TO BE COMPUTED ACCORDING TO WHETHER MODE = 1 OR MODE = 2.

MATRIX A IS THE GIVEN MATRIX OF A SYSTEM OF M LINEAR EQUATIONS IN N

UNKNOWNS. MATRIX W IS A GIVEN DIAGONAL MATRIX OF WEIGHTS WITH ALL

DIAGONAL ELEMENTS NONNEGATIVE. LET H = (SQRT(W))*A.

MODE = 1 INDICATES THAT IF Nl.LT.N THE ORIGINAL MATRIX (M BY N)

IS TO BE REPLACED BY A SMALLER MATRIX (M BY Nl ) WHOSE

COLUMNS ARE LINEARLY INDEPENDENT. AND A SOLUTION IS TO BE

SOUGHT FOR THE SMALLER SYSTEM OF EQUATIONS. THUS N - Nl

COLUMNS OF THE ORIGINAL MATRIX H ARE DELETED. AND COEFFICIENTS

CORRESPONDING TO THESE N - Nl DELETED COLUMNS WILL BE SET

EQUAL TO ZERO. (MODE 1 CORRESPONDS TO SUBROUTINE L2A.

)

MODE = 2 INDICATES THAT A SOLUTION IS SOUGHT FOR A LEAST SQUARES

PROBLEM HAVING N ELEMENTS IN THE SOLUTION VECTOR. IN ORDER

DRV 0010

DRV 0020

DRV 0030

DRV 04

DRV 0050

DRV 0060

DRV 07

DRV 0080

DRV 0090

DRV 0100

DRV 0110

DRV 0120

DRV 0130

DRV 0140

DRV 0150

DRV 0160

DRV 0170

DRV 0180

DRV 0190

DRV 0200

DRV 0210

DRV 0220

DRV 0230

DRV 0240

DRV 0250

DRV 0260

DRV 0270

DRV 0280

DRV 0290

SDRV 0300

DRV 0310

DRV 32

DRV 330

DRV 34

DRV 0350

DRV 360

DRV 37

DRV 38

EDRV 0390

DRV 40

SDRV 041

DRV 0420

EDRV 0430

DRV 44

DRV 45

DRV 0460

DRV 0470

DRV 0480

DRV 0490

DRV 0500

DRV 0510

DRV 520

DRV 53

DRV 0540

DRV 0550

DRV 0560

DRV 0570

27
TO OBTAIN A UNIQUE SOLUTION. THE CONDITION THAT THE

SOLUTION VECTOR BE OF MINIMAL EUCLIDEAN NORM IS IMPOSED.

(MODE 2 CORRESPONDS TO SUBROUTINE L2B.

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

THE ARGUMENT LISTS IN THE CALLING SEQUENCES FOR THE FIVE SUBROUTINES

CALLED BY THIS MAIN PROGRAM ARE ALL ALIKE. OF THE 21 PARAMETERS

APPEARING IN THESE ARGUMENT LISTS. ONE IS THE INPUT PARAMETER MODE

ALREADY DESCRIBED. 13 ARE DIMENSIONED ARRAYS. AND 7 ARE DIMENSIONING

PARAMETERS. THE DIMENSIONING PARAMETERS ARE USED TO TRANSMIT

DIMENSIONS OF ARRAYS TO SUBROUTINES. FIVE ADDITIONAL INPUT QUANTITIESDRV

(M. N. Ml. L AND TOD NEEDED IN THE CALLING SEQUENCES FOR L2A AND L2B DRV

WILL BE GIVEN VALUES IN THE TESTING SUBROUTINES. SEE SUBROUTINE TEST5DRV

FOR COMMENTS ON THESE PARAMETERS. DRV

DRV

DRV

DRV

DRV

DRV

DRV

DESCRIPTION OF DIMENSIONED ARRAYS —

A TWO-DIMENSIONAL ARRAY OF SIZE (MM.NNP1). ON ENTRY TO L2A OR

L2B. THIS ARRAY CONTAINS THE GIVEN MATRIX OF A SYSTEM OF M

LINEAR EQUATIONS IN N UNKNOWNS. WHERE THE FIRST Ml EQUATIONS

ARE TO BE SATISFIED EXACTLY.

B TWO-DIMENSIONAL ARRAY OF SIZE (MM.LL). ON ENTRY TO L2A OR L2B. DRV

B CONTAINS THE L GIVEN RIGHT-HAND SIDES (VECTORS OF

OBSERVATIONS).

W ONE-DIMENSIONAL ARRAY OF LENGTH (MM). ON ENTRY TO L2A OR L2B.

W CONTAINS THE DIAGONAL ELEMENTS OF A GIVEN DIAGONAL MATRIX OF DRV

WEIGHTS. ALL NONNEGATIVE. (THE FIRST Ml ELEMENTS OF W ARE SET DRV

EQUAL TO 1.0 BY THE PROGRAM WHEN Ml IS GREATER THAN ZERO.) ON DRV

EXIT FROM L2A OR L2B. THE ORIGINAL ELEMENTS OF W HAVE BEEN

REPLACED BY THEIR SQUARE ROOTS.

X TWO-DIMENSIONAL ARRAY OF SIZE (NN.LL).

L2B. X CONTAINS THE SOLUTION VECTORS.

RES TWO-DIMENSIONAL ARRAY OF SIZE (MM.LL).

L2B. RES CONTAINS THE RESIDUAL VECTORS.

R TWO-DIMENSIONAL ARRAY OF SIZE (NN.NN).

CONTAINS THE LOWER TRIANGULAR PORTION OF THE SYMMETRIC

UNSCALED COVARIANCE MATRIX.

Q TWO-DIMENSIONAL ARRAY OF SIZE (MM.NN). USED INTERNALLY ONLY

IN SUBROUTINE L2A AND RELATED SUBROUTINES.

QR TWO-DIMENSIONAL ARRAY OF SIZE (MMPNN.NN). ON EXIT FROM L2B.

QR CONTAINS THE LOWER TRIANGULAR PORTION OF THE SYMMETRIC

UNSCALED COVARIANCE MATRIX.

IPIVOT ONE-DIMENSIONAL ARRAY OF LENGTH (NN). ON EXIT FROM L2A OR

L2B. THIS ARRAY RECORDS THE ORDER IN WHICH THE COLUMNS OF H

WERE SELECTED BY THE PIVOTING SCHEME IN THE COURSE OF THE

ORTHOGONAL DECOMPOSITION.

ON EXIT FROM L2A OR

ON EXIT FROM L2A OR

ON EXIT FROM L2A. R

FAIL

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

ONE-DIMENSIONAL ARRAY OF LENGTH (MNL). IF L2A IS TO BE CALLED. DRV

MNL MUST EQUAL AT LEAST 4*(M+N)+2*L. IF L2B IS TO BE CALLED. DRV

MNL MUST EQUAL AT LEAST 6*(M+N)+2*L. THIS VECTOR IS USED (1) DRV

FOR INTERNAL WORK SPACE AND (2) FOR RETURNING INFORMATION ON DRV

THE BEHAVIOR OF THE ITERATIVE REFINEMENT PROCEDURE. DRV

(A) NUMIT IS THE NUMBER OF ITERATIONS CARRIED OUT DURING THE DRV

ITERATIVE REFINEMENT IN ATTEMPTING TO OBTAIN A SOLUTION DRV

FOR THE K-TH RIGHT-HAND SIDE. ON EXIT FROM L2A OR L2B. DRV

C(K) = +NUMIT IF THE SOLUTION CONVERGED. AND DRV

C(K) = -NUMIT IF THE SOLUTION FAILED TO CONVERGE. DRV

(B) DIGITX GIVES AN ESTIMATE OF THE NUMBER OF CORRECT DIGITS DRV

IN THE INITIAL SOLUTION OF THE COEFFICIENTS FOR THE K-TH DRV

RIGHT-HAND SIDE. ON EXIT FROM L2A OR L2B. C(K+L) = DIGITX.DRV

ONE-DIMENSIONAL ARRAY OF LENGTH (LL). USED IN SUBROUTINE DRV

580

0590

60

0610

0620

0630

0640

0650

66

0670

0680

0690

0700

0710

0720

0730

0740

0750

0760

0770

0780

0790

0800

0810

0820

0830

0840

0850

0860

0870

0880

0890

0900

0910

0920

0930

0940

0950

0960

0970

0980

0990

1000

1010

1020

1030

1040

1050

1060

1070

1080

1090

1100

1110

1120

1130

1140

1150

28
OUTPUT IN CONNECTION WITH INFORMATION ON THE SUCCESS OR

FAILURE OF THE ITERATIVE REFINEMENT PROCEDURE.

SDX ONE-DIMENSIONAL ARRAY OF LENGTH (NN). USED FOR STANDARD

DEVIATIONS OF COEFFICIENTS IN SUBROUTINE OUTPUT. ' SF ONE-DIMENSIONAL ARRAY OF LENGTH (NNPLL). USED FOR SCALE

FACTORS IN SUBROUTINES TEST4. TESTS AND SCALE.

DIMENSIONING PARAMETERS —

MM MUST SATISFY MM.GE.M.

NN MUST SATISFY NN.GE.N.

LL MUST SATISFY LL.GE.L.

MMPNN MUST SATISFY MMPNN. GE.M+N.

NNP1 MUST SATISFY NNPI.GE.N+1.

NNPLL MUST SATISFY NNPLL. GE .N+L.

MNL MUST SATISFY MNL.GE.K*( M+N ) +2*Lt WHERE K EQUALS 4 IF

SUBROUTINE L2A IS TO BE CALLED* AND K EQUALS 6 IF SUBROUTINE L2B IS TO BE CALLED.

DIMENSIONED ARRAYS ARRANGED BY TYPE —

INTEGER IPIVOT(NN)

REAL A(MMtNNPl) t B(MM.LL). C(MNL)

REAL Q(MM.NN), R(NN.NN). QR(MMPNN.NN)

REAL RES(MM.LL). SDX(NN), SF(NNPLL)t W(MM). X(NN.LL)

LOGICAL FAIL(LL)

INTEGER IPIV0T(20)

REAL AO0.21). BOO. 4), C(308»

REAL QO0.20). R(20.20). QR(50.20)

REAL RESO0.4), SDX<20)« SF(24), WOO). X(20.4)

LOGICAL FAIL<4)

IN THE FOLLOWING DATA STATEMENT. NR IS THE CARD READER DEVICE

AND NW IS THE PRINTER DEVICE NUMBER.

DATA NR.NW / 5. 6 /

DATA MM / 30 /

DATA NN / 20 /

DATA LL / 4 /

DATA MMPNN / SO /

DATA NNP1 / 21 /

DATA NNPLL / 24 /

DATA MNL /308 /

READ (NR. 80000) KIND. MODE

IF (KIND.LT.l .OR. KIND.GT.5) WRITE (NW. 90000)

IF (KIND.LT.l .OR. KIND.GT.5) STOP

IF (M0DE.NE.2) MODE = 1 WRITE (NW. 90010) KIND. KIND .MODE

GO TO (10.20.30.40.50). KIND

TEST1 — PROBLEMS DUE TO HASTINGS.

10 CALL TEST1 (A. B. C. FAIL. IPIVOT. LL» MM. MMPNN. MNL. MODE.

* NN. NNPLL. NNP1, Q. QR . R. RES. SDX. SF. W* X)

GO TO 60

TEST2 — PROBLEMS DUE TO LAWSON AND HANSON.

DRV 1160

DRV 1170

DRV 1180

DRV 1190

DRV 1200

DRV 1210

DRV 122

DRV 1230

DRV 1240

DRV 1250

DRV 1260

DRV 1270

DRV 1280

DRV 1290

DRV 1300

DRV 1310

DRV 132

DRV 1330

DRV 1340

DRV 1350

DRV 1360

DRV 1370

DRV 1380

DRV 1390

DRV 1400

DRV 1410

DRV 1420

DRV 1430

DRV 1440

DRV 1450

DRV 1460

DRV 1470

DRV 1480

DRV 1490

DRV 1500

DRV 1510

DRV 1520

DRV 1530

DRV 1540

DRV 1550

DRV 1560

DRV 1570

DRV 1580

DRV 1590

DRV 1600

DRV 1610

DRV 162

DRV 163

DRV 1640

DRV 1650

DRV 1660

DRV 1670

DRV 168

DRV 1690

DRV 1700

DRV 1710

DRV 1720

DRV 173

29
20 CALL TEST2 (Ac B. C» FAIL. IPIVOTe LL * MM. MMPNN. MNL. MODE*

* NN, NNPLL, NNPl, Q, QR , R, RES* SDX, SF, W. X)

GO TO 60

C

C TEST3 — PROBLEMS DUE TO HANSON.

C

30 CALL TEST3 (A. B, C. FAIL* IPIVOT. LL. MM* MMPNNt MNL. MODEt

* NN. NNPLL* NNPl, Q, QR, R, RES, SDX, SF, W, X)

GO TO 60

C

C TEST4 PROBLEMS DUE TO BJORCK AND WAMPLER.

C

40 CALL TEST4 (Ac Bt Ct FAIL* IPIVOTe LL. MM. MMPNN. MNL. MODE.

* NN. NNPLL. NNPl. Q. QR. R. RES. SDX. SF. W. X)

GO TO 60

C

C TEST5 - GENERAL PROBLEMS.

C

50 CALL TESTS (A, B. C. FAIL. IPIVOTe LL . MM. MMPNN. MNL. MODE.

* NN. NNPLL. NNPl, Q, QR , R, RES, SDX, SF, W, X)

C

60 STOP

C

80000 FORMAT (215)

C

90000 FORMAT (60H0*** ERROR — KIND MUST EQUAL 1. 2. 3e 4 OR 5.

*IND =. 15)

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

DRV

HERE KDRV

DRV

90010 FORMAT (1H1.6HKIND =.I2/10H CALL TEST.I1.34H TO OBTAIN LEAST SQUARDRV

*ES SOLUTIONS/7H MODE =,12) DRV

C DRV

END DRV

1740

1750

1760

1770

1780

1790

1800

1810

182

1830

184

1850

1860

1870

1880

1890

1900

1910

1920

1930

1940

1950

196

197

1980

1990

2000

2010

2020

2030

2040

2050

SUBROUTINE TEST1 (A. B. C. FA I * NN. NNPLL. NNPl, Q, QR , R, RE

SUBROUTINE TEST1 CALLS SUBROUTINE

(MATRICES A AND B) FOR TESTING LEA

IN ORTHO IS BASED ON EXAMPLE 2 IN TEST DATA FOR STATISTICAL ALGOR

BY W. K. HASTINGS, JOURNAL OF T ASSOCIATION* VOL. 67, 1972. PP.

HASTINGS DID NOT GIVE FORMULAS FOR

PROBLEMS HAVING LINEAR CONSTRAINTS

ARE INCLUDED IN ORDER TO OBTAIN DA

THIS SUBROUTINE CALLS A MACHINE-DE

WHICH PRINTS THE NUMBER OF SECONDS

RUN BEGAN. TIME IS PRINTED IN SEC

OF A SECOND.

THE SEQUENCE OF INPUT CARDS (OR CA

1. CARD IN (15) FORMAT GIVING THE

PROBLEMS TO BE SOLVED.

2. FOR EACH OF THE NPROB PROBLEMS.

VALUES OF M, N AND Ml.

L* IPIVOT, LL. MM. MMPMN. MNL. MODE

S* SDX, SF, W, X)

ORTHO TO GENERATE INPUT DATA

ST SQUARES PROGRAMS. THE ALGORITHM

ITHMS LEAST SQUARES AND ANOVA,

HE AMERICAN STATISTICAL

874-879.

THE COEFFICIENTS IN THE CASE OF

• BUT SOME PROBLEMS OF THIS TYPE

TA ON TIMING.

PENDENT EXTERNAL ROUTINE (CLOCKT)

ELAPSED SINCE EXECUTION OF THE

ONDS. TO THE CLOSEST THOUSANDTH

RD IMAGES) FOR THIS SUBROUTINE IS - VALUE OF NPROB, THE NUMBER OF

A CARD IN (315) FORMAT GIVING THE

»TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

TS1

-TS1

TS1

TS1

TS1

TS1

0010

0020

0030

0040

0050

0060

0070

0080

0090

0100

0110

0120

0130

0140

0150

0160

0170

0180

0190

0200

0210

22

0230

0240

30
c

cc

c

LOGICAL FAIL(LL)

IN THE FOLLOWING DATA

AND NW IS THE PRINTER

C M TOTAL NUMBER OF EQUATIONS. M MUST BE EVEN. TS1 0250

C N NUMBER OF UNKNOWN COEFFICIENTS. TS1 0260

C MI NUMBER OF LINEAR CONSTRAINTS ( 0.LE.M1 .LE «M AND Ml.LE.N). TS1 0270

C TS1 0280

INTEGER IPIVOT(NN) TSl 0290

REAL A(MM.NNPl). B(MM.LL)t C(MNL) TSI 0300

REAL Q(MMtNN)» R(NN«NN)« QR(MMPNN«NN) TSl 0310

REAL RES(MM«LL>t SDX(NN). SFC NNPLL) • W(MM) t X(NN.LL) TSl 0320

REAL AC. TOL TSl 0330

TSl 03*0

TSl 0350

STATEMENTt NR IS THE CARD READER DEVICE TSl 0360

DEVICE NUMBER. TSl 0370

TSl 0380

DATA NRtNW /5t6/ TSl 0390

C TSl 0400

CALL CLOCKT(AC) TSl 0410

C TSl 0420

READ (NR, 80000) NPROB TSl 0430

IPROB = TSl 0440

10 IPROB = IPROB + 1 TSl 0450

WRITE (NW. 90000) IPROB TSl 0460

C TSl 0470

CALL CLOCKT(AC) TSl 0480

C TSl 0490

READ (NR»80010) M.N. MX TSl 0500

L = 1 TSl 0510

TOL = 0.0 TSl 0520

IWGHT = 1 TSl 0530

C TSl 0540

C IWGHT = 1 DENOTES THAT ALL WEIGHTS ARE EQUAL TO 1.0. TSl 0550

C TSl 0560

WRITE <NW.90010) TSl 0570

WRITE (NW.90020) M» N» Ml » L» IWGHT. MODE* TOL TSl 0580

C TSl 0590

C SET PARAMETERS WHICH ALLOCATE VECTOR C (TO BE USED AS WORK AREA IN TSl 0600

C SUBROUTINE ORTHO). TSl 0610

C TSl 0620

TSl 0630

TSl 0640

TSl 0650

TSl 0660

TSl 0670

TSl 0680

C TSl 0690

C C(K6) IS A VECTOR OF LENGTH N. TSl 0700

C TSl 0710

CALL CLOCKTCAO TSl 0720

C TSl 0730

IF (MODE.EQ.l) CALL ORTHO ( M. N» MM, MM t Qc Bt C(K1). C(K2). TSl 0740

* C(K3), C(K4). C(K5). C(K6>. A) TSl 0750

C TSl 0760

IF (M0DE.EQ.2) CALL ORTHO (M. N. MM. MMPNN* QR» B. C(K1), C(K2). TSl 0770

* C(K3)t C(K4). C(K5). C(K6). A) TSl 0780

C TSl 0790

CALL CLOCKT(AC) TSl 0800

C TSl 0810

NNN = M/2 TSl 0820

31

Kl = 1 K2 = Kl + N + 1 K3 = K2 N

K4 = K3 N

K5 = K4 + N

K6 = K5 N

) IS A VECTOR OF LENGTH
C M MUST BE EVEN.

C

IF (2+NNN .NE. M) GO TO 90

GO TO (20.50) . MODE

20 DO 40 I = 1«M

DO 30 J=1.N

A( It J) = Q( It J)

30 CONTINUE

W(I ) = 1.0

40 CONTINUE

C

CALL CLOCKT(AC)

C

C MATRIX A© MATRIX B AND VECTOR OF WEIGHTS COULD BE PRINTED HERE.

C

CALL L2A (Mt Nt Mlt Lt At Bt Wt TOLt MMt NN.

* Nit IPIVOTt Xt RESt Rt Qt Ct IFAULT)

C

CALL CLOCKT(AC)

C

GO TO 80

50 DO 70 I=ltM

DO 60 J=ltN

A(ItJ) = QR(I.J)

60 CONTINUE

W(I ) = 1.0

70 CONTINUE

C

CALL CLOCKT(AC)

C

C MATRIX At MATRIX B AND VECTOR OF WEIGHTS COULD BE PRINTED HERE.

C

CALL L2B (M, N. Ml. L. A. B. W. TOL. MM. NN* MMPNN.

* Nl, IPIVOTt Xt RESt QRt Ct IFAULT)

C

CALL CLOCKT(AC)

C

80 CALL OUTPUT ( Bt Ct IFAULTt IPIVOTt Lt M. MODEt Ml g Nt Nit QRt

* RESt Wt Xt LLt MM. MMPNN. MNL • NN. FAIL. SDX)

C

90 CALL CLOCKT(AC)

C

IF ( IPROB.LT.NPROB) GO TO 10

RETURN

C

80000 FORMAT (15)

80010 FORMAT (315)

C

90000 FORMAT ( 1H1 » 1 1 5( 1 H* ) . 4X. 7HPR0BLEM. I 4)

90010 FORMAT ( 1 HO . 3X. 1 HM» 4X e 1 HN 9 3X, 2HM1 9 4X, 1 HL. 5X » SHI WGHT » 6X . 4HMOD

* 3HTOL)

90020 FORMAT ( 4I5t 2 1 10. Gl 5.8)

C

END

SUBROUTINE TEST2 (At Bt Ct FA ILt IPIVOTt LLt MMt MMPMNt MNLt

TS1 0830

TS1 0840

TS1 0850

TS1 0860

TS1 0870

TS1 0880

TS1 089

TS1 0900

TS1 0910

TS1 0920

TS1 0930

TS1 0940

TS1 0950

TS1 0960

TS1 0970

TS1 0980

TS1 0990

TS1 1000

TS1 1010

TS1 1020

TS1 1030

TS1 1040

TS1 1050

TS1 1060

TS1 1070

TS1 1080

TS1 1090

TS1 1100

TS1 1110

TS1 1120

TS1 1130

TS1 1140

TS1 1150

TSl 1160

TS1 1170

TSl 1180

TSl 1190

TSl 1200

Rt TSl 1210

TSl 122

TSl 1230

TSl 1240

TSl 1250

TSl 1260

TSl 1270

TSl 1280

TSl 1290

TSl 130

TSl 1310

TSl 1320

• 7X. TSl 1330

TSl 1340

TSl 1350

TSl 1360

TSl 1370

MODE • TS2 0010

32
* NN. NNPLLi NNP1. Q, QR * R, RES* SDX, SF. W, X) TS2 0020

C TS2 0030

C PORTIONS OF THIS SUBROUTINE ARE TAKEN FROM PROGRAM PR0G2. TS2 0040

C CHARLES L. LAWSON AND RICHARD J. HANSON. TS2 0050

C SOLVING LEAST SQUARES PROBLEMS. COPYRIGHT 1974. PAGES 281-282, TS2 0060

C REPRINTED BY PERMISSION OF PRENTICE-HALL. INC., ENGLEWOOD CLIFFS, N.J.TS2 0070

C TS2 0080

C SUBROUTINE TEST2 CALLS FUNCTION GEN TO GENERATE LEAST SQUARES DATA TS2 0090

C MATRICES A AND B. TS2 0100

C TS2 0110

C FOUR SETS OF PROBLEMS ARE GENERATED WITH 18 PROBLEMS IN EACH SET, TS2 0120

C VALUES OF M AND N ARE SELECTED FROM THE SET (1, 2, 3, 6, 7, 8). TS2 0130

C THE TOLERANCE PARAMETER. TOL. IS SET EITHER TO 0,0 OR 0.5. TS2 0140

C TS2 0150

C SET 1. DATA MATRICES A AND B ARE INTEGERS BETWEEN -500 AND +500. TS2 0160

C SOME OF THE A-MATRICES ARE RANK-DEFICIENT. TS2 0170

C TOL = 0.0. TS2 0180

C SET 2. DATA OF SET 1 MODIFIED BY ADDING NOISE TO THE INTEGERS. TS2 0190

C TOL = 0.0. TS2 0200

C SET 3. DATA AS IN SET 1. TS2 0210

C TOL = 0.5. TS2 0220

C SET 4. DATA AS IN SET 2. TS2 0230

C TOL = 0.5. TS2 0240

C TS2 0250

C THIS SUBROUTINE REQUIRES NO INPUT FROM THE USER. TS2 0260

C TS2 0270

INTEGER IPIVOT(NN) TS2 0280

REAL A(MM.NNPl). B(MM.LL). C(MNL) TS2 0290

REAL Q(MM.NN). R(NN.NN). QR(MMPNN.NN) TS2 0300

REAL RES(MM.LL). SDX(NN), SF(NNPLL), W ( MM ) , X(NN.LL) TS2 0310

REAL ANOISE, ANORM, DUMMY, GEN, TAU, TOL TS2 0320

LOGICAL FAIL(LL) TS2 0330

C TS2 0340

C IN THE FOLLOWING DATA STATEMENT. NW IS THE PRINTER DEVICE NUMBER. TS2 0350

C TS2 0360

DATA NW / 6 / TS2 0370

C TS2 0380

IPROB = TS2 0390

DO 90 IT0L=1,2 TS2 0400

DO 80 N0ISE=1.2 TS2 0410

ANORM = 500.0 TS2 0420

ANOISE = 0.0 TS2 0430

TAU = 0.5 TS2 0440

IF (NOISE.EQ.l) GO TO 10 TS2 0450

ANOISE = l.E-4 TS2 0460

TAU = ANORM * ANOISE * 10.0 TS2 0470

10 CONTINUE TS2 0480

C TS2 0490

C INITIALIZE THE DATA GENERATION FUNCTION. TS2 0500

C TS2 0510

DUMMY = GEN(-l.O) TS2 0520

WRITE (NW, 90000) ANOISE. ANORM. TAU TS2 0530

DO 70 MN1=1.6»5 TS2 0540

MN2 = MN1 + 2 TS2 0550

DO 60 M=MN1.MN2 TS2 0560

DO 50 N=MNl,MN2 TS2 0570

C TS2 0580

C GENERATE DATA FOR MATRICES A AND B AND SET WEIGHTS EQUAL TO UNITY. TS2 0590

33
20

30

DO 30 1=1.

DO 20 J=l.N

A(ItJ) = GEN(ANOISE)

CONTINUE

B(I.l) = GEN(ANOISE)

W(I» = 1.0

CONTINUE

Ml =

L = 1 IWGHT = 1 TOL = 0,0

IF (ITOL.EQ.2) TOL = TAU

IPROB = IPROB 1 (NW. 90010) IPROB

(NW. 900 20)

(NW. 90030)

WRITE

WRITE

WRITE M»N.M1bL.I WGMT . MODE* TOL

Cc

c PRINT At B AND W. WRITE

DO 40

40

(NW. 90040)

1=1.

WRITE (NW, 90050)

CONTINUE

( A(I,J),J=1,N).(B( I*K)*K=1*L)*W( I)

EQ.l ) CALL

IPIVOT. X t

IF (M0DE.EQ.2) CALL

NNt MMPNN* Nl* IPI\

L2A (M, N, Ml*

RES. R» Q* c*

L2B (M. N. Ml.

OT. Xt RESi> QR

L* Aj B»

I FAULT

W e TOL«

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

MM* TS2

TS2

TS2

MM* TS2

TS2

TS2

Ml* N* TS2

FAIL. SDX)TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

TS2

RELATIVE NOISE LEVEL OF THE GENERATED DATA WILLTS2

* BE.G15.8/33H0THE MATRIX NORM IS APPROXIMATELY* Gl 5. 8/43H0THE ABS0LTS2

*UTE PSEUDORANK TOLERANCE* TAU* IS.G15.8) TS2

90010 FORMAT ( 1H1 , 1 15 ( 1 H*e 4X .7HPR0BLEM . 14 I TS2

90020 FORMAT < I HO* 3X* 1HM. 4X* 1 HN* 3X. 2HM1 * 4X. 1HL* 5X* 5HI WGHT* 6X*4HM0DE* 7X» TS2

* 3HTOL) TS2

90030 FORMAT ( 415.2110. G15.8) TS2

FORMAT (41H0MATRIX A* MATRIX B AND VECTOR OF WEIGHTS/) TS2

FORMAT (1X.8G15.8) TS2

TS2

END TS2

IF (MODE.

NN* Nl*

L*

» C< A* B* W*

IFAULT)

TOL(

CALL OUTPUT (Be

Nl* QR* R. RES*

C*

W*

IFAULT* IPIVOT. L. M«

X* LL* MM* MMPNN. MNL

MODE*

NN.

50 CONTINUE

60 CONTINUE

70 CONTINUE

80 CONTINUE

90 CONTINUE

RETURN

C

90000 FORMAT (1H0.54HT

90040

90050

C

0600

0610

0620

0630

0640

0650

0660

670

0680

0690

0700

0710

0720

0730

0740

0750

0760

0770

0780

0790

0800

0810

0820

0830

0840

0850

0860

0870

0880

0890

0900

0910

0920

930

0940

0950

0960

0970

0980

0990

1000

1010

1020

1030

1040

1050

1060

1070

1080

1090

1100

1110

SUBROUTINE TEST3 (A* B* C* FAIL* IPIVOT* LL * MM* MMPNN* MNL* M0DE.TS3 0010

* NN* NNPLL* NNP1. Q, QR • R. RES. SDX. SF. W, X) TS3 0020

C TS3 0030

C SUBROUTINE TEST3 GENERATES LEAST SQUARES TEST PROBLEMS CONSTRUCTED TS3 0040

34
FROM HADAMARD MATRICES. THE METHOD OF CONSTRUCTING THESE PROBLEMS IS TS3 0050

DUE TO RICHARD J. HANSON. THIS SUBROUTINE CALLS SUBROUTINE HADMAR TS3 0060

WHICH FURNISHES THE HADAMARD MATRICES. TS3 0070

TS3 0080

THE SEQUENCE OF INPUT CARDS (OR CARD IMAGES) FOR THIS SUBROUTINE IS TS3 0090

1. CARD IN (15) FORMAT GIVING THE VALUE OF NK3« A PARAMETER WHICH TS3 0100

SPECIFIES HOW MANY SETS OF PROBLEMS ARE TO BE CONSTRUCTED AND TS3 0110

SOLVED. EACH SET WILL CONSIST OF 10 PROBLEMS USING THE FOLLOWING TS3 0120

COMBINATIONS OF M AND N. TS3 0130

(M, N) = (2,2). (4,2), ( 4, 4 ) • (8,2). (8*4), (8,8), (16.2), (16,4), TS3 0140

(16,8), (16,16). TS3 0150

2. FOR EACH OF THE NK3 SETS OF PROBLEMS, A CARD IN (15) FORMAT GIVING TS3 0160

A VALUE OF K3. K3, A NONNEGATIVE INTEGER, IS TO BE USED AS A POWERTS3 0170

OF 2 WHICH WILL DETERMINE THE CONDITION NUMBER OF MATRIX A. TS3 0180

(CONDITION NUMBER OF A EQUALS 2**K3.) TS3 0190

3. FOR EACH OF THE NK3 SETS OF PROBLEMS, 6 CARDS IN (1413) FORMAT TS3 0200

GIVING VALUES OF IRHO(I), I=1,...,MMN, WHERE MMN = M - N (FOR TS3 0210

PROBLEMS HAVING M.GT.N). THESE ARE TO BE FURNISHED IN THE TS3 0220

FOLLOWING SEQUENCE TS3 0230

CARD 1 M = 4 N = 2 MMN = 2 TS3 0240

CARD 2 M = 8 N = 2 MMN = 6 TS3 0250

CARD 3 M = 8 N = 4 MMN = 4 TS3 0260

CARD 4 M = 16 N = 2 MMN = 14 TS3 0270

CARD 5 M = 16 N = 4 MMN = 12 TS3 0280

CARD 6 M = 16 N = 8 MMN = 8 TS3 0290

THE ELEMENTS OF IRHO ARE TO BE INTEGERS WHICH ARE UNIFORMLY TS3 0300

DISTRIBUTED ON THE CLOSED INTERVAL (0, BNORM), WHERE TS3 0310

BNORM = N*SQRT(M). TS3 0320

TS3 0330

INTEGER IPIVOT(NN) TS3 0340

INTEGER IRHO(16) TS3 0350

REAL A(MMeNNPl), B(MM,LL)» C(MNL) TS3 0360

REAL Q(MM.NN), R(NN.NN), QR(MMPNN.NN) TS3 0370

REAL RES(MMdLL), SDX(NN), SF(NNPLL), W(MM), X(NN,LL) TS3 0380

REAL BHAT(16)» RH0(16)» S(16*16), TOL« U(16,16), V(16,16) TS3 0390

REAL Z(16, 16) TS3 0400

LOGICAL FAIL(LL) TS3 0410

TS3 0420

IN THE FOLLOWING DATA STATEMENT, NR IS THE CARD READER DEVICE TS3 0430

AND NW IS THE PRINTER DEVICE NUMBER. TS3 0440

TS3 0450

DATA NR, NW / 5, 6 / TS3 0460

TS3 0470

COMPUTE AND PRINT HADAMARD MATRICES OF ORDER 2, 4, 8 AND 16. TS3 0480

TS3 0490

N = 2 TS3 0500

WRITE (NW, 90000) TS3 0510

TS3 0520

10 CALL HADMAR (N, V) TS3 0530

TS3 0540

WRITE (NW, 90010) N TS3 0550

DO 20 1=1, N TS3 0560

WRITE (NW, 90020) ( V ( I , J ) • J= 1 c N ) TS3 0570

20 CONTINUE TS3 0580

N = 2*N TS3 0590

IF (N.LE.16) GO TO 10 TS3 0600

TS3 0610

READ (NR08OOOO) NK3 TS3 0620

35
WRITE (NW. 90030) NK3

IPROB = KA = 30 KA = KA + 1 READ (NR, 80000) K3

KK = 2**K3

K3 CONTROLS VALUE OF S(N.N).

S(N.N) MAY BE SET EQUAL TO THE FOLLOWING POWERS OF 2 —

2**12 = 4096

2**16 = 65536

2**20 = 1048576

2**24 = 16777216

2**25 = 33554432

M = 2

40 N = 2

50 DO 70 I = 1.M

W(I) = 1.0

DO 60 J = 1«N

S( It J) =0.0

60 CONTINUE

70 CONTINUE

S(l, 1 ) = 1.0

S(N«N) = FLOAT(KK)

IF (N.EQ.2) GO TO 90

IF N.GT.2 SET S(Itl) EQUAL TO NUMBERS BETWEEN 1.0 AND S(N.N)* FOR

l.LT.I.LT.N. THE S(I.I) ARE EITHER EQUALLY SPACED OR ALMOST EQUALLY

SPACED. HERE THE SCI. I) ARE SET ARBITRARILY. THESE VALUES COULD BE

CHOSEN RANDOMLY.

NM1 = N - 1 IFACT = (KK-1 l/CN-1

DO 80 1=2. NM1

Sdol) = 1.0 + FLOAT< IFACT* < 1-1 >

80 CONTINUE

CALL SUBROUTINE HADMAR TO OBTAIN V (N BY N).

90 CALL HADMAR ( N. V)

MULTIPLY MATRIX S (M BY N) BY V (N BY N) TO OBTAIN MATRIX Z (M BY N).

DO 120 I=1.M

DO 110 J=1.N

Z< I.J) = 0.0

DO 100 K=1.N

Zfl.J) = Z(I.J) + S( I,K)*V(K.J)

100 CONTINUE

110 CONTINUE

120 CONTINUE

CALL SUBROUTINE HADMAR TO OBTAIN U (M BY M).

CALL HADMAR (M. U)

MULTIPLY MATRIX U (M BY M) BY Z (M BY N) TO OBTAIN MATRIX A (M BY N),

36

TS3 0630

TS3 0640

TS3 0650

TS3 0660

TS3 0670

TS3 0680

TS3 0690

TS3 0700

TS3 0710

TS3 0720

TS3 0730

TS3 0740

TS3 750

TS3 0760

TS3 0770

TS3 0780

TS3 0790

TS3 0800

TS3 0810

TS3 0820

TS3 0830

TS3 0840

TS3 850

TS3 0860

TS3 0870

TS3 0880

TS3 0890

TS3 0900

TS3 0910

TS3 0920

TS3 0930

TS3 0940

TS3 0950

TS3 0960

TS3 0970

TS3 0980

TS3 0990

TS3 1000

TS3 1010

TS3 1020

TS3 1030

TS3 1040

TS3 1050

TS3 1060

TS3 1070

TS3 1080

TS3 1090

TS3 1100

TS3 1110

TS3 1120

TS3 1130

TS3 1140

TS3 1150

TS3 1160

TS3 1170

TS3 1180

TS3 1190

TS3 1200
00 150 1=1.

DO 140 J = 1.N

A(ItJ) = 0.0

DO 130 K=1.M

Ad, J) = A(ItJ) + U< I.K)*Z(K. J)

130 CONTINUE

140 CONTINUE

150 CONTINUE

DO 160 I=1.M

BHAT(I) = FLOAT(N)

RHO(I) = 0.0

160 CONTINUE

READ IRHO EXCEPT WHEN N = M.

THE DATA USED HERE AS IRHO WERE TAKEN FROM —

(1) A MILLION RANDOM DIGITS BY THE RAND CORPORATION. THE FREE

GLENCOE. ILLINOIS. 1955.

(2) TABLES OF RANDOM PERMUTAT I ONS * L. E. MOSES AND R. V. OAKFO

STANFORD UNIVERSITY PRESS. STANFORD. CALIFORNIA. 1963.

AS AN ALTERNATIVE TO READING IN DATA FOR IRHO. ONE COULD CALL

SUBROUTINE TO GENERATE RANDOM INTEGERS.

IF (N.EQ.M) GO TO 180

MMN = M - N

READ (NR. 80010) ( IRHO( I ) * I =1 . MMN)

DO 170 I=1«MMN

IPN = I + N

RHO(IPN) = FLOAT( IRHO( I)

)

170 CONTINUE

MULTIPLY MATRIX U (M BY M) BY RHO ( M BY 1) TO OBTAIN R ( M BY 1

180 DO 210 1=1.

DO 200 J=l.

R( It J) = 0.0

DO 190 K=1.M

R(I« J) = R( I, J)

190 CONTINUE

200 CONTINUE

210 CONTINUE

FOR ALPHA = 0.0» SET B =

U( I.K

)

*RHO(K)

BHAT.

DO 220 1=1.

B(I.l) = BHAT(I)

220 CONTINUE

IF (N.EQ.M) GO TO 240

COMPUTE B(1.2). B(I.3)» B(I»4) FOR ALPHA = 0.01. 1.0* 100.0

M.GT.N.

DO 230 I=1.M

B(I«2) = BHAT(I) + 0.01 * R(Itl)

B(I,3) = BHAT(I) 1.00 * R(I.l)

B(Ie4) = BHAT(I) + 100. * R(I.1>

TS3 1210

TS3 1220

TS3 1230

TS3 1240

TS3 1250

TS3 1260

TS3 1270

TS3 1280

TS3 1290

TS3 130

TS3 1310

TS3 1320

TS3 1330

TS3 1340

TS3 1350

TS3 1360

PRESS. TS3 1370

TS3 1380

RD« TS3 1390

TS3 1400

A TS3 1410

TS3 1420

TS3 1430

TS3 1440

TS3 1450

TS3 1460

TS3 1470

TS3 1480

TS3 1490

TS3 1500

TS3 1510

TS3 1520

TS3 1530

). TS3 1540

TS3 1550

TS3 1560

TS3 1570

TS3 1580

TS3 1590

TS3 1600

TS3 1610

TS3 1620

TS3 1630

TS3 1640

TS3 1650

TS3 1660

TS3 1670

TS3 1680

TS3 169

TS3 1700

TS3 1710

EN TS3 1720

TS3 1730

TS3 1740

TS3 1750

TS3 1760

TS3 1770

TS3 1780

37
230

240

250

260

CONTINUE

IPROB = IPROB + 1 WRITE (NW, 90040) IPROB

IF (N.EQ.M) L = 1 IF (N.LT.M) L = 4

Ml = IWGHT = 1 TOL = 0.0

WRITE (NW, 90050)

WRITE (NW,90060) M,

t Ml , L. IWGHT, MODE. TOL

IF (N.EQ.M) GO TO 250

WRITE (NW«90070) ( I RHO( I ) • 1=1 • MMN

WRITE (NWo90080) K3,KK

DO 260 1=1.

RHO( I ) = S( It I)

CONTINUE

WRITE

WRITE

(NW 8 90090)

(NW,90100

(RHO( I ), 1=1, N)

PRINT A, B AND W.

270

WRITE (NW, 90110)

DO 270 1=1,

WRITE (NW, 90100)

CONTINUE

(A(I,J),J=1,N),(B(I,K),K=1,L).W(I)

IF (MODE.EQ.l) CALL

Nl, IPIVOT, X. RES.

L2A (M, N, Ml, L,

R, Q. C. IFAULT)

A, B, W, TOL, MM, NN<

IF (M0DE.EQ.2) CALL L2B (M,

* Nl, IPIVOT, X, RES, QR, C,

CALL

* RES,

OUTPUT

W, X,

(B.

LL,

C.

MM,

IFAULT.

MMPNN,

N. Ml. L. A. B* W, TOL, MM, NN, MMPNN,

IFAULT)

IPIVOT. L. M. MODE. Ml. N, Nl, QR, R,

MNL, NN, FAIL, SDX)

280

C

80000

80010

C

90000

90010

90020

90030

90040

90050

IF (N.EQ.M) GO TO 280

N = 2*N

GO TO 50

M = 2*M

IF (M.LE.16) GO TO 40

IF (KA.LT.NK3) GO TO 30

RETURN

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

* 3HTOL)

90060 FORMAT

90070 FORMAT

90080 FORMAT

90090 FORMAT

(15)

(1413)

(1H1)

(25H0HADAMARD MATRIX OF ORDER, 13/)

( 16F6.0)

(1H0//23H NUMBER OF K3 VALUES IS, 15)

( 1H1, 1 15

1H*)»4X,7HPR0BLEM, 14)

( 1H0.3X.1HM.4X.1HN.3X.2HM1 . 4X . 1 HL. 5X . 5HI WGHT .6X . 4HM0DE * 7X,

(415, 2110, G15.8)

(5H0IRHO.24I5)

(13H0K3 AND 2**K3, I 5, I 1 ) (23H0DIAGONAL ELEMENTS OF S/

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

TS3

1790

1800

1810

1820

1830

1840

1850

1860

1870

1880

1890

1900

1910

1920

1930

1940

1950

1960

1970

1980

1990

2000

2010

2020

2030

2040

2050

2060

2070

2080

2090

2100

2110

2120

2130

2140

2150

2160

2170

2180

2190

2200

2210

2220

2230

2240

2250

2260

2270

2280

2290

2300

2310

2320

2330

2340

2350

2360

38
90100 FORMAT (1X.8G15.8) TS3 2370

90110 FORMAT (41H0MATRIX A. MATRIX B AND VECTOR OF WEIGHTS/) TS3 2380

C TS3 2390

END TS3 2400

SUBROUTINE TEST4 (A* B. Ct FAIL. IPIVOT. LL f MM. MMPNN. MNL. M0DE.TS4 0010

* NNt NNPLLt NNP1* Q, QR , R, RES, SDX, SF • W, X) TS4 0020

C TS4 0030

C SUBROUTINE TEST4 GENERATES DATA FOR WEIGHTED POLYNOMIALS. THE METHOD TS4 0040

C IS BASED ON — TS4 0050

C 1. COMMENT ON THE ITERATIVE REFINEMENT OF LEAST-SQUARES SOLUTIONS. TS4 0060

C BY AKE BJORCK. JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION. TS4 0070

C VOL. 73. 1978. PP • 161-166. TS4 0080

C 2. A REPORT ON THE ACCURACY OF SOME WIDELY USED LEAST SQUARES TS4 0090

C COMPUTER PROGRAMS. BY ROY H. WAMPLER. JOURNAL OF THE AMERICAN TS4 0100

C STATISTICAL ASSOCIATION. VOL. 65. 1970. PP • 549-565. TS4 0110

C TS4 0120

C THE SEQUENCE OF INPUT CARDS (OR CARD IMAGES) FOR THIS SUBROUTINE IS —TS4 0130

C 1. CARD IN (15) FORMAT GIVING VALUE OF THE SCALING PARAMETER ISCALE. TS4 0140

C ISCALE = MEANS THAT DATA ARE NOT TO BE SCALED. TS4 0150

C ISCALE = 1 MEANS DATA MATRIX A IS TO BE SCALED. BUT MATRIX B TS4 0160

C IS NOT TO BE SCALED. SCALE FACTORS WILL BE AUTOMATICALLY TS4 0170

C COMPUTED FROM THE DATA. THE J-TH SCALE FACTOR. SF(J). IS TS4 0180

C BASED ON THE EUCLIDEAN NORM OF THE J-TH COLUMN OF A

TS4 0190

C ISCALE = 2 MEANS DATA MATRICES A AND B ARE TO BE SCALED. TS4 0200

C MATRIX A IS SCALED IN THE SAME MANNER AS WHEN ISCALE = 1. TS4 0210

C AND MATRIX E IS SCALED IN A SIMILAR MANNER. TS4 0220

C ISCALE = 3 MEANS DATA MATRICES A AND B ARE TO BE SCALED. THE TS4 0230

C SCALE FACTORS ARE FURNISHED BY THE USER. THEY ARE READ IN TS4 0240

C AS ITEM 4. DESCRIBED BELOW. ONLY WHEN ISCALE = 3. TS4 0250

C 2. CARD IN (215) FORMAT GIVING VALUES OF MIN AND MAX. TS4 0260

C MIN AND MAX WILL BE USED IN DETERMINING WEIGHTS TO BE APPLIED TO TS4 0270

C THE FIRST. MIDDLE AND LAST OBSERVAT I IONS. THESE OBSERVATIONS ARE TS4 0280

C GIVEN WEIGHT W = 2**(IP-1) WHERE IP = MI N. MI N+ 1 . M IN+2 • • • • • MAX. TS4 0290

C 3. CARDS IN (F3.0.F9.0.F7.0) FORMAT GIVING VALUES FOR AA(I*2). TS4 0300

C Y2(I) AND DEL(I). I =1 . 2. • . • • 21 . FROM THESE QUANTITIES AND TS4 0310

C FROM THE WEIGHTS W. THE MATRICES A AND B ARE CONSTRUCTED. TS4 0320

C 4. CARD<S) GIVING SCALE FACTORS IN (8F10.0) FORMAT. TO BE READ ONLY TS4 0330

C IF ISCALE = 3. TS4 0340

C TS4 0350

INTEGER IPIVOT(NN) TS4 0360

REAL A(MM.NNPl). B(MM.LL). C(MNL) TS4 0370

REAL Q(MM.NN). R(NN.NN). QR(MMPNN.NN) TS4 0380

REAL RES(MM.LL). SDX(NN). SF( NNPLL ) . W ( MM ) . X(NN.LL) TS4 0390

REAL AA(21«6). DEL(21). T» TOL. U. Y2(21) TS4 0400

LOGICAL FAIL(LL) TS4 0410

C TS4 0420

C IN THE FOLLOWING DATA STATEMENT. NR IS THE CARD READER DEVICE TS4 0430

C AND NW IS THE PRINTER DEVICE NUMBER. TS4 0440

C TS4 0450

DATA NR, NW / 5. 6 / TS4 0460

DATA IONE. ITWO. ITHREE / 1. 2. 3 / TS4 0470

C TS4 0480

READ (NR.80000) ISCALE TS4 0490

READ (NR.80010) MIN.MAX TS4 0500

IF (M0DE.NE.2) MODE = 1 TS4 0510

IPROB = TS4 0520

39
M = 21

N = 6

Ml =

L = 1 IWGHT = 2

TOL = 0.0

NM1 = N - 1 WRITE (NW e 90000)

WRITE (NWt90010) M, N ,M1 ,L I WGHT .MODE. TOL t M I N, MAX

WRITE (NW» 90020) ISCALE

DO 20 I=1»M

READ (NFU80020) AA( I . 2) . Y2( I ) .DEL ( I

W( I ) = 1.0

AA( It 1 ) = 1.0

DO 10 J=2«NM1

JP1 = J + 1 AAUoJPl) = AA(ItJ) * AA(I»2)

10 CONTINUE

20 CONTINUE

C

NPL = N + L

IF (ISCALE. NE. 3) GO TO 40

READ (NR, 80030) ( SF ( J ) . J=l .NPL ) DO 30 J=l e NPL

IF (SF( J).EQ.0.0) SF(J) = 1.0

30 CONTINUE

C

40 DO 80 IP = MINC MAX

DO 60 1=1.

IF (I.EQ.l .OR. I.EQ.ll .OR. I.EQ.21) W(I) = 2.**(IP-1)

DO 50 J=leN

A(IeJ) = AA(I.J) * W( I)

50 CONTINUE

B( Is 1) = Y2( I) * W( I) T = (2.**(IP-D) * (1.0/W(D)

T = T * DEL( I)

B( It 1 ) = B( It 1 ) + T

W( I) = 1.0

60 CONTINUE

IPROB = IPROB + 1 C

C PRINT A, Be We IP AND PROBLEM NUMBER.

C

WRITE (NW, 90030) IP, IPROB

DO 70 I=ltM

WRITE (NWt90040) ( A ( It J ) t J = 1 1 N ) . ( B( I • K ) t K = l ,L ) . W( I

70 CONTINUE

C

IF ( ISCALE. GT.O) CALL SCALE (ISCALE, ITWO, Me Ne Lt Wt MM. NNe

* NNe At Bt Re RESe SF, Ut Xi IFAULT)

C

IF ( ISCALE. GT.O) WRITE (NW, 90050)

IF ( ISCALE. GT.O ) WRITE (NWt90040) ( SF ( J ) * J= It NPL

C

IF (MODE.EQ.l) CALL L2A (M, N, Ml, L, A, B, W, TOL, MM, NN,

* Nit IPIVOTo Xt RESe Re Qe Ce IFAULT)

C

IF (M0DE.EQ.2) CALL L2B (M, N, Ml, L, A, B, W, TOL, MM, NN,

40

TS4 0530

TS4 0540

TS4 0550

TS4 560

TS4 0570

TS4 0580

TS4 590

TS4 0600

TS4 0610

TS4 0620

TS4 0630

TS4 0640

TS4 0650

TS4 0660

TS4 0670

TS4 0680

TS4 0690

TS4 0700

TS4 0710

TS4 0720

TS4 0730

TS4 0740

TS4 750

TS4 0760

TS4 0770

TS4 0780

TS4 0790

TS4 0800

TS4 0810

TS4 0820

TS4 0830

TS4 0840

TS4 0850

TS4 0860

TS4 0870

TS4 0880

TS4 0890

TS4 0900

TS4 0910

TS4 0920

TS4 0930

TS4 0940

TS4 0950

TS4 0960

TS4 0970

TS4 0980

TS4 0990

TS4 1000

TS4 1010

TS4 1020

TS4 1030

TS4 104

TS4 1050

TS4 1060

TS4 1070

TS4 1080

TS4 1090

TS4 1100
* MMPNN. Nl« IPIVOT. X. RES. QR. C. IFAULT) TS4 1110

C TS4 1120

IF (ISCALE.GT.O .AND. MODE.EQ.l) CALL SCALE (ISCALEo ITHREE. TS4 1130

* M, N. L c W. MM* NN. NN, A, Be R. RES* SF * U. X. IFAULT) TS4 1140

C TS4 1150

IF (ISCALE.GT.O .AND. M0DE.EQ.2) CALL SCALE (ISCALE* ITHREE* TS4 1160

* M* N 9 Lb W* MM* MMPNN. NN« A* B* QR* RES* SF * U, X, IFAULT) TS4 1170

C TS4 1180

CALL OUTPUT (B* C, IFAULT* IPIVOT. L* M» MODE* Ml* N* Nl* QR* R. TS4 1190

* RES* W. X* LL» MM, MMPNN* MNL* NN* FAIL. SDX

)

TS4 1200

C TS4 1210

80 CONTINUE TS4 1220

RETURN TS4 1230

C TS4 1240

80000 FORMAT (15) TS4 1250

80010 FORMAT (215) TS4 1260

80020 FORMAT ( F3 . 0* F9. 0* F7 .0

)

TS4 1270

80030 FORMAT (3F10.0) TS4 1280

C TS4 1290

90000 FORMAT ( 1 HO. 3X. 1 HM. 4X. 1 HN. 3X. 2HM1 » 4X» 1 HL* 5X » 5HI WGHT. 6X. 4HM0DE. 7X. TS4 1300

* 3HT0L*7X»3HMIN.2X* 3HMAX) TS4 1310

90010 FORMAT ( 4 15 . 2 1 1 . G15. 8 , 2 I 5

)

TS4 1320

90020 FORMAT (9H0ISCALE =.12) TS4 1330

90030 FORMAT (41H1MATRIX A* MATRIX B AND VECTOR OF WE IGHTS, 5X ,4H IP =,I4,TS4 134

* 30X,7HPROBLEM,I4/) TS4 1350

90040 FORMAT (IX.8G15.8) TS4 1360

90050 FORMAT (14H0SCALE FACTORS/) TS4 1370

C TS4 1380

END TS4 1390

SUBROUTINE TEST5 (A. B* C* FAIL* IPIVOT* LL » MM, MMPNN. MNL* M0DE.TS5 0010

* NN, NNPLL, NNP1. Q, QR . R, RES, SDX, SF, W, X) TS5 0020

C TS5 0030

C SUBROUTINE TESTS READS AND PRINTS INPUT DATA FOR LEAST SQUARES TS5 0040

C PROBLEMS. CALLS SUBROUTINE SCALE TO SCALE THE DATA (IF SCALING IS TS5 0050

C REQUESTED), CALLS EITHER SUBROUTINE L2A OR L2B TO COMPUTE SOLUTIONS, TS5 0060

C AND CALLS SUBROUTINE OUTPUT TO PRINT COMPUTED RESULTS. TS5 0070

C TS5 0080

C THE USER INDICATES THAT DATA ARE TO BE SCALED BY GIVING A POSITIVE TS5 0090

C VALUE (1* 2 OR 3) TO THE PARAMETER ISCALE. THE PURPOSE OF SCALING TS5 0100

C IS TO MITIGATE, IF POSSIBLE, ROUNDING ERROR PROBLEMS WHICH CAN OCCUR TS5 0110

C IN CONNECTION WITH SOLVING ILL-CONDITIONED SYSTEMS OF EQUATIONS. TS5 0120

C TS5 0130

C SEE MAIN PROGRAM FOR A DESCRIPTION OF THE 21 PARAMETERS APPEARING IN TS5 0140

C THE ARGUMENT LIST OF SUBROUTINE TESTS. TS5 0150

C TS5 0160

C THE SEQUENCE OF INPUT CARDS (OR CARD IMAGES) FOR THIS SUBROUTINE IS —TS5 0170

C 1. CARD IN (15) FORMAT GIVING VALUE OF THE SCALING PARAMETER ISCALE. TS5 0180

C THIS VALUE OF ISCALE PERTAINS TO THE COMPLETE SET OF PROBLEMS WHICHTS5 0190

C FOLLOWS. IT IS READ IN ONLY ONCE FOR ANY SET OF PROBLEMS. TS5 0200

C ISCALE = MEANS THAT DATA ARE NOT TO BE SCALED. TS5 0210

C ISCALE = 1 MEANS DATA MATRIX A IS TO BE SCALED, BUT MATRIX B TS5 0220

C IS NOT TO BE SCALED. SCALE FACTORS WILL BE AUTOMATICALLY TS5 0230

C COMPUTED FROM THE DATA. THE J-TH SCALE FACTOR, SF ( J ) * IS TS5 0240

C BASED ON THE EUCLIDEAN NORM OF THE J-TH COLUMN OF A

TS5 0250

C ISCALE = 2 MEANS DATA MATRICES A AND B ARE TO BE SCALED. TS5 0260

C MATRIX A IS SCALED IN THE SAME MANNER AS WHEN ISCALE = 1* TS5 0270

41
AND MATRIX B IS SCALED IN A SIMILAR MANNER.

ISCALE = 3 MEANS DATA MATRICES A AND B ARE TO BE SCALED. THE

SCALE FACTORS ARE FURNISHED BY THE USER. THEY ARE READ IN

AS ITEM 6. DESCRIBED BELOW. ONLY WHEN ISCALE = 3.

PROBLEM HEADING CARD. IN (80A1) FORMAT.

PARAMETER CARD IN ( 615. 5X. F 1 0. ) FORMAT. GIVING VALUES OF THE

PARAMETERS M» N. Ml. L. ITYPE. IWGHT. TOL.

M

N

Ml

L

ITYPE

IWGHT

TOTAL NUMBER OF EQUATIONS.

NUMBER OF UNKNOWN COEFFICIENTS.

NUMBER OF LINEAR CONSTRAINTS ( 0.LE.M1 .LE.M AND Ml.LE.N).

NUMBER OF RIGHT-HAND SIDES (VECTORS OF OBSERVATIONS).

PARAMETER WHICH SPECIFIES WHETHER OR NOT DATA FOR A POLYNOMIAL TYPE FIT ARE TO BE READ IN.

ITYPE = 1 INDICATES POLYNOMIAL TYPE.

ITYPE = 2 INDICATES NON-POLYNOMIAL TYPE.

PARAMETER WHICH SPECIFIES WHETHER OR NOT WEIGHTS ARE TO

BE READ IN.

IWGHT = 1 INDICATES WEIGHTS ARE NOT TO BE READ IN. (THE

PROGRAM SETS ALL WEIGHTS EQUAL TO 1.0.)

INDICATES WEIGHTS ARE TO BE READ IN.

USED IN DETERMINING THE RANK OF MATRIX H.

(SQRT(W) )*A.

ZERO. THE TOLERANCE USED IN THE

SUBROUTINE WILL BE BASED ON MACHINE

6.

IWGHT = 2

TOL PARAMETER

WHERE H = NOTE —

( 1 ) IF TOL EQUALS

DECOMPOSITION

PRECISION.

(2) IF TOL IS GREATER THAN ZERO. THIS VALUE OF TOL WILL BE

USED IN SETTING AN ABSOLUTE TOLERANCE FOR COMPARISON

WITH DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX OBTAINED

IN THE DECOMPOSITION SUBROUTINE. THE VALUE OF TOL CAN

BE BASED ON KNOWLEDGE CONCERNING THE ACCURACY OF THE

DATA.

CARD GIVING FORMAT OF THE DATA CARDS (CONTAINING A. B AND

POSSIBLY W) WHICH FOLLOW. THIS FORMAT CARD IS IN (80A1)

FORMAT.

DATA CARDS FOR THE ARRAYS A. B AND POSSIBLY W. THERE ARE FOUR

POSSIBLE CONFIGURATIONS FOR THE DATA, DEPENDING ON THE VALUES

OF ITYPE AND IWGHT. (FOR POLYNOMIAL FITS. THE FIRST POWER OF A IS READ IN AND HIGHER POWERS ARE COMPUTED BY THE PROGRAM WHEN

N.GT.2.) THE FOUR CONFIGURATIONS ARE ILLUSTRATED BELOW BY

SHOWING WHAT THE CARD (OR CARDS) FOR THE I-TH ROW OF OATA

CONTAINS.

A. ITYPE =1* IWGHT =1.

POLYNOMIAL TYPE FIT. EQUAL WEIGHTS. NOT TO BE READ IN.

A(1.2) B(I.l) B(1.2) ... B(I.L)

B. ITYPE = 1, IWGHT = 2.

POLYNCMIAL TYPE FIT. UNEQUAL WEIGHTS. TO BE READ IN.

A(I«2) B(I.l) B(1.2) ... B(I.L) W(I)

C. ITYPE = 2. IWGHT =1.

NON-POLYNOMIAL TYPE FIT. EQUAL WEIGHTS.

MX. 1) A(I,2) ... A(I.N) B(I,1) B(1.2) ...

D. ITYPE = 2. IWGHT = 2.

NON-POLYNOMIAL TYPE FIT. UNEQUAL WEIGHTS.

A(I.l) A(1.2) ... A(I.N) B(I.l) B(1.2) ... *** OMIT THIS ITEM UNLESS ISCALE = 3. *** WHEN ISCALE = 3. SCALE FACTORS ARE TO BE READ IN.

IN (8F10.0) FORMAT. A TOTAL OF N + L SCALE FACTORS ARE

TO BE FURNISHED. THE FIRST N PERTAINING TO THE N COLUMNS OF

NOT TO BE

• B(I.L)

READ IN,

TO BE READ IN« B(I.L) W(I)

TS5

TS5

TS5

TSS

TS5

TS5

TS5

TSS

TS5

TS5

TS5

TSS

TS5

TSS

TS5

TS5

TS5

TS5

TSS

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TSS

TS5

TS5

TS5

TS5

TS5

TS5

TSS

TS5

TSS

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TSS

TS5

TS5

TSS

TS5

TSS

TS5

TS5

TS5

TS5

TS5

0280

0290

0300

0310

0320

0330

0340

350

0360

0370

0380

0390

0400

0410

0420

0430

0440

0450

0460

0470

0480

0490

0500

0510

0520

0530

0540

0550

0560

0570

0580

0590

0600

0610

0620

0630

0640

0650

660

0670

0680

0690

0700

0710

0720

0730

0740

0750

0760

0770

0780

0790

0800

0810

0820

0830

0840

0850

42
C MATRIX A AND THE LAST L PERTAINING TO THE L COLUMNS OF B. TS5 0860

C EACH SCALE FACTOR IS USED TO MULTIPLY THE CORRESPONDING TS5 0870

C COLUMN OF A OR B BEFORE THE LEAST SQUARES FIT IS PEFORMED. TS5 0880

C 7. CARD IN (15) FORMAT GIVING VALUE OF THE PARAMETER IFDONE. IF TS5 0890

C THE PROBLEM AT HAND IS TO BE FOLLOWED BY ANOTHER PROBLEM, TS5 0900

C IFDONE = 1. OTHERWISE, IFDONE EQUALS ANY INTEGER EXCEPT 1. TS5 0910

C IF IFDONE = 1. GO JO ITEM 2 ABOVE FOR THE NEXT HEADING CARD. TS5 0920

C TS5 0930

INTEGER IPIVOT(NN) TS5 0940

INTEGER IFMT(80), IHEAD(80) TS5 0950

REAL A(MM»NNP1), B(MMt LL)t C(MNL) TS5 0960

REAL Q(MM,NN), R(NN.NN), QR(MMPNN.NN) TS5 0970

REAL RES(MM,LL), SDX(NN). SF(NNPLL). W ( MM ) * X(NN,LL) TS5 0980

REAL TOL, U TS5 0990

LOGICAL FAIL<LL) TS5 1000

C TS5 1010

C IN THE FOLLOWING DATA STATEMENT, NR IS THE CARD READER DEVICE TS5 1020

C AND NW IS THE PRINTER DEVICE NUMBER, TS5 1030

C TS5 1040

DATA NR, NW / 5, 6 / TS5 1050

C TS5 1060

DATA IONE, ITWO, ITHREE / 1, 2, 3 / TS5 1070

C TS5 1080

READ (NR, 80000) ISCALE TS5 1090

IF (ISCALE. LT.O .OR. ISCALE. GT. 3) ISCALE = TS5 1100

WRITE (NW, 90000) ISCALE TS5 1110

IPROB = TS5 1120

10 READ (NR,80010) IHEAD TS5 1130

IPROB = IPROB + 1 TS5 1140

IFAULT = TS5 1150

U = 0.0 TS5 1160

WRITE (NW, 90010) IPROB TS5 1170

WRITE (NW,90020) IHEAD TS5 1180

READ (NR,80020) M, N, Ml , L, I TYPE , I WGHT, TOL TS5 1190

WRITE (NW,90030) TS5 1200

WRITE (NW, 90040) M, N, Ml , L, ITYPE, IWGHT, MODE, I SCALE, TOL TS5 1210

READ (NR,80010) IFMT TS5 1220

WRITE (NW, 90050) IFMT TS5 1230

GO TO (20,100), ITYPE TS5 1240

C TS5 1250

C TYPE 1. POLYNOMIAL FIT. TS5 1260

20 GO TO (30,50), IWGHT TS5 1270

C A. EQUAL WEIGHTS, NOT TO BE READ IN. TS5 1280

30 DO 40 I=1,M TS5 1290

READ (NRelFMT) A ( I , 2 ) , ( B( I, K ) , K= 1 , L

)

TS5 1300

W(I) = 1.0 TS5 1310

40 CONTINUE TS5 1320

GO TO 70 TS5 1330

C B. UNEQUAL WEIGHTS, TO BE READ IN. TS5 1340

50 DO 60 1=1,

TS5 1350

READ (NR,IFMT) A( I , 2 ) . ( B ( I, K ) , K= 1 ,L ) , W( I ) TS5 1360

60 CONTINUE TS5 1370

C TS5 1380

C CALL SUBROUTINE SCALE TO COMPUTE MEAN OF VECTOR A(I,2), IF DATA ARE TS5 1390

C TO BE SCALED. TS5 1400

C TS5 1410

70 IF ( ISCALE. GT.O .AND.N.GT.l) CALL SCALE (ISCALE, IONE, M, N, TS5 1420

* L, W, MM. NN, NN, A, B, R, RES, SF , U, X, IFAULT) TS5 1430

43
IF DATA ARE SCALED* PRINT U. IN SCALED POLYNOMIAL-TYPE PROBLEMS. U

EQUALS THE MEAN OF THE INDEPENDENT VARIABLE (VECTOR A(I,2)).

U IS SUBTRACTED FROM EACH ELEMENT OF A(1.2) BEFORE POWERS OF

THIS VECTOR ARE GENERATED.

IF ( ISCALE.GT.O) WRITE (NW. 90060) U

DO 90 1 = 1,

A(I. 1) = 1.0

IF (N.EQ.l ) GO TO 90

A{ 1,2) = A( 1,2) - U

IF (N.EQ.2) GO TO 90

DO 80 J=3.N

A(I.J) = AC I.2)**(J-1)

80 CONTINUE

90 CONTINUE

GO TO 150

TYPE 2. NON-POLYNOMIAL FIT.

100 GO TO (110,130). IWGHT

C. EQUAL WEIGHTS. NOT TO BE READ IN.

110 DO 120 1=1.

READ (NR.IFMT) ( A< I, J ) • J=l • N ) . ( B( I, K) * K=l ,L

W(I) = 1.0

120 CONTINUE

GO TO 150

D. UNEQUAL WEIGHTS, TO BE READ IN.

130 DO 140 1=1.

READ (NR.IFMT) ( A( I. J ) . J=l . N ) » ( B( I. K ) , K=l ,L) . W( I

140 CONTINUE

150 IF (M1.EQ.0 .OR. IWGHT. EQ . 1 ) GO TO 170

INSURE THAT WEIGHTS EQUAL 1.0 FOR THE FIRST Ml EQUATIONS WHEN Ml

IS GREATER THAN ZERO.

DO 160 1=1. Ml

W( I ) = 1.0

160 CONTINUE

PRINT A, B AND W.

170 WRITE (NW, 90070)

KZ = DO 180 1=1,

IF (W(I).EQ.O.O) KZ = KZ + 1 WRITE (NW.90080) ( A( I • J ) . J= 1 • N) , ( B( I . K ) . K=l 9 L ) , W ( I

180 CONTINUE

NPL = N + L IF (ISCALE.NE.3) GO TO 200

READ (NR, 80030) ( SF( J ) , J=l .NPL

DO 190 J=1.NPL

IF (SF(J) .EQ.O.O) SF(J) = 1.0

190 CONTINUE

CALL SUBROUTINE SCALE TO SCALE DATA MATRICES A AND B, IF REQUESTED.

SCALE FACTORS WILL BE COMPUTED IF ISCALE EQUALS 1 OR 2.

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

TS5

1440

1450

1460

1470

1480

1490

1500

1510

1520

1530

1540

1550

1560

1570

1580

1590

1600

1610

1620

1630

1640

1650

1660

1670

1680

1690

1700

1710

1720

1730

1740

1750

1760

1770

1780

1790

1800

1810

1820

1830

1840

1850

1860

1870

1880

1890

1900

1910

1920

1930

194

1950

1960

1970

1980

1990

2000

2010

44
200 IF (ISCALE.GT.O) CALL SCALE (ISCALE, ITWO, M, N. L. W, MM* TS5 2020

* NN, NNt A* B. Rt RES* SFt U. X* IFAULT) TS5 2030

C TS5 2040

IF (IFAULT. EQ.4) GO TO 240 TS5 2050

C TS5 2060

C PRINT SCALE FACTORS (SF) IF DATA WERE SCALED. TS5 2070

C TS5 2080

IF (ISCALE.GT.O) WRITE (NW, 90090) TS5 2090

IF (ISCALE.GT.O) WRITE (NW, 90080) ( SF( J ) . J= 1 .NPL) TS5 2100

GO TO (210*220)* MODE TS5 2110

C TS5 2120

C CALL SUBROUTINE L2A OR L2B TO COMPUTE LEAST SQUARES SOLUTIONS. TS5 2130

C TS5 2140

210 CALL L2A (M* N* Ml* L. A* B. W* TOL* MM* NN* TS5 2150

* Nl* IPIVOT* X* RES. R. Q. C* IFAULT) TS5 2160

C TS5 2170

GO TO 230 TS5 2180

C TS5 2190

220 CALL L2B (M. N. Ml. L. A. B* W* TOL. MM, NN. MMPNN* TS5 2200

* Nl. IPIVOT* X* RES* QR* C* IFAULT) TS5 2210

C TS5 2220

230 IF (IFAULT. GE.l .AND. IFAULT. LE. 4) GO TO 240 TS5 2230

C TS5 2240

C IF DATA WERE SCALED. CALL SUBROUTINE SCALE TO ADJUST COMPUTED RESULTS TS5 2250

C (COEFFICIENTS* RESIDUALS AND COVARIANCE MATRIX) FOR SCALING. TS5 2260

C TS5 2270

IF (ISCALE.GT.O .AND. MODE.EQ.l) CALL SCALE (ISCALE. ITHREE. TS5 2280

* M. N* L* W* MM* NN* NN* A* B* R* RES* SF* U* X* IFAULT) TS5 2290

C TS5 2300

IF (ISCALE.GT.O .AND. MODE.EQ.2) CALL SCALE (ISCALE* ITHREE. TS5 2310

* M. N» L* W* MM* MMPNN* NN* A* B* OR* RES* SF* U* X* IFAULT) TS5 2320

C TS5 2330

C CALL SUBROUTINE OUTPUT TO PRINT COMPUTED RESULTS. TS5 2340

C TS5 2350

240 CALL OUTPUT (B. C. IFAULT. IPIVOT. L. M. MODE, Ml. N, Nl , QR, R, TS5 2360

* RES* W* X* LL* MM* MMPNN* MNL* NN* FAIL* SDX

)

TS5 2370

C TS5 2380

READ (NR*80000) IFDONE TS5 2390

IF ( IFDONE. EQ.l ) GO TO 10 TS5 2400

C TS5 2410

RETURN TS5 2420

C TS5 2430

80000 FORMAT (15) TS5 2440

80010 FORMAT (80A1) TS5 2450

80020 FORMAT ( 615* 5X.F10.0

)

TS5 2460

80030 FORMAT (8F10.0) TS5 2470

C TS5 2480

90000 FORMAT (9H0ISCALE =* 15) TS5 2490

90010 FORMAT ( 1H1 , 1 13( 1 H*) ,4X,7HPR0BLEM , 14) TS5 2500

90020 FORMAT (1H0*80A1) TS5 2510

90030 FORMAT ( 1H0* 3X» 1HM* 4X. 1 HN* 3X. 2HM1 * 4X* 1 HL* 5X» 5HITYPE* 5X. SHI WGHT.6X , TS5 2520

* 4HM0DE*4X*6HISCALE,7X,3HT0L) TS5 2530

90040 FORMAT ( 41 5* 4 1 1 0* Gl 5.8

)

TS5 2540

90050 FORMAT (8H0FORMAT *80A1) TSS 2550

90060 FORMAT (28H0U = MEAN OF VECTOR A(I,2) =,G15.8/72H U IS SUBTRACTED TS5 2560

FROM EACH ELEMENT OF A(I*2) IN CONNECTION WITH SCALING.) TS5 2570

90070 FORMAT (41H0MATRIX A* MATRIX 8 AND VECTOR OF WEIGHTS/) TS5 2580

90080 FORMAT (1X*8G15.8) TS5 2590

45
90090 FORMAT (14H0SCALE FACTORS/)

C

END

SUBROUTINE OUTPUT

* OR. R« RES* W* X, (B» C* IFAULT, IPIVOT. Lo M» MODE* Ml,

LL, MM. MMPNN, MNL, NN . FAIL* SOX)

N» Nl.

C

c

c

ccc

cc

cc

cc

cc

ccc

ccc

c

c

cc

c

c

c

cc

c

c

ccc

c

c

cccc

c

cc

c

c

c

cc

ccc

SUBROUTINE OUTPUT PRINTS THE SOLUTIONS OF LEAST SQUARES PROBLEMS

WHICH WERE OBTAINED FROM EITHER SUBROUTINE L2A OR SUBROUTINE L2B.

SEE MAIN PROGRAM FOR

INPUT VARIABLES

B TWO-DIMENSIONAL ARRAY OF SIZE (MM.LL).

FURTHER DETAILS.

C VECTOR OF LENGTH MNL. SEE MAIN PROGRAM.

IFAULT FAULT INDICATOR WHICH IS ZERO IF NO ERRORS WERE ENCOUNTERED

AND POSITIVE IF ERRORS WERE DETECTED OR IF EVIDENCE OF SEVERE

ILL-CONDITIONING WAS FOUND. DIAGNOSTIC MESSAGES ARE PRINTED

FROM SUBROUTINE ERROR. IF IFAULT IS SET EQUAL TO 1» 2* 3* 4*

5 6 OR 7, EXECUTION IS TERMINATED. EXECUTION CONTINUES WHEN

IFAULT IS SET EQUAL TO 8. 9 OR 10 PROVIDED THAT A SOLUTION

WAS OBTAINED FOR AT LEAST ONE RIGHT-HAND SIDE. THE VALUE OF

IFAULT IS USED TO INDICATE THE FOLLOWING —

12

3

4

5

6

7

8

9

10

OR L).

M+M.GT. MMPNN.

11 =

IPIVOT

L

MMODE

Ml

N

Nl

QR

R

RES

W

XLL

MM

MMPNN

MNL

NN

NO ERRORS ENCOUNTERED.

BAD INPUT PARAMETER <M« N

BAD INPUT PARAMETER (Ml).

BAD DIMENSION. EITHER M.GT.MM, N.GT.NN OR

AT LEAST ONE WEIGHT IS NEGATIVE.

EITHER MATRIX H OR MATRIX OF CONSTRAINTS EQUALS ZERO.

CONSTRAINTS ARE LINEARLY DEPENDENT.

ALL SOLUTIONS FAILED TO CONVERGE.

SOLUTION FAILED TO CONVERGE FOR AT LEAST ONE RIGHT-HAND

SIDE.

LARGE NUMBER OF ITERATIONS REQUIRED FOR CONVERGENCE.

ESTIMATED NUMBER OF DIGITS IN INITIAL SOLUTION OF

COEFFICIENTS IS SMALL.

DIAGONAL ELEMENT OF COVARIANCE MATRIX WAS COMPUTED TO BE

NEGATIVE OWING TO ROUNDING ERROR.

VECTOR OF LENGTH NN. SEE MAIN PROGRAM.

NUMBER OF RIGHT-HAND SIDES (VECTORS OF OBSERVATIONS).

TOTAL NUMBER OF EQUATIONS.

PARAMETER WHICH INDICATES WHETHER SOLUTION WAS OBTAINED FROM

L2A OR L2B. SEE MAIN PROGRAM FOR FURTHER DETAILS.

NUMBER OF LINEAR CONSTRAINTS.

NUMBER OF UNKNOWN COEFFICIENTS.

COMPUTED RANK OF THE SYSTEM OF EQUATIONS.

TWO-DIMENSIONAL ARRAY OF SIZE (MMPNN, NN). SEE MAIN PROGRAM.

TWO-DIMENSIONAL ARRAY OF SIZE (NN*NN). SEE MAIN PROGRAM.

TWO-DIMENSIONAL ARRAY OF SIZE (MM.LL). SEE MAIN PROGRAM.

VECTOR OF LENGTH MM. SEE MAIN PROGRAM.

TWO-DIMENSIONAL ARRAY OF SIZE (NNtLL). SEE MAIN PROGRAM.

DIMENSIONING PARAMETER. SEE MAIN PROGRAM.

SEE MAIN PROGRAM.

SEE MAIN PROGRAM.

SEE MAIN PROGRAM.

SEE MAIN PROGRAM.

DIMENSIONING PARAMETER.

DIMENSIONING PARAMETER.

DIMENSIONING PARAMETER.

DIMENSIONING PARAMETER.

INTERNAL VARIABLES —

TS5 2600

TS5 2610

TS5 2620

OUT 0010

OUT 0020

OUT 0030

OUT 0040

OUT 0050

OUT 0060

OUT 0070

OUT 0080

OUT 0090

OUT 0100

OUT 0110

OUT 0120

OUT 0130

OUT 0140

OUT 0150

OUT 0160

OUT 0170

OUT 0180

OUT 0190

OUT 0200

OUT 0210

OUT 22

OUT 0230

OUT 0240

OUT 0250

OUT 0260

OUT 0270

OUT 0280

OUT 0290

OUT 300

OUT 0310

OUT 0320

OUT 0330

OUT 0340

OUT 0350

OUT 0360

OUT 370

OUT 0380

OUT 390

OUT 0400

OUT 0410

OUT 0420

OUT 0430

OUT 044

OUT 0450

OUT 0460

OUT 0470

OUT 0480

OUT 0490

OUT 0500

OUT 0510

OUT 0520

OUT 0530

46
FAIL VECTOR OF LENGTH LL.

FAILED TO CONVERGE.

CONVERGED.

SDX VECTOR OF LENGTH NN.

IF FAIL(K) = .TRUE. THE K-TH SOLUTION

IF FAIL(K) = .FALSE. THE K-TH SOLUTION

STANDARD DEVIATIONS OF COEFFICIENTS (X)

INTEGER IPIVOT ( NN

)

REAL B(MM c LL)e C(MNL)* QR ( MMPNN» NN ) • R(NN*NN). RES(MM»LL)

REAL W(MM), X(NN*LD* SDX(NN)

REAL SD. SNORM. SS. WSQt Z DOUBLE PRECISION SUM

LOGICAL FAIL(LL)

IN THE FOLLOWING DATA STATEMENT. NW IS THE PRINTER DEVICE NUMBER.

DATA NW / 6 /

KZ = DO 10 I=1*M

IF ( W(I).EQ.O.O) KZ = KZ + 1 10 CONTINUE

PRINT COMPUTED RESULTS.

WRITE (NW*90000)

WRITE (NW* 90010) MODE.IFAULT

IF (IFAULT.GE.l .AND. IFAULT.LE.4) GO TO 190

WRITE (NW. 90020) Nl

IF (N1.EQ.0) GO TO 190

IF (N1.LT.M1) GO TO 190

WRITE (NW. 90030)

WRITE (NW«90040) ( IP IVOT( J ) , J = 1 . N 1

)

IF (Nl.EQ.N) GO TO 20

N1P1 = Nl + 1 WRITE (NW*90050)

WRITE (NW«90040) ( IPIVOT* J ) * J=N1P1 • N)

20 NDF = M - Nl - KZ

WRITE (NW«90060) KZ.NDF

WRITE (NW. 90070)

DO 40 K=1»L

K2 = L + K IF (C(K) .LT.0.0) GO TO 30

FAIL(K) = .FALSE.

WRITE (NW. 90080) K*C(K).C(K2)

GO TO 40

30 FAIL(K) = .TRUE.

C(K) = -C(K)

WRITE (NW. 90090) K.C(K)*C(K2)

40 CONTINUE

IF (IFAULT.EQ.7) GO TO 190

DO 150 K=l .L

COMPUTE SUM OF SQUARED RESIDUALS. NORM OF RESIDUALS. RESIDUAL

STANDARD DEVIATION* STANDARD DEVIATIONS OF COEFFICIENTS* AND

PREDICTED VALUES. PRINT THESE QUANTITIES, TOGETHER WITH

COEFFICIENTS. OBSERVED VALUES. RESIDUALS AND WEIGHTS. WEIGHTS ARE

RESTORED TO THEIR ORIGINAL VALUES BEFORE BEING PRINTED. (WITHIN

SUBROUTINE L2A OR L2B. THE ORIGINAL WEIGHTS WERE REPLACED BY THEIR

SQUARE ROOTS.)

OUT 0540

OUT 0550

OUT 0560

OUT 0570

OUT 0580

OUT 0590

OUT 0600

OUT 0610

OUT 0620

OUT 630

OUT 0640

OUT 0650

OUT 0660

OUT 0670

OUT 0680

OUT 0690

OUT 0700

OUT 0710

OUT 0720

OUT 0730

OUT 0740

OUT 0750

OUT 0760

OUT 0770

OUT 0780

OUT 0790

OUT 0800

OUT 0810

OUT 0820

OUT 0830

OUT 0840

OUT 0850

OUT 0860

OUT 0870

OUT 0880

OUT 0890

OUT 0900

OUT 0910

OUT 0920

OUT 0930

OUT 0940

OUT 0950

OUT 0960

OUT 0970

OUT 0980

OUT 0990

OUT 1000

OUT 1010

OUT 102

OUT 1030

OUT 1040

OUT 1050

OUT 1060

OUT 1070

OUT 1080

OUT 1090

OUT 1100

OUT 1110

47
C OUT 1120

IF (FAIL(K)) GO TO 150 OUT 1130

WRITE (NW, 90100) K OUT 1140

SS = 0.0 OUT 1150

SUM =0.0 OUT 1160

IF (M.EQ.M1) GO TO 110 OUT 1170

M1P1 = Ml + 1 OUT 1180

DO 50 I=M1P1,M OUT 1190

SUM = SUM + DBLE(RES( ItK) )**2 * DBLE ( W( I ) **2) OUT 1200

50 CONTINUE OUT 1210

SS = SUM OUT 1220

IF ( M.LE.N1 +KZ) GO TO 110 OUT 1230

IF (MODE.EQ.2 .AND. Nl.LT.N) GO TO 110 OUT 1240

SD = SQRT(SS/FLOAT(NDF)

OUT 1250

IF (MODE.EQ.2) GO TO 70 OUT 1260

DO 60 J=1«N OUT 1270

IF (R< J« J ).LT.0.0) R(J.J) = 0.0 OUT 1280

SDX(J) = SD*SQRT(R( J, J)

OUT 1290

60 CONTINUE OUT 1300

GO TO 90 OUT 1310

70 DO 80 J=1«N OUT 1320

IF <QR( J» J) .LT.O.O) QR(J.J) = 0.0 OUT 1330

SDX(J) = SD*SQRT(QR( J, J)) OUT 1340

80 CONTINUE OUT 1350

90 WRITE (NW«90110) OUT 1360

DO 100 J=1.N OUT 1370

WRITE (NW.90120) J, X( J, K) t SDX ( JOUT 1380

100 CONTINUE OUT 1390

GO TO 130 OUT 1400

110 WRITE (NW 9 90130) OUT 1410

DO 120 J=ltN OUT 1420

WRITE (NWe90120) J»X(J«K) OUT 1430

120 CONTINUE OUT 1440

130 WRITE <NW*90140) OUT 1450

DO 140 I=l e M OUT 1460

Z = B(I»K) - RES(I.K) OUT 1470

V.SQ = WCX)*W(I) OUT 1480

WRITE (NW S 90120) I • B( I • K) • Z. RES ( I « K ) • WSQ OUT 1490

140 CONTINUE OUT 1500

SNORM = SQRT(SS) OUT 1510

WRITE (NW«90150) SStSNORM OUT 1520

IF ((MODE.EQ.2 .AND. Nl.LT.N) .OR. ( M.EQ. Nl +KZ) ) GO TO 150 OUT 1530

WRITE (NW.90160) SD OUT 1540

150 CONTINUE OUT 1550

C OUT 1560

C PRINT LOWER TRIANGULAR PORTION OF SYMMETRIC UNSCALED COVARIANCE OUT 1570

C MATRIX. OUT 1580

C OUT 1590

IF (MODE.EQ.2 .AND. Nl.LT.N) GO TO 190 OUT 1600

WRITE (NW«90170) OUT 1610

IF (MODE.EQ.2) GO TO 170 OUT 1620

DO 160 I=1 9 N OUT 1630

WRITE (NWc90180) ( R ( I * J ) e J= 1 • I

OUT 1640

160 CONTINUE OUT 1650

GO TO 190 OUT 1660

170 DO 180 1=1 »N OUT 1670

WRITE (NW, 90180) ( QR( I , J ) , J= 1 . I ) OUT 1680

180 CONTINUE OUT 1690

48
190 RETURN

C

C FORMAT STA

C

90000 FORMAT

90010 FORMAT

90020 FORMAT

900 30 FORMAT

(17H0C0MPJTED RESULTS)

(7H0MODE =» I4e5Xe8HIFAULT =, 14

)

(44H0N1 = COMPUTED RANK OF SYSTEM OF EQUATIONS

(87H0COLUMNS OF H = <SQRT(W))*A WERE SELECTED

*NG SCHEME IN THE FOLLOWING ORDER/)

90040 FORMAT (3014)

90050 FORMAT (98H0THE FOLLOWING COLUMNS OF H ARE LINEARLY D

* MODE It THEY DID NOT ENTER THE REGRESS ION. /30H IF MO

TERED LAST*/)

90060 FORMAT (25H0NUMBER OF ZERO WEIGHTS =e I 3e 5X, 1 7HDEG . OF

*)

90070 FORMAT ( 1H0« 1 5Xe 9HREP0RT ONe 1 2 Xe 9HNUMBER 0F»4X«27HEST

* OF C0RRECT/13H B-VECTOR NO. • 3Xe 1 1 HCONVERGENCE. 10X * 10

X.26HDIGITS IN INITIAL SOLUTION/)

90080 FORMAT ( I 7» 9Xp 9HC0NVERGED. 1 4X 9 F4 . 0. I 0X» Gl 5. 8

)

90090 FORMAT ( I 7e 9Xe 1 8HFAI LED TO CONVERGEe 5X t F4 • • 1 OX , Gl 5 .8

90100 FORMAT ( 26H0S0LUT ION FOR B-VECTOR NO., 13)

90110 FORMAT ( 1 HO. 2 7X. 1 8HSTANDARD DEV I AT I0N/5X. 1H J, 4X S 1 4HC0

* 5Xel7H0F COEFFICIENT J)/)

90120 FORMAT ( 16 . Gl 7. 8 , 3G21 • 8

)

90130 FORMAT ( lHOe 4X. 1 H J • 4X, 1 4HC0EFF I C IENT( J ) /)

90140 FORMAT ( 1 H0« 4X, 1 HI t 4X t 1 1HOBSERVED ( I ) * 9X t 1 2HPRED ICTEDC

*IDUAL( I)«12X,9HWEIGHT( I)/)

90150 FORMAT (30H0SUM OF SQUARED RESIDUALS = • G 1 5 .8/1 X . 1

7

DUALS, 11X,1H=,G15.8>

90160 FORMAT (30H RESIDUAL STANDARD DEVIATION =«G15.8)

90170 FORMAT ( 27H0UNSCALED COVARIANCE MATRIX/)

90180 FORMAT (1X,8G15.8)

C

END

OUT 1700

OUT 1710

OUT 172

OUT 1730

OUT 1740

OUT 1750

=.14) OUT 1760

BY THE PIVOTIOUT 1770

OUT 1780

OUT 1790

EPENDENT. IFOUT 1800

DE 2, THEY ENOUT 1810

OUT 1820

FREEDOM =,I30UT 1830

OUT 1840

I MATED NUMBEROUT 1850

HITERATI0NS,30UT 1860

OUT 1870

OUT 188

) OUT 1890

OUT 1900

EFFICIENTC J). OUT 1910

OUT 192

OUT 1930

OUT 1940

I ) ,10X,11HRES0UT 1950

OUT 1960

HNORM OF RESIOUT 1970

OUT 1980

OUT 1990

OUT 2000

OUT 2010

OUT 2020

OUT 2030

SUBROUTINE ORTHO ( Ne IRe MM, MMM, X, Ye A. Be Ce De Ee BETAe Z)

SUBROUTINE ORTHO

PROGRAMS. THE AL

TEST DATA FOR

BY W. K. HASTI

ASSOCIATION, V FIRST A SET OF OR

THEN A SET OF X-V

THESE ARE USED AS

FORMULAS FOR THE

QUANTITIES.

N AS USED HERE CO

NN AS USED HERE C

X AND Y IN THIS S

IN SUBROUTINE TES

INPUT VARIABLES -

N LENGTH OF

IR NUMBER OF

GENERATES DATA (X AND Y) FOR TESTING LEAST SQUARES

GORITHM IS BASED ON EXAMPLE 2 IN

STATISTICAL ALGORITHMS -- LEAST SQUARES AND ANOVA,

NGSe JOURNAL OF THE AMERICAN STATISTICAL

OL. 67, 1972, PP. 874-879.

THONORMAL VECTORS IS CONSTRUCTED IN THE ARRAY Z.

ECTORS AND A Y-VECTOR ARE GENERATED FROM Z, AND

INPUT FOR LEAST SQUARES PROBLEMS. HASTINGS GAVE

LEAST SQUARES SOLUTIONS AND CERTAIN RELATED

RRESPONDS TO LOWER CASE N IN HASTINGS.

ORRESPONDS TO UPPER CASE N IN HASTINGS.

UBROUTINE CORRESPOND TO A AND Be RESPECT I VELYe

T2 (THE CALLING PROGRAM).

THE X-. Y- AND Z-VECTORS. N MUST BE EVEN.

X-VECTORS TO BE GENERATED.

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

0010

0020

0030

0040

0050

0060

0070

0080

0090

0100

0110

0120

0130

0140

0150

0160

0170

0180

0190

0200

0210

0220

49
MM DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN

ARRAYS Y AND Z. N.LE.MM.

MMM DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN

ARRAY X. N.LE.MMM.

OUTPUT VARIABLES —

X TWO-DIMENSICNAL ARRAY OF SIZE (MMM,IR).

Y TWO-DIMENSIONAL ARRAY OF SIZE (MM,1).

INTERNAL VARIABLES —

AB

C

D

EBETA

Z

(WORK-SPACE) OF LENGTH N + CWORK-SPACE) OF LENGTH N.

VECTOR

VECTOR

VECTOR (WORK-SPACE)

VECTOR (WORK-SPACE)

VECTOR (WORK-SPACE)

OF LENGTH N.

OF LENGTH N.

OF LENGTH N.

VECTOR (WORK-SPACE) OF LENGTH N.

TWO-DIMENSIONAL ARRAY OF SIZE (MM.IR+1).

OF Z ARE ORTHONORMAL.

THE IR + 1 COLUMNS

CONSTANTS A(K). B(J) AND C(J) REQUIRED BY THIS ALGORITHM CAN BE SET

BY THE USER. HERE WE USE THE FOLLOWING VALUES —

A(K) = K*K (K=1,2,...,IR+1

B(J) = 1000 ( J=l,2. •••« IR-l

B(IR) = 1 C(J) = 1 (J=1.2..... IR)

REAL A(l), B<1), BETA(l). C(l), DO), E<1), X(MMM.l). Y(MM.l)

REAL Z(MM. 1

REAL ALPHA. ENN. ES# PI. SSREG, SSTOT. SUMYY, T

IN THE FOLLOWING DATA STATEMENT. NW IS THE PRINTER DEVICE NUMBER.

DATA NW / 6 /

PI = 4.0 * ATAN(l.O)

M = IR + 1 NN = N/2

WRITE (NW.90000)

WRITE (NW.90010)

IF (2*NN .NE. N)

IF (2*NN .NE. N)

IF (IR.LT.2)

IF (IR.LT.2)

N.NN. IR.M

WRITE (NW, 90020)

RETURN

WRITE (NW, 90030)

RETURN

10

20

ENN = FLOAT(NN)

DO 40 J=1.M

IF (J.GT.NN+1) GO TO 20

IS = J - 1 ES = IS

ALPHA = 1 .OZSQRT(ENN)

IF (IS.EQ.O .OR. IS.EQ

DO 10 1=1,

T = I Z(I.J) = ALPHA * COS(T*PI*ES/ENN)

CONTINUE

GO TO 40 IS=J-NN-1

ES = IS

NN) ALPHA = 1.0/SQRT(2.0*ENM)

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

0230

0240

0250

0260

0270

0280

0290

0300

0310

0320

330

340

350

360

370

380

0390

0400

0410

0420

0430

0440

0450

0460

0470

0480

0490

0500

0510

0520

0530

540

0550

560

0570

0580

0590

0600

0610

0620

630

0640

0650

0660

0670

0680

0690

0700

0710

0720

0730

0740

0750

0760

0770

0780

079

0800

50
ALPHA = 1 .0/SQRT(ENN)

DO 30 I=UN

T = I Z(I,J) = ALPHA * SIN(T*PI*ES/ENN)

30 CONTINUE

40 CONTINUE

DO 50 J=1«IR

B(J) = 1000.0

C( J) = 1 .0

50 CONTINUE

B(IR) = 1.0

TO CONSTRUCT AN ILL-CONDITIONED SYSTEM, ONE OR MORE OF THE B VALUES

SHOULD BE SMALL RELATIVE TO THE OTHER B VALUES.

FORM X-VECTORS.

DO 80 1=1 S N DO 70 J=l 9 IR

X( I, J) = 0.0

DO 60 K=1,J

X(IeJ) = X(I,J) + B(K)*Z(I,K)

60 CONTINUE

X( I. J) = C( J)*X(I, J>

70 CONTINUE

80 CONTINUE

FORM Y-VECTOR.

DO 90 K=1#M

A(K) = FLOAT(K*K)

90 CONTINUE

SUMYY = 0.0

DO 120 1 = 1, N DO 1 10 L=l, 1 Y(UL) = 0.0

DO 100 K=1,M

Y<I,L) = Y(I,L) + A(K)*Z

100 CONTINUE

SUMYY = SUMYY + Y(I,L)*Y(I

110 CONTINUE

120 CONTINUE

c

c COMPUTE D, E AND BETA.

c

L)

DO 140 J=1,IR

D(J) = 0.0

E(J) = 0.0

DO 130 K=1,J

D( J) = D( J) + A<K)*B(K)

E(J) = E(J) + B(K)*B(K)

130 CONTINUE

140 CONTINUE

IR1 = IR - 1 BETA(l) = (1.0/C(1)) * (D(l)/E(l) - ( D ( 2 ) -D { 1 ) ) /( E( 2) -E( 1 ) )

)

BETA(IR) = (1.0/C(IR)) * ( < D< I R ) -D < IR 1 ) )/ < E ( IR) -E < IR 1 > > ) IF (IR1.LT.2) GO TO 160

DO 150 J=2,IR1

ORT 0810

ORT 0820

ORT 0830

ORT 84

ORT 0850

ORT 0860

ORT 87

ORT 0880

ORT 0890

ORT 0900

ORT 091

ORT 0920

ORT 0930

ORT 94

ORT 0950

ORT 0960

ORT 0970

OPT 0980

ORT 0990

ORT 1000

ORT 1010

ORT 1020

ORT 1030

ORT 1040

ORT 1050

ORT 1060

ORT 1070

ORT 1080

ORT 1090

ORT 1100

ORT 1110

ORT 1120

ORT 1 130

ORT 1140

ORT 1 150

ORT 1 160

ORT 1 170

ORT 1180

ORT 1190

ORT 1200

ORT 1210

ORT 122

ORT 1230

ORT 1240

ORT 1250

ORT 1260

ORT 1270

ORT 1280

ORT 1290

ORT 1300

ORT 1310

ORT 1320

ORT 1330

ORT 1340

ORT 1350

ORT 1360

ORT 1370

ORT 1380

51
BETA(J) = (1,0/C( J) )*< (D(J)-D(J-l ))/(E« J)-E(J-1 ) ) -

* (DC J + l )-DC J))/(E(J+ 1 )-E(J) ) ) 150 CONTINUE

160 SSREG a 0.0

DO 170 J=ltIR

SSREG = SSREG C( J) *D( J )*BETA< J

170 CONTINUE

SSTOT =0.0

DO 180 K = l*M

SSTOT = SSTOT A(K)*A(K)

180 CONTINUE

PRINT REGRESSION SUM

Y- SQUARED.

WRITE (NW, 90040)

WRITE (NW t 90050)

WRITE (NW. 90060)

RETURN

OF SQUARES. TOTAL SUM OF SQUARES. AND SUM OF

(BETA( J), J=l, IR)

SSREG* SSTOT. SUMYY

C

90000

90010

90020

90030

90040

90050

90060

C

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

END

(20H0 N

(415)

(48H0***

(40H0***

(46H0BETA

(IX.8G15.8)

(20H0SSREG.

NN IR M)

ERROR — N (LENGTH OF VECTOR) IS NOT

ERROR — IR CANNOT BE LESS THAN 2./)

COEFFICIENTS CALCULATED FROM C* D AND

SSTOT. SUMYY.5X.3G15.8)

EVEN./)

E./)

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

ORT

1390

1400

1410

142

1430

1440

1450

1460

1470

1480

1490

1500

1510

1520

1530

1540

1550

1560

1570

1580

1590

160

1610

1620

1630

1640

1650

1660

1670

FUNCTION GEN(ANOISE)

C CHARLES L. LAWSON AND RICHARD J. HANSON.

C SOLVING LEAST SQUARES PROBLEMS. COPYRIGHT 1974. PAGE 311.

C REPRINTED BY PERMISSION OF PRENTICE-HALL* INC.. ENGLEWOOD CLIFFS*

C GENERATE NUMBERS FOR CONSTRUCTION OF TEST CASES.

IF (ANOISE) 10.30.20

10 MI=891

MJ=457

1=5

J=7

AJ = 0.

GEN=0.

RETURN

C

C THE SEQUENCE OF VALUES OF J IS BOUNDED BETWEEN 1 AND 996

C IF INITIAL J = 1.2*3.4,5.6.7*8. OR 9. THE PERIOD IS 332

20 J=J*MJ

J=J-997*( J/997)

AJ=J-498

C THE SEQUENCE OF VALUES OF I IS BOUNDED BETWEEN 1 AND 999

C IF INITIAL I = 1.2.3.6.7. OR 9. THE PERIOD WILL BE 50

C IF INITIAL I = 4 OR 8 THE PERIOD WILL BE 25

C IF INITIAL 1=5 THE PERIOD WILL BE 10

30 I=I*MI

I=I-1000*( 1/1000)

AI=I-500

GEN=AI+AJ* ANOISE

GEN0O100

GEN

GEN

,J.GEN

GEN00400

GEN00500

GEN00600

GEN00700

GEN00800

GEN00900

GEN01000

GEN01 100

GEN01200

GEN01300

GEN01400

GEN01500

GEN01600

GEN01700

GEN01800

GEN01900

GEN02000

GEN02100

GEN02200

GEN02300

GEN02400

GEN02500

GEN02600

52
cc

cc

cccc

c

RETURN

END

SUBROUTINE HADMAR (N» H)

SUBROUTINE HADMAR GENERATES HADAMARD MATRICES OF ORDER 2» 4. 8 AND 16

THESE WILL BE USED IN THE CONSTRUCTION OF LEAST SQUARES TEST PROBLEMS

INPUT VARIABLE

N ORDER OF THE DESIRED HADAMARD MATRIX,

OUTPUT VARIABLE —

H HADAMARD MATRIX OF ORDER N.

INTEGER M2(2»2>i N4(4»4) t N8<8»8)» N16(16.16)

REAL H(16i 16)

IN THE FOLLOWING DATA STATEMENTt NN IS THE PRINTER DEVICE NUMBER.

DATA NW / 6 /

DEFINE N2. AN HADAMARD MATRIX OF ORDER 2.

N2( It 1 ) = 1 N2(1.2) = 1 N2<2.1) = 1 N2(2«2) = -1

IF (N.EQ.2) GO TO 200

DEFINE N4. AN HADAMARD MATRIX OF ORDER 4.

DO 20 1=1,4

DO 10 J=1.4

N4(I,J) = 1 10 CONTINUE

20 CONTINUE

N4<2,2> = -1 N4(2.4) = -1 N4<3,3) = -1 N4(3,4) = -1 N4(4,2) = -]

N4(4.3) = -1 IF (N.EQ.4) GO TO 400

COMPUTE N8 AS THE TENSOR PRODUCT OF N2 ANO N4.

DO 40 I = l»4

DO 30 J = 1.4

N8(I»J) = N4(I.J)

N8CI+4.J) = N4(I«J)

N8(I,J+4) = N4(ItJ)

N8U+4.J+4) = -N4(I.J)

30 CONTINUE

40 CONTINUE

IF (N.EQ.8) GO TO 800

COMPUTE N16 AS THE TENSOR PRODUCT OF N4 AND N4.

GEN02700

GEN02 800

HAD 0010

HAD 0020

• HAD 0030

• HAD 0040

HAD 0050

HAD 0060

HAD 0070

HAD 0080

HAD 0090

HAD 0100

HAD 0110

HAD 0120

HAD 0130

HAD 0140

HAD 0150

HAD 0160

HAD 0170

HAD 0180

HAD 0190

HAD 0200

HAD 0210

HAD 0220

HAD 0230

HAD 0240

HAD 0250

HAD 0260

HAD 0270

HAD 0280

HAD 0290

HAD 0300

HAD 0310

HAD 0320

HAD 0330

HAD 0340

HAD 0350

HAD 0360

HAD 0370

HAD 0380

HAD 0390

HAD 0400

HAD 0410

HAD 0420

HAD 0430

HAD 0440

HAD 0450

HAD 0460

HAD 0470

HAD 0480

HAO 0490

HAD 0500

HAD 0510

HAD 52

HAD 0530

HAD 0540

53
DO 60 1=1,4

DO 50 J=l e 4

N16(I« J) = N4(I • J)

N16( I« J+4) = N4( I .J)

N16(I» J+8) = N4( I >J)

N16( I,J+12) = N4( I iJ)

N16(I+4, J+4) = N4(I .J)

N16( 1+4, J+8) = -N4( I • J)

N16( I+4*J+12) = -N4( I .J)

N16( I+8« J+8) = N4( I,.J)

N16( I+8« J+12) = -N4( I .J)

N16( 1+12, J+12) = N4( I >J)

50 CONTINUE

60 CONTINUE

DO 80 1=1,16

DO 70 J=I«16

N16C J, I ) = N16(I, J)

70 CONTINUE

80 CONTINUE

IF (N.EQ.16) GO TO 1600

WRITE (NW, 90000)

N

RETURN

200 DO 220 I=1»N

DO 210 J=1.N

H(I,J) =N2(I,J)

210 CONTINUE

220 CONTINUE

RETURN

400 DO 420 1=1,

DO 410 J=1

N H(I,J) = N4( I , J)

410 CONTINUE

420 CONTINUE

RETURN

800 DO 820 1=1,

DO 810 J=ltN

H(I,J) = N8(I,J)

810 CONTINUE

820 CONTINUE

RETURN

1600 DO 1620 1=1 ,N

DO 1610 J=1,N

H(I, J) = N16( I , J)

1610 CONTINUE

1620 CONTINUE

RETURN

C

90000 FORMAT (52H0*** * 15/)

ERROR

—

N MUST EQUAL 2, 4« 8 OR 16.

END

HAD 0550

HAD 0560

HAD 0570

HAD 0580

HAD 0590

HAD 0600

HAD 610

HAD 0620

HAD 0630

HAD 0640

HAD 0650

HAD 66

HAD 0670

HAD 0680

HAD 69

HAD 0700

HAD 0710

HAD 0720

HAD 0730

HAD 0740

HAD 0750

HAD 0760

HAD 0770

HAD 0780

HAD 0790

HAD 0800

HAD 0810

HAD 0820

HAD 0830

HAD 0840

HAD 0850

HAD 0860

HAD 0870

HAD 0880

HAD 0890

HAD 0900

HAD 0910

HAD 0920

HAD 0930

HAD 0940

HAD 0950

HAD 0960

HAD 0970

HAD 0980

HAD 0990

HAD 1000

HAD 1010

HAD 1020

HAD 1030

HAD 1040

HAD 1050

HAD 1060

HERE

N : = c HAD 1070

HAD 1080

HAD 1090

HAD 1100

54
SUBROUTINE SCALE ( ISCALE. NCALL, M, N, Li * R. RES* SF. U. X* IFAULT)

W. MM, MMNN. NN. A. B, SCL

SCL

SCL

SUBROUTINE SCALE SCALES THE DATA (A AND B) FOR LEAST SQUARES PROBLEMS SCL

IN ORDER TO MITIGATEt IF POSSIBLE* ROUNDING ERROR PROBLEMS WHICH CAN SCL

OCCUR IN CONNECTION WITH SOLVING ILL-CONDITIONED SYSTEMS OF EQUATIONS. SCL

SCL

IN THE CASE OF POLYNOMIAL TYPE PROBLEMS. THE MEAN OF THE INDEPENDENT SCL

VARIABLE (VECTOR A(I,2>) IS COMPUTED SO THAT IT CAN BE SUBTRACTED FROMSCL

EACH ELEMENT OF A(I,2> BEFORE POWERS OF THIS VECTOR ARE GENERATED. SCL

SCL

AFTER A SOLUTION IS OBTAINED FOR A SCALED PROBLEM. THE COEFFICIENTS SCL

(X)« RESIDUALS (RES). AND COVARIANCE MATRIX (R) MUST BE ADJUSTED TO SCL

ACCOUNT FOR SCALING. SCL

SCL

INPUT VARIABLES — SCL

ISCALE PARAMETER WHICH INDICATES THE TYPE OF SCALING TO BE DONE. SCL

ISCALE = I MEANS DATA MATRIX A IS TO BE SCALED, BUT MATRIX B SCL

IS NOT TO BE SCALED. SCALE FACTORS WILL BE AUTOMATICALLY SCL

COMPUTED FROM THE DATA. THE J-TH SCALE FACTOR. SF(J), IS SCL

BASED ON THE EUCLIDEAN NORM OF THE J-TH COLUMN OF A. SCL

ISCALE = 2 MEANS DATA MATRICES A AND B ARE TO BE SCALED. SCL

MATRIX A IS SCALED IN THE SAME MANNER AS WHEN ISCALE =1. SCL

AND MATRIX B IS SCALED IN A SIMILAR MANNER. SCL

ISCALE = 3 MEANS DATA MATRICES A AND B ARE TO BE SCALED. THE SCL

SCALE FACTORS ARE FURNISHED BY THE USER AND ARE ENTERED SCL

AS INPUT DATA IN THE CALLING PROGRAM. SCL

NCALL PARAMETER INDICATING WHICH OPERATIONS ARE TO BE PERFORMED. SCL

NCALL = 1 MEANS THAT THE MEAN OF THE VECTOR A(I»2) IS TO BE SCL

COMPUTED IN POLYNOMIAL TYPE PROBLEMS. SCL

NCALL = 2 MEANS — SCL

(A) SCALE FACTORS ARE TO BE COMPUTED IF ISCALE = 1 OR 2 . SCL

(B) MATRIX A IS TO BE SCALED. SCL

(C) MATRIX B IS TO BE SCALED IF ISCALE = 2 OR 3

.

SCL

NCALL = 3 MEANS THAT COEFFICIENTS (X). RESIDUALS (RES). AND SCL

COVARIANCE MATRIX (R) ARE TO BE ADJUSTED TO ACCOUNT FOR SCL

SCALING. THE MATRIX OF OBSERVATIONS (Bl WILL ALSO BE SCL

ADJUSTED SO THAT IT IS RESTORED TO ITS ORIGINAL FORM (BEFORESCL

MN

L

WMM

MMNN

NN

SCALING).

TOTAL NUMBER OF EQUATIONS.

NUMBER OF UNKNOWN COEFFICIENTS.

NUMBER OF VECTORS OF OBSERVATIONS.

VECTOR OF WEIGHTS (OF LENGTH M).

DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN

ARRAYS A, B AND RES.

DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN

ARRAY R.

DIMENSIONING PARAMETER SPECIFYING MAXIMUM NUMBER OF ROWS IN

ARRAY X.

INPUT AND OUTPUT VARIABLES —

B

RRES

SF

C U

TWO-DIMENSIONAL ARRAY OF SIZE (MM.N). SEE MAIN PROGRAM FOR

FURTHER DETAILS.

TWO-DIMENSIONAL ARRAY OF SIZE (MM.L). SEE MAIN PROGRAM.

TWO-DIMENSIONAL ARRAY OF SIZE (MMNN.N). SEE MAIN PROGRAM.

TWO-DIMENSIONAL ARRAY OF SIZE (MM.L). SEE MAIN PROGRAM.

VECTOR OF SCALE FACTORS (OF LENGTH N + L).

IN POLYNOMIAL TYPE PROBLEMS U IS THE MEAN OF THE VECTOR

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

SCL

0010

0020

0030

0040

0050

0060

0070

0080

0090

0100

0110

0120

0130

0140

0150

0160

0170

0180

0190

0200

0210

0220

0230

0240

0250

260

0270

0280

0290

0300

0310

0320

0330

340

0350

360

0370

380

390

0400

0410

0420

0430

0440

0450

0460

0470

0480

0490

0500

510

0520

0530

0540

550

0560

0570

580

55
C A(I»2). IN NON-POLYNOMIAL TYPE PROBLEMS. U = 0. SCL 0590

C X TWO-DIMENSICNAL ARRAY OF SIZE (NN.L). SEE MAIN PROGRAM. SCL 0600

C SCL 0610

C OUTPUT VARIABLE SCL 0620

C IFAULT FAULT INDICATOR WHICH ON EXIT EQUALS EITHER OR 4. SCL 0630

C IFAULT = 4 MEANS THAT AT LEAST ONE WEIGHT WAS FOUND TO BE SCL 0640

C NEGATIVE. (ALL WEIGHTS SHOULD BE NONNEGATI VE.

SCL 0650

C IFAULT = MEANS ALL WEIGHTS ARE NONNEGATI VE. SCL 0660

C SCL 0670

REAL A(MM.N). B(MM,L). R(MMNN,N). RES(MM,L)t SF ( 1 ) , U, W(M) SCL 0680

REAL X(NN.L) SCL 0690

DOUBLE PRECISION SUM SCL 0700

C SCL 0710

IF (NCALL.EQ.2) GO TO 20 SCL 0720

IF (NCALL.EQ.3) GO TO 120 SCL 0730

C SCL 0740

C NCALL = 1. SCL 0750

C SCL 0760

SUM = 0.0 SCL 0770

DO 10 1=1. M SCL 0780

SUM = SUM + DBLE(A(I.2)) SCL 0790

10 CONTINUE SCL 0800

U = SUM SCL 0810

U = U / FLOAT(M) SCL 0820

RETURN SCL 0830

C SCL 0840

C NCALL = 2. SCL 0850

C SCL 860

20 IFAULT = SCL 0870

DO 70 J = 1»N SCL 0880

IF ( ISCALE.EQ.3) GO TO 50 SCL 0890

SUM =0.0 SCL 0900

SF(J) =1.0 SCL 0910

DO 40 1=1.SCL 0920

IF (W( I ).GE.0.0) GO TO 30 SCL 0930

IFAULT = 4 SCL 0940

CALL ERROR (IFAULT. I, W( I ) ) SCL 0950

30 SUM = SUM + DBLE( A(I. J))*DBLE( A(I. J))*DBLE(W( I >) SCL 0960

40 CONTINUE SCL 0970

IF ( IFAULT. EQ. 4) RETURN SCL 0980

IF (SUM. EQ. 0.0) GO TO 70 SCL 0990

SF(J) = 1.0/DSQRT(SUM) SCL 1000

50 DO 60 1=1.

SCL 1010

A(I,J) = A(I.J)*SF(J) SCL 1020

60 CONTINUE SCL 1030

70 CONTINUE SCL 1040

C SCL 1050

DO 110 K=1«L SCL 1060

NPK = N + K SCL 1070

IF (ISCALE.EQ.3) GO TO 90 SCL 1080

SUM = 0.0 SCL 1090

SF(NPK) = 1.0 SCL 1100

IF (ISCALE.EQ.l ) GO TO 110 SCL 1110

DO 80 I=1«M SCL 1120

SUM = SUM + DBLE(B( I.K) )*DBLE(B( I.K))*DBLE( W( I )) SCL 1130

80 CONTINUE SCL 1140

IF (SUM. EQ. 0.0) GO TO 110 SCL 1150

SF(NPK) = 1.0/DSQRT(SUM) SCL 1160

56
90 DO 100 I=liM SCL 1170

B(I.K) = B( I,K)*SF(NPK) SCL 1180

100 CONTINUE SCL 1190

110 CONTINUE SCL 1200

RETURN SCL 1210

C SCL 1220

C NCALL = 3. SCL 1230

C SCL 1240

120 DO 170 K=1#L SCL 1250

NPK = N + K SCL 1260

DO 130 J=1.N SCL 1270

X(J#K) = (X( J,K)*SF( J ) ) / SF(NPK) SCL 1280

130 CONTINUE SCL 1290

DO 140 1=1 .

M

SCL 1300

RES(I.K) = RES(ItK) / SF(NPK) SCL 1310

B(I.K) = B(I.K) / SF(NPK) SCL 1320

140 CONTINUE SCL 1330

IF (K.GT.l) GO TO 170 SCL 1340

DO 160 1=1.

N

SCL 1350

DO 150 J=ltl SCL 1360

R(I.J) = Rdi J)*SF( I )*SF( J) SCL 1370

150 CONTINUE SCL 1380

160 CONTINUE SCL 1390

170 CONTINUE SCL 1400

IF (U.EQ.0.0) RETURN SCL 1410

C SCL 1420

C IN SCALED POLYNOMIAL PROBLEMS. ADJUST COEFFICIENTS (X) AND UNSCALED SCL 1430

C COVARIANCE MATRIX FOR SUBTRACTION OF MEAN (U). SCL 1440

C REFERENCE — SCL 1450

C G. A. F. SEBER. LINEAR REGRESSION ANALYSIS (1977), THEOREM SCL 1460

C 1.4 AND COROLLARIES. PAGES 10-11. SCL 1470

C SCL 1480

C SCL 1 490

C FILL OUT THE SYMMETRIC COVARIANCE MATRIX SO THAT ALL ELEMENTS ARE SCL 1500

C PRESENT. SCL 1510

C SCL 1520

DO 190 1=1.

N

SCL 1530

DO 180 J = I«N SCL 1540

R(I.J) = R(J.I) SCL 1550

180 CONTINUE SCL 1560

190 CONTINUE SCL 1570

DO 260 1=1,

N

SCL 1580

SF( I ) = 1.0 SCL 1590

I PI = I + 1 SCL 160

IF (IP1.GT.N) GO TO 210 SCL 1610

DO 200 J=IP1.N SCL 1620

SF(J) = -DBLE(FLOAT(J-l ))/DBLE( FLOAT( J-I ) ) * DBLE ( SF ( J-l )

)

SCL 163

* * DBLE(U) SCL 1640

200 CONTINUE SCL 1650

210 DO 230 K=1«L SCL 1660

SUM = 0.0 SCL 1670

DO 220 J=I.N SCL 1680

SUM = SUM + DBLE(SF(J)> * DBLE(X(J,K)) SCL 1690

220 CONTINUE SCL 1700

X< I.K) = SUM SCL 171

230 CONTINUE SCL 1720

DO 250 J=I.N SCL 1730

SUM = 0.0 SCL 1740

57
DO 240 K=I.N

SUM = SUM OBLE(SF(K))*DBLE«RCK.J)) sc. \ yln 240 CONTINUE *„ lit

R<I ' J > ' SUM Icl \7al 250 CONTINUE ~ l L

&0

260 CONTINUE Z.ZT 1790

c SCL 1800

SF(J) = 1.0 __. IP1 = J 1 SCL 183 ° IF (XP1.GT.NJ GO TO 280 J?" J®*®

00 270 I = IP1.N ll^ \H° SF(I) = -DBLE<FLOAT(I-l>)/DBLE<FLOATC I-J>> * DBLE < SF< I-l ) > SCL 1870

270 CONTINUE f^ JJfJ

280 DO 300 I-l. J "J 00«

00 290 K=J.N SCL l91 ° " ** •»•" SCL 1920 SUM = SUM OBLE(SF(KJ|*DBLE(R<I.KI) SCL 1 <530 290 CONTINUE fl, J 9

f°

300 CONTINUE rij" J*J°

310 CONTINUE ??" 196°

END !J

L 1980

SCL 1990

APPENDIX B. LISTING OF DATA

DATA1A

DATA1B

DATA3A

DATA3B

DATA3C

DATA3D

DATA3E

DATA4

DATA5A

DATA5B

DATA5C

DATA5D

DATA5E

58
***** 0ATA1A *****

10 2

10 3

10 4

10 6

10 8

10 8 4

20 2 20 5

20 10

20 15

20 15 5 30 5

30 10

30 15

30 20

30 20 10

40 5 40 10

40 20

40 30

40 30 10

50 10

50 20

50 30

50 40

50 40 10

100 10

100 20

100 40

100 60

100 80

100 80 10

150 10

150 20

150 40

150 60

150 60

150 80 10

200 10

200 20

200 40

200 60

200 80

200 80 10

'*** DATA1B *****

16

10 8

10 8 4

20 15

20 15 5

30 20

59
30 20 10

40 30

40 30 10

50 40

50 40 10

100 80

100 80 10

150 80

150 80 10

200 80

200 80 10

***** TEST.DATA3A *****

1 4

5 4 1 2 4

2 9 7 4

1 2 4 1 6 2 6 5 6 2 3

5 10 9 14 14 5 6 3 2 3 12

6 25 19 27 28

***** TEST.DATA3B *****

2 3

3 4 2 5 4

3 1 1 1 6 3 5 3 4 1 2 8 5 8 7

4 1 6 12 4 10 9 6 11 1 5 16

29 26 3 4 6 15 8 32

***** TEST.0ATA3C *****

3

2 1 5 3 3

8 10 5 8

8 7 5 8 7 7 4 3 1 4

11 8 12 13 14 10 4 16 6 8 12 1 1 32 16 4 7 24 24 31 15

***** TEST.DATA3D *****

4 2

1 4 1 1 10 7 2 10

6 3 7 2 2 2 3 4 6 1 7

14 13 9 12 3 12 1 11 6 7 6 1 1 15 28 16 12 14 9 16

***** TEST.DATA3E *****

2 1 12 5

6 11 7 5

60
745521785435

16 3 15 187390 12 84

18 30 29 22 14 21 2

***** TEST.DATA4 *****

0.• 100000. 759.

1.. 111111. -2048.

2.> 124992. 2048.

3<• 142753. -2048.

4<> 165984. 2523.

5.• 196875. -2048.

6.> 238336. 20 48.

7,i 294117. -2048.

8,> 368928. 1838.

9.> 468559. -2048.

10.. 600000. 20 48.

11., 771561. -2048.

12., 992992. 1838.

13.. 1275603. -2048.

14.> 1632384. 2048.

15.> 2078125. -2048.

16.. 2629536. 2523.

17.» 330 5367. -2048.

18.> 4126528. 2048.

19., 5116209. -2048.

20.> 6300000* 759.

***** TEST.DATA5A *****

(1) WAMPLERt

21 6

(F2.0.2F8.0)

J. A MER. STAT. ASSN. 1970*

1 ,549. 5TH

0.

OEG. POLYNOMIALS. EQUAL WEIGHTS.

1 760

1 6 -2042

2 63 2111

3 364 -1684

4 1365 3686

5 3906 1858

6 9331 11379

7 19608 17560

8 37449 39287

9 66430 64382

10 111 1 11 113159

11 177156 175108

12 271453 273291

13 402234 40 0186

14 579195 581243

15 813616 811568

16 1118481 1 121004

17 1503598 1506550

18 2000719 2002767

19 2613660 2611612

20 3368421

1FIRST C

3369180

(2) )EGREE POLYNOMIAL. UNEQUAL WEIGHTS.

61
6 2 112 1 0.

(3F3.0)

1. 2. 2.

2. 2. 1.

3. 5. 1.

4. 4. 1.

5. 7. 1.

6. 7. 2.

1 (3) DATA WITH UNEQUAL WEIGHTS. TWO COLUMNS LINEARLY DEPENDENT,

7 6 2 2 2 1 0.

<2F3.0 e F3.1.F3.0tF4.2tF3.0.F5.1tF4.0tF2.0)

5 2 .25 2 13.0 130 2

2

1

1

1

1

1

1

1

(4)

123

2

1

1

1EXA

12

(8F3.0)

11

1 -1

2

3

1 -3

1

,25

,00

,00

,00

,00

,00

170

182

88

17,0

18.2

8.8

-3.0 -30

2.8 28

2.1 21

MPLE WITH WEIGHTS AND CONSTRAINTS.

6 3 12 2 1

(5)

11

-1

-1

1

1

1

-1

-1

1

-1

-1

1

-1

-1

11

1

14401211 INVERSE OF HILBERT MATRIX OF ORDER 4. M = 4» N = 4. Ml - . 0.

(5F7.0)

16.

-120.

240.

-140.

1 (6) INVERSE OF HILBERT MATRIX OF ORDER 4,

4

(5F7.(

16.

-120.

240.

-140.

1 (7) BUSINGER-GOLUBt 6512211 NUM. MATH. 1965o P.269t INVERSE OF HILBERT MATRIXt ORDER 6.

0.

(5F10.0,

-120. 240. -140.

1200. -2700. 1680.

-2700. 6480. -4200.

1680. -4200. 2800.

0)

-120« -140.

1680.

240.

1200. -2700.

-2700. 6480. -4200.

1680. -4200. 2800.

-4.

60.

-180.

140.

1 -4.

60.

-180.

140.

M = 4, N = 4* Ml

0.

= 4,

-6

33

5 12 • 10X. 2F10.0)

36. -630.

30. 14700.

60. -88200.

3360.

-88200.

-7560.

21 1680.

7560.

-220500.

564480. -1411200. 1512000.

463. 463

13860. -17820

97020. 93555

62
-7560, 211680

7560. -2 20 500

-2772. 83160

1 (8) EXAMPLE WITH X

1 1 (2F3.0)

1. 0.

1 (9) ALBERT* REGRESS

3 4

(7F4.0)

1 • 1 • -1. 0.

0. 1 •

-1411200. 3628

1512000. -3969

-582120. 1552

= (HENCE XNORM 12 10

800. -3969000.

000. 4410000.

320. -1746360.

-258720.

291060.

-1 16424.

-261800.

286288.

-118944.

0). TOL = -1.

-1 ON ENTRY TO L2A OR L2B«

ION AND THE MOORE

3 2 12 -PENROSE INVERSE* 1972*

0.

P. 63.

1.

0. 1 1. 1

1 (10) FIFTH DEGREE P 216

(F8.0,3F9.0. 2F10.0*

1000000. 0.

1000000. 100000.

1000000. 200000.

1000000.

1000000.

1000000.

1000000.

1000000.

1000000.

1000000.

300000.

400000.

500000.

600000.

700000.

800000.

900000.

1000000. 1000000.

1000000. 1100000.

1000000. 1200000.

1000000. 1300000.

1000000. 1400000.

1000000. 1500000.

1000000. 1600000.

1000000. 1700000.

1000000. 1800000.

1000000. 1900000.

1000000. 2000000.

1 (11) LAWSON-HANSON*

8 6 4

(7F6.0)

155. 105. -445.

355. 305. -245.

-445. -495. -45.

-245. -295. 155.

-45. -95. 355.

155. 105. -445.

355. 305. -245.

-445. -495. -45.

1 (12) LAWSON-HANSON*

8 6 4

(7F6.0)

155. 105. -445.

355. 305. -245.

• 0.

• 1.

• 0. OLYNOMI

1 2

F13. 2*F

10000

40000

90000

160000

250000

360000

490000

640000

810000

000000

210000

440000

690000

960000

250000

560000

890000

240000

610000

000000

0.

0.

1.

AL WIT

2 10.0)

10

80

270

640

1250

2160

3430

5120

7290

10000

13310

17280

21970

27440

33750

40960

49130

58320

68590

80000

H HE

1

0.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

00.

AVY WEIGHTS* MATRIX A SCALED.

0.

IFAULT=11

0. 0. 100000.1853 16777216

100. 1. -8277497.00

1600. 32. 8513600.00

8100. 243. -8245855.00

25600. 1024. 10500192.00

62500. 3125. -8191733.00

129600. 7776. 8626944.00

240100. 16607. -8094491.00

409600. 32768. 7897376.00

656100. 59049. -7920049.00

1000000. 100000. 600000.50 16777216

1464100* 161051. -7617047.00

2073600. 248832. 8521440.00

2856100. 371293. -7113005.00

3841600. 537824. 10020992.00

5062500. 759375. -6310483.00

6553600. 1048576. 12963744.00

8352100. 1419857. -5083241.00

10497600. 1889568. 12515136.00

13032100. 2476099. -3272399.00

16000000. 3200000. 6300000.1853 16777216

SOLVING LEAST SO. 12 11 UARES PROBLEMS. 1974, SET I EX.16.

0.

-495. -45. -95. -245.

-295. 155. 105. -295.

-95. 355. 305. 155.

105. -445. -495. 105.

305. -245. -295. -445.

-495. -45. -95. -495.

-295. 155. 105. -45.

-95. 355. 305. -95.

IFAULT=10

SOLVING LEAST SQ 12 12

-495. -45. -95.

-295. 155. 105.

HARES PROBLEMS* P. 252. SET 1* EX.16* T0L=.5.

0.5

-245.

-295.

63
-445. -495.

-245. -295.

-45. -95.

155. 105.

355. 305.

-445. -495.

1 (13) BJORCK8 6

(6F12.0)

20160.

945.

-952560.

-40320.

1 1430720.

456120.

-58212000.

-2236080.

149688000.

5599440.

-20 4324120.

-7495488.

141261120.

5105100.

-38918880.

-13899 60.

1 (14) LAWSON6 8

(9F6.0)

-245. -295.

355. 305.

-45. -95.

-45,

155,

-95,

105,

355,

-445,

305.

-495.

355. 305. -245. -295.

-445. -495. -45. -95.

-245. -295,

-45. -95,

155,

355,

105.

305.

GOLUB* BIT 1967* P. 322.

2 3 2 11

-92400.

945.

4656960.

-40320.

-58212000.

3256120.

304920000.

-136080.

-800415000.

7279440.

1109908800.-

-6095488.

-776936160.

6305100.

21621 £000.

-339960.

221760

8400945

-11642400

4159680

149688000

3256120

-800415000

-136080

2134 440000

7279440

2996753760

-6095488

2118916800

6305100

-594594000

-339960

155.

105.

-445.

-495.

-45.

-95.

HILBERT MATRIX INVERSE. ORDER 8. IFAULT=8,

0.

-288288. 192192. -51480.

15567552. -10594584. 2882880.

-204324120. 141261120. -38918880.

1109908800. -776936160. 216216000.

-2996753760. 2118916800. -594594000.

4249 941696.-30 30051024. 8 56215360.

-3030051024. 2175421248. -618377760.

856215360. -618377760. 176679360.

-445. -495.

155. 105.

-245. -295.

1 (15) EXAMPLE

6 3

(4F2.0) 1111

2 2 2 1 10 1 12 4 1 13 9 1 14 9 1

HANSONt SOLVING LEAST SQ

6 12 11

155. 105. -445. -495.

-245. -295. 155. 105.

355. 305. -245. -295.

-45. -95. 355. 305.

-445. -495. -45. -95.

155. 105. -445. -495.

UARES PROBLEMSt 1974. SET 1. EX.12. IFAULT=7

0.

-45,

-445.

155.

-245,

355,

-45,

-95.

-495.

105.

-295.

305.

-95.

355.

305.

-245.

-295.

155.

105.

1 .

11 , 4,

9,

16,

WITH SINGULAR MATRIX OF

3 12 11 CONSTRAINTS.

0.

Ml = 3t Nl = 2, IFAULT=6

(16) EXAMPLE

3 2 (3F2.0)

1

1

1 1 (17) EXAMPLE

2 1 (2F3.0.F4.0)

1.1. 0.

WITH MATRIX A EQUAL TO 12 11 ZERO (HENCE RANK EQUALS ZERO).

0.

IFAULT=5

WITH ZERO AND NEGATIVE 12 2 1 WEIGHTS.

0,

IFAULT=4

64
1. 1. -1.

1 (18) EXAMPLE WHERE N EXCEEDS THE CORRESPONDING DIMENSION LIMIT.

1 21 I 2 1 2 0.

(22F2.0) 123456789012345678901 1 1 (19) 112 EXAMPLE WHERE Ml EXCEEDS M AND N. 12 11 0.

(2F3.0)

1. 1.

1 (20) EXAMPLE WITH M = 0. 10 12 11 0.

(2F3.0)

1* 1* 1 •

IFAULT=3

IFAULT=2

IFAULT=1

***** TEST.OATA5B *****

(21 ) FI

8

(3F3.0)

1* 2*

RST DEGREE POLYNOMIAL. POSITIVE AND ZERO WEIGHTS.

2 112 0.

COMPARE EXAMPLE (2).

2.

5.

4.

7.

7.

9.

6.

1 (22) EX

14

(8F3.0)

11

1 -1

1

1

1

1

2.

1.

1.

1.

1.

2.

0.

0.

AMPLE WITH WEIGHTS AND CONSTRAINTS. COMPARE EXAMPLE (4).

6 3 12 2 0.

1

1

-1

-1

11

1

-1

-1

-1

-1

-1

-1

1

-1

-1

o

(23) EX

4

(4F3.0)

1.0.

0. 1.

0. 0.

0. 0.

AMPLE WITH ZERO WEIGHTS. WHERE RANK A = RANK H

2 12 2 0.

1. 1. 2. 1. 3. 0. 4. 0,

= 2, H = (SQRT(W) )*A.

65
1 (24) EXAMPLE WITH ZERO WEIGHTS* WHERE RANK A = 2* RANK H = 1. H = (SQRT(W))*A.

4 2 12 2 0.

(4F3.0)

1. 0. 1* 1.

0. 1. 2. 0.

0. 0* 3. 0*

0* 0. 4* 1*

***** TEST.DATA5C *****

(3) DATA 7602221 WITH UNEQUAL WEIGHTS* TWO COLUMNS LINEARLY DEPENDENT.

0.

(2F3.0*F3.1*F3.0*F4.2*F3.0*F5.1*F4.0*F2.0)

1 1 • 5 2 •25 2 13.0 130 2

1 2 • 5 2 •25 3 17.0 170 2

3 • 3 •00 3 18.2 182 1 2 • 1 •00 1 8.8 88 1 1 • -3 •00 -3.0 -30 1 1 • •00 2.8 28 1

1 (9) ALB

• 1 •00 2.1 21 1

ERT* REGRESSION AND THE MOORE--PENROSE INVERSE* 1972* P.

3 4 3 2 1 2 0.

(7F4.0)

1. • 1 • 1. 1 • 0* 0.

0. 1 • -1 • 0* • 1 • 0*

1. 1

1 (12) LA

• 0. 1. 0. 0. 1.

WSON-HANSON* SOLVING LEAST SQUARES PROBLEMS* P. 252* SET

8 6 4 12 1 2 0.5

(7F6.0)

155. 105. -445. -495. -45. -95. -245.

355. 3 05. -245. -295. 155. 105. -295.

-445. -495. -45. -95. 355. 305. 155.

-245. -295. 155. 105. -445. -495. 105.

-45. -95. 355. 305. -245. -295. -445.

155. 105. -445. -495. -45. -95. -495.

355. 305. -245. -295. 155. 105. -45.

-445.

1 (14) LA

-495. -45. -95. 355. 305. -95.

WSCN-HANSON, SOLVING LEAST SQUARES PROBLEMS* 1974* SET 1 6 8 6 12 1 1 0.

(9F6.0)

-245. -295. 155. 105. -445. -495. -45. -95. 355. 1.

355. 305. -245. -295. 155. 105 . -445. -495. 305. 1.

-45. -95. 355. 305. -245. -295. 155. 105. -245. 1.

-445. -495. -45. -95. 355. 305. -245. -295. -295. 4.

155. 105. -445. -495. -45. -95. 355. 305. 155. 9.

-245.

1 (18) EX

-295. 155. 105. -445. -495. -45. -95. 105. 16.

AMPLE WHERE N EXCEEDS THE CORRESPONDING DIMENSION LIMIT.

1 21 12 1 2 0.

(22F2.0)

12 3 4 5 6 789012345 6 7 8 9 11

63.

1* EX.16* T0L=.5.

* EX.12. IFAULT=7

IFAULT=3

66
***** TEST.DATA5D *****

(25) BENNETT AND FRANKLIN. 1954i PP. 379-385. MISSING DATA.

20 8 1I 2 2 0.

(10F4.0)

1. 1. 0. 0. 0.. 1. 0. 0., 28.

1* 1 • 0. 0. 0.> 0. 1. 0.. 20.

1. 1. 0. 0. 0,> 0. 0. 1.. 11 . 1. 1 . 0. 0. 0.. -1. -1. -1.. 10.

1. 0. 1. 0. 0,. 1. 0. 0.> 0. 0.

1. 0. 1 • 0. 0.. 0. 1. 0.> 16.

1. 0. 10. 0.> 0. 0. 1.> 15.

1. 0. 1 • 0. 0. -1. -1. -1 . 8.

1 • • • 1 . 0., 1 • 0. 0.> 29.

1. 0. 0. 1. 0.. 0. 1. 0.. 13.

1. 0. 0. 1 • 0.> 0. 0. ll> 16.

1. 0. 0. 1 . 0., -1. -1. -1.. 0. 0.

1. 0. 0. 0. 1.i 1 • 0. 0.> 27.

1. 0. 0. 0. 1 . » 0. 1. 0.> 10.

1. 0. 0. 0. 1 . i 0. 0. 1.> 18.

1. 0. 0. 0. 1 , . -1. -1. -1.. 11.

1 • -1 . -1

.

-1 • -1 i 1. 0. 0.> 28.

1. -1 . -1

.

-1 • -1 . 0. 1. 0.. 11.

1 . -1 . -1

.

-1. -1.• 0. 0. ll. 15.

1. -1 . -1.

1 (26) NBS TEO

-1 . -1,. -1. -1. -I.. 10.

1NICAL CJOTE 552. TWOWAY ANOVA. PAGES 134-135 — WITH ZERO WEIGHTS

12 6 11 2 2 0.

(8F3.0)

1 1 1 8 1 110 1 2 110 1 1 110-1 -1 -1 3 10 11 4 10 1 1 10 10 1 3 1 10 1-1 -1 -1 1 1

1 -1 -I 1 2

1 -I -1 1 1 1 1-1-1 1 1

1 -1 -1 -1

1 -1 -1 4 1

1 (27) 6TH DEGREE POMrNOMIAL FITTEiD TO 5TH DEGREE DATA — SCALING IS NOT USED.

21 7 o ;> 1 1 0.

(F2.0.2F8.0)

1 760

1 6 -2042

2 63 21 11

3 364 -1684

4 1365 3888

5 3906 1858

6 93 31 1 1379

7 19608 17560

8 37449 39287

9 66430 64382

10 111111 113159

67
11 177156 175108

12 271453 273291

13 402234 400186

14 579195 581243

15 813616 811568

16 1 118481 1 121004

17 1508598 1506550

18 2000719 2002767

19 2613660 2611612

20 3368421 3369180

1 (28) WAMPLERt J. AMER. STAT. ASSN. 1970t P.549 f 5TH OEG. POLYNOMIALt EQUAL WEIGHTS.

21 6 1 1 1 1 0.

(F3.0. 8X.F8.0)

100 1 760

101 6 -2042

102 63 2111

103 364 -1664

104 1365 3888

105 3906 1858

106 9331 11379

107 19608 17560

108 37449 39287

109 66430 64382

110 111111 113159

111 177156 175108

112 271453 273291

113 402234 400186

114 579195 581243

115 813616 811568

116 1118481 1121004

117 1508598 1506550

118 2000719 2002767

119 2613660 2611612

120 3368421 3369180

***** TEST.DATA5E *****

(27) 6TH OEG. POLYNOMIAL FITTED TO 5TH DEG. DATA. SCALING USEO ( I SCALE = 3).

21 7 2 1 1 0.

(F2.0.2F8.0)

1 760

1 6 -2042

2 63 21 11

3 364 -1684

4 1365 3886

5 3906 1858

6 9331 11379

7 19608 17560

8 37449 39287

9 66430 64382

10 1111 11 113159

11 177156 175108

12 271453 273291

13 402234 400186

14 57919 5 581243

68
15 813616 811568

16 1118481 1121004

17 1508598 1506550

18 2000719 2002767

19 2613660 2611612

20 3368421 3369180

1000000.0 100000.00 10000.000 1000.0000 100.00000 10.000000 1.0000000 1.0000000

1.0000000

APPENDIX C. SAMPLE OUTPUT FROM FIVE PROBLEMS USING TEST5

Fig. 2. Example 9 from DATA5C, MODE = 1, ISCALE = 0.

Fig. 3. Example 9 from DATA5C, MODE = 2, ISCALE = 0.

Fig. A. Example 22 from DATA5B, MODE = 1, ISCALE = 0.

Fig. 5. Example 28 from DATA5D, MODE = 1, ISCALE = 0.

Fig. 6. Example 28 from DATA5D, MODE = 1, ISCALE = 2.

The five examples given in this appendix illustrate various

features of the L2A and L2B algorithms. Input parameters and input

data (matrix A, matrix B, and vector of weights) are printed. If

data are to be scaled before solutions are computed, scale factors

are printed. Computed results include the computed rank of the

system of equations; the sequence in which columns were selected by

the pivoting scheme; information on the behavior of the iterative

refinement procedure; coefficients and their standard deviations;

observations, predicted values and residuals; sum of squared residuals,

norm of residuals, and residual standard deviation; and unsealed

covariance matrix.

69
oooooooo o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o • • • P4 o o o o o o o o o o o o o o o o o o o o o o • o o o o o o o o o o o o o o o o

a.LUQao

13z

1-1 3o

_l

-Jo

LL

UJI

UJ

a o o o UJ o o o I • UJ O • o o u

cm -1 *4 • • (0

r» < a u o «* CO z4 o 1- HI o o o o co o o o > a o o o H

UJ UJ <-t o o o a > Q o o o z o o o o CM UJ tM Z o• oo o• I

1- UJ —4 • w* CO II V o CD a 10 z CO z Q

UJ

K *4 1- o UJ

aX

X o o o 1-4 1- i CJ> 13 o o o H o UJ 3 1-4 o o o < UJ a M UJ o o o 3 -1 Q 3B o o o o UJ o o o o UJ CO Z IL o o o O • • o u. UJ

UJ 1-4 <-4 • o Cf X UJ CM a 1 LU H a o o z

UJ 3 O »- o 1- < z i-t UJ C0 * < > o o o > ** o o o II CO »» z Q o o o 3 o z o o o l- u. w 1-4 < o o o J o H

CO _i ro o o o 3 cr V) CD o o o < * o UJ o • • u. z to a X • -4 *4 CO C4 < *^ o t—

ta

UJ ^4 o **

a _) II

a 2 o•

i-< 3

CO

a

UJX »

<t z UJ 1-

»- u. o o o a w« 3 IL

a. Z * N • o o o a o

UJ «** < o o o a z

CD o o o UJ o (0

-l 1-

X o o o (- II o z < < t—

o o o 3 z 2 ro z a o o o 0. UJ II 3 ^ a. i- • o • z o _l o o < W4 • 1-4 o o •4 o *w» u. z u z z <J zo

t-4 CO

CO

UJa

u>

UJa

UJX

a.

LUKZ

LUOz>WX

LUoozKZ

LUQZ

LUa

LUQ>

-Ja<

LUZ Q Z

UJ t-i H< CO Z K

h-4 1-4 K C5

CO **

lua

CO I o o o o o o o o o o o o o o

(J\ Q\ Q\

C7> On C7> N S N N N I'- ve »0 <0 CM CM CM

GO <0 CD • • • N N N

HH * ro

'»-' . ro ro

-1 ro ro < ro ro 3 ro ro Q ro ro 4 ro ro 10 ro ro

LU • •

or

a Li. z z **> LU o o o t

UJ

MM ^^

a

a 1- • • • 1- 1-

a. LU < <-4 CM *-4 < z **

CD

a 1-4 LU CM CM

O o M s ro

u. z UJ > i-^ 10 10 o o ** <o ro o 3 i- __^ UJ u <t <t o o a \0 ro Z •4 Q 1-4 o O o o LU ^ ro i • LL * tf o o Y- o ro i O Q U. *4 -4 o o U vo ro i LU a. LU s N o o H4 \0 ro i Q < o * * o o Q »0 ro i o CJ • • • • LU • • z a. 1 < u. a. H o

<-4 C0

O

UJo •

II Zo z

UJ

a

LU

Q

LU

Q

LU oz ^»

CO C3 o O

O -J »- 1- a a a a a *-* I cr LU LU LU LU o H ^ o o > > > > >- Z >0 ro o o H4 o o < -4 a z Z z z u LU <0 ro o o <mt o o < LU LU o o o o LU 1-4 <0 ro o o Q o o < 2 a o u u u > u

*4

<0

«0

ro

ro oo oo

LU> o o c o o < O CD LL v0 ro o o a o o c a LL o ro o o LU o o < UJ • a UJ «0 ro o o CO • o < N oz o

u. o

(J • •1 • • CDo

«4 •

u.o o zo

a. i- »rt CM ro 1-4 LU u 1- -> *-4 CM ro * 1-4 -4 cm

r CD LU 3 Z > _J 3 1 o Z CO CO

70
XM

LU

O

O

O o o o o e o o o o o o o o o o o o o

i-I

UJ e o o o o o o © o o o o o o o o o o o o o co

no •H c •o >*-l 4-1 •H c o 3

CO 3 iH 3 O sr O

M-4 co « T3 u CO c 0) 3

•H CO

u

H < CO < 60 Q X CO c

•H

g •H e h E 3 o 4-1 3 •H M CO i-H CO IH e o 4J

c_> 43 o> <r O

n>

X . c

i-H CN •H &

en

a -o

n) QJ c 01 X X. cfl 1-4 w H M oe

"4-1 60 • • o •H CN rH

0) 0) • W x> U

&oQ CO •H O o

Pm s 4-1 <

•-4 ro « ro w ro ro ro

-1 ro ro ro < l*> 10 ro 3 ro ro ro Q ro ro ro M ro ro ro

</> ro ro ro

UJ • • • a l M * ro * w ro ro ro J ro ro ro < ro ro ro 3 ro ro ro a ro ro ro

»* ro ro ro to ro ro ro W • • • a l l oooooooo Z ^ z • o "» o -> M ** w s^ K K 1- H o o •* N N < Z ** ro N N < Z #% * r- N o o ro cm tM M UJ CM CM o o M ro N ro ro CM CM IM UJ (M CM o o 1-4 •* ro N ro CM CM o o ro o 1 ro m ) ro ro ) ro n i ro N | in in

11' • • 11 O > M If) if) o o w< ro <0 ro ro

O o > M If) If) o o Nrf ro ro <0 ro o O o o If) UJ u * * o o Q ro <0 ro ro If) If) UJ U <t < o o Q ro ro <0 ro 10 10 o o ro a W4 o e o o UJ ro >0 ro ro ro ro a « o o o o UJ ro ro <0 ro ro ro o o N u. 4* « o o H ro O ro ro r- N u. * • o o H ro ro <0 ro N N o o r> a u. •* «* o o U ro o ro ro r» r> a U. w* •< o o U ro ro *> (0 N N o o to a UJ s N o o M ro >0 ro ro in in a. UJ N N o o •-• ro ro <0 ro if) 10 • • • < o * » o o a ro <o ro • • • < o « * o o Q ro ro <o • • • oZ u • • • • UJ

or •1 • • az u • • • •

ti • • • < IL a < U. a

II II ii zoM

CM•

(A

O II II ii zoM

ro•V) o II II II Za

X

Ml

<o

<0 o o

e> o o o l/> H O (A K o V) H ce « o o J < z M -1 < z *» -l < K <0 o o < *4 ^ < M n < 1-4 < <0 o o 3 > a %* 3 > a •^^ 3 > Z <0 o o

UJ o H «.

a UJ o H •»

a UJ »0 o o h4 o K Z ro <0 o o •-i o o o M o i- Z ro ro o e w o o e » O UJ • • •

JJ trt 1 X _l< 3 3

jj

a

x -> « to 3 UJ

U Ul ro <0 o o w o o o «/>

u UJ ro ro o o *•» o o o V)

u o UJ M ro <0 o o Q o o o UJ to o UJ M ro ro o o Q o o o UJ </> a z

IX > o ro >0 o o UJ o o o a. -I

a > <J ro ro o o UJ o o o a -Ja < < 1 M ro >0 o o > o o o < < 1 IM ro ro o o > o o o < < M Q to U. ro <0 o o a o o o a 3 a to u. ro ro o e a o o o o 3 a a z u. ro >0 o o UJ o o o Ul O z u. ro ro o e UJ o o o Ul O z < < a. UJ fO NO o o (/> o » o a. * < a. UJ ro ro o e w o o • a. HH < > o ro o o K o o • • • • CD • w4 • < 0) »- o o • • • • CO • • «4 < (rt i- o •o 10 o o </) u. u I o 3 UJ V) u. u o 3 UJ w u <o ro o o 3 a O a O ex <o ro o o

/> _l z (/)

J z co -1 o •o ro o o

Ul < o U. < u u. < UJ <0 ro o o L O 3 3 M IL o 3 1-4 u. o 3 j <0 ro o o e> K -> »» CM ro •* 1-4 «-4 CJ ro o O K -» w* CM ro •* 1-4 *4 CM ro o a < v0 ro o o e a 3 O

1-4 3 £ w 3 z 1-4

u • • • •

c/)

J £ a tf) -I z IX CO V) 1 UJ 8 3 o UJ O 3 o UJ z /) Z a CO Z a U) W) z a 3

71
o o o o o o o o o o o o

ro

*0 oo

_i o O O

y- oooo o o o o o o o o o o o o o o o o o o o o o o o • o o o o o o o o o o o o o o o o o o o o o o • o o o o o o © o o o

a.wQ

<ro

e>zo

_lJo

u.

UII

I- zo

l-l COCO

ID

a.(J

111

IDX

ex. ujz

UJazoQ

>-UJX • O o o UJ • o o o 2 •> a oo oo oo

UJI

•4

•> UJ o • o o o UJ cm J t-M • • if) a n < o

0*

u o 2

w* CO

l-l z

u. • 1- M

UJ o o o o

10 o o o > a o o o M •

UJ UJ <VJ o o o a 1- > a o o o z z o o o o CM UJ UJ M 2 • oo o• I oz

UJ f>4 • -4 UJ CO II V a o CD UJ a 10 a z CO z Q

UJ

H «-*H o UJ > a I X e o o tH 1- J i O o o o o 1- o ct UJ 3 1—1 o o o < UJ < a i-i UJ o o o 3 J UJ o 3 o o o a UJ z o e o o UJ 10 M 2 u.o o• o• oo IL UJ

-1

UJ -* »* • o a UJ I UJ CM a. 1 UJ a •

1-

a

>- o

i- o 2UJ 3 < 1-

tf)

O iu »- < X < Z i^ UJ if) * _l < > o o o > -^ u. o o o II if) 0m, O Q z a o o o 3 UJ o z o o o h- IL %* if) a

l-H < o o o -1 o H Z UJ

</> _jm o o o D

a 2 t10 CD o o o < * o 3 Z

UJ o • • IL z V) _l UJ a X • «-t *-t if) hH < ** o o HH t- a cj >

UJ •H o **

a _l II UJ

a 2 o• < if)

Q

UJ

X

O

X z H- • « 2 UJ K h*

t- u. o o o a CM o u. 3 » a Z * r*- © o o o a o O CM

ID «— < o o o Q 2 _l

CD o o o UJ o CO CM

J UJ

J iX « o o h- II

u z O

Q < < l-H o o o 3 2 u. o 2 ro 2 ct o o o a UJ II O 2

*^

a i- • o • 2 a _J *" UJ

a o < —4 • -) o o —« o I u. S* u. 2 u 2 z o 1- —

Ko Z

UJ o

or l-l a. H o D u O

U. CO o

-J or <

UJ l-lO

O 0> CO H N O O 2 HH O N N D z >0 S N z l-l o O v0 o CM CM a z N 00 CO UJ W • • • - N N N < CO 2 K

l-l i-i r* »-

V)

UJ

a

II 2o CO Q u. z

UJ o o

UJ

M

a. IX

H • • •

u. UJ

CO < CM CM <-<

u. 2 UJ o 3z

h- •o

UJQ 5

coCO

l-l ?c

•Ho

cuCO Phoo

r* CJ o

Ulu •

II z z a

a

a o o UJ UJ LU UJ z **

CO o C5 C5 o » —i H 1- a Ct a a ct *** o I ct UJ LU UJ LU o H 1 o o > > > > \- Z o N

N o M a z z z z u LU o *0 <o o

UJ UJ o o o o LU •* o >0 \0

O 3 a u u o <J >1 U oo

<0

O o o o CO IL o >o <o o

ce IL o >0 vO o

UJ • ct UJ CM vO CM CM N oz o

IL ou • • t •

LUO cto zo a 1- *>4 CM ro l-l UJ CJ »-

n t-i CM ro «f CO UJ D 2 > -1 O 1 O z CD CO _,

72
m IU

•-i Hi3 o o o o o o o o o o o o o o o o o o o o o

** <t ro •*

N— in ro ro

-I ro ro ro < ro ro ro 3 D P) 10 Q ro ro ro

»-• ro ro ro

</) ro ro ro

UJ • • • a 1

«« «•

N

M MO * ro cm w <0 ro ro ro o a »o ro ro ro in

UJ >o ro ro ro ro

t- >o ro ro ro

N

u <0 ro ro ro

N

»-H <o ro ro ro in Q <0 ro ro • •

UJ • • • a 1 ao

UJ>a

UJco

CDO o o o o o o o o o o o o o o o o o O O O O • o o

—I • »

>-> -* <\j ro

II II

C/> <3a

if) UJ CO

a.

J

cu>C

•Ho

T33

(1) CO

(X)

CVJ ozaoo

UJ>I (0

a.o

u.zo ou

N4 o o o — o o o K o o o I o o o

15 o o o

hH o o o

UJ o o o 2 • • • M ro <* ro w ro ro ro J ro to ro < ro ro ro 3 ro ro ro Q ro ro ro

"- ro ro ro

If) ro ro to

UJ • • •

or 1

** ro

N

t—i

ro S ro ro cm

w ro o ro ro o Q ro <o ro ro tn

UJ ro »0 ro ro ro h- ro *0 ro ro s U ro O ro ro N

>— ro >n ro ro

m

Q ro >0 ro • •

UJ • • • a 1 a

•* ro ro

ro ro ro

ro ro «h ro ro **

ro ro ro ro ro in

ro ro o

ro ro rol -*

— cm ro ^ <-> o o o o o o o o o o o o o o o o o o o o o o • o • —I t — cm ro

II II

to

-I<3Q

V)

UJ CO

a. j<3a

to

UJa

LLO2aoz

roozaa

i-o

UJ>I CO

(Xo

U.zo3

_lo

co

o>

tou

o>>c

•Ho

XI

Hi§3OuKZ

LU

U.

U.UJOu

ai 0* c

•H

H-l o o o •U V o o o 3 H o o o O X o o o u

<J o o o X3 M o o o CO UJ o o o S • • • (^

w4 «*4 w*

0> U3

itOo

COuo ^ o

H-« <t ro <t 01 W ro ro ro >

-1 ro ro ro c < o

ro ro ro 3 ro ro ro •H a ro ro ro U M ro ro ro 3

CO ro ro ro

r-\ UJ • • • O

CO

a l 1 0)

ai

S-J ,C

4-1

cu X.H

** •* r«- <N <

** * ro s ro CM w ro ro >0 ro o O •H a ro ro G ro in

UJ ro ro

O ro ro 4J i- ro ro o ro N n) u ro ro c ro r» 00 E

*h ro ro •a ro in 3 o ro ro •a • •

•H

CO «*

UJ • • • 3 X

ora * CO oI n in o o

•C £ O '0

vO nO

eg <£> • •I o o o o o o o o o o o o o o o o o o o o o o o • • • —

-« cvj ro •*m <-< cm ro

II II

CO

-I<3Q

V)

UJ V)

Q

&0 •H

CU

4-1o Eo QJ u CO <u u

OJ OS >c

cu •H .H O &

•a3

a) 0> X CO W (X

01 J3

0)>

•HM

73
oooooooooooooo

oooooooooooooo ••••••••••••oo o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o e o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o

OPIN

awo

a.o o o o o o o o o o o o o o o o 1- o o o o o o o o o o o o o o -1 u z o o o o o o o o o o o o o o J UJ o • o o o o o • o o • o o o o o a •-

9* • • • • • 94 • • 94 • • • • IL

a »- o 1 1 O 3 o UJ O _l

J

C

X o o o \- IL CO

1- o o o o o o o o o o o o o o o o o z J o o o o o o o o o o o o o o o M a. < • o o o o o o o o o o o o o o o UJ <-<

*•** • o o o o o o o o o o o o o o UJ

m i- * o o o o o o o o o o o o o o z Z M

*» o o o o o o o o o o o o o o UJ 3 Z o o o o o o o o o o o o o o X z •-

UJ lil o • o o o o o o o • o o • • o <J J J 94 • • • • • • • 9* • • 94 94 • 10 o z a < 1 1 1 UJ •* z o C3 »- < CO z < (0 X M 94 Z K

01 1- >->

M o o o o o o o o o o o o o o o o »- C3 ai o o o o o o o o o o o o o o > (0 M a o o o o o o o o o o o o o o 94 UJ Q < UJ f* o o o o o o o o o o o o o o a II a o o o o o o o o o o o o o o o Z o o o o o o o o o o o o o o o <o UJ z O z o o o o o o o o o o o o o o I o co u •94 o o o o •1 o o o o •1 o o •1 II > Q

UJ

UJ

IL Z o a • CD a a i- co CO IL UJ < K CO z o oo or Z 1- C\J K o UJ IL Z UJ M I X o o o o o o o o o o o o o o »-) h- o o t- < C5 C5 o o o o o o o o o o o o o o 1- u Z "• a 3 N4 o o o o o o o o o o o o o o < UJ • K M UJ o o o o o o o o o o o o o o 3 -1 19 co 3 o o o o o o o o o o o o o o o UJ UJ z o o o o o o o o o o o o o o UJ CO a o u. o o o o o o o o o o o o o o u o • • o o • o o • o o • • o o u.o

UJ

a UJ CM

a 1 1 UJ z a o z 3 < > i-CJ o UJ

I- <

CM

UJ CO M UJ co » V

H- > o o o o o o o o o o o o o o > 0m II z z X o o o o o o o o o o o o o o II CO — 10 O UJ o o o o o o o o o o o o o o o o 3 CO o

•-< z o o o o o o o o o o o o o o »- IL

w H- t- a. UJ < o o o o o o o o o o o o o o _J o 1- X a uj 3 _J 94 o o o o o o o o o o o o o o 3 a <o o a >

CD o o o o o o o o o o o o o o < * o •H

a z X • • • • o o o • • • o o o • IL z co UJ UJ O K X CO 94 < «* 3 a o

1-1 t-+ 1 ha CM 2 z

ro o

-J o

II oa

UJ • < CO UJ

X UJ • J n Z UJ K 4 N o a u. o o o o o o o o o o o o o o a 94 3 IL z z z >o CO o o o o o o o o o o o o o o o a o IL < »»* < o o o o o o o o o o o o o o o z o or X o o o o o o o o o o o o o o LU o to 1*1 o

UJ fX o o o o o o o o o o o o o o 1- II o z a. \- < fr-l o o o o o o o o o o o o o o D z UJ CJ «* z •* z a o o o o o o o o o o o o o o a ID II D CO UJ cm r4 a i- o o o o o • o z o -J - z >

CM o

u. <z • • • • • 9« • oo oz

94z oo z 1 CO

74
Io

oooooooooooooo oooooooooooooo

oooooooooooooo oooooooooooooo

oooooooooooooo • .......... .00

in

CD

ino

CD

C\J inN

ro

«-« CM

«•

n

n ~ o m

co o

in

w

m

© * cvj • •1 00000*S0-I<ll/1NI1)* oooaw^NcvJsroooooin OOO^IOONOfflCOfOWONH ooon(\H>p)N*noNrtio OOO*rtfflNH(M0NOIfl<t oooNiniflo>inNro — no© oooa>oro»oo-»mincvjcvjro 00000 • * S. »n • • • co in 1 1 O

(M •*

N <£> CD

<t <t

O

IP CDO

CO roN

O in <* O ~4 in

<\J ~* CD z — o ->

CVJ

CO 00

•* *

<* It

10 in

o> *> o <* o ro <* in n r- 00 tO o

00 ro co m •* vo N 00 N o 1 I OOOrtDCOOO OOOh*u1N* o o o * O CVJ 0> O O O vO CVJ O vO o o o in 00 -* s O O O CVJ CVJ •* o 0©00»<>OrO-«COrO-»rON • • • — minwmcvjnincvj-* vOP)(VJ»»»»»«»«»» l 1 1 1 11 o •* -< ro N m ro in cvj —

-< ro o cjv s

cvj in

n cvj vo 00

00 o ro co

N •*

cvj cvj ro cvj o <j>

oj ro cvj o> 00 o

vON CVJ

N Cvj vO O CO 00

ro ro * • o >0

*0 • •

II II II I I

0000

0000 o o o o 0000 0000

Q

UJ

a.HI>zou oz

a.o

t-u

UJ>I CDao

a.za z »*

UJ CVJ M 00

Cvl 1—t

«4

U. in

u.UJ .

^N

u * o> 00 cvj ro o cjv * Cjv •* n o> o ^ in

-1 Cjv

N

N

O

r>. o — o cvj in 00 n ro co * 00 CM Cvj Oi J 4 ' '

in

-I<O oooooooooooooo

oooooooooooooo

'O ro cvj

-« cvj ro •* in *o « -i(\in*ifi\ONoo(fOH cvj ro <i

V)

UJ CO

a. j< Q O

UJQ

a<Da

in

u.o

in >

UJQQa<az M <

<n

h

uj in a

-i<Q

10

UJa. <

UJ

<_> z<a<>ouQ

LU

-I<U

<nz3 N

O -* •« cvj 1* O O

4) «—t

o> 00 vO

<e >n CVJ ro CD

in «o in

m

O

00

N <* ^*O

—i

>o

N

CVJ <t tf

N 00 Cjv CO in ^ in CM ro 00 o 1 I vO vO

<t <*

in in

CO 00 o o o o N N * .I o — - o I I - Cjv O Cvj vo inr- •* o —

CVJ CJV CO 0>

in cvj

— cvj o 1 I CO vO ro 00 m o - o

vO S * «

co in

ro 00 • . 1 c

03C •H O 4J •H 3

4-1 i-H

CO

O > 03 U

0) 0)

03 4=

X> 4-1 O

60C

0) c

4-1 •H c CO H 4-1 CO J2 M O

4-1 03

C c •H

CJa

S-i

M

CO

O

01

C c 60

•H •H iH

01

01 l-l u CO

CO ^

03 03 C 4-1 O X.

•H 60

4-1

H

co 01 3 &

O*

CO

01 01

01N

MX c

4-1 01>

4-1 •H

03 60

l-l •H * 4H 03G

01 J= •H • H 4J CO

CO c >

l-l •H

i-H 01 4-1 03 cfl w XI > O ^j

0> S O 03 & X

00 4-1

O

a

•H 4-1 0> 03 03 03 3 CO 01 rH X. n 4-1 PQ 0> in x: u < Eh H 4-4 <Q •a

03 0) E 4-1 4-1 O J3 3 M 60 D,

14-1 •H 1

0)O

CM s u

CSI rH QJ

0) CO

u

^H 3 CO

I" a*

OJ 03

CO c i-H X 3 CO Pd 3 C XI

01 •H • > 03 <r •H 01 60 M

60 0) 4-1 •H M 3

Pm CO

X

75
ooooooooooooooooooooo

ooooooooooooooooooooo

ooooooooooooooooooooo oooooooooooooooo • • • • • OOOOOOOOOO • • • • • • <t O N CM O oooooo • • • • o*B — ©ro©o©© — © O • • • • •0>ONNlftO0>(84'<0otftNO«*

<0O»*\000C0«-'N0i-tf-<NNOC0-<<-<ir>O>£>r0 N<Mcvi^ro-i-<«-fO<o-<^(vj*moo — — cvjcvirn I I

V)HZ

13M©2

-1 o < o D -1 o O o o

UJ t- oo • o

-1 o < • 4zoz UJ o

>- -J J < o u a ©

N4 •o

UJa

UJ —t z o

h- o © z •

<j> «

in• © a X

tH

X • o o o 3 M s l-l UJ

o> 36 •4

u. • o z© UJ 44 a © a o < > i- • i- o

1- i—»

UJ < > K© a • <» z a o <

UJ j •4 • z CO CD < It • o> X

-»

X •4

— o CO C£

a 5 o H a o <

UJ • z

_i ro a Z <0 u.

m z w < <3 i-< t—i

** Z •4 z a. 00 CM (X iCM o <

N^ u. z + + + + + + + + + + + + + + + + + + + + + o o © ** O* © © f- w4 a> o CM s CM © CM s o N N O o o o •* CM -4 in -4 00 ro o CO #- If) •* N •« © N ro o o «-4 CO N in 00 CM in CM CM — in <* ro 44 © <t * © © CM o o o CM © CM n in ro © in o ro «• <* ro ro * N ro ro o 44 •* Q> © © 00 CM o> oo o in CM CM © •* o CM N © © o in o in *4 N ro o © ro r-t ao © <t CM f-4 o a> © © © o o t-4 f4 CM CM ro * <t in © © N © 0> o

CM w4

CM CM

CM

CM ro *CM CM o a o a> 0> O* a <7» a> <7> o o> CT> <* (? 0> 0> a> O' o & o o o o O O o O o o o o o o O O o o O o o + * + + + + + + + + + + + + + + + + + + + o o CM 00 $

CM

O o o © o <* * © CM ro •* CM

© CM

O o <* <M 00 © N © 0> «-« o o a ro O © a N N 0> O o o ro o in O S & 00 00 o N «-4 N © o ro © N ro o o © <* in 00 in <* N <t in «"4 O in « 0* o © © N ro © o o CM in o in CM o o 44 * CO ro o © O o ro © © ro o <* 00 CM © r* © 44 © *4 -0 44) N ro © « w* N ro O N o o o r* — CM CM ro ro <t * in in © © r- © © o o o

CM CM

ooooooooooooooooooooo o — ©N«*©©roc\io*o»*©N«*©©rocMO»o ooocm©<\i — -* — cMorocMa>**r-o — ro©o orocMN©©oOf>-oo©a©©©©©0'*o 00-»N*N-lfi(M()- N*CM-«OO-Jr0©ffl OrO©WCM©OCM©0>rO©0'*fflCM©0<t©CM oooo — -<»-cvJCMCMroro«t'**©©©©©f«-

ooooooooooooooooooooo o —1

<t o ©

©

© o <* PMt o fH <* & ©

©

© 0> tf «4 o o o o o wt <M ro <t © © o CM H* © 0» CM © © CM © o o N <t © © o CM •* © © w* ro 10 N (7> CM <t © 0> «* <t o o o o o •"« »4 •-• ~* •M CM CM CM CM CM ro ro ro ro * <t

ooooooooooooooooooooo ooooooooooooooooooooo

ooooooooooooooooooooo

0-<CMrO'*©©N©C>0-4CMrO<f©©N©0>0 OOOOOOOOOO — -< — -<f-^*-<rt-4-(CM

ooooooooooooooooooooo

ooooooooooooooooooooo

a

UJa

a.a

e>zM3o

_l

_lo

u.

UJX

UJz

UJX

<J©oz

•4Ho>Mn

ro UJi M © ©z O o UJ

1-4

K

H

O < UJ D J O UJ

UJ

©

IL UJ o a

UJ z 3 o UJ

1- < © * > 4*1

ii

© 3

t- u. **

-i o (- o or < ^ o

u. z ©

©

M < W

Ha

J ll 3 a © UJ X

UJ

H

a. w*

D LL a O

Q z

UJ o ©

»- II

u z D z a Ul II D z a -I o o *N o u z z o zo

1-4 ©©m

a.o

UJ

a.

UJx

a:UJHZ

tu

I-oz>

UJI

I-©oozZ©QZ©a©a>

or<©z©< i-© i <

ILO©zzJo

<J

13Z © o

-I

_lo

LL©I

I- zoa©©a

u.

u.o©QXo©3o

a.©N

u.o

a.©©zoz o© ( aaoo©o©© > zza© < < ( Z »

V)©©o

a.©mzDZZo

I-ao

a.

_

© c a

76
_l<Q—

LU

or z ~ o -> o o o o o o o o o o o o o o o o o o o o o o o o

rooI — «* N o in o s - to n •* ro SOJ)

<o ro o

inn oo • — N

0* • • Q

UJh-U

»-( Q

HIaa ozaoHU

HI>I 00ao

ILzo z

HI

i-i u—

Ii. h

HIOu o o o o o o o o o o o o o o o o o o

CVJ oI ro oo in ON* mo- <0 oo ro o oo in • in o o O N o • a H N • I

ooooooooooooooooooooo ooooooooooooooooooooo ooooooooooooooooooooo ooooooooooooooooooooo

Nonoo>niflo-**nior) ONONOniMOOIJlH o^- • • n m oo •

00 o • fO 00

— CM

— • • in cvj ro • •

ro in o> cm

a <o o>

vo

n oo a> o oo in

O ro

— *

>o in— <*

cm 0*> ro

— o , SO(D0i«

in>oONCorocMco

— rorom«f©o BiM0l,l*-inNOO<0<fNO

rocMCMino*co

— ini>roi>-ino* • • ro - - o

00 vO

O

in to

a

roro

— o — cMinin*

— cMinooo^oooovo—

n—

cm - <t i i i

CVJ Ov

—

—

— *4-NhP)HH- p4 f4 I I I I I I N

N I

ro

c\iov<oo*romo~<*rON OONONOr)0>00*0 • • n in oo • • • • in ^

coco • • • in o h

m • • cJNroov — cviroovooov —

fOcvjcvjcvjscoincvjo^N niocm — oincM<oinoom

ro ro ro

—

© cm * ro o >0 o

cm

— <roco

——

—

— ^* I I I I I I I

o>-**cvjino>oo>cvj^ NI"l*ONt««ti • • • • • in o ao o <o cm o oo cm — a> en * in n

ovfoininmomrocvjcvi rocrinNocovo^^o J-OSOOOM1* in — ocvi<7>rovoomo — ro in n o* — — cvjcvjro o o o o ooooooooo ooooooooo ooooooooo o o o o o • • • • • • • • • ov o s cm in ©

CM— •tfC0C0r»->O<O00—— CM— CMtf)— >OCM

<*

— aooomromcviroroinio. o —— cvjoo

O

— vO OO 00— N<7>^—

N

N

© 00

—

— tf>o

ooooooooo oooo**«*» • • • <* O N

ovoo

— orooooinvo

t> oo * vo o in >

(M

CMCM

— fO——— fOvO—

— I I o

00 vO - - a - vO

\0 ro

cm <* in co

— -cvjcvjio

— cvj to * m vo

—— cvj ro * m«oNooo>o

— cvjro^mvONooo>o—

— -^^^

—_— - — (vjrvj +o in ro

cvj

m in o • •

in vo ro

oo vo in oo oo ro

cvi « * * in in • <o —

ii ii ii

z> > a lu - Q

V)uj in

a

cc _j

a < < O Z - < W IHl tfl a

Q

HI

tr<a

to

OO

C

oo

C 3 c •H O

•H iH <4-l W Cfl 3 O CO W •H O 4-1 < m 3

<C

OX

H .c •H < 4-1 M O •H 4-1 S <fl e B oM • VO

4-1 o X

OO II CN <-t w CM 3 0) & CO 4-1 (0 M X >4-l w C O

tfl ^ • • c CO

in iH tou QJ • W X>

ooQ 01 •H O x; o

F*4 £ 4-1 JJ

<0

1-4Iin

in

ooNN

in

inCM - «t — —4 I I CM >0

<*

N

O h- 09 O

vO 00

in

v

cm «o - in • •I

00 o CM o — — 1 CVJ N 1

<ta

N

N s in

CM 00 o

CM <o

—

tf It ro

CO

a

—

ro >o

m • •l • o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o XH

or<2

HIUZ<—a<>o

<_> Q

HI_J<O

10Z

OOOOO

OOOOO o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o

77
ooooooooocooooooooooo ooooooooooooooooooooo

ooooooooooooooooooooo oooooooooooooooo • • • • • OOOOOOOOOO • • • • • • <t O N CM O oooooo • • • •Qi<oxioniooifl<0'*flo O • • • • iOOSNI0OCM0**O»Ni0« • w — *coooN>oooco — ^cvi-.wm=-<£)w — o> o^rtoocoiflniONnninnoHHiMooxio oo-<\0(Doo<-'SO| ^, <-<sso<D-<-imovoro n n n h n h - ~< ro <o - — w<tinoo-<-<cvjcjro I I

in a. o o • o o o o o o o o o o o o o o o o o o o o o in o

<0 o o o o o o o o o o o o o o o o o o o o o CO z

t- • o o o o o o o o o o o o o o o o o o o • in

M

X o • • • o o o o o o o o o o o o o • • • o o 3

u» o o CO

N • • • o o o o o o o • • • r- CO o o o o tIM o •* >0 o O in <* • o o o o o • <± in >0 o o <* o a J u

UJ o o s CON CM CM ro • o o o • ro CM CM s oo

N o o N -1 UJ 3 o o CVJ >0 N <mt o <t CM • o • CM <t o «M r- NO CM a o in o or w4 m ro w* s ro *4 CM fO *4 • f-4 ro CM 9* ro r» f* ro in «* • u. a. _l O • 1 1 1 1 1 1 1 1 1 1 o < o o ID u 3 -I o z «* X o O o w o 1- IL UJ t- o _J 1 o o < o o o o o o o o o o o o o o o o o o o o o N z • o u o o o o o o o o o o o o o o o o o o o o o * M a

_l o <n o o o o o o o o o o o o o o o o o o o o o 4-4 UJ < • • o o o o o o o o o o o o o o o o o o o • CO UJ CO M X o • • • • o o o o o o o o o o o • • • • o o z z z K o M vO •4 <0 • • o o o o o o o • • vO *4 o •H o in UJ 3 o >H o o a o a m <0 • • o o o • • «0 in 0> o 0» «0 o <* X z z UJ (VJ 3 o in o « CM CM m *< «o • o • vO *M in CM CM <t o in o in o

>- _l -» «o «* CM »* vO CM CO p* ~t • (-4 ** 00 CM -0 •M CM <t o «M • V)

a

J < Z UJ o u O V H a MH

roo z

fl <z • U 1 H H o UJ o o o o o o o o o o o o o o o o o o o o o ro o in 1-

UJ z o o o o o o o o o o o o o o o o o o o o o <* > f* co Q z o o o o o o o o o o o o o o o o o o o o o *4 M UJ

UJ

— o o o o o o o o o o o o o o o o o o o o o o CM

a II

X

Q (J • o o o o o o o o o o o o o o o o o o o • N

>- o o • • • • • o o o o o o o o o • • • • • o CM <o UJ z

in z z o o CM to vO

m • • o o o o o • • in <0 ro CM o o o X

a

-• o CM «-4 <* V* CM <f

N • • o • • N «* CM •-4 <t rt CM o in 1-

a IL » ~4 N in ro CM •4 o CM 00 H • ** CO CM >0 »H CM ro m r- —< • UJ o

ON *» 1 1 1 1 1 1 1 1 1 1 II > UJ

•* eg GO a.

a

in • CM 10 IL UJ • •* CO o z a CO a 1- -« o w 1- 1 o UJ IL z I o < X o o o o o o o o o o o o o o o o o o o o o CO H-» 1- o 3

•> 15 o o o o o o o o o o o o o o o o o o o o o o o in tu z o 3 o u. * o o o o o o o o o o o o o o o o o o o o o •* < UJ • s M o o UJ o o o o o o o o o o o o o o o o o o o o o -0 3 _) w

0< • 2 o o o o o o o o o o o o o o o o o o o o o CM o UJ UJ H o H • o o o o o o o o o o o o o o o o o o o • * UJ 10 o

•* z IL o • • • • • • o o o o o o o • • • • • • o < • V* UJ o o #* <t a vO in >0 • • • o • • • o in <o o< «• t-4 o <t u. UJ z z P4 co vD <* n CM <"4 a> <* PM] • *4 <* <J< w* CM ro <* >o 00 *4 • o a

10 UJ «4 UJ

a UJ CO 0. -J o z 3 < > II UJ y- *4 o UJ o • I- u o K •t »- H4 ** X UJ 1 (0 * < CM u > o o o o o o o o o o o e o o o o o o o o o 00 > #4» II z

h- « < o o o o o o o o o o o o o o o o o o o o o o II V)

— •* o

<0

M UJ

a o o o o o o o o o o o o o o o o o o o o o <t 3 CO * ** W z o o o o o o o o o o o o o o o o o o o o o r* r- u. «* 1- 1- a o < z < o o o o o o o o o o o o o o o o o o o o o ro -1 o K X a. U -J - • o o o o o o o o o o o o o o o o © o o o o o o 3 or ro o o z 00 a a. CO • o o o o o o o o o o o o o o o o o o o • <0 < i£ o M4 0. < U. o u. o o ro u. z (0 UJ UJ • a K X o CO N <0 in <* ro CM *H • r* CM ro <* in <0 N CO o *4 • V) I"H < w 3 a

->

XU

Q 1-4 1 1 1 1 1 1 1 1 1 1 1-a tO - o 00 UJ UJ a J II o • z a > 1- K to 3 a a. a. o u < a V) UJ X UJ UJ • U. < Z o UJ H U) N j ro o a o o o o o o o o o o o o o o o o o o o o o K 0> a. V* 3 u. a z o u. t- w o o o o o o o o o o o o o o o o o o o o o u 00 a o IL s w z CD < o o o o o o o o o o o o o o o o o o o o o < N a z o < < 3 o o o o o o o o o o o o o o o o o o o o o U. f-4 UJ o 10 CM * K UJ CO X o o o o o o o o o o o o o o o o o o o o o CM 1- II u z a. < z tmt o o o o o o o o o o o o o o o o o o o o o UJ CO 3 z UJ

**» z -i Z V) a. o o o o o o o o o o o o o o o o o o o o o -1 **

a UJ II 3 CD

CO cm

a II -» t- < CM z Q _l - z

CM o < o • o o V4> o 3

u. 3 3 z (0

u z z o z

78
X

O

HI3

ooooooooooqoooooooooo ooooooooooooooooooooo

ooooooooooooooooooooo ooooooooooooooooooooo ooooooooooooooooooooo

00 o X

c •H •H

•H u M

CO CO 4-1 3 e COS

* 4-) Q 3 \C m CO CD ?s -C O

X < 4-1 iH 1 O •H

IS

tsi CD B CD tOo

ou

X!

4J

14-1 CN

<4-l in 00 II o • eg CM vO

CNJ w .*

CD hJ C 01 •

.H < CO XI t C/l M

o M

c

cu x:

4J

4-1

B

ino

I

ino

l CO •. 3 en in

• oo O in N vO c

•H

i-H

<4-i

CO 10 S

&o Q CO •H cm in •H O o in *

fc JS CO < ooO

CM

0> P> CM

<7> O O

O 0> O

- 0100*10000'

** roinoioo»<»*N

J coo>o>o©ooo»

< o>o*o»ooooo>

Q • N N CO I") CO CO N

m CO-*<ttft\l<S-<»-ct

to inoooinooo^

UJ N(Mt\JtMNCMWCM — t\|C\lt\(-»

a. i I i i i i

N N

O O

oxuiOriOiOo

hcjiioi,)o> iomC}> ot>ooooo>*

• ••••••o 00 N N 10 O N 00 • «*«**CM*<**0* oooinooom «NN(\|(\|N(\|N

I I I I

roo

o in

i CO e CO N

CO ro m N *4 o

O •o N

V4 o ro • in N

ro »* CM

to • •

z «-» o -> M <*> O

uz

a

t-V

UJ>

I CDa

o

u.z

o

r-Jo

in

>

UJo

Q

a

<o

z

<

rZ

UJ +o

>0

in s

ro • cm

_ CM 0> .4

<j> cm in •

cm ro -• o

0> O •* CO N

N N o m • •* M « <f |<) M H

CD - •* O •

<* CO ro « — CM

oo roCM

a

UJ

o

UJa

a

N <* <o oo <0 S

ro ro •* * 0* »*

(j> cd O ro O o (400010 rtlDQ O • • O <J> O • * '-0 O O • 4 O O

CM 00 0* CO CO a O r>

a* o> o* a> 0>

o o o o o oCO o o o

o> • • •

• H CO ffl in in co oi h

in n

o o

>0 ro * _. NCDoi-ininnoi'** oo*oj--<j(^H«©coorO(o

ro>0'*'*-*r>-<-'CMO»ro«^oo»<<o

»* in >o ro — rocji-«ro>«»« — <M'*inco«-4-<cMCMn

t-Z

UJ

1-4 u

l-l u.u.

UJu

o 0> - o

f +

>0 CM O CM

(J. o

Q> ro o o

o o

g\ <J<

• •

I

o ro ro o

• N ro o

*0 o o o

o • o o

ro — o

o o

•» >o a

a> o ij> 0> o> .<t

o

o

o

I I

a

UJ>

a.

UJm

m

o

ooooooooooooooooooooo oooooooooooooo>o> oooooooo • • • oooo • • • • c co •*

• •••aoNcMinoo>„ — 4j-C0C0N'OCDC0 — —eg — c\jin-« — cocoinroinc\jiOfoinroo-<»<cM HiflcOH)rtS()i^H •> N CO *4 «4 cM(M ,-, rO'-<-<«"«ro\o««^CM

I I

MONNO ONNOIOOmO COiOOlON^H

>0 CM — C7> 00«<0

in o <o ro

in o »* -* cm <vj ro

» N rO 4 10 (0 • *Nr0*l0'0NC00>O«Nn<'ini0N 00 0> O H

— - CM <M

• <J\ (J\ O -

o o

CM CO

<t • •

in o o

in •* <o

io «h ro

co o> cm

ii ll ii zo

M

v> r-

_J <

< *

O UJ M O

(0

uj m q

a j a

< <

O 3 Q uoz a p-< < < in h

3 uj <n o

en

u.

u. o

o

z

z a o o

<n z

<

z

UJu

z

<a

<

>

a

uQ

UJ

_J<u

enz

3

O

I ro <0 N CM & l> •* iO • NO* « cm in n

oo • ro in — in m o

o> to o >o

ro ro

I

o

ro

10 • «0 * N

ct« in « ro o

ro o • < •» O o ro ro cm — in r» co ao o o ro • *

CM 00 0> vO "0

-. -. g> o •

CM — 0> P4 »< O

+ +

CO CM N

ro m —

CM .0 0> « CM CO N ro 0> N <0 N

m CM •* • • •

I

0>

I

- 0> *

vO 00 ro c*

I

(7> ro m

ro co • in (M -<

79 U. S. GOVERNMENT PRINTING OFFICE : 1980 311-046/168
NBS-114A (REV. 9-78)

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO. NBS TN 1126

2. Gov't Accession No.

4. TITLE AND SUBTITLE

Problems Used in Testing the Efficiency and Accuracy of

the Modified Gram-Schmidt Least Squares Algorithm

3. Recipient's Accession; No.

5. Publication Date August 1980

6. Performing Organization Code

[tali

kdoli

7. AUTHOR(S)

Roy H. Wampler

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, DC 20234

19. m$(s<£ffmk/mti Unit No.

;!i;!i

11. Contract/Grant No.

id

!;" flV

flit).

k

la

hi

pneei

Mil

12. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street, city, state, ZIP)

Same as item 9.

13. Type of Report & Period Covered Final

15. SUPPLEMENTARY NOTES

| | Document describes a computer program; SF-185, FIPS Software Summary, is attached.

16. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant bibliography or

literature survey, mention it here.)

In preparing "Algorithm 544: L2A and L2B, Weighted Least Squares Solutions

by Modified Gram-Schmidt with Iterative Refinement" for publication in ACM Transactions on Mathematical Software (Vol. 5, 1979) , the Fortran computer program was extensively

tested. This note describes the various types of problems which were used to explore

the efficiency and accuracy of this algorithm. The Fortran subprograms which

performed the various tests are listed in an appendix. Also listed are the data used

in executing the testing routines as well as typical output from several different

types of problems. Among the testing routines is one which is suitable for handling

general linear least squares problems. Here, the user has the option of scaling

his raw data in order to mitigate the effects of ill-conditioning.

wis

to *

Hit!

:-.•

nil!

w

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of the first key word unless a proper name;

separated by semicolons)

Algorithms; curve fitting; least squares; modified Gram-Schmidt; pseudoinverse;

regression; statisics; test problems; test results.

18. AVAILABILITY Unlimited

I I For Official Distribution. Do Not Release to NTIS

JQC] Order From Sup. of Doc, U.S. Government Printing Office, Washington, DC

20402.

Order From National Technical Information Service (NTIS), Springfiel-d, VA. 22161

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF

PRINTED PAGES

83

22. Price

$4.00

USCOMM-DC

.*
NBS TECHNICAL PUBLICATIONS

PERIODICALS

JURNAL OF RESEARCH—The Journal of Research of the

titional Bureau of Standards reports NBS research and develop-

"|nt in those disciplines of the physical and engineering sciences in

v"ich the Bureau is active. These include physics, chemistry,

a;ineering, mathematics, and computer sciences. Papers cover a

"Triad range of subjects, with major emphasis on measurement

rthodology and the basic technology underlying standardization. ^o included from time to time are survey articles on topics -^sely related to the Bureau's technical and scientific programs.

$ a special service to subscribers each issue contains complete

;itions to all recent Bureau publications in both NBS and non- US media. Issued six times a year. Annual subscription: domestic

71); foreign $16.25. Single copy, S3 domestic; $3.75 foreign. OTE: The Journal was formerly published in two sections: Secn A "Physics and Chemistry" and Section B "Mathematical

ences." ENSIONS/NBS—This monthly magazine is published to in- m scientists, engineers, business and industry leaders, teachers,

""Vdents, and consumers of the latest advances in science and

t hnology, with primary emphasis on work at NBS. The magazine

hlights and reviews such issues as energy research, fire protech, building technology, metric conversion, pollution abatement,

•ijilth and safety, and consumer product performance. In addi-

|m, it reports the results of Bureau programs in measurement

sndards and techniques, properties of matter and materials,

ej>ineering standards and services, instrumentation, and

ttomatic data processing. Annual subscription: domestic $11;

eign $13.75. NONPERIODICALS

nographs—Major contributions to the technical literature on

tfious subjects related to the Bureau's scientific and technical acuities. indbooks—Recommended codes of engineering and industrial

! dctice (including safety codes) developed in cooperation with in- vested industries, professional organizations, and regulatory

dies.

ecial Publications—Include proceedings of conferences sponred by NBS, NBS annual reports, and other special publications

propriate to this grouping such as wall charts, pocket cards, and

jliographies.

iplied Mathematics Series—Mathematical tables, manuals, and

idies of special interest to physicists, engineers, chemists,

ologists, mathematicians, computer programmers, and others gaged in scientific and technical work.

itional Standard Reference Data Series—Provides quantitative

|ta on the physical and chemical properties of materials, comlled from the world's literature and critically evaluated,

jeveloped under a worldwide program coordinated by NBS under

ie authority of the National Standard Data Act (Public Law

1-396).

NOTE: The principal publication outlet for the foregoing data is the Journal of Physical and Chemical Reference Data (JPCRD)

published quarterly for NBS by the American Chemical Society

(ACS) and the American Institute of Physics (AIP). Subscriptions,

reprints, and supplements available from ACS, 1 155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical information developed at the Bureau on building materials, components,

systems, and whole structures. The series presents research results,

test methods, and performance criteria related to the structural and

environmental functions and the durability and safety characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in their treatment of a subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures

published by the Department of Commerce in Part 10, Title 15, of

the Code of Federal Regulations. The standards establish nationally recognized requirements for products, and provide all concerned interests with a basis for common understanding of the

characteristics of the products. NBS administers this program as a supplement to the activities of the private sector standardizing

organizations. Consumer Information Series—Practical information, based on NBS research and experience, covering areas of interest to the con- sumer. Easily understandable language and illustrations provide

useful background knowledge for shopping in today's technological marketplace. Order the above NBS publications from: Superintendent of Documents, Government Printing Office. Washington. DC 20402. Order the following NBS publications—F1PS and NBSIR's—from

the National Technical Information Services, Springfield. VA 22161

.

Federal Information Processing Standards Publications (FIPS

PUB)—Publications in this series collectively constitute the Federal Information Processing Standards Register. The Register

serves as the official source of information in the Federal Govern- ment regarding standards issued by NBS pursuant to the Federal Property and Administrative Services Act of 1949 as amended,

Public Law 89-306 (79 Stat. 1127), and as implemented by Executive Order 11717(38 FR 12315, dated May 11, 1973) and Part 6 of Title 15 CFR (Code of Federal Regulations).

NBS Interagency Reports (NBSIR)—A special series of interim or

final reports on work performed by NBS for outside sponsors

(both government and non-government). In general, initial dis- tribution is handled by the sponsor; public distribution is by the National Technical Information Services, Springfield, VA 22161,

in paper copy or microfiche form.

BIBLIOGRAPHIC SUBSCRIPTION SERVICES

he following current-awareness and literature-survey bibliographies

Ire issued periodically by the Bureau:

ryogenic Data Center Current Awareness Service. A literature sur-

!sy issued biweekly. Annual subscription: domestic $25; foreign

30.

iquefied Natural Gas. A literature survey issued quarterly. Annual

inscription: $20.

Superconducting Devices and Materials. A literature survey issued

quarterly. Annual subscription: $30. Please send subscription or- ders and remittances for the preceding bibliographic services to the National Bureau of Standards, Cryogenic Data Center (736) Boulder, CO 80303.
U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington. D.C. 20234

OFFICIAL BUSINESS

Penalty for Private Use. S300

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE

COM-21S

SPECIAL FOURTH-CLASS RATE

BOOK