5631runlisting.txt

ian@ian-HP-Stream-Notebook-PC-11:~/563dir$ ls
563  5631.f  5631run  563.f  563.sh  563.txt  fort.7  xx00  xx01
ian@ian-HP-Stream-Notebook-PC-11:~/563dir$ ./5631run < 563.txt
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   6
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   0
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
     1.000     1.000
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     1.000     2.000
    -1.000     0.000     0.000
     0.000    -1.000     0.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    2
      RESIDUAL L1 NORM =   3.0000000E+00
0***** CYCLE NUMBER     1
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    1
      RESIDUAL L1 NORM =   3.0000000E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    4
      RESIDUAL L1 NORM =   3.0000000E+00
0***** CYCLE NUMBER     3
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    1
      RESIDUAL L1 NORM =   3.0000000E+00
0***** CYCLE NUMBER     4
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         1    4
      RESIDUAL L1 NORM =   3.0000000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  3.0000000E+00
0X-VECTOR
  1.0000000E+00  1.0000000E+00
0ACTIVE EQUATIONS/CONSTRAINTS
    1    4
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   5
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   2
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
     3.000    -1.000
0ARRAYS  E  AND  F
 
     1.000     1.000     1.000
     1.000     2.000     1.000
     1.000     3.000     2.000
     1.000     4.000     3.000
     1.000     5.000     2.000
     1.000     1.000     2.000
    -1.000    -1.000    -3.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       3.0000000E+00 -1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2
      RESIDUAL L1 NORM =   1.1000000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   1.1000000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.0000000E+00
      X-VECTOR...
       2.0000000E+00  0.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    5
      RESIDUAL L1 NORM =   3.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   3.0000000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  3.0000000E+00
0X-VECTOR
  2.0000000E+00  0.0000000E+00
0ACTIVE EQUATIONS/CONSTRAINTS
    6    5
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   0
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.6875000E+00
      X-VECTOR...
       3.7500000E-01  6.2500000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  1.0000000E+00
0X-VECTOR
  5.0000000E-01  5.0000000E-01
0ACTIVE EQUATIONS/CONSTRAINTS
    3
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   2
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     0.000    -1.000
     0.000     1.000    -1.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   1.8250000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.1250000E+00
      X-VECTOR...
       3.7500000E-01  6.2500000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   1.0000000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  1.0000000E+00
0X-VECTOR
  5.0000000E-01  5.0000000E-01
0ACTIVE EQUATIONS/CONSTRAINTS
    3
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   2
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     0.000     3.000
     0.000     1.000     3.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   2.6250000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.2500000E+00
      X-VECTOR...
       7.5000000E-01  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         2
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   5.2500000E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   2.5000000E-01
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    1
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   5.0000000E+00
0***** CYCLE NUMBER     3
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         1    2
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   4.1250000E+00
0***** CYCLE NUMBER     4
      STEP TAKEN =   2.0000000E+00
      X-VECTOR...
       3.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    4
      RESIDUAL L1 NORM =   5.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   2.6250000E+00
0***** CYCLE NUMBER     5
      STEP TAKEN =   2.0000000E+00
      X-VECTOR...
       3.0000000E+00  3.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    5
      RESIDUAL L1 NORM =   9.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   1.1250000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  9.0000000E+00
0X-VECTOR
  3.0000000E+00  3.0000000E+00
0ACTIVE EQUATIONS/CONSTRAINTS
    4    5
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   1
0NUMBER OF INEQUALITY CONSTRAINTS =   0
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     1.000     1.500
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   2.1750000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.2083334E+00
      X-VECTOR...
       6.2500000E-01  8.7500000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         4
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   1.0000000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  1.0000000E+00
0X-VECTOR
  7.5000000E-01  7.5000000E-01
0ACTIVE EQUATIONS/CONSTRAINTS
    4
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   1
0NUMBER OF INEQUALITY CONSTRAINTS =   0
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     1.000     3.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   2.3250000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.2500000E+00
      X-VECTOR...
       7.5000000E-01  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         2
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   2.2500000E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   2.5000000E-01
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    1
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   2.0000000E+00
0***** CYCLE NUMBER     3
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    1
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   1.1250000E+00
0***** CYCLE NUMBER     4
      STEP TAKEN =   1.0000000E+00
      X-VECTOR...
       1.0000000E+00  2.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         1    4
      RESIDUAL L1 NORM =   3.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   3.7500000E-01
0TERMINATION CODE =   1
0L1 VALUE.     =  3.0000000E+00
0X-VECTOR
  1.0000000E+00  2.0000000E+00
0ACTIVE EQUATIONS/CONSTRAINTS
    1    4
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   1
0NUMBER OF INEQUALITY CONSTRAINTS =   2
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     1.000     6.000
     1.000     0.000     2.000
     0.000     1.000     2.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   3.6000000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.0000000E+00
      X-VECTOR...
       1.0000000E+00  1.2500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         1
      RESIDUAL L1 NORM =   1.5000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   7.0000000E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.0000000E+00  0.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         3    1
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   8.1250000E+00
0***** CYCLE NUMBER     3
      STEP TAKEN =   5.0000000E+00
      X-VECTOR...
       1.0000000E+00  5.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         1    4
      RESIDUAL L1 NORM =   9.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   2.1250000E+00
0***** CYCLE NUMBER     4
      STEP TAKEN =   1.0000000E+00
      X-VECTOR...
       2.0000000E+00  4.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    5
      RESIDUAL L1 NORM =   9.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   1.1250000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  9.0000000E+00
0X-VECTOR
  2.0000000E+00  4.0000000E+00
0ACTIVE EQUATIONS/CONSTRAINTS
    4    5
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   1
0NUMBER OF INEQUALITY CONSTRAINTS =   2
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     1.000     2.000
     1.000     0.000     1.000
     0.000     1.000     1.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   3.0000000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   9.6875000E-01
      X-VECTOR...
       8.7500000E-01  1.1250000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         4
      RESIDUAL L1 NORM =   1.2500000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   1.3750000E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   8.3333336E-02
      X-VECTOR...
       1.0000000E+00  1.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    6
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   1.0000000E+00
0TERMINATION CODE =   1
0L1 VALUE.     =  1.0000000E+00
0X-VECTOR
  1.0000000E+00  1.0000000E+00
0ACTIVE EQUATIONS/CONSTRAINTS
    4    6
1
0NUMBER OF VARIABLES =   2
0NUMBER OF EQUATIONS =   3
0NUMBER OF EQUALITY CONSTRAINTS =   1
0NUMBER OF INEQUALITY CONSTRAINTS =   2
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
    -3.000    -2.750
0ARRAYS  E  AND  F
 
     1.000     0.000     1.000
     0.000     1.000     1.000
     1.000     1.000     1.000
     1.000     1.000     2.000
    -1.000     0.000     0.000
     0.000    -1.000     0.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -3.0000000E+00 -2.7500000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4500000E+01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   2.2250000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   1.1250000E+00
      X-VECTOR...
       3.7500000E-01  6.2500000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   3.0000000E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   2.1250000E+00
0***** CYCLE NUMBER     3
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   2.0156250E+00
0***** CYCLE NUMBER     4
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   1.95E-03
      PENALTY FUNCTION VALUE =   2.0019531E+00
0***** CYCLE NUMBER     5
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   2.44E-04
      PENALTY FUNCTION VALUE =   2.0002441E+00
0***** CYCLE NUMBER     6
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   3.05E-05
      PENALTY FUNCTION VALUE =   2.0000305E+00
0***** CYCLE NUMBER     7
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   3.81E-06
      PENALTY FUNCTION VALUE =   2.0000038E+00
0***** CYCLE NUMBER     8
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   4.77E-07
      PENALTY FUNCTION VALUE =   2.0000005E+00
0***** CYCLE NUMBER     9
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       5.0000000E-01  5.0000000E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         3
      RESIDUAL L1 NORM =   1.0000000E+00
      PENALTY PARAMETER =   5.96E-08
      PENALTY FUNCTION VALUE =   2.0000000E+00
0TERMINATION CODE =   2
0L1 VALUE.     =  1.0000000E+00
0X-VECTOR
  5.0000000E-01  5.0000000E-01
0ACTIVE EQUATIONS/CONSTRAINTS
    3
1
0NUMBER OF VARIABLES =   5
0NUMBER OF EQUATIONS =  15
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   0
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
     0.000     0.000     0.000     0.000     0.000
0ARRAYS  E  AND  F
 
     5.000     3.000     4.000    12.000     4.000     7.000
     9.000     7.000     3.000    19.000    13.000     4.000
     6.000     6.000     0.000    12.000    12.000     2.000
     9.000     9.000     7.000    25.000    11.000     7.000
     3.000     0.000     1.000     4.000     2.000     7.000
     8.000     1.000     8.000    17.000     1.000     7.000
     1.000     9.000     8.000    18.000     2.000     3.000
     0.000     9.000     3.000    12.000     6.000     3.000
     3.000     1.000     1.000     5.000     3.000     5.000
     6.000     7.000     6.000    19.000     7.000     1.000
     6.000     1.000     9.000    16.000    -2.000     4.000
     0.000     4.000     8.000    12.000    -4.000     1.000
     0.000     5.000     7.000    12.000    -2.000     6.000
     7.000     3.000     2.000    12.000     8.000     6.000
     5.000     4.000     9.000    18.000     0.000     0.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
        15
      RESIDUAL L1 NORM =   6.3000000E+01
0***** CYCLE NUMBER     1
      STEP TAKEN =   3.7207752E-03
      X-VECTOR...
       3.8859557E-02  8.5410953E-02 -1.0269672E-01
  2.1573827E-02  2.2696729E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
        15    2
      RESIDUAL L1 NORM =   4.9712200E+01
0***** CYCLE NUMBER     2
      STEP TAKEN =   2.0141857E-02
      X-VECTOR...
      -4.9746182E-02  1.7868021E-01 -1.0321485E-01
  2.5719091E-02  2.3214892E-01
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    3
      ACTIVE EQUATIONS/CONSTRAINTS...
        15    2    8
      RESIDUAL L1 NORM =   4.8888325E+01
0***** CYCLE NUMBER     3
      STEP TAKEN =   1.8686117E+13
      X-VECTOR...
      -8.9805350E+06 -8.9805360E+06  1.4967559E+07
 -2.9935115E+06  1.1974047E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    4
      ACTIVE EQUATIONS/CONSTRAINTS...
        15    2    8    9
      RESIDUAL L1 NORM =   4.4000000E+01
0***** CYCLE NUMBER     4
      STEP TAKEN =   1.7592186E+13
      X-VECTOR...
      -1.1765278E+07 -1.1765278E+07  1.2182816E+07
 -2.0876900E+05  1.1974047E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    4
      ACTIVE EQUATIONS/CONSTRAINTS...
        15    2    8    9
      RESIDUAL L1 NORM =   5.5000000E+01
0***** CYCLE NUMBER     5
      STEP TAKEN =   1.6202532E+00
      X-VECTOR...
      -1.1611258E+07 -1.1611258E+07  1.1926117E+07
 -1.5742909E+05  1.1768687E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    4
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    8    9   11
      RESIDUAL L1 NORM =   8.4250000E+01
0***** CYCLE NUMBER     6
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -1.1611258E+07 -1.1611258E+07  1.1926117E+07
 -1.5742909E+05  1.1768687E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.0025000E+02
0***** CYCLE NUMBER     7
      STEP TAKEN =   1.6080402E-01
      X-VECTOR...
      -1.1611258E+07 -1.1611258E+07  1.1926117E+07
 -1.5742908E+05  1.1768688E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
        11
      RESIDUAL L1 NORM =   7.3500000E+01
0***** CYCLE NUMBER     8
      STEP TAKEN =   1.2944984E-01
      X-VECTOR...
      -1.1611258E+07 -1.1611258E+07  1.1926117E+07
 -1.5742917E+05  1.1768688E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
        11    8
      RESIDUAL L1 NORM =   6.5000000E+01
0***** CYCLE NUMBER     9
      STEP TAKEN =   2.6823521E-01
      X-VECTOR...
      -1.1611258E+07 -1.1611258E+07  1.1926117E+07
 -1.5742917E+05  1.1768688E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    3
      ACTIVE EQUATIONS/CONSTRAINTS...
        11    8   12
      RESIDUAL L1 NORM =   6.4000000E+01
0***** CYCLE NUMBER    10
      STEP TAKEN =   1.1728124E+13
      X-VECTOR...
      -1.3131956E+07 -1.3131956E+07  9.9108610E+06
  1.6105481E+06  1.1521409E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    4
      ACTIVE EQUATIONS/CONSTRAINTS...
        11    8   12    1
      RESIDUAL L1 NORM =   1.3100000E+02
0***** CYCLE NUMBER    11
      STEP TAKEN =   6.0000000E+00
      X-VECTOR...
      -7.7962315E+06 -7.7962310E+06  1.0179860E+06
  3.3891232E+06  4.4071085E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    4
      ACTIVE EQUATIONS/CONSTRAINTS...
        11    8   12    2
      RESIDUAL L1 NORM =   7.9500000E+01
0***** CYCLE NUMBER    12
      STEP TAKEN =   4.1048433E+13
      X-VECTOR...
      -3.1803530E+06 -3.1803520E+06 -6.6751455E+06
  4.9277495E+06 -1.7473975E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
        11    8   12    2    3
      RESIDUAL L1 NORM =   5.9000000E+01
0***** CYCLE NUMBER    13
      STEP TAKEN =   1.1000000E+01
      X-VECTOR...
      -1.8462325E+05 -1.8462375E+05 -3.4982340E+07
  1.7583482E+07 -1.7398858E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
        11    8    2    3    4
      RESIDUAL L1 NORM =   1.0900000E+02
0***** CYCLE NUMBER    14
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
      -1.8462325E+05 -1.8462375E+05 -3.4982340E+07
  1.7583482E+07 -1.7398858E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   1.4100000E+02
0***** CYCLE NUMBER    15
      STEP TAKEN =   3.5310346E-01
      X-VECTOR...
      -1.8462323E+05 -1.8462364E+05 -3.4982340E+07
  1.7583482E+07 -1.7398858E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         2
      RESIDUAL L1 NORM =   1.3700000E+02
0***** CYCLE NUMBER    16
      STEP TAKEN =   1.6559337E+00
      X-VECTOR...
      -1.8462339E+05 -1.8462334E+05 -3.4982340E+07
  1.7583482E+07 -1.7398858E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    7
      RESIDUAL L1 NORM =   1.1800000E+02
0***** CYCLE NUMBER    17
      STEP TAKEN =   4.5978971E+00
      X-VECTOR...
      -1.8462428E+05 -1.8462230E+05 -3.4982340E+07
  1.7583482E+07 -1.7398858E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    3
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    7    6
      RESIDUAL L1 NORM =   1.3900000E+02
0***** CYCLE NUMBER    18
      STEP TAKEN =   8.7569994E+14
      X-VECTOR...
       1.1756572E+07  1.1756571E+07 -3.3480304E+07
  1.0861867E+07 -2.2618436E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    7    6   13   15
      RESIDUAL L1 NORM =   9.7000000E+01
0***** CYCLE NUMBER    19
      STEP TAKEN =   1.0000000E+00
      X-VECTOR...
       1.0904221E+07  1.0904220E+07 -3.4332656E+07
  1.1714218E+07 -2.2618436E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         2    7    6   15   12
      RESIDUAL L1 NORM =   1.3800000E+02
0***** CYCLE NUMBER    20
      STEP TAKEN =   2.0000000E+01
      X-VECTOR...
       1.4713142E+07  1.4713140E+07 -1.9968536E+07
  2.6276970E+06 -1.7340836E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   15   12    4
      RESIDUAL L1 NORM =   1.5500000E+02
0***** CYCLE NUMBER    21
      STEP TAKEN =   1.0000000E+00
      X-VECTOR...
       1.6588984E+07  1.6588981E+07 -2.1447394E+07
  2.4292055E+06 -1.9018186E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   12    4   11
      RESIDUAL L1 NORM =   9.6000000E+01
0***** CYCLE NUMBER    22
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       1.6588984E+07  1.6588981E+07 -2.1447394E+07
  2.4292055E+06 -1.9018186E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   12   11    9
      RESIDUAL L1 NORM =   1.0500000E+02
0***** CYCLE NUMBER    23
      STEP TAKEN =   2.2500000E+00
      X-VECTOR...
       1.9311392E+07  1.9311388E+07 -1.4595018E+07
 -2.3581870E+06 -1.6953202E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   12   11    4
      RESIDUAL L1 NORM =   1.1700000E+02
0***** CYCLE NUMBER    24
      STEP TAKEN =   2.0000000E+00
      X-VECTOR...
       1.8706412E+07  1.8706408E+07 -1.6117768E+07
 -1.2943220E+06 -1.7412088E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   12   11    3
      RESIDUAL L1 NORM =   1.2600000E+02
0***** CYCLE NUMBER    25
      STEP TAKEN =   2.0000000E+00
      X-VECTOR...
       1.9695956E+07  1.9695952E+07 -2.2511222E+07
  1.4076335E+06 -2.1103588E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   12    3   15
      RESIDUAL L1 NORM =   1.7600000E+02
0***** CYCLE NUMBER    26
      STEP TAKEN =   1.7000000E+01
      X-VECTOR...
       2.2515204E+07  2.2515200E+07 -4.0726408E+07
  9.1056040E+06 -3.1620806E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   12    3   10
      RESIDUAL L1 NORM =   1.0700000E+02
0***** CYCLE NUMBER    27
      STEP TAKEN =   1.7777779E+01
      X-VECTOR...
       4.1678680E+06  4.1678660E+06 -2.3235802E+07
  9.5339680E+06 -1.3701834E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6    3   10   15
      RESIDUAL L1 NORM =   5.4000000E+01
0***** CYCLE NUMBER    28
      STEP TAKEN =   9.0000000E+00
      X-VECTOR...
      -5.6515560E+06 -5.6515570E+06 -1.7870968E+07
  1.1761262E+07 -6.1097050E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   10   15    4
      RESIDUAL L1 NORM =   1.0700000E+02
0***** CYCLE NUMBER    29
      STEP TAKEN =   2.2000000E+01
      X-VECTOR...
       3.1413264E+07  3.1413256E+07 -3.5736768E+07
  2.1617570E+06 -3.3575020E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   15    4    1
      RESIDUAL L1 NORM =   2.6300000E+02
0***** CYCLE NUMBER    30
      STEP TAKEN =   8.7500000E+00
      X-VECTOR...
       1.2114470E+07  1.2114466E+07 -2.6434456E+07
  7.1599975E+06 -1.9274466E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         7    6   15    4    8
      RESIDUAL L1 NORM =   2.2700000E+02
0***** CYCLE NUMBER    31
      STEP TAKEN =   4.4571430E+01
      X-VECTOR...
      -1.1558740E+06 -1.1558790E+06 -3.0041384E+07
  1.5598632E+07 -1.4442754E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   15    4    8   12
      RESIDUAL L1 NORM =   1.7100000E+02
0***** CYCLE NUMBER    32
      STEP TAKEN =   7.7241378E+00
      X-VECTOR...
       3.1208325E+06  3.1208285E+06 -2.5385128E+07
  1.1132150E+07 -1.4252981E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    4    8   12    2
      RESIDUAL L1 NORM =   1.0500000E+02
0***** CYCLE NUMBER    33
      STEP TAKEN =   8.0000000E+00
      X-VECTOR...
       1.2699942E+07  1.2699938E+07 -2.0665710E+07
  3.9828860E+06 -1.6682827E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    8   12    2   13
      RESIDUAL L1 NORM =   7.2000000E+01
0***** CYCLE NUMBER    34
      STEP TAKEN =   1.7500000E+00
      X-VECTOR...
       9.3276760E+06  9.3276720E+06 -2.2327152E+07
  6.4997405E+06 -1.5827415E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    8   12    2    4
      RESIDUAL L1 NORM =   1.0900000E+02
0***** CYCLE NUMBER    35
      STEP TAKEN =   1.3333334E+00
      X-VECTOR...
       1.4067246E+07  1.4067242E+07 -1.9564868E+07
  2.7488135E+06 -1.6816058E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    8   12    4    1
      RESIDUAL L1 NORM =   1.2800000E+02
0***** CYCLE NUMBER    36
      STEP TAKEN =   4.0000000E+00
      X-VECTOR...
       1.2236546E+07  1.2236542E+07 -1.9295850E+07
  3.5296545E+06 -1.5766198E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    4    1    2
      RESIDUAL L1 NORM =   1.3700000E+02
0***** CYCLE NUMBER    37
      STEP TAKEN =   6.0000002E-01
      X-VECTOR...
       1.1320098E+07  1.1320095E+07 -1.8395500E+07
  3.5377035E+06 -1.4857799E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    4    2    8
      RESIDUAL L1 NORM =   1.6100000E+02
0***** CYCLE NUMBER    38
      STEP TAKEN =   3.6000000E+01
      X-VECTOR...
       3.2882516E+07  3.2882524E+07  1.8823108E+07
 -2.5852812E+07 -7.0297050E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2    8   15
      RESIDUAL L1 NORM =   2.9800000E+02
0***** CYCLE NUMBER    39
      STEP TAKEN =   7.7241378E+00
      X-VECTOR...
       3.6016244E+07  3.6016252E+07  2.4232196E+07
 -3.0124220E+07 -5.8920250E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2    8   11
      RESIDUAL L1 NORM =   3.9400000E+02
0***** CYCLE NUMBER    40
      STEP TAKEN =   1.8000000E+01
      X-VECTOR...
       4.9500484E+07  4.9500492E+07  2.5445120E+07
 -3.7472800E+07 -1.2027685E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2   11   14
      RESIDUAL L1 NORM =   5.7900000E+02
0***** CYCLE NUMBER    41
      STEP TAKEN =   4.8500000E+01
      X-VECTOR...
       8.5835536E+07  8.5835544E+07  2.8713498E+07
 -5.7274520E+07 -2.8561026E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2   11    4
      RESIDUAL L1 NORM =   5.2000000E+02
0***** CYCLE NUMBER    42
      STEP TAKEN =   1.8933333E+02
      X-VECTOR...
       1.4913504E+07  1.4913520E+07  2.2333982E+07
 -1.8623744E+07  3.7102360E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2   11    7
      RESIDUAL L1 NORM =   4.2400000E+02
0***** CYCLE NUMBER    43
      STEP TAKEN =   1.2857142E+01
      X-VECTOR...
       2.7756528E+07  2.7756542E+07  2.3489226E+07
 -2.5622878E+07 -2.1336545E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2   11    4
      RESIDUAL L1 NORM =   3.7750000E+02
0***** CYCLE NUMBER    44
      STEP TAKEN =   2.1333334E+01
      X-VECTOR...
       3.5747744E+07  3.5747756E+07  2.4208044E+07
 -2.9977896E+07 -5.7698530E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6   12    2   11    3
      RESIDUAL L1 NORM =   2.5100000E+02
0***** CYCLE NUMBER    45
      STEP TAKEN =   3.5764706E+01
      X-VECTOR...
       3.7641788E+07  3.7641792E+07  2.0850236E+07
 -2.9246012E+07 -8.3957760E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11    3   10
      RESIDUAL L1 NORM =   3.6300000E+02
0***** CYCLE NUMBER    46
      STEP TAKEN =   2.0765957E+01
      X-VECTOR...
       3.9792680E+07  3.9792676E+07  1.7037078E+07
 -2.8414878E+07 -1.1377797E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11    3    7
      RESIDUAL L1 NORM =   2.4300000E+02
0***** CYCLE NUMBER    47
      STEP TAKEN =   1.8909090E+01
      X-VECTOR...
       3.9091104E+07  3.9091104E+07  1.8280848E+07
 -2.8685976E+07 -1.0405126E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11    3   15
      RESIDUAL L1 NORM =   1.4900000E+02
0***** CYCLE NUMBER    48
      STEP TAKEN =   5.0000000E+00
      X-VECTOR...
       3.8526136E+07  3.8526136E+07  1.9282442E+07
 -2.8904288E+07 -9.6218450E+06
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11    3    7
      RESIDUAL L1 NORM =   1.8300000E+02
0***** CYCLE NUMBER    49
      STEP TAKEN =   2.7636364E+01
      X-VECTOR...
       3.9551516E+07  3.9551512E+07  1.7464626E+07
 -2.8508068E+07 -1.1043441E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11    3   12
      RESIDUAL L1 NORM =   3.6900000E+02
0***** CYCLE NUMBER    50
      STEP TAKEN =   2.8000000E+01
      X-VECTOR...
       5.3057968E+07  5.3057968E+07  1.9570274E+07
 -3.6314120E+07 -1.6743847E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11   12    8
      RESIDUAL L1 NORM =   2.3400000E+02
0***** CYCLE NUMBER    51
      STEP TAKEN =   1.1200000E+01
      X-VECTOR...
       5.8701496E+07  5.8701496E+07  2.0450096E+07
 -3.9575796E+07 -1.9125700E+07
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         6    2   11   12    8
      RESIDUAL L1 NORM =   2.3400000E+02
0TERMINATION CODE =   3
0L1 VALUE.     =  2.3400000E+02
0X-VECTOR
  5.8701496E+07  5.8701496E+07  2.0450096E+07 -3.9575796E+07 -1.9125700E+07
0ACTIVE EQUATIONS/CONSTRAINTS
    6    2   11   12    8
1
0NUMBER OF VARIABLES =   7
0NUMBER OF EQUATIONS =   9
0NUMBER OF EQUALITY CONSTRAINTS =   0
0NUMBER OF INEQUALITY CONSTRAINTS =   5
0NUMBER OF ITERATIONS =  50
0STARTING POINT 
 
     0.000     0.000     0.000     0.000     0.000     0.000     0.000
0ARRAYS  E  AND  F
 
     8.000    32.000     8.000     0.000     0.000     0.000     0.000     2.000
     1.000    23.000    23.000     1.000     0.000     0.000     0.000     1.000
     0.000     8.000    32.000     8.000     0.000     0.000     0.000     0.000
     0.000     1.000    23.000    23.000     1.000     0.000     0.000     0.000
     0.000     0.000     8.000    32.000     8.000     0.000     0.000     0.000
     0.000     0.000     1.000    23.000    23.000     1.000     0.000     0.000
     0.000     0.000     0.000     8.000    32.000     8.000     0.000     0.000
     0.000     0.000     0.000     1.000    23.000    23.000     1.000     1.000
     0.000     0.000     0.000     0.000     8.000    32.000     8.000     2.000
     1.000    -2.000     1.000     0.000     0.000     0.000     0.000     0.000
     0.000     1.000    -2.000     1.000     0.000     0.000     0.000     0.000
     0.000     0.000     1.000    -2.000     1.000     0.000     0.000     0.000
     0.000     0.000     0.000     1.000    -2.000     1.000     0.000     0.000
     0.000     0.000     0.000     0.000     1.000    -2.000     1.000     0.000
0***** CYCLE NUMBER     0
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         5    3    4    6   10   12   14
      RESIDUAL L1 NORM =   6.0000000E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   6.0000000E+00
0***** CYCLE NUMBER     1
      STEP TAKEN =   2.0338982E-02
      X-VECTOR...
      -4.2372881E-03 -1.8361585E-03  5.6497194E-04
 -4.2372892E-04 -1.4124293E-03  4.1666668E-02
  8.4745765E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         3    4    6   10   12   14    9
      RESIDUAL L1 NORM =   3.4372883E+00
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   3.4406781E+00
0***** CYCLE NUMBER     2
      STEP TAKEN =   2.9787239E-01
      X-VECTOR...
       8.4219866E-02  4.1666668E-02 -8.8652340E-04
 -8.8652462E-04 -8.8652648E-04  4.1666668E-02
  8.4219858E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    6   10   12   14    9    1
      RESIDUAL L1 NORM =   6.8085110E-01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   6.8085110E-01
0***** CYCLE NUMBER     3
      STEP TAKEN =   6.9930041E-03
      X-VECTOR...
       8.8842981E-02  4.0501166E-02 -8.4763573E-04
 -8.7412569E-04 -9.0061611E-04  4.1666668E-02
  8.4233947E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    6   12   14    9    1    2
      RESIDUAL L1 NORM =   6.5034968E-01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   6.5034968E-01
0***** CYCLE NUMBER     4
      STEP TAKEN =   6.8965647E-03
      X-VECTOR...
       8.8793106E-02  4.0517241E-02 -8.6206751E-04
 -8.6206879E-04 -8.6207048E-04  4.0517241E-02
  8.8793114E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         4    6   12    9    1    2    8
      RESIDUAL L1 NORM =   6.2068963E-01
      PENALTY PARAMETER =   1.00E+00
      PENALTY FUNCTION VALUE =   6.2068963E-01
0***** CYCLE NUMBER     5
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       8.8793084E-02  4.0517237E-02 -8.6206850E-04
 -8.6206745E-04 -8.6206821E-04  4.0517241E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    0
      RESIDUAL L1 NORM =   6.2069046E-01
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   7.7586308E-02
0***** CYCLE NUMBER     6
      STEP TAKEN =   3.0179567E-10
      X-VECTOR...
       8.8793084E-02  4.0517237E-02 -8.6206885E-04
 -8.6206920E-04 -8.6207088E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         8    9
      RESIDUAL L1 NORM =   6.2069005E-01
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   7.7586256E-02
0***** CYCLE NUMBER     7
      STEP TAKEN =   4.1711112E-10
      X-VECTOR...
       8.8793084E-02  4.0517241E-02 -8.6206681E-04
 -8.6206675E-04 -8.6207094E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
         8    9    1    2    6
      RESIDUAL L1 NORM =   6.2068969E-01
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   7.7586219E-02
0***** CYCLE NUMBER     8
      STEP TAKEN =   1.4658029E-08
      X-VECTOR...
       8.8793069E-02  4.0517244E-02 -8.6207001E-04
 -8.6206803E-04 -8.6206949E-04  4.0517237E-02
  8.8793136E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    6
      ACTIVE EQUATIONS/CONSTRAINTS...
         8    9    1    2    6    4
      RESIDUAL L1 NORM =   6.2068957E-01
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   7.7586196E-02
0***** CYCLE NUMBER     9
      STEP TAKEN =   1.6770223E-02
      X-VECTOR...
       8.5235119E-02  4.1666664E-02 -1.9017928E-03
  1.7288525E-04 -1.9017924E-03  4.1666657E-02
  8.5235186E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         8    9    1    2    6    4   10
      RESIDUAL L1 NORM =   5.7261437E-01
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   7.5726151E-02
0***** CYCLE NUMBER    10
      STEP TAKEN =   2.4586922E-07
      X-VECTOR...
       8.5235134E-02  4.1666664E-02 -1.9018046E-03
  1.7289752E-04 -1.9018047E-03  4.1666672E-02
  8.5235141E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         8    9    1    6    4   10   14
      RESIDUAL L1 NORM =   5.7261395E-01
      PENALTY PARAMETER =   1.25E-01
      PENALTY FUNCTION VALUE =   7.5726151E-02
0***** CYCLE NUMBER    11
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       8.5235134E-02  4.1666660E-02 -1.9017933E-03
  1.7288793E-04 -1.9017961E-03  4.1666660E-02
  8.5235119E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2
      RESIDUAL L1 NORM =   5.7261473E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.3096470E-02
0***** CYCLE NUMBER    12
      STEP TAKEN =   8.3061750E-04
      X-VECTOR...
       8.5326426E-02  4.1676797E-02 -1.8527203E-03
 -1.2802519E-03 -7.0778339E-04  4.2315580E-02
  8.5338943E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    3
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2   12
      RESIDUAL L1 NORM =   7.4095511E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.1577424E-02
0***** CYCLE NUMBER    13
      STEP TAKEN =   1.6961407E-03
      X-VECTOR...
       8.5125126E-02  4.0889658E-02 -1.0618004E-03
 -1.1658865E-03 -1.2699724E-03  4.1666668E-02
  8.4603310E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    4
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2   12    9
      RESIDUAL L1 NORM =   6.8176502E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.0652578E-02
0***** CYCLE NUMBER    14
      STEP TAKEN =   3.5539768E-03
      X-VECTOR...
       8.5495397E-02  4.0576484E-02 -7.7448040E-04
 -9.4151666E-04 -1.1085527E-03  4.1666668E-02
  8.4441893E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    5
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2   12    9    4
      RESIDUAL L1 NORM =   6.7419827E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.0534348E-02
0***** CYCLE NUMBER    15
      STEP TAKEN =   9.7747566E-03
      X-VECTOR...
       8.6555094E-02  4.0604118E-02 -8.5107231E-04
 -8.7522122E-04 -8.9936989E-04  4.1666668E-02
  8.4232710E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    6
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2   12    9    4    6
      RESIDUAL L1 NORM =   6.6620064E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.0409385E-02
0***** CYCLE NUMBER    16
      STEP TAKEN =   2.1177571E-02
      X-VECTOR...
       8.8842981E-02  4.0501162E-02 -8.4763713E-04
 -8.7412592E-04 -9.0061454E-04  4.1666668E-02
  8.4233955E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2   12    9    4    6    1
      RESIDUAL L1 NORM =   6.5034968E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.0161714E-02
0***** CYCLE NUMBER    17
      STEP TAKEN =   2.8037418E-02
      X-VECTOR...
       9.0059482E-02  4.0109031E-02 -4.9560983E-04
 -1.1682247E-03 -1.8408394E-03  4.1666668E-02
  8.5174181E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
        14    2   12    9    4    1    8
      RESIDUAL L1 NORM =   6.4485985E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   1.0075935E-02
0***** CYCLE NUMBER    18
      STEP TAKEN =   6.8965587E-03
      X-VECTOR...
       8.8793114E-02  4.0517237E-02 -8.6206896E-04
 -8.6206908E-04 -8.6206885E-04  4.0517241E-02
  8.8793114E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         2   12    9    4    1    8    6
      RESIDUAL L1 NORM =   6.2068963E-01
      PENALTY PARAMETER =   1.56E-02
      PENALTY FUNCTION VALUE =   9.6982755E-03
0***** CYCLE NUMBER    19
      STEP TAKEN =   0.0000000E+00
      X-VECTOR...
       8.8793069E-02  4.0517241E-02 -8.6206966E-04
 -8.6206850E-04 -8.6206943E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    1
      ACTIVE EQUATIONS/CONSTRAINTS...
         9
      RESIDUAL L1 NORM =   6.2069011E-01
      PENALTY PARAMETER =   1.95E-03
      PENALTY FUNCTION VALUE =   1.2122875E-03
0***** CYCLE NUMBER    20
      STEP TAKEN =   3.6189013E-10
      X-VECTOR...
       8.8793069E-02  4.0517241E-02 -8.6206925E-04
 -8.6206920E-04 -8.6206908E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    2
      ACTIVE EQUATIONS/CONSTRAINTS...
         9   12
      RESIDUAL L1 NORM =   6.2069005E-01
      PENALTY PARAMETER =   1.95E-03
      PENALTY FUNCTION VALUE =   1.2122852E-03
0***** CYCLE NUMBER    21
      STEP TAKEN =   3.5670693E-09
      X-VECTOR...
       8.8793069E-02  4.0517241E-02 -8.6206896E-04
 -8.6206896E-04 -8.6206885E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    3
      ACTIVE EQUATIONS/CONSTRAINTS...
         9   12    6
      RESIDUAL L1 NORM =   6.2068999E-01
      PENALTY PARAMETER =   1.95E-03
      PENALTY FUNCTION VALUE =   1.2122851E-03
0***** CYCLE NUMBER    22
      STEP TAKEN =   8.4548120E-08
      X-VECTOR...
       8.8793069E-02  4.0517248E-02 -8.6207106E-04
 -8.6206960E-04 -8.6206809E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    6
      ACTIVE EQUATIONS/CONSTRAINTS...
         9   12    6    8    2    1
      RESIDUAL L1 NORM =   6.2068969E-01
      PENALTY PARAMETER =   1.95E-03
      PENALTY FUNCTION VALUE =   1.2122846E-03
0***** CYCLE NUMBER    23
      STEP TAKEN =   7.1594093E-07
      X-VECTOR...
       8.8793077E-02  4.0517244E-02 -8.6206914E-04
 -8.6206902E-04 -8.6206879E-04  4.0517237E-02
  8.8793129E-02
      NUMBER OF ACTIVE EQUATIONS/CONSTRAINTS    7
      ACTIVE EQUATIONS/CONSTRAINTS...
         9   12    6    8    2    1    4
      RESIDUAL L1 NORM =   6.2068963E-01
      PENALTY PARAMETER =   1.95E-03
      PENALTY FUNCTION VALUE =   1.2122844E-03
0TERMINATION CODE =   1
0L1 VALUE.     =  6.2068993E-01
0X-VECTOR
  8.8793069E-02  4.0517248E-02 -8.6206850E-04 -8.6206966E-04 -8.6206594E-04
  4.0517237E-02  8.8793121E-02
0ACTIVE EQUATIONS/CONSTRAINTS
    9   12    6    8    2    1    4
1
ian@ian-HP-Stream-Notebook-PC-11:~/563dir$ cat 5631.f
      INTEGER I,IFL,INDX(50),MXS,NACT,NEQC,NEQNS,NIQC
      INTEGER NVARS,NXINPT,NXOUTP,NXPRMP
      LOGICAL PSW
      REAL E(500),EL1N,F(50),RES(50),X(10),W(570)
C
C     ***************
C     A SIMPLE DRIVER PROGRAM FOR   CL1.
C
C     NXINPT,NXOUTP,NXPRMP   ARE UNIT NUMBERS FOR
C     INPUT, OUTPUT, AND PROMPTING RESPECTIVELY.
C
C     PROBLEM DIMENSIONS ARE READ BY  DATFL1.
C     DATA ARRAYS ARE READ BY  DATFL2.
C     SEVERAL SUCCESSIVE PROBLEMS CAN BE SOLVED.
C
C     A TYPICAL INPUT SEQUENCE IS AS FOLLOWS...
C
C       NEW                           (NEW=1)
C       NVARS NEQNS NEQC NIQC
C       MXS
C       PSW                           (1 FOR PRINTING)
C       E(1,1)   THROUGH    E(N,1)
C       F(1)
C       .
C       .
C       .
C       E(1,K)   THROUGH    E(N,K)
C       F(K)
C       E(1,K+1) THROUGH    E(N,K+1)
C       F(K+1)
C       .
C       .
C       .
C       E(1,L)   THROUGH    E(N,L)
C       F(L)
C       E(1,L+1) THROUGH    E(N,L+1)
C       F(L+1)
C       .
C       .
C       .
C       E(1,M)   THROUGH    E(N,M)
C       F(M)
C       X(1)     THROUGH    X(N)
C       NEW                           (1 FOR MORE, 0 TO STOP)
C
C     WHERE
C       K STANDS FOR NEQNS,
C       L STANDS FOR NEQNS+NEQC,
C       M STANDS FOR NEQNS+NEQC+NIQC.
C
C     FOR FURTHER INFORMATION, CONSULT THE
C     COMMENTS IN  DATFL1,DATFL2, OR RUN THE
C     PROGRAM IN INTERACTIVE MODE AND RESPOND
C     TO THE PROMPTS PUT OUT ON UNIT  NXPRMP.
C     ***************
C
      DATA NXINPT /5/ , NXOUTP /6/ , NXPRMP /7/
C
   10 CONTINUE
C
C     ***************
C     READ DATA
C     ***************
C
      WRITE (NXOUTP,1000)
      CALL DATFL1(NEQNS,NEQC,NIQC,NVARS,MXS,PSW,NXINPT,NXOUTP,NXPRMP)
      CALL DATFL2(NEQNS,NEQC,NIQC,NVARS,E,X,F,NXINPT,NXOUTP,NXPRMP)
C
C     ***************
C     ENTER  CL1
C     ***************
C
      CALL CL1(NEQNS,NEQC,NIQC,NVARS,NACT,IFL,MXS,PSW,
     *         E,NVARS,X,F,EL1N,RES,INDX,W)
C
C     ***************
C     PRINT RESULTS, AND
C     GO BACK FOR NEW PROBLEM
C     ***************
C
      WRITE (NXOUTP,1010) IFL,EL1N
      WRITE (NXOUTP,1020)  (X(I),I=1,NVARS)
      IF (NACT .LE. 0)  GO TO 20
        WRITE (NXOUTP,1030)
        WRITE (NXOUTP,1040)  (INDX(I),I=1,NACT)
   20 CONTINUE
      GO TO 10
C
C     ***************
C     FORMATS
C     ***************
C
 1000 FORMAT(1H1)
 1010 FORMAT(21H0TERMINATION CODE =   ,I2/
     *       16H0L1 VALUE.     = ,1PE15.7/
     *       9H0X-VECTOR )
 1020 FORMAT(5(1PE15.7))
 1030 FORMAT(29H0ACTIVE EQUATIONS/CONSTRAINTS)
 1040 FORMAT(10I5)
      END
      SUBROUTINE DATFL1(NEQNS,NEQC,NIQC,NVARS,MXS,PSW,
     *                  NXINPT,NXOUTP,NXPRMP)
C
      INTEGER K,MXS,NEQNS,NEQC,NIQC,NVARS
      INTEGER NXINPT,NXOUTP,NXPRMP
      LOGICAL PSW
C
C     ***************
C     READ  1  TO EXECUTE A PROBLEM.
C     READ  0  TO STOP.
C     FORMAT  I1
C     ***************
C
      WRITE (NXPRMP,1000)
      READ (NXINPT,1010)  K
      IF (K .EQ. 0)  STOP
C
C     ***************
C     READ  NVARS,NEQNS,NEQC,NIQC.
C     FORMAT  4I5
C     ***************
C
      WRITE (NXPRMP,1020)
      READ (NXINPT,1030) NVARS,NEQNS,NEQC,NIQC
      WRITE (NXOUTP,1040) NVARS,NEQNS,NEQC,NIQC
C
C     ***************
C     READ  MXS.
C     FORMAT  I5
C     ***************
C
      WRITE (NXPRMP,1050)
      READ (NXINPT,1030) MXS
      WRITE (NXOUTP,1060) MXS
C
C     ***************
C     READ  1  TO GET INTERMEDIATE PRINTING.
C     READ  0  TO GET NONE.
C     FORMAT  I1
C     ***************
C
      WRITE (NXPRMP,1070)
      READ (NXINPT,1010) K
      PSW = .FALSE.
      IF (K .EQ. 1)  PSW = .TRUE.
      RETURN
C
C     ***************
C     FORMATS
C     ***************
C
 1000 FORMAT(27H0NEW PROBLEM (1=YES,0=STOP))
 1010 FORMAT(I1)
 1020 FORMAT(20H0NUMBER OF VARIABLES/
     *        11X,9HEQUATIONS/11X,20HEQUALITY CONSTRAINTS/
     *        11X,28HINEQUALITY CONSTRAINTS (4I5))
 1030 FORMAT(10I5)
 1040 FORMAT(23H0NUMBER OF VARIABLES = ,I3/
     *       23H0NUMBER OF EQUATIONS = ,I3/
     *       34H0NUMBER OF EQUALITY CONSTRAINTS = ,I3/
     *       36H0NUMBER OF INEQUALITY CONSTRAINTS = ,I3)
 1050 FORMAT(26H0NUMBER OF ITERATIONS (I5))
 1060 FORMAT(24H0NUMBER OF ITERATIONS =  ,I3)
 1070 FORMAT(35H0INTERMEDIATE PRINTING (1=YES,0=NO))
      END
      SUBROUTINE DATFL2(NEQNS,NEQC,NIQC,NVARS,E,X,F,
     *                  NXINPT,NXOUTP,NXPRMP)
C
      INTEGER I,IX,J,NEQNS,NEQC,NIQC
      INTEGER NVARS,NXINPT,NXOUTP,NXPRMP
      REAL E(NVARS,1),F(1),X(1)
C
      IF (NEQNS .LE. 0)  GO TO 20
C
C     ***************
C     IF  NEQNS .GT. 0, READ IN
C     EQUATION COEFFICIENTS
C     FOLLOWED BY RIGHT HAND SIDES.
C     FORMAT  8(F8.3,2X)
C     ***************
C
        WRITE (NXPRMP,1000)
        DO 10 I=1,NEQNS
          WRITE (NXPRMP,1010)  NVARS,I
          READ (NXINPT,1020)  (E(J,I),J=1,NVARS)
          WRITE (NXPRMP,1030)
          READ (NXINPT,1020) F(I)
   10   CONTINUE
   20 CONTINUE
      IF (NEQC .LE. 0)  GO TO 40
C
C     ***************
C     IF  NEQC .GT. 0, READ IN
C     EQUALITY-CONSTRAINT COEFFICIENTS
C     FOLLOWED BY RIGHT HAND SIDES.
C     FORMAT  8(F8.3,2X)
C     ***************
C
        WRITE (NXPRMP,1040)
        DO 30 I=1,NEQC
          IX = I + NEQNS
          WRITE (NXPRMP,1050)  NVARS,I
          READ (NXINPT,1020)  (E(J,IX),J=1,NVARS)
          WRITE (NXPRMP,1030)
          READ (NXINPT,1020)  F(IX)
   30   CONTINUE
   40 CONTINUE
      IF (NIQC .LE. 0)  GO TO 60
C
C     ***************
C     IF  NIQC .GT. 0, READ IN
C     INEQUALITY-CONSTRAINT COEFFICIENTS
C     FOLLOWED BY RIGHT HAND SIDES.
C     FORMAT  8(F8.3,2X)
C     ***************
C
        WRITE (NXPRMP,1060)
        DO 50 I=1,NIQC
          IX = I + NEQNS + NEQC
          WRITE (NXPRMP,1070)  NVARS,I
          READ (NXINPT,1020)  (E(J,IX),J=1,NVARS)
          WRITE (NXPRMP,1030)
          READ (NXINPT,1020)  F(IX)
   50   CONTINUE
   60 CONTINUE
C
C     ***************
C     READ IN STARTING VALUES FOR  X.
C     FORMAT  8(F8.3,2X)
C     ***************
C
      WRITE (NXPRMP,1080)
      READ (NXINPT,1020)  (X(I),I=1,NVARS)
      WRITE (NXOUTP,1090)
C
C     ***************
C     ECHO ALL ARRAYS
C     ***************
C
      WRITE (NXOUTP,1100)  (X(I),I=1,NVARS)
      WRITE (NXOUTP,1110)
      IX = NEQNS + NEQC + NIQC
      DO 70 I=1,IX
        WRITE (NXOUTP,1100)  (E(J,I),J=1,NVARS),F(I)
   70 CONTINUE
      RETURN
C
C     ***************
C     FORMATS
C     ***************
C
 1000 FORMAT(18H0SPECIFY EQUATIONS)
 1010 FORMAT(1H ,I3,28H  COEFFICIENTS FOR EQUATION ,
     *       I3,13H (8(F8.3,2X)))
 1020 FORMAT(8(F8.3,2X))
 1030 FORMAT(24H RIGHT-HAND SIDE  (F8.3))
 1040 FORMAT(29H0SPECIFY EQUALITY CONSTRAINTS )
 1050 FORMAT(1H ,I3,
     *       39H  COEFFICIENTS FOR EQUALITY CONSTRAINT  ,
     *       I3,10H (F8.3,2X))
 1060 FORMAT(31H0SPECIFY INEQUALITY CONSTRAINTS )
 1070 FORMAT(1H ,I3,
     *       41H  COEFFICIENTS FOR INEQUALITY CONSTRAINT ,
     *       I3,10H (F8.3,2X))
 1080 FORMAT(37H0SPECIFY STARTING POINT  (8(F8.3,2X)))
 1090 FORMAT(16H0STARTING POINT /1H )
 1100 FORMAT(8(2X,F8.3))
 1110 FORMAT(18H0ARRAYS  E  AND  F/1H )
      END
      SUBROUTINE CL1(NEQNS,NEQC,NIQC,NVARS,NACT,IFL,MXS,PSW,
     *               E,NER,X,F,EL1N,RES,INDX,W)
C
      INTEGER IFL,INDX(1),MXS,NACT,NEQC,NEQNS,NER,NIQC,NVARS
      LOGICAL PSW
      REAL E(NER,1),EL1N,F(1),RES(1),W(1),X(1)
C
C
C     **********************************************************
C     A PROGRAM FOR THE SOLUTION IN THE   L1  SENSE
C     OF A LINEAR-EQUATION SYSTEM
C     (WITH OR WITHOUT LINEAR CONSTRAINTS).
C     RICHARD H. BARTELS AND ANDREW R. CONN.
C     LATEST UPDATE .... 13 APRIL, 1980.
C
C     DEVELOPMENT OF THIS PROGRAM WAS SUPPORTED
C     IN PART BY U. S. NSF GRANT DCR75-07817
C     AND BY FUNDS FROM THE NATIONAL BUREAU OF STANDARDS,
C     AND IN PART BY CANADIAN NRC GRANT A8639.
C
C     +++++ PARAMETERS +++++
C     ----------------------------------------------------------
C                           INPUT
C     NAME  TYPE  SUBSCRPT  OUTPUT        DESCRIPTION
C                           SCRATCH
C     ..........................................................
C     NEQNS INT.    NONE      IN      NUMBER OF EQUATIONS
C                                     (MAY BE ZERO)
C
C     NEQC  INT.    NONE      IN      NUMBER OF EQUALITY
C                                     CONSTRAINTS
C                                     (MAY BE ZERO)
C
C     NIQC  INT.    NONE      IN      NUMBER OF INEQUALITY
C                                     CONSTRAINTS
C                                     (MAY BE ZERO)
C
C     NVARS INT.    NONE      IN      NUMBER OF VARIABLES
C
C     NACT  INT.    NONE      OUT     NUMBER OF ACTIVE
C                                     EQUATIONS/CONSTRAINTS
C                                     AT TERMINATION
C                                     (IF ANY, THEIR ASSOCIATED
C                                      COLUMN POSITIONS IN  E  WILL
C                                      BE LISTED IN  INDX(1)
C                                      THROUGH  INDX(NACT) )
C
C     IFL   INT.    NONE      OUT     TERMINATION CODE
C                                     (SEE BELOW)
C
C     MXS   INT.    NONE      IN      MAXIMUM NUMBER OF STEPS
C                                     ALLOWED
C
C     PSW   LOGIC.  NONE      IN      PRINT SWITCH
C                                     (SEE BELOW)
C
C     E     REAL     2        IN      EQUATION/CONSTRAINT MATRIX
C                                     THE FIRST  NEQNS  COLUMNS
C                                     (SEE NOTE BELOW) SPECIFY
C                                     EQUATIONS, THE REMAINING
C                                     COLUMNS (IF ANY) SPECIFY
C                                     CONSTRAINTS.
C
C     NER   INT.    NONE      IN      ROW DIMENSION OF E
C
C     X     REAL     1        IN      STARTING VALUES FOR THE
C                                     UNKNOWNS (USE ZEROS IF NO
C                                     GUESS IS AVAILABLE)
C                             OUT     TERMINATION VALUES FOR
C                                     THE UNKNOWNS
C
C     F     REAL     1        IN      EQUATION/CONSTRAINT
C                                     RIGHT-HAND SIDES
C
C     EL1N  REAL    NONE      OUT     L1 NORM OF EQUATION
C                                     RESIDUALS AT TERMINATION
C
C     RES   REAL     1        OUT     EQUATION/CONSTRAINT
C                                     RESIDUALS AT TERMINATION
C
C     INDX  INT.     1        OUT     INDEX VECTOR USED TO RECORD
C                                     THE ORDER IN WHICH THE COLUMNS
C                                     OF  E  ARE BEING PROCESSED
C
C     W     REAL     1        SCR.    WORKING STORAGE
C     ----------------------------------------------------------
C
C     +++++ PURPOSE +++++
C     ----------------------------------------------------------
C     THIS SUBROUTINE SOLVES THE   NEQNS BY NVARS
C     SYSTEM OF EQUATIONS
C
C                       (A-TRANSPOSE) * X   ==   B
C
C     SUBJECT TO THE  NEQC   CONSTRAINTS
C
C                       (G-TRANSPOSE) * X  .EQ.  H
C
C     AND THE  NIQC  INEQUALITY CONSTRAINTS
C
C                       (C-TRANSPOSE) * X  .GE.  D
C
C     FOR THE UNKNOWNS  X(1),...,X(NVARS).
C
C     THE PROBLEM MUST BE WELL-POSED, NONTRIVIAL
C     AND OVERDETERMINED IN THE SENSE THAT
C
C                          NVARS .GE. 1
C                          NEQNS .GE. 0
C                          NEQC  .GE. 0
C                          NIQC  .GE. 0
C               NEQNS+NEQC+NIQC  .GE. NVARS.
C
C     FURTHER, NO COLUMN OF  A, G  OR  C  SHOULD BE ZERO.
C     IF THESE CONDITIONS ARE NOT MET, THE PROGRAM
C     WILL TERMINATE WITHOUT PERFORMING ANY SUBSTANTIVE
C     COMPUTATIONS.
C
C     A POINT  X  IS A SOLUTION IF IT MINIMIZES THE EQUATION
C     RESIDUALS FROM AMONG ALL POINTS WHICH SATISFY THE
C     CONSTRAINTS.  AT ANY (NONDEGENERATE) SOLUTION
C     THERE WILL BE  NACT  EQUATIONS AND CONSTRAINTS
C     WHOSE RESIDUALS
C
C          (A(I)-TRANSPOSE) * X - B(I)
C
C          (G(I)-TRANSPOSE) * X - H(I)
C
C     AND
C
C          (C(I)-TRANSPOSE) * X - D(I)
C
C     ARE ZERO.
C
C     THE COLUMNS OF  (A,G,C)  CORRESPONDING TO THE ZERO RESIDUALS
C     ARE REFERRED TO AS  ACTIVE COLUMNS  THROUGHOUT THIS LISTING.
C     THE NUMBERS OF THE ACTIVE COLUMNS ARE MAINTAINED AS THE
C     ENTRIES  1,...,NACT  OF THE ARRAY  INDX.
C
C     A SOLUTION  X  IS FOUND BY MINIMIZING A PIECEWISE
C     LINEAR PENALTY FUNCTION FORMED FROM THE  L1
C     NORM OF THE EQUATION RESIDUALS AND THE SUM OF THE
C     INFEASIBILITIES IN THE CONSTRAINTS.
C     THE MINIMIZATION PROCEEDS IN A STEP-BY-STEP
C     FASHION, TERMINATING AFTER A FINITE NUMBER OF STEPS.
C
C     NOTE THAT  A, G  AND  C  APPEAR TRANSPOSED IN THE
C     PROBLEM FORMULATION.  HENCE IT IS THE COLUMNS OF  (A,G,C)
C     WHICH DEFINE THE EQUATIONS AND CONSTRAINTS RESPECTIVELY.
C
C     THE ARRAY  E  IS A COMPOSITE OF   A, G AND C
C     AND  F  IS A COMPOSITE OF  B, H  AND  D.
C     E  SHOULD CONTAIN  A  AS ITS FIRST  NEQNS  COLUMNS.
C     IT SHOULD CONTAIN  G  AS ITS NEXT  NEQC  COLUMNS AND
C     CONTAIN  C  AS ITS REMAINING  NIQC  COLUMNS.
C     SIMILARLY  F  SHOULD CONTAIN  B  AS ITS FIRST
C     NEQNS  COMPONENTS,  H  AS ITS NEXT  NEQC  COMPONENTS
C     AND  D  AS ITS LAST  NIQC  COMPONENTS.
C     ----------------------------------------------------------
C
C     +++++ ARRAYS +++++
C     ----------------------------------------------------------
C     E  IS TO BE DIMENSIONED AT LEAST    N  BY  M,
C     X                       AT LEAST    N,
C     F                       AT LEAST    M,
C     RES                     AT LEAST    M,
C     INDX                    AT LEAST    M,
C     W                       AT LEAST    ((3*N*N+11*N+2)/2) + (2*M).
C
C                                         WHERE  N = NVARS  AND
C                                         M = NEQNS+NEQC+NIQC
C     ----------------------------------------------------------
C
C     +++++ INITIALIZATION +++++
C     ----------------------------------------------------------
C     THE USER MUST INITIALIZE
C
C          NEQNS,NEQC,NIQC,NVARS,MXS,PSW,E,NER,X,F .
C
C     THE FOLLOWING ARE SET BY  CL1
C     AND DO NOT REQUIRE INITIALIZATION
C
C          NACT,INDX,RES .
C
C     THE ARRAY  W  IS USED AS SCRATCH SPACE.
C     ----------------------------------------------------------
C
C     +++++ TERMINATION CODES AND INTERMEDIATE PRINTING +++++
C     ----------------------------------------------------------
C     MXS  SETS A LIMIT ON THE NUMBER OF MINIMIZATION STEPS TO BE
C     TAKEN.
C
C     UPON TERMINATION  IFL  WILL BE SET ACCORDING TO
C     THE FOLLOWING CODE ...
C
C             IFL = 1 .... SUCCESSFUL TERMINATION.
C
C             IFL = 2 .... UNSUCCESSFUL TERMINATION.
C                          CONSTRAINTS CANNOT BE SATISFIED.
C                          PROBLEM IS INFEASIBLE.
C
C             IFL = 3 .... LIMIT IMPOSED BY  MXS  REACHED
C                          WITHOUT FINDING A SOLUTION.
C
C             IFL = 4 .... PROGRAM ABORTED.
C                          NUMERICAL DIFFICULTIES
C                          DUE TO ILL-CONDITIONING.
C
C             IFL = 5 .... NEQNS, NVARS, NEQC AND/OR
C                          NIQC  HAVE IMPROPER VALUES
C                          OR  E  CONTAINS A ZERO COLUMN.
C
C     IN ALL CASES THE OUTPUT PARAMETERS  X,EL1N AND RES
C     WILL CONTAIN THE VALUES WHICH THEY REACHED AT TERMINATION.
C
C     INTERMEDIATE PRINTING WILL BE TURNED OFF IF  PSW = .FALSE. .
C     ON THE OTHER HAND,  DETAILS OF EACH MINIMIZATION CYCLE
C     WILL BE PRINTED IF  PSW  IS SET TO  .TRUE.
C     ----------------------------------------------------------
C
C     +++++ REMARKS AND USER CAUTIONS +++++
C     ----------------------------------------------------------
C       1.  BEYOND SOME PRECAUTIONARY STEPS TAKEN IN
C           CERTAIN DIVISIONS, NO SPECIAL
C           OVERFLOW/UNDERFLOW PROTECTION IS PROVIDED.
C       2.  ALL TOLERANCES FOR CHECKING ZEROS AND LINEAR
C           DEPENDENCIES ARE DETERMINED FROM THE QUANTITY
C           EPS  WHICH APPEARS IN DATA DECLARATIONS IN  CL1
C           AND SEVERAL OF ITS SUBROUTINES.   EPS  CAN BE SET TO THE
C           LEAST POSITIVE NUMBER SATISFYING  (1.0 + EPS) .GT. 1.0
C           IN THE PRECISION OF ARITHMETIC BEING USED.  WITH THIS
C           SETTING,  CL1  USES AN EXTREMELY STRICT ZERO TOLERANCE.
C           FOR A MORE FORGIVING VERSION OF  CL1,  EPS MAY BE REMOVED
C           FROM THE DATA DEFINITIONS AND INCLUDED AS A USER-SPECIFIED
C           ZERO-TESTING PARAMETER IN THE ARGUMENT LIST.  IN SUCH AN
C           EVENT, IF THE PROBLEM DATA IS GIVEN TO  NDIG  SIGNIFICANT
C           DIGITS, THEN  10.0**(-NDIG)  IS A REASONABLE CHOICE
C           FOR THE VALUE OF  EPS.
C           OVERFLOW CHECKING PRIOR TO DIVISION IS
C           DONE USING THE QUANTITY  BIG,  ALSO SPECIFIED AS DATA.
C           BIG  SHOULD BE THE LARGEST REPRESENTABLE FLOATING
C           POINT NUMBER.
C       3.  THIS IS A SINGLE PRECISION VERSION OF  CL1.
C           TO CHANGE THIS CODE INTO DOUBLE PRECISION ...
C           A.  CHANGE ALL OCCURRENCES OF  - REAL -
C               DECLARATIONS TO  - DOUBLE PRECISION -
C           B.  CHANGE ALL OCCURRENCES OF  - SYSTEM ROUTINES -
C               (AS LISTED IN THE HEADING OF EACH SUBROUTINE)
C               TO THEIR CORRESPONDING DOUBLE PRECISION VERSIONS
C           C.  CHANGE ALL OCCURRENCES OF THE STRINGS  E+  AND  E-
C               TO  D+  AND  D-  RESPECTIVELY
C           D.  CHANGE ALL BASIC LINEAR ALGEBRA ROUTINES  (BLAS)
C               TO THEIR DOUBLE-PRECISION EQUIVALENTS
C           E.  BOTH  EPS  AND  BIG  WILL HAVE TO BE CHANGED
C               TO THEIR DOUBLE PRECISION EQUIVALENTS
C           F.  THE REFERENCES TO  - IFIX -  IN SUBROUTINES
C               - RESID -  AND  - GETV -  MUST BE CHANGED
C               FROM THE FORM
C                    IFIX(FLOAT(K)*UNIF(...))
C               TO THE FORM
C                    IFIX(FLOAT(K)*SNGL(UNIF(...)))
C           G.  THE REFERENCE TO  - FLOAT -  IN SUBROUTINE
C               - RESID -  MUST BE CHANGED FROM THE FORM
C                    FLOAT(...)
C               TO THE FORM
C                    DBLE(FLOAT(...))
C           H.  REMOVE THESE COMMENT CARDS  (3., 3.A.-3.H.).
C     ----------------------------------------------------------
C
C     +++++ PORTABILITY NOTE +++++
C     ----------------------------------------------------------
C     THE INTRINSIC FORTRAN FUNCTIONS WHICH HAVE BEEN USED
C     ARE NOT DECLARED IN TYPE STATEMENTS IN  CL1  OR ANY
C     OF ITS SUBROUTINES.  A VERY FEW FORTRAN COMPILERS DO
C     REQUIRE SUCH TYPING, AND THEY WILL FAIL QUITE EXPLICITLY
C     BECAUSE OF OUR OMISSION.  ON THE OTHER HAND, SOME FEW
C     OTHER FORTRAN COMPILERS DO NOT EXPECT SUCH TYPING, AND
C     THEIR MODE OF FAILURE IS MORE SUBTILE.  THESE COMPILERS
C     MAY ATTACH THE NAMES OF THE INTRINSIC FUNCTIONS TO THE
C     LIST OF SUBROUTINES TO BE OBTAINED FROM EXTERNAL SOURCES,
C     CAUSING A FAULT IN LIBRARY SEARCHING DURING PRE-EXECUTION
C     LOADING.   SINCE THERE IS NO WAY TO WIN THE GAME, WE
C     HAVE OPTED TO LOSE IT, IN THOSE FEW CASES WHERE WE MUST,
C     IN THE MOST IMMEDIATE AND OBVIOUS FASHION.
C     NOTE -- INTRINSIC FUNCTIONS ARE LISTED AMONG THE
C                  - SYSTEM ROUTINES -
C              AT THE START OF EACH SUBPROGRAM.
C
C     THE QUANTITIES  EPS  AND  BIG  WHICH ARE TO BE FOUND
C     IN DATA STATEMENTS ARE MACHINE DEPENDENT.  SEE  2.  UNDER
C     + REMARKS AND USER CAUTIONS +.
C
C     MANY OF THE COMPUTATIONS IN THIS PROGRAM ARE
C     DEFINED IN TERMS OF THE  BASIC LINEAR ALGEBRA
C     SUBROUTINES  (BLAS).  REFER TO THE APPENDIX
C     OF THIS CODE FOR FURTHER INFORMATION.
C     **********************************************************
C
      INTEGER DDX,GRDX,ICYC,IADDC,IDELC
      INTEGER PX,PTEX,RRX,TOPX,ZZX
      REAL ALPHA,AMAG,CGMAG,PEN,PENPAR
C
C     ////////////////  BEGIN PROGRAM  /////////////////////////
C
      CALL SETUP
     *                    (NEQNS,NEQC,NIQC,NVARS,DDX,GRDX,PX,PTEX,RRX,
     *                     TOPX,ZZX,ICYC,IFL,E,NER,AMAG,CGMAG,PENPAR)
C
C
   10   CONTINUE
C
          CALL NEWPEN
     *                    (IADDC,IDELC,NACT,NEQNS,NEQC,
     *                     NIQC,NVARS,IFL,E,NER,X,F,RES,
     *                     W(PTEX),ALPHA,PENPAR,INDX)
C
   20     CONTINUE
C
            CALL UPDATE
     *                    (IADDC,IDELC,NACT,NEQNS,NEQC,NIQC,NVARS,
     *                     ICYC,IFL,MXS,E,NER,X,F,RES,
     *                     W(GRDX),EL1N,PEN,PENPAR,
     *                     INDX,W(ZZX),NVARS,W(DDX),W(RRX),
     *                     W(TOPX))
            CALL MONIT
     *                    (NACT,NEQC,NIQC,NVARS,ICYC,PSW,X,
     *                     ALPHA,EL1N,PEN,PENPAR,INDX)
            CALL FINDP
     *                    (IDELC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL,E,NER,
     *                     X,F,RES,W(GRDX),W(PX),EL1N,AMAG,CGMAG,
     *                     PENPAR,INDX,W(ZZX),NVARS,W(DDX),
     *                     W(RRX),W(TOPX))
            CALL STEP
     *                    (IADDC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL,E,NER,
     *                     X,RES,W(GRDX),W(PX),W(PTEX),ALPHA,PENPAR,
     *                     INDX,W(TOPX))
C
            IF (IFL .EQ. 0)  GO TO 20
C
          IF (IFL .EQ. 2
     *           .AND. (CGMAG + PENPAR*AMAG) .NE. CGMAG)  GO TO 10
C
      RETURN
      END
C
C     ---------------
C     SECOND LEVEL SUBROUTINES --
C          SETUP,NEWPEN,UPDATE,MONIT,FINDP,STEP,REFINE
C     ---------------
C
      SUBROUTINE SETUP
     *                    (NEQNS,NEQC,NIQC,NVARS,DDX,GRDX,PX,PTEX,RRX,
     *                     TOPX,ZZX,ICYC,IFL,E,NER,AMAG,CGMAG,PENPAR)
C
      INTEGER DDX,GRDX,ICYC,IFL,NEQC,NEQNS,NER
      INTEGER NIQC,NVARS,PX,PTEX,RRX,TOPX,ZZX
      REAL AMAG,CGMAG,E(NER,1),PENPAR
C
C     ***************
C     CL1  VERSION.
C
C     SET UP THE PROGRAM
C     PARAMETERS AND INDICES.
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C     +++++++++++++++
C
      INTEGER I,J,NCOLS,NQNP1
      REAL OCT,TMP,ZERO
C
      DATA OCT /8.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
C     ***************
C     CHECK VALIDITY OF PROBLEM DIMENSIONS
C     ***************
C
      NCOLS = NEQNS + NEQC + NIQC
      IF (NVARS .GE. 1
     *       .AND. NEQC .GE. 0
     *        .AND. NIQC .GE. 0
     *         .AND. NEQNS .GE. 0
     *          .AND. NCOLS .GE. NVARS
     *           .AND. NER .GE. NVARS)  GO TO 10
        IFL = 5
        GO TO 100
   10 CONTINUE
C
C     ***************
C     SET UP INDICES FOR THE TEMPORARY STORAGE VECTOR  W.
C     ***************
C
      NQNP1 = NEQNS + 1
      GRDX = 1
      PX = GRDX + NVARS
      PTEX = PX + NVARS
      DDX = PTEX + NCOLS
      RRX = DDX + NVARS
      ZZX = RRX + (((NVARS + 1)*(NVARS + 2))/2)
      TOPX = ZZX + NVARS*NVARS
C
C     ***************
C     AMAG  IS A ROUGH ESTIMATE OF THE NORM OF  A.
C     CGMAG  IS A ROUGH ESTIMATE OF THE NORM OF  (G,C).
C     TOGETHER THEY ARE USED TO DETERMINE WHEN THE
C     PENALTY PARAMETER IS TOO SMALL AND WHEN THE
C     RESTRICTED GRADIENT IS ZERO.
C     ***************
C
      AMAG = ZERO
      CGMAG = ZERO
      IF (1 .GT. NEQNS)  GO TO 50
        DO 40 J=1,NEQNS
          TMP = ZERO
          DO 20 I=1,NVARS
            TMP = TMP + ABS(E(I,J))
  20     CONTINUE
          IF (TMP .GT. ZERO)  GO TO 30
            IFL = 5
            GO TO 100
   30     CONTINUE
          IF (TMP .GT. AMAG)  AMAG = TMP
   40   CONTINUE
   50 CONTINUE
      IF (NQNP1 .GT. NCOLS)  GO TO 90
        DO 80 J=NQNP1,NCOLS
          TMP = ZERO
          DO 60 I=1,NVARS
            TMP = TMP + ABS(E(I,J))
  60     CONTINUE
          IF (TMP .GT. ZERO)  GO TO 70
            IFL = 5
            GO TO 100
   70     CONTINUE
          IF (TMP .GT. CGMAG)  CGMAG = TMP
   80   CONTINUE
   90 CONTINUE
C
C     ***************
C     INITIALIZE  IFL,ICYC,PENPAR
C     ***************
C
      IFL = 2
      ICYC = -1
      PENPAR = OCT
  100 CONTINUE
      RETURN
      END
      SUBROUTINE NEWPEN
     *                    (IADDC,IDELC,NACT,NEQNS,NEQC,NIQC,
     *                     NVARS,IFL,E,NER,X,F,RES,PTE,
     *                     ALPHA,PENPAR,INDX)
C
      INTEGER IADDC,IDELC,IFL,INDX(1),NACT
      INTEGER NEQC,NEQNS,NER,NIQC,NVARS
      REAL ALPHA,E(NER,1),F(1),PENPAR,PTE(1),RES(1),X(1)
C
C     ***************
C     CL1  VERSION.
C
C     BEGIN A ROUND OF MINIMIZATION STEPS
C     WITH A NEW PENALTY PARAMETER VALUE.
C     ***************
C
C     +++++++++++++++
C     BLAS  SDOT
C     +++++++++++++++
C
      INTEGER I,NCOLS
      REAL OCT,ONE,ZERO
C
      REAL SDOT
C
      DATA ZERO /0.0E+00/
      DATA ONE /1.0E+00/
      DATA OCT /8.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
C     ***************
C     SET PENALTY PARAMETER VALUE.
C     ERASE RECORD OF ACTIVE EQUATION/CONSTRAINTS.
C     ***************
C
      IF (IFL .NE. 2)  GO TO 20
        NCOLS = NEQNS + NEQC + NIQC
        IFL = 0
        NACT = 0
        IADDC = 0
        IDELC = 0
        ALPHA = ZERO
        PENPAR = PENPAR/OCT
C
C     ***************
C     INITIALIZE  INDX,RES,PTE,INDX
C     ***************
C
        DO 10 I=1,NCOLS
          RES(I) = SDOT(NVARS,E(1,I),1,X,1) - F(I)
          PTE(I) = ZERO
          INDX(I) = I
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE UPDATE
     *                     (IADDC,IDELC,NACT,NEQNS,NEQC,NIQC,NVARS,
     *                     ICYC,IFL,MXS,E,NER,X,F,RES,GRD,
     *                     EL1N,PEN,PENPAR,INDX,ZZ,
     *                     NZZR,DD,RR,W)
C
      INTEGER IADDC,IDELC,ICYC,IFL,INDX(1),MXS
      INTEGER NACT,NEQC,NEQNS,NER,NIQC,NVARS,NZZR
      REAL DD(1),E(NER,1),EL1N,F(1),GRD(1),PEN,PENPAR
      REAL RES(1),RR(1),W(1),X(1),ZZ(NZZR,1)
C
C     ***************
C     CL1  VERSION.
C
C     PREPARATION FOR NEXT MINIMIZATION STEP.
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C     +++++++++++++++
C
      INTEGER NALLQ,NCOLS
      REAL ONE,ZERO
C
      DATA ONE /1.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
C     ***************
C     DETERMINE THE ACTIVE EQUATIONS AND ACTIVE
C     CONSTRAINTS.  COMPUTE RESIDUALS AND FUNCTION VALUE.
C     UPDATE THE  Z*D*R  DECOMPOSITION.
C     ***************
C
      NALLQ = NEQNS + NEQC
      NCOLS = NALLQ + NIQC
      IF (IFL .NE. 0)  GO TO 20
        ICYC = ICYC + 1
        IF (ICYC .LE. MXS)  GO TO 10
          IFL = 3
          GO TO 20
   10   CONTINUE
          CALL DELCOL(IADDC,IDELC,NACT,NVARS,ZZ,NZZR,DD,RR,INDX)
          CALL RESID(IADDC,NACT,NCOLS,NVARS,E,NER,X,F,RES,INDX)
          CALL ADDCOL(IADDC,IDELC,NACT,NVARS,ZZ,NZZR,DD,RR,
     *                E,NER,INDX,W)
          CALL OBJECT(IADDC,NACT,NEQNS,NALLQ,NCOLS,NVARS,E,NER,
     *                RES,GRD,EL1N,PEN,PENPAR,INDX)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE MONIT
     *                    (NACT,NEQC,NIQC,NVARS,ICYC,PSW,X,
     *                     ALPHA,EL1N,PEN,PENPAR,INDX)
C
      INTEGER ICYC,INDX(1),NACT,NEQC,NIQC,NVARS
      LOGICAL PSW
      REAL ALPHA,EL1N,PEN,PENPAR,X(1)
C
C     ***************
C     CL1  VERSION.
C
C     PRINT OUT INFORMATION AT
C     EACH STEP IF  PSW = .TRUE.
C
C     TO CHANGE THE OUTPUT MEDIUM,
C     CHANGE THE  DATA  DECLARATION
C     APPEARING BELOW.
C     ***************
C
      INTEGER I,NXOUTP
C
      DATA NXOUTP /6/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (.NOT. PSW)  GO TO 20
        WRITE (NXOUTP,6000) ICYC,ALPHA
        WRITE (NXOUTP,6010) (X(I),I=1,NVARS)
        WRITE (NXOUTP,6020)  NACT
        IF (NACT .LE. 0)  GO TO 10
          WRITE (NXOUTP,6030)
          WRITE (NXOUTP,6040) (INDX(I),I=1,NACT)
   10   CONTINUE
        WRITE (NXOUTP,6050) EL1N
        IF (NEQC + NIQC .GT. 0)  WRITE (NXOUTP,6060) PENPAR,PEN
   20 CONTINUE
      RETURN
C
 6000 FORMAT(20H0***** CYCLE NUMBER ,I5/
     *       6X,13HSTEP TAKEN = ,1PE15.7/
     *       6X,11HX-VECTOR...)
 6010 FORMAT(5X,3(1PE15.7))
 6020 FORMAT(6X,38HNUMBER OF ACTIVE EQUATIONS/CONSTRAINTS,I5)
 6030 FORMAT(6X,31HACTIVE EQUATIONS/CONSTRAINTS...)
 6040 FORMAT(5X,7I5)
 6050 FORMAT(6X,19HRESIDUAL L1 NORM = ,1PE15.7)
 6060 FORMAT(6X,20HPENALTY PARAMETER = ,1PE10.2/
     *       6X,25HPENALTY FUNCTION VALUE = ,1PE15.7)
      END
      SUBROUTINE FINDP
     *                    (IDELC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL,
     *                     E,NER,X,F,RES,GRD,P,EL1N,AMAG,CGMAG,
     *                     PENPAR,INDX,ZZ,NZZR,DD,RR,W)
C
      INTEGER IDELC,IFL,INDX(1),NACT,NEQC,NEQNS,NER,NIQC,NVARS,NZZR
      REAL AMAG,CGMAG,DD(1),E(NER,1),EL1N,F(1),GRD(1),P(1),PENPAR
      REAL RES(1),RR(1),W(1),X(1),ZZ(NZZR,1)
C
C     ***************
C     CL1  VERSION.
C
C     DETERMINE DESCENT DIRECTION  P
C     (OR DISCOVER OPTIMALITY)
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C
C     BLAS  SASUM,SCOPY,SSCAL
C
C     EPS  IS THE SMALLEST POSITIVE NUMBER WHICH
C     SATISFIES   (1.0 + EPS) .GT. 1.0   IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING,
C      EPS  CAN BE SET TO A USER-SPECIFIED TOLERANCE.)
C     +++++++++++++++
C
      INTEGER COEFX,I,IX,NALLQ,NALQP1,NCOLS,NQNP1,TOPX
      LOGICAL FAIL
      REAL EPS,GRDNRM,ONE,PNRM,PROD,TEST,ZERO
C
      REAL SASUM
C
C     +++++++++++++++
C     THE FOLLOWING DECLARATIONS ARE NECESSARY
C     FOR PORTABILITY WHEN  SCOPY  IS USED, AS
C     IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE
C     (ONE=UNITY  AND  ZERO=ZIP  IN THIS CASE).
C     +++++++++++++++
C
      REAL UNITY(1),ZIP(1)
      EQUIVALENCE (ONE,UNITY(1)),(ZERO,ZIP(1))
C
      DATA EPS /0.0E-31/
      DATA ONE /1.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IDELC = 0
      IF (IFL .NE. 0)  GO TO 90
        NALLQ = NEQNS + NEQC
        NALQP1 = NALLQ + 1
        NCOLS = NALLQ + NIQC
        NQNP1 = NEQNS + 1
        COEFX = 1
        TOPX = COEFX + NVARS
C
C     ***************
C     PROJECT THE NEGATIVE OF THE RESTRICTED GRADIENT
C     ONTO THE ORTHOGONAL COMPLEMENT OF THE SPACE
C     SPANNED BY THE ACTIVE COLUMNS.
C     ***************
C
        CALL ZDRPOC(NVARS,NACT,ZZ,NZZR,DD,GRD,P,FAIL)
        IF (.NOT.FAIL)  GO TO 5
          IFL = 4
          GO TO 90
    5   CONTINUE
        CALL SSCAL(NVARS,-ONE,P,1)
        PNRM = SASUM(NVARS,P,1)
        GRDNRM = SASUM(NVARS,GRD,1)
C
C     ***************
C     IF THE PROJECTION IS NOT ZERO,
C     IT WILL SERVE AS A DESCENT DIRECTION.
C
C     OTHERWISE FIND THE REPRESENTATION OF
C     THE RESTRICTED GRADIENT AS A LINEAR
C     COMBINATION OF THE ACTIVE COLUMNS.
C     THE COEFFICIENTS OF THE LINEAR COMBINATION
C     ARE TO BE STORED IN THE ARRAY  COEF
C     (THAT IS, IN  W(COEFX),...,W(COEFX+NACT-1)).
C     ***************
C
        IF (PNRM .GT. EPS*(AMAG*PENPAR + CGMAG))  GO TO 90
        IF (NACT .EQ. 0)  GO TO 50
          CALL ZDRGNV(NVARS,NACT,ZZ,NZZR,RR,GRD,W(COEFX),FAIL)
          IF (.NOT. FAIL)  GO TO 10
            IFL = 4
            GO TO 90
   10     CONTINUE
C
C     ***************
C     CONVERT THE COEFFICIENTS OF THE LINEAR
C     COMBINATION INTO A DESCENT DIRECTION  P ,
C     OR DETERMINE OPTIMALITY.
C
C     IF THE OPTIMALITY TEST IS NOT SATISFIED,
C     GETV  WILL INDICATE AN EQUATION/CONSTRAINT
C     TO BE DELETED FROM ACTIVITY BY THE VALUE
C     OF  IDELC.  FOR OPTIMALITY,  IDELC=0.
C     ***************
C
          CALL GETV(IDELC,NACT,NVARS,NEQNS,NALLQ,E,NER,GRD,
     *              W(COEFX),PENPAR,INDX)
          PNRM = ZERO
          IF (IDELC .EQ. 0)  GO TO 20
            CALL ZDRGIT(NVARS,NACT,ZZ,NZZR,RR,W(COEFX),P,FAIL,W(TOPX))
            IF (.NOT. FAIL)  PNRM = SASUM(NVARS,P,1)
            IF (.NOT. FAIL)  GO TO 20
              IFL = 4
              GO TO 90
   20   CONTINUE
C
C     ***************
C     IF A DESCENT DIRECTION  P  COULD HAVE BEEN FOUND,
C     IT HAS BEEN OBTAINED BY THIS POINT IN THE PROGRAM.
C
C     CHECK FOR OPTIMALITY.
C
C     PNRM  HAS BEEN SET EXACTLY ZERO
C     AFTER THE CALL TO SUBROUTINE  GETV
C     IF THE OPTIMALITY CONDITIONS ARE SATISFIED.
C     THE CHECK BELOW HAS BEEN MADE SOMEWHAT
C     COMPLICATED TO ALLOW FOR THE RARE EVENT THAT
C     THE RESTRICTED GRADIENT IS ZERO AND NO
C     COLUMNS ARE ACTIVE,  OR THAT THE  L1  NORM OF
C               (A-TRANSPOSE) * X - F
C     IS COMPUTATIONALLY ZERO.
C     (THE CALL TO THE SUBROUTINE  REFINE
C      MAY BE OMITTED, IF DESIRED.)
C     ***************
C
        IF (PNRM .LE. EPS*(AMAG*PENPAR + CGMAG))  GO TO 50
        DO 40 I=1,NEQNS
          TEST = ABS(F(I))
          DO 30 IX=1,NVARS
            PROD = ABS(E(IX,I)*X(IX))
            IF (PROD .GT. TEST)  TEST = PROD
   30     CONTINUE
          IF (ABS(RES(I)) .GT. EPS*TEST)  GO TO 90
   40   CONTINUE
   50   CONTINUE
          IFL = 1
          CALL REFINE(NACT,NEQNS,NCOLS,NVARS,IFL,E,NER,
     *                X,F,EL1N,RES,INDX,ZZ,NZZR,RR,W)
          IF (IFL .NE. 1)  GO TO 90
C
C     ***************
C     IF THE PROBLEM HAS CONSTRAINTS,
C     CHECK FEASIBILITY.
C     ***************
C
          IF (NQNP1 .GT. NCOLS)  GO TO 90
            DO 80 I=NQNP1,NCOLS
              TEST = ABS(F(I))
              DO 60 IX=1,NVARS
                PROD = ABS(E(IX,I)*X(IX))
                IF (PROD .GT. TEST)  TEST = PROD
   60         CONTINUE
              TEST = EPS*TEST
              IF (I .GT. NALLQ)  GO TO 70
                IF (ABS(RES(I)) .LE. TEST)  GO TO 80
                  IFL = 2
                  GO TO 90
   70         CONTINUE
                IF (RES(I) .GE. (-TEST))  GO TO 80
                  IFL = 2
                  GO TO 90
   80       CONTINUE
   90 CONTINUE
      RETURN
      END
      SUBROUTINE STEP
     *                   (IADDC,NACT,NEQNS,NEQC,NIQC,NVARS,IFL,
     *                    E,NER,X,RES,GRD,P,PTE,ALPHA,PENPAR,INDX,ALF)
C
      INTEGER IADDC,IFL,INDX(1),NACT,NEQC,NEQNS,NER,NIQC,NVARS
      REAL ALPHA,ALF(1),E(NER,1),GRD(1),P(1)
      REAL PENPAR,PTE(1),RES(1),X(1)
C
C     ***************
C     CL1  VERSION.
C
C     PIECEWISE LINEAR LINE SEARCH.
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES ABS,SIGN
C
C     BLAS  SASUM,SAXPY,SDOT
C
C     EPS  IS THE SMALLEST POSITIVE NUMBER WHICH
C     SATISFIES   (1.0 + EPS) .GT. 1.0   IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING,
C      EPS  CAN BE SET TO A USER-SPECIFIED TOLERANCE.)
C
C     BIG  IS THE LARGEST POSITIVE NUMBER
C     WHICH CAN BE REPRESENTED IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     +++++++++++++++
C
      INTEGER I,IIN,IX,JX,NACTP1,NALLQ,NCOLS,NUM,NUMNAC
      REAL BIG,DEN,EPS,GRDNRM,ONE,PNRM,PTG
      REAL RATIO,RESID,TMP,TWO,ZERO
C
      REAL SASUM,SDOT
C
      DATA BIG /1.0E+31/
      DATA EPS /0.0E-31/
      DATA ONE /1.0E+00/
      DATA TWO /2.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
C     ***************
C     THIS ROUTINE DETERMINES ALL OF THE RATIOS  ALF
C     OF THE FORM
C        -RES(I)/((E(.,I)-TRANSP)*P),
C              FOR  I = K+1,...,MPL
C     WHICH ARE NONNEGATIVE AND HENCE INDICATE DISTANCES
C     FROM THE POINT  X  TO BREAKPOINTS WHICH WILL
C     BE ENCOUNTERED IN TRAVEL ALONG DIRECTION  P.
C     THE INDEX VECTOR  INDX  IS REARRANGED SO THAT
C     ITS  K+1  THROUGH  NUM  COMPONENTS CORRESPOND TO
C     THESE NONNEGATIVE RATIOS.
C     THE RESULTS ARE HEAPED SO THAT THE  ALF  VALUES CAN
C     BE INSPECTED IN ORDER FROM SMALLEST TO LARGEST.
C     THE BREAKPOINT  ALPHA  GIVING THE MINIMUM OBJECTIVE
C     FUNCTION VALUE IS FOUND, AND  X  IS
C     ADJUSTED TO  X + ALPHA*P .
C
C     THE INNER PRODUCTS  (E(.,I)-TRANSPOSE)*P  ARE SAVED
C     FOR LATER USE IN UPDATING THE RESIDUAL VALUES.
C     ***************
C
      ALPHA = ZERO
      IF (IFL .NE. 0)  GO TO 90
        NALLQ = NEQNS + NEQC
        NCOLS = NALLQ + NIQC
        NACTP1 = NACT + 1
        NUM = 0
        IF (1 .GT. NACT)  GO TO 20
          DO 10 I=1,NACT
            IX = INDX(I)
            PTE(IX) = SDOT(NVARS,E(1,IX),1,P,1)
   10     CONTINUE
   20   CONTINUE
        IF (NACTP1 .LE. NCOLS)  GO TO 30
          IFL = 1
          GO TO 90
   30   CONTINUE
          DO 50 I=NACTP1,NCOLS
            IX = INDX(I)
            RESID = RES(IX)
            DEN = SDOT(NVARS,E(1,IX),1,P,1)
            PTE(IX) = DEN
            IF (SIGN(ONE,RESID) .EQ. SIGN(ONE,DEN)
     *                    .AND. RESID .NE. ZERO)  GO TO 50
              RESID = ABS(RESID)
              DEN = ABS(DEN)
              IF (DEN .GE. ONE)  GO TO 40
                IF (RESID .GE. DEN*BIG)  GO TO 50
   40         CONTINUE
                RATIO = RESID/DEN
                NUM = NUM + 1
                NUMNAC = NUM + NACT
                JX = INDX(NUMNAC)
                INDX(NUMNAC) = IX
                INDX(I) = JX
                ALF(NUM) = RATIO
   50     CONTINUE
          IF (NUM .GT. 0)  GO TO 60
            IFL = 2
            GO TO 90
   60     CONTINUE
C
C     ***************
C     HEAP THE POSITIVE RATIOS
C     ***************
C
            CALL DKHEAP(.TRUE.,NUM,INDX(NACTP1),ALF)
C
C     ***************
C     TRAVEL ALONG  P  UNTIL NO FURTHER DECREASE IN THE
C     PENALTY FUNCTION IS POSSIBLE
C     ***************
C
            IIN = NUM
            PTG = SDOT(NVARS,GRD,1,P,1)
            PNRM = SASUM(NVARS,P,1)
            GRDNRM = SASUM(NVARS,GRD,1)
            DO 70 I=1,NUM
              IX = INDX(NACTP1)
              TMP = -SIGN(ONE,RES(IX))
              IF (IX .LE. NALLQ)  TMP = TMP*TWO
              IF (IX .LE. NEQNS)  TMP = TMP*PENPAR
              PTG = PTG + TMP*PTE(IX)
              IF (PTG .GE. (-EPS*GRDNRM*PNRM))  GO TO 80
                CALL DKHEAP(.FALSE.,IIN,INDX(NACTP1),ALF)
   70       CONTINUE
            IFL = 2
            GO TO 90
   80       CONTINUE
            IADDC = NACTP1
C
C     ***************
C     ADJUST  X  TO  X + ALPHA*P
C     ***************
C
            ALPHA = ALF(1)
            CALL SAXPY(NVARS,ALPHA,P,1,X,1)
   90 CONTINUE
      RETURN
      END
      SUBROUTINE REFINE(NACT,NEQNS,NCOLS,NVARS,IFL,E,NER,
     *                  X,F,EL1N,RES,INDX,ZZ,NZZR,RR,W)
C
      INTEGER IFL,INDX(1),NACT,NCOLS,NEQNS,NER,NVARS,NZZR
      REAL E(NER,1),EL1N,F(1),RES(1),RR(1),W(1),X(1),ZZ(NZZR,1)
C
C     ***************
C     A ROUTINE FOR REFINING THE SOLUTION
C     PRODUCED BY  CL1.
C
C     (THIS ROUTINE MAY BE OMITTED IF DESIRED.)
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C
C     BLAS  SDOT
C     +++++++++++++++
C
      INTEGER I,IX
      LOGICAL FAIL
      REAL TMP,ZERO
C
      REAL SDOT
C
      DATA ZERO /0.0E+00/
C
C     /////////////// BEGIN PROGRAM ///////////////
C
      IF (NACT .EQ. 0)  GO TO 40
        DO 10 I=1,NACT
          IX = INDX(I)
          RES(I) = F(IX)
   10   CONTINUE
        CALL ZDRGIT(NVARS,NACT,ZZ,NZZR,RR,RES,X,FAIL,W)
        IF (.NOT. FAIL)  GO TO 20
          IFL = 4
          GO TO 40
   20   CONTINUE
        EL1N = ZERO
        DO 30 I=1,NCOLS
          TMP = SDOT(NVARS,E(1,I),1,X,1) - F(I)
          RES(I) = TMP
          IF (I .LE. NEQNS)  EL1N = EL1N + ABS(TMP)
   30   CONTINUE
   40 CONTINUE
      RETURN
      END
C
C     ---------------
C     THIRD LEVEL SUBROUTINES --
C          DELCOL,RESID,ADDCOL,OBJECT,GETV
C     ---------------
C
      SUBROUTINE DELCOL(IADDC,IDELC,NACT,NROW,ZZ,NZZR,DD,RR,INDX)
C
      INTEGER INDX(1),NZZR,NACT,IDELC,IADDC,NROW
      REAL DD(1),RR(1),ZZ(NZZR,1)
C
C     ***************
C     CL1  VERSION OF  DELCOL.
C
C     THIS ROUTINE ADMINISTERS THE DELETION OF THE COLUMN
C     INDICATED BY THE VALUE OF IDELC
C     FROM AN  NROW BY NACT   Z*D*R   DECOMPOSITION.
C     NOTE THAT THE VALUE OF IDELC
C     IS THE NUMBER OF A COLUMN IN THE DECOMPOSITION
C     RATHER THAN A NUMBER WHICH REFERS TO
C     A COLUMN IN THE MATRIX  E.
C     (THE  E-COLUMN  NUMBERS CORRESPONDING TO
C     THE COLUMNS OF THE FACTORIZATION ARE TO BE
C     FOUND IN   INDX(1),...,INDX(NACT) .
C     THE CONTENTS OF   INDX(NACT+1),...,INDX(IADDC)
C     INDICATE COLUMNS OF  E  WHICH ARE SLATED FOR
C     ADDITION TO THE DECOMPOSITION.)
C     THE VECTOR  INDX   IS REARRANGED BY
C     PERMUTING THE ELEMENT WHICH CORRESPONDS TO
C     THE DELETION OUT TO THE   IADDC-TH  POSITION.
C     NACT  AND  IADDC  ARE DECREASED ACCORDINGLY.
C     ***************
C
      INTEGER I,IDLP1,IXDLC
      LOGICAL FAIL
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (IDELC .EQ. 0)  GO TO 20
        IDLP1 = IDELC + 1
        IXDLC = INDX(IDELC)
        DO 10 I=IDLP1,IADDC
          INDX(I - 1) = INDX(I)
   10   CONTINUE
        INDX(IADDC) = IXDLC
        IADDC = IADDC - 1
        CALL ZDRCOU(NROW,NACT,ZZ,NZZR,DD,RR,IDELC,FAIL)
        IDELC = IXDLC
   20 CONTINUE
      RETURN
      END
      SUBROUTINE RESID(IADDC,NACT,NCOLS,NVARS,E,NER,X,F,RES,INDX)
C
      INTEGER IADDC,INDX(1),NACT,NCOLS,NER,NVARS
      REAL E(NER,1),F(1),RES(1),X(1)
C
C
C     ***************
C     COMPUTE THE RESIDUALS
C          (E(.,IX)-TRANSP)*X - F(IX)  .
C     THE RESIDUALS ARE STORED IN THE ARRAY  RES.
C     INDX  IS REARRANGED SO THAT THE ZERO RESIDUALS
C     CORRESPOND TO  INDX(1),...,INDX(IADDC)  .
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS,IFIX,FLOAT,SQRT
C
C     BLAS  SDOT
C
C     EPS  IS THE SMALLEST POSITIVE NUMBER WHICH
C     SATISFIES   (1.0 + EPS) .GT. 1.0   IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING,
C      EPS  CAN BE SET TO A USER-SPECIFIED TOLERANCE.)
C     +++++++++++++++
C
      INTEGER I,IADP1,IDUMMY,IRAND,IX,J,NACTP1
      REAL EPS,PROD,TEMP,TEST,TOL,ZERO
C
      REAL SDOT,UNIF01
C
      DATA EPS /0.0E-32/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      TOL = EPS*SQRT(FLOAT(NVARS))
      NACTP1 = NACT + 1
      IF (1 .GT. IADDC)  GO TO 20
C
C     ***************
C     ZERO OUT ALL RESIDUALS KNOWN TO BE ZERO.
C     ***************
C
        DO 10 I=1,IADDC
          IX = INDX(I)
          RES(IX) = ZERO
   10   CONTINUE
   20 CONTINUE
C
C     ***************
C     COMPUTE THE REMAINING RESIDUALS.
C     DETECT ANY MORE RESIDUALS WHICH
C     ARE COMPUTATIONALLY ZERO, AND
C     SET THEM EXACTLY ZERO.  THEIR
C     ASSOCIATED INDICES ARE PERMUTED
C     SO THAT THEY ARE STORED IN
C     INDX(NACT+1),...,NACT(IADDC).
C
C     (A FAIRLY TIGHT ZERO CHECK IS USED.
C     IT IS FAR LESS EXPENSIVE IN RUNNING
C     TIME TO NEGLECT AN EXTRA ZERO
C     RESIDUAL THAN TO ACCEPT IT AND RISK
C     INVOKING THE ANTI-CYCLING
C     MECHANISMS IN THE PROGRAM.
C     THE ACCURACY OF THE SOLUTION AS
C     FINALLY DETERMINED IS NOT AFFECTED.)
C     ***************
C
      IADP1 = IADDC + 1
      IF (IADP1 .GT. NCOLS)  GO TO 50
        DO 40 I=IADP1,NCOLS
          IX = INDX(I)
          TEMP = SDOT(NVARS,E(1,IX),1,X,1) - F(IX)
          TEST = ABS(F(IX))
          DO 25 J=1,NVARS
            PROD = ABS(E(J,IX)*X(J))
            IF (PROD .GT. TEST)  TEST = PROD
   25     CONTINUE
          TEST = TOL*TEST
          IF (ABS(TEMP) .LE. TEST)  GO TO 30
            RES(IX) = TEMP
            GO TO 40
   30     CONTINUE
            IADDC = IADDC + 1
            INDX(I) = INDX(IADDC)
            INDX(IADDC) = IX
            RES(IX) = ZERO
   40   CONTINUE
   50 CONTINUE
C
C     ***************
C     IF ANY NEW ZERO RESIDUALS HAVE
C     BEEN FOUND, RANDOMIZE THEIR
C     ORDERING AS AN ANTI-CYCLING
C     DEVICE FOR  ADDCOL.
C     ***************
C
      IF (IADDC .LE. NACTP1)  GO TO 70
        DO 60 I=NACTP1,IADDC
          IRAND = I + IFIX(FLOAT(IADDC - I + 1)*UNIF01(0,IDUMMY))
          IX = INDX(IRAND)
          INDX(IRAND) = INDX(I)
          INDX(I) = IX
   60   CONTINUE
   70 CONTINUE
      RETURN
      END
      SUBROUTINE ADDCOL(IADDC,IDELC,NACT,NVARS,ZZ,NZZR,DD,RR,
     *                  E,NER,INDX,W)
C
      INTEGER IADDC,IDELC,INDX(1),NACT,NER,NVARS,NZZR
      REAL DD(1),E(NER,1),RR(1)
      REAL W(1),ZZ(NZZR,1)
C
C     ***************
C     CL1 VERSION OF ADDCOL.
C
C     THIS ROUTINE ADMINISTERS THE ADJUSTMENT OF THE
C     Z*D*R   DECOMPOSITION FOR ANY NEW ZERO RESIDUALS.
C     THE DATA CORRESPONDING TO THE ZERO RESIDUALS IS INDEXED
C     IN  INDX(NACT+1),...,INDX(IADDC).
C     ***************
C
C     +++++++++++++++
C     BLAS  SASUM
C
C     EPS  IS THE SMALLEST POSITIVE NUMBER WHICH
C     SATISFIES   (1.0 + EPS) .GT. 1.0   IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING,
C      EPS  CAN BE SET TO A USER-SPECIFIED TOLERANCE.)
C     +++++++++++++++
C
      INTEGER I,ISTRT,IX,NACTP1,TOPX
      LOGICAL FAIL
      REAL COLNRM,EPS,PRJNRM
C
      REAL SASUM
C
      DATA EPS /0.0E-31/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      TOPX = NVARS + 1
      ISTRT = NACT + 1
      IF (ISTRT .GT. IADDC)  GO TO 20
C
C     ***************
C     CANDIDATES FOR ADDITION TO THE  Z*D*R
C     FACTORIZATION ARE INSPECTED IN RANDOM
C     ORDER TO HINDER CYCLING.
C     THE RANDOMIZATION WAS CARRIED OUT BY  RESID.
C
C     IF A CANDIDATE HAS JUST BEEN RELEASED
C     FROM THE FACTORIZATION OR IS DEPENDENT UPON THE
C     COLUMNS IN THE FACTORIZATION,
C     THEN IT IS OMITTED FROM ADDITION.
C
C     UPON EXIT, INDICES OF SUCH OMITTED
C     COLUMNS ARE TO BE FOUND IN
C          INDX(NACT+1),...,INDX(IADDC) .
C     ***************
C
        DO 10 I=ISTRT,IADDC
          NACTP1 = NACT + 1
          IX = INDX(I)
          CALL ZDRPOC(NVARS,NACT,ZZ,NZZR,DD,E(1,IX),W,FAIL)
          COLNRM = SASUM(NVARS,E(1,IX),1)
          PRJNRM = SASUM(NVARS,W,1)
          IF (PRJNRM .LE. EPS*COLNRM .OR. IX .EQ. IDELC)  GO TO 10
            INDX(I) = INDX(NACTP1)
            INDX(NACTP1) = IX
            CALL ZDRCIN(NVARS,NACT,ZZ,NZZR,DD,RR,E(1,IX),FAIL,W)
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE OBJECT(IADDC,NACT,NEQNS,NALLQ,NCOLS,NVARS,
     *                  E,NER,RES,GRD,EL1N,PEN,PENPAR,INDX)
C
      INTEGER IADDC,INDX(1),NACT,NALLQ,NEQNS,NCOLS,NER,NVARS
      REAL E(NER,1),EL1N,GRD(1),PEN,PENPAR,RES(1)
C
C     ***************
C     CL1 VERSION OF OBJECT.
C
C     THIS ROUTINE ADMINISTERS THE EVALUATION OF THE
C     PENALTY (OBJECTIVE) FUNCTION GIVEN THE EQUATION
C     AND CONSTRAINT RESIDUALS.  IT ALSO COMPUTES THE
C     RESTRICTED GRADIENT OF THE FUNCTION.
C
C     COLUMNS WHICH ARE NOT IN THE  Z*D*R FACTORIZATION
C     BUT WHICH ARE ASSOCIATED WITH ZERO RESIDUALS MUST
C     BE INCLUDED IN THE RESTRICTED GRADIENT WITH RANDOM
C     SIGNS AS AN ANTI-CYCLING DEVICE.
C     THE INDICES OF THESE COLUMNS ARE TO BE
C     FOUND IN  INDX(NACT+1),...,INDX(IADDC)
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS,SIGN
C
C     BLAS  SAXPY,SCOPY
C     +++++++++++++++
C
      INTEGER I,IDUMMY,IX,NACTP1
      REAL HALF,ONE,SGN,TMP,ZERO
C
      REAL UNIF01
C
C     +++++++++++++++
C     THE FOLLOWING DECLARATIONS ARE NECESSARY
C     FOR PORTABILITY WHEN  SCOPY  IS USED, AS
C     IT IS BELOW, TO FILL ARRAYS WITH A SINGLE
C     VALUE  (ZERO=ZIP  IN THIS CASE).
C     +++++++++++++++
C
      REAL ZIP(1)
      EQUIVALENCE (ZERO,ZIP(1))
C
      DATA HALF /0.5E+00/
      DATA ONE /1.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      NACTP1 = NACT + 1
      EL1N = ZERO
      PEN = ZERO
      CALL SCOPY(NVARS,ZIP,0,GRD,1)
      IF (NACTP1 .GT. NCOLS)  GO TO 30
        DO 20 I=NACTP1,NCOLS
          IX = INDX(I)
          TMP = RES(IX)
          SGN = SIGN(ONE,TMP)
          IF (I .LE. IADDC)  SGN = SIGN(ONE,HALF - UNIF01(0,IDUMMY))
          IF (IX .GT. NALLQ .AND. SGN .GT. ZERO)  GO TO 20
            TMP = ABS(TMP)
            IF (IX .GT. NEQNS)  GO TO 10
              EL1N = EL1N + TMP
              TMP = TMP*PENPAR
              SGN = SGN*PENPAR
   10       CONTINUE
            PEN = PEN + TMP
            CALL SAXPY(NVARS,SGN,E(1,IX),1,GRD,1)
   20   CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE GETV(IDELC,NACT,NVARS,NEQNS,NALLQ,E,NER,
     *                GRD,COEF,PENPAR,INDX)
C
      INTEGER IDELC,INDX(1),NACT,NALLQ,NEQNS,NER,NVARS
      REAL COEF(1),E(NER,1),GRD(1),PENPAR
C
C     ***************
C     CL1  VERSION.
C
C     SET UP THE RIGHT-HAND-SIDE VECTOR
C     (AND STORE IN THE ARRAY  COEF)
C     FOR THE LINEAR PROBLEM WHICH DETERMINES
C     A DESCENT DIRECTION  P  IN THE CASE WHERE
C     THE PROJECTION OF THE RESTRICTED GRADIENT IS ZERO.
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS,FLOAT,IFIX,SIGN
C
C     BLAS  SAXPY
C
C     EPS  IS THE SMALLEST POSITIVE NUMBER WHICH
C     SATISFIES   (1.0 + EPS) .GT. 1.0   IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     (ALTERNATIVELY, FOR LESS STRICT ZERO CHECKING,
C      EPS  CAN BE SET TO A USER-SPECIFIED TOLERANCE.)
C     +++++++++++++++
C
      INTEGER I,IDUMMY,IRAND,IX
      REAL CF,EPS,ONE,OPE,S,TMP,TMPSAV,ZERO
C
      REAL UNIF01
C
      DATA EPS /0.0E-31/
      DATA ONE /1.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
C     ***************
C     FIND THE MOST OUT-OF-KILTER
C     COEFFICIENT.  BEGIN INSPECTING
C     THE COEFFICIENTS AT A RANDOM INDEX
C     TO HINDER CYCLING.  SET  COEF
C     TO ZERO ON THE FLY.
C     ***************
C
      OPE = ONE + EPS
      IDELC = 0
      TMPSAV = ZERO
      IF (1 .GT. NACT)  GO TO 40
        IRAND = IFIX(FLOAT(NACT)*UNIF01(0,IDUMMY))
        DO 30 I=1,NACT
          IRAND = IRAND + 1
          IF (IRAND .GT. NACT)  IRAND = 1
          IX = INDX(IRAND)
          CF = COEF(IRAND)
          COEF(IRAND) = ZERO
          IF (IX .GT. NALLQ)  GO TO 10
            IF (IX .LE. NEQNS)  CF = CF/PENPAR
            TMP = OPE - ABS(CF)
            GO TO 20
   10     CONTINUE
            TMP = CF + EPS
   20    CONTINUE
         IF (TMP .GE. TMPSAV)  GO TO 30
           IDELC = IRAND
           S = SIGN(ONE,CF)
           TMPSAV = TMP
   30   CONTINUE
C
C     ***************
C     IF NO COEFFICIENTS ARE OUT OF KILTER,
C     THEN RETURN.  OTHERWISE SET A
C     VALUE IN AN APPROPRIATE COMPONENT
C     (INDICATED BY  IDELC)  OF  COEF
C     AND ADJUST THE RESTRICTED GRADIENT
C     IF NECESSARY.
C     ***************
C
        IF (IDELC .EQ. 0)  GO TO 40
          COEF(IDELC) = ONE
          IX = INDX(IDELC)
          IF (IX .GT. NALLQ)  GO TO 40
            COEF(IDELC) = -S
            TMP = -S
            IF (IX .LE. NEQNS)  TMP = TMP*PENPAR
            CALL SAXPY(NVARS,TMP,E(1,IX),1,GRD,1)
   40 CONTINUE
      RETURN
      END
C
C     ---------------
C     FOURTH LEVEL SUBROUTINES --
C               DKHEAP,UNIF01,ZDRCIN,ZDRCOU,
C               ZDRGIT,ZDRGNV,ZDRPOC
C     ---------------
C
      SUBROUTINE DKHEAP(MAKE,IR,INDX,ARAY)
C
      INTEGER INDX(1),IR
      LOGICAL MAKE
      REAL ARAY(1)
C
C     ***************
C     AN ADAPTATION OF D. E. KNUTH,S HEAPING
C     ROUTINES (SEE VOLUME 3 OF
C          THE ART OF COMPUTER PROGRAMMING  ).
C     IF  MAKE  IS  .TRUE.,  THE FULL HEAP BUILDING
C     PROCESS IS CARRIED OUT ON
C          ARAY(1),...,ARAY(IR) ,
C     AND THE VALUE OF  IR  IS UNCHANGED.
C     IF  MAKE  IS  .FALSE.,  ONE STEP OF THE SORTING
C     PROCESS IS CARRIED OUT TO PROVIDE THE NEXT
C     ELEMENT OF  ARAY  IN ORDER,  AND THE VARIABLE
C     IR  IS DECREASED BY ONE.  THE INTERRUPTION OF THE
C     SORTING PHASE IS BUILT IN VIA THE FLAG  ONCE.
C     INDX  IS AN INDEX VECTOR ASSOCIATED WITH
C     ARAY  WHICH MUST BE REARRANGED IN PARALLEL
C     WITH IT.
C     ***************
C
      INTEGER I,IL,IT,J
      LOGICAL ONCE
      REAL T
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (IR .GT. 1)  GO TO 5
        IF (.NOT.MAKE)  IR = 0
        RETURN
    5 CONTINUE
C
C     ***************
C     TEST WHETHER OR NOT THE INITIAL
C     HEAP IS TO BE BUILT
C     ***************
C
      IL = 1
      IF (MAKE)  IL = (IR/2) + 1
      ONCE = .FALSE.
C
C     ***************
C     THE LOOP BEGINS HERE
C     ***************
C
   10 CONTINUE
      IF (IL .GT. 1)  GO TO 20
C
C     ***************
C     THE SORTING PHASE USES THIS BRANCH
C     ***************
C
        IF (MAKE .OR. ONCE)  RETURN
        ONCE = .TRUE.
        IT = INDX(IR)
        T = ARAY(IR)
        INDX(IR) = INDX(1)
        ARAY(IR) = ARAY(1)
        IR = IR - 1
        IF (IR .GT. 1)  GO TO 30
          INDX(1) = IT
          ARAY(1) = T
          RETURN
   20 CONTINUE
C
C     ***************
C     THE HEAP-BUILDING PHASE USES THIS BRANCH
C     ***************
C
        IL = IL - 1
        IT = INDX(IL)
        T = ARAY(IL)
   30 CONTINUE
C
C     ***************
C     THE REMAINING STATEMENTS ARE COMMON
C     TO BOTH PHASES AND EMBODY THE
C     HEAP-RECTIFYING (SIFTING) SECTION
C     ***************
C
      J = IL
   40 CONTINUE
      I = J
      J = 2*J
      IF (J - IR)  50, 60, 70
   50 CONTINUE
      IF (ARAY(J) .LE. ARAY(J + 1))  GO TO 60
        J = J + 1
   60 CONTINUE
      IF (T .LE. ARAY(J))  GO TO 70
        INDX(I) = INDX(J)
        ARAY(I) = ARAY(J)
        GO TO 40
   70 CONTINUE
      INDX(I) = IT
      ARAY(I) = T
      GO TO 10
      END
      REAL FUNCTION UNIF01(ISEED,IX)
C
      INTEGER ISEED,IX,IX0
C
      DATA IX0/2/
C
C     +++++++++++++++
C     SYSTEM ROUTINES  FLOAT,MOD
C     +++++++++++++++
C
C     --------------------------------------------------------------
C     --------------------------------------------------------------
C
C     *****PURPOSE-
C     THIS FUNCTION RETURNS A PSEUDO-RANDOM NUMBER DISTRIBUTED
C     UNIFORMLY IN THE INTERVAL (0,1).
C
C     *****PARAMETER DESCRIPTION-
C     ON INPUT-
C
C     ISEED,  IF IT IS NONZERO MODULO 9973, BECOMES THE
C          NEW SEED, I.E. IT REPLACES THE INTERNALLY STORED
C          VALUE OF IX0.  ON MACHINES WHERE FORTRAN VARIABLES
C          RETAIN THEIR VALUES BETWEEN CALLS, THE INTERNALLY
C          STORED VALUE IF IX0 IS THE VALUE ASSIGNED TO  IX  IN
C          THE PREVIOUS INVOCATION OF  UNIF01.  OTHERWISE -- AND
C          IN THE FIRST CALL TO  UNIF01 --  IX0=2.
C
C     ON OUTPUT-
C
C     IX IS THE NEXT INTEGER IN A PSEUDO-RANDOM SEQUENCE OF
C          INTEGERS BETWEEN  1  AND  9972  AND IS GENERATED FROM ITS
C          PREDECESSOR  IX0  (I.E.  FROM  ISEED,  IF  ISEED  IS NONZERO
C          MODULO 9973).  IX  IS THE VALUE WHICH  ISEED  SHOULD HAVE
C          IN THE NEXT INVOCATION OF  UNIF01  TO GET THE NEXT
C          PSEUDO-RANDOM NUMBER.  THE CALLER WILL OFTEN PASS THE
C          SAME VARIABLE FOR  ISEED  AS FOR  IX,
C          E.G.  X = UNIF01(IX,IX).
C
C     *****APPLICATION AND USAGE RESTRICTIONS-
C     UNIF01  SHOULD ONLY BE USED WHEN PORTABILITY IS IMPORTANT AND A
C     COURSE RANDOM NUMBER GENERATOR SUFFICES.  APPLICATIONS REQUIRING
C     A FINE, HIGH PRECISON GENERATOR SHOULD USE ONE WITH A MUCH
C     LARGER MODULUS.
C
C     *****ALGORITHM NOTES-
C     UNIF01 WILL RUN ON ANY MACHINE HAVING AT LEAST 20 BITS OF AC-
C     CURACY FOR FIXED-POINT ARITHMITIC.  IT IS BASED ON A GENERATOR
C     RECOMMENDED IN (3), WHICH PASSES THE SPECTRAL TEST WITH FLYING
C     COLORS -- SEE (1) AND (2).
C
C     REFERENCES-
C     (1) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL
C     RANDOM-NUMBER GENERATORS-  AN EMPIRICAL VIEW,
C     MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV.
C
C     (2) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2
C     (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS.
C
C     (3) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER
C     GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18,
C     PP. 586-593.
C
C     *****GENERAL-
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989.
C
C     PERMISSION FOR THE USE OF  UNIF01  IN  CL1  WAS
C     GENEROUSLY GIVEN BY  V. KLEMA  AND  D. HOAGLIN.
C
C     --------------------------------------------------------------
C     --------------------------------------------------------------
C
      IF (ISEED .EQ. 0) GO TO 10
      IX = MOD(ISEED, 9973)
      IF (IX .NE. 0) IX0 = IX
C  ***
C  IN ORDER THAT ALL FIXED-POINT CALCULATIONS REQUIRE ONLY 20 BIT
C  ARITHMETIC, WE USE TWO CALLS TO  MOD  TO COMPUTE
C  IX0 = MOD(3432*IX0, 9973).
C  ***
 10   IX0 = MOD(52*MOD(66*IX0, 9973), 9973)
      IX = IX0
      UNIF01 = FLOAT(IX0)/9973.0E+00
      RETURN
      END
      SUBROUTINE ZDRCIN(N,K,ZZ,NZZR,DD,RR,COL,FAIL,W)
C
      INTEGER K,N,NZZR
      LOGICAL FAIL
      REAL COL(1),DD(1),RR(1),W(1),ZZ(NZZR,1)
C
C     ***************
C     PREPARED BY RICHARD BARTELS
C     THE UNIVERSITY OF WATERLOO
C     COMPUTER SCIENCE DEPARTMENT
C     LATEST UPDATE .... 30 NOVEMBER, 1979.
C
C     GIVEN THE FACTORIZATION
C
C          ZZ*DD*RR
C
C     OF SOME  N BY K  MATRIX
C
C       (0 .LE. K .LT. N)
C         (N .GE. 1),
C
C     WHERE
C
C          (ZZ-TRANSP)*(ZZ) = (DD-INV),
C          DD  IS DIAGONAL AND NONSINGULAR,
C     AND
C          RR  HAS ZEROS BELOW THE DIAGONAL,
C
C     AND GIVEN A  (K+1)TH  COLUMN
C     TO BE ADDEDTO THE ORIGINAL MATRIX,
C     THIS PROGRAM UPDATES  ZZ,DD AND RR.
C
C     THE VALUE OF  K  IS INCREASED BY ONE.
C
C     W  IS A SCRATCH ARRAY.
C
C     USE IS MADE OF ROUTINES FROM THE LIBRARY
C     OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS).
C
C     PARAMETERS...
C
C                     INPUT/
C       NAME   TYPE   OUTPUT/   SUB-    DESCRIPTION
C                     SCRATCH  SCRIPTS
C       -------------------------------------------
C       N      INT.      I              NUMBER OF ROWS
C
C       K      INT.     I/O             NUMBER OF COLUMNS
C
C       ZZ     REAL     I/O       2     SCALED ORTHOGONAL
C                                       MATRIX
C
C       NZZR   INT.      I              ROW DIMENSION OF ZZ
C
C       DD     REAL     I/O       1     DIAGONAL SCALING
C                                       MATRIX (DIAGONAL
C                                       ELEMENTS ONLY)
C
C       RR     REAL     I/O       1     RIGHT-TRIANGULAR
C                                       MATRIX IN COMPACT FORM.
C
C       COL    REAL      I        1     COLUMN TO BE
C                                       ADDED TO  RR
C
C       FAIL   LOG.      O             .TRUE.  IF  K,N
C                                       ARE IMPROPER
C
C       W      REAL     SCR       1     WORKSPACE
C       -------------------------------------------
C
C     THE  I-TH  SEGMENT OF THE ARRAY  RR  IS  N-I+2 SPACES
C     LONG AND CONTAINS  1  WORK SPACE FOLLOWED BY THE
C     K-I+1  ELEMENTS OF ROW  I  FOLLOWED BY  N-K
C     SCRATCH SPACES.
C     ***************
C
C     +++++++++++++++
C     BLAS  SCOPY,SDOT,SROTM,SROTMG
C     +++++++++++++++
C
      INTEGER I,J,JDEL,KP1,KP2
      REAL DI,ONE,PARAM(5),WI,ZERO
C
      REAL SDOT
C
C     +++++++++++++++
C     THE FOLLOWING DECLARATIONS ARE NECESSARY
C     FOR PORTABILITY WHEN  SCOPY  IS USED, AS
C     IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE
C     (ONE=UNITY  AND  ZERO=ZIP  IN THIS CASE).
C     +++++++++++++++
C
      REAL UNITY(1),ZIP(1)
      EQUIVALENCE (ONE,UNITY(1)),(ZERO,ZIP(1))
C
      DATA ONE /1.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (K .GE. 0 .AND. K .LT. N .AND. N .LE. NZZR)  GO TO 5
        FAIL = .TRUE.
        RETURN
    5 CONTINUE
      IF (K .GT. 0)  GO TO 20
C
C     ***************
C     FOR THE SPECIAL CASE THAT THE
C     FACTORIZATION WAS VACUOUS,
C     RESET THE ARRAYS  ZZ AND DD
C     TO REPRESENT THE  IDENTITY.
C     ***************
C
        K = 0
        CALL SCOPY(N,UNITY,0,DD,1)
        CALL SCOPY(N*N,ZIP,0,ZZ,1)
        DO 10 I=1,N
          ZZ(I,I) = ONE
   10   CONTINUE
   20 CONTINUE
      KP1 = K + 1
      KP2 = K + 2
C
C     ***************
C     TRANSFORM THE INCOMING COLUMN,
C     AND STORE THE RESULT IN  W.
C     ***************
C
      DO 30 I=1,N
        W(I) = SDOT(N,ZZ(1,I),1,COL,1)
   30 CONTINUE
C
C     ***************
C     ZERO OUT THE SPIKE WHICH WOULD RESULT FROM
C     STORING  W  IN  RR.   UPDATE  ZZ  AND  DD.
C     ***************
C
      IF (KP2 .GT. N)  GO TO 50
        DO 40 I=KP2,N
          DI = DD(I)
          WI = W(I)
          CALL SROTMG(DD(KP1),DI,W(KP1),WI,PARAM)
          W(I) = WI
          DD(I) = DI
          CALL SROTM(N,ZZ(1,KP1),1,ZZ(1,I),1,PARAM)
   40   CONTINUE
   50 CONTINUE
C
C     ***************
C     STORE THE NEW COLUMN, WHICH IS STILL
C     IN THE ARRAY  W,  INTO  RR.
C     ***************
C
      J = KP2
      JDEL = N
      DO 60 I=1,KP1
        RR(J) = W(I)
        J = J + JDEL
        JDEL = JDEL - 1
   60 CONTINUE
      K = KP1
      FAIL = .FALSE.
      RETURN
      END
      SUBROUTINE ZDRCOU(N,K,ZZ,NZZR,DD,RR,IC,FAIL)
C
      INTEGER IC,K,N,NZZR
      LOGICAL FAIL
      REAL DD(1),RR(1),ZZ(NZZR,1)
C
C     ***************
C     PREPARED BY RICHARD BARTELS
C     THE UNIVERSITY OF WATERLOO
C     COMPUTER SCIENCE DEPARTMENT
C     LATEST UPDATE .... 30 NOVEMBER, 1979.
C
C     GIVEN THE FACTORIZATION
C
C          ZZ*DD*RR
C
C     OF SOME  N BY K  MATRIX
C
C       (1 .LE. K .LE. N)
C          (N .GE. 1),
C
C     WHERE
C
C          (ZZ-TRANSP)*(ZZ) = (DD-INV),
C          DD  IS DIAGONAL AND NONSINGULAR,
C     AND
C          RR  HAS ZEROS BELOW THE DIAGONAL,
C
C     AND GIVEN THE INDEX  IC  OF A COLUMN
C     TO BE REMOVED  (1 .LE. IC .LE. K),
C     THIS PROGRAM UPDATES  ZZ,DD AND RR .
C
C     THE VALUE OF  K  IS DECREASED BY ONE, AND
C     THE COLUMN ORDERING IN  RR  IS CHANGED.
C
C     USE IS MADE OF ROUTINES FROM THE LIBRARY
C     OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS).
C
C     PARAMETERS...
C
C       NAME   TYPE    INPUT/   SUB-    DESCRIPTION
C                      OUTPUT  SCRIPTS
C       -------------------------------------------
C       N      INT.      I              NUMBER OF ROWS
C
C       K      INT.     I/O             NUMBER OF COLUMNS
C
C       ZZ     REAL     I/O       2     SCALED ORTHOGONAL
C                                       MATRIX
C
C       NZZR   INT.      I              ROW DIMENSION OF ZZ
C
C       DD     REAL     I/O       1     DIAGONAL SCALING
C                                       MATRIX (DIAGONAL
C                                       ELEMENTS ONLY)
C
C       RR     REAL     I/O       1     RIGHT-TRIANGULAR
C                                       MATRIX IN COMPACT FORM.
C
C       IC     INT.      I              INDEX OF COLUMN
C                                       TO BE REMOVED
C
C       FAIL   LOG.      O              .TRUE.  IF  K,N,IC
C                                       ARE IMPROPER
C       -------------------------------------------
C
C     THE  I-TH  SEGMENT OF THE ARRAY  RR  IS  N-I+2 SPACES
C     LONG AND CONTAINS  1  WORK SPACE FOLLOWED BY THE
C     K-I+1  ELEMENTS OF ROW  I  FOLLOWED BY  N-K
C     SCRATCH SPACES.
C     ***************
C
C     +++++++++++++++
C     BLAS  SROTM,SROTMG
C     +++++++++++++++
C
      INTEGER I,IM1,J,JEND,JINC,JSTRT,KM1,LSTRT
      REAL DI,PARAM(5),RJ
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (K .GE. 1 .AND. K .LE. N .AND. N .LE. NZZR)  GO TO 5
        FAIL = .TRUE.
        RETURN
    5 CONTINUE
      KM1 = K - 1
C
C     ***************
C     SPECIAL CASES ARE HANDLED FIRST.
C     1.  K=1 AND THE FACTORIZATION BECOMES NULL.
C     2.  IC=K AND THE UPDATING IS TRIVIAL.
C     ***************
C
      IF (K .GT. 1)  GO TO 10
        K = 0
        FAIL = .FALSE.
        RETURN
   10 CONTINUE
      IF (IC .LT. K)  GO TO 20
        K = KM1
        FAIL = .FALSE.
        RETURN
   20 CONTINUE
C
C     ***************
C     GENERAL UPDATING STEP.
C     THE COLUMN TO BE DELETED MUST BE PERMUTED
C     TO THE RIGHT, AND SUBDIAGONAL ELEMENTS
C     WHICH RESULT IN  RR  HAVE TO BE
C     TRANSFORMED TO ZERO.
C     ***************
C
      JSTRT = IC + 1
      JEND = K
      JINC = N
      DO 50 I=1,K
C
C     ***************
C     PERMUTATION OF THE  I-TH  ROW OF RR.
C     ***************
C
        DO 30 J=JSTRT,JEND
          RR(J) = RR(J + 1)
   30   CONTINUE
        IF (I .LE. IC)  GO TO 40
C
C     ***************
C     TRANSFORMATION OF THE CURRENT AND LAST
C     ROWS  (I AND I-1)  OF RR  AS WELL AS
C     CORRESPONDING CHANGES TO  ZZ AND DD.
C
C     THE EXTRA VARIABLES  DI  AND  RJ
C     ARE USED TO AVOID AN ERROR MESSAGE
C     FROM THE  PFORT VERIFIER, AND THEY
C     MAY BE REMOVED, IF DESIRED, SO THAT
C     THE CALL TO  SROTMG  WOULD BE
C
C     CALL SROTMG(DD(IM1),DD(I),RR(LSTRT),RR(JSTRT),PARAM)
C
C     ***************
C
          IM1 = I - 1
          DI = DD(I)
          RJ = RR(JSTRT)
          CALL SROTMG(DD(IM1),DI,RR(LSTRT),RJ,PARAM)
          RR(JSTRT) = RJ
          DD(I) = DI
          CALL SROTM(JEND - JSTRT + 1,RR(LSTRT + 1),1,
     *               RR(JSTRT + 1),1,PARAM)
          CALL SROTM(N,ZZ(1,IM1),1,ZZ(1,I),1,PARAM)
          JSTRT = JSTRT + 1
   40   CONTINUE
C
C     ***************
C     INDEX UPDATING
C     ***************
C
        LSTRT = JSTRT
        JSTRT = JSTRT + JINC
        JEND = JEND + JINC
        JINC = JINC - 1
   50 CONTINUE
      K = KM1
      FAIL = .FALSE.
      RETURN
      END
      SUBROUTINE ZDRGIT(N,K,ZZ,NZZR,RR,GV,SOL,FAIL,W)
C
      INTEGER K,N,NZZR
      LOGICAL FAIL
      REAL GV(1),RR(1),W(1),SOL(1),ZZ(NZZR,1)
C
C     ***************
C     PREPARED BY RICHARD BARTELS
C     THE UNIVERSITY OF WATERLOO
C     COMPUTER SCIENCE DEPARTMENT
C     LATEST UPDATE .... 30 NOVEMBER, 1979.
C
C     GIVEN THE FACTORIZATION
C
C          ZZ*DD*RR
C
C     OF SOME  N BY K  MATRIX
C
C       (1 .LE. K .LE. N)
C         (N .GE. 1),
C
C     WHERE
C
C          (ZZ-TRANSP)*(ZZ) = (DD-INV),
C          DD  IS DIAGONAL AND NONSINGULAR,
C     AND
C          RR  HAS ZEROS BELOW THE DIAGONAL,
C
C     AND GIVEN AN ARBITRARY VECTOR  GV  OF
C     APPROPRIATE DIMENSION, THIS ROUTINE FINDS THE
C     VECTOR  SOL  SATISFYING THE UNDERDETERMINED SYSTEM
C
C          (ZZ*DD*RR-TRANSP.)*(SOL) = (GV).
C
C     THAT IS,
C
C          (SOL) = ((ZZ*DD*RR)-GEN.INV.-TRANSP.)*(GV).
C
C     THE ARRAY  DD  IS NOT NEEDED BY  ZDRGIT.
C
C     USE IS MADE OF ROUTINES FROM THE LIBRARY
C     OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS).
C
C     W  IS A SCRATCH ARRAY.
C
C     PARAMETERS...
C
C                     INPUT/
C       NAME   TYPE   OUTPUT/   SUB-    DESCRIPTION
C                     SCRATCH  SCRIPTS
C       -------------------------------------------
C       N      INT.      I              NUMBER OF ROWS
C
C       K      INT.     I/O             NUMBER OF COLUMNS
C
C       ZZ     REAL     I/O       2     SCALED ORTHOGONAL
C                                       MATRIX
C
C       NZZR   INT.      I              ROW DIMENSION OF ZZ
C
C       RR     REAL     I/O       1     RIGHT-TRIANGULAR
C                                       MATRIX IN COMPACT FORM.
C
C       GV     REAL      I        1     GIVEN VECTOR
C
C       SOL    REAL      O        1     SOLUTION
C
C       FAIL   LOG.      O              .TRUE. IF  N,K
C                                       ARE IMPROPER, OR IF
C                                       RR  IS SINGULAR
C
C       W      REAL     SCR       1     WORKSPACE
C       -------------------------------------------
C
C     THE  I-TH  SEGMENT OF THE ARRAY  RR  IS  N-I+2 SPACES
C     LONG AND CONTAINS  1  WORK SPACE FOLLOWED BY THE
C     K-I+1  ELEMENTS OF ROW  I  FOLLOWED BY  N-K
C     SCRATCH SPACES.
C
C     IF  GV  AND  SOL  ARE DIMENSIONED TO THE
C     MAXIMUM OF  N  AND  K , THEN THE SAME
C     STORAGE ARRAY MAY BE USED FOR BOTH OF
C     THESE VECTORS.
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C
C     BLAS  SAXPY,SCOPY
C
C     BIG  IS THE LARGEST POSITIVE NUMBER
C     WHICH CAN BE REPRESENTED IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     +++++++++++++++
C
      INTEGER I,J,JDEL
      REAL BIG,WI,ONE,RRJ,ZERO
C
C     +++++++++++++++
C     THE FOLLOWING DECLARATIONS ARE NECESSARY
C     FOR PORTABILITY WHEN  SCOPY  IS USED, AS
C     IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE
C     (ZERO=ZIP  IN THIS CASE).
C     +++++++++++++++
C
      REAL ZIP(1)
      EQUIVALENCE (ZERO,ZIP(1))
C
      DATA BIG /1.0E+31/
      DATA ONE /1.0E+00/
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (K .GE. 1 .AND. K .LE. N .AND. N .LE. NZZR)  GO TO 10
        FAIL = .TRUE.
        RETURN
   10 CONTINUE
C
C     ***************
C     FIRST SOLVE  (RR-TRANSP.)*(W) = (GV)
C     ***************
C
      CALL SCOPY(K,GV,1,W,1)
      J = 2
      JDEL = N + 1
      DO 30 I=1,K
        RRJ = RR(J)
        WI = W(I)
        IF (ABS(RRJ) .GE. ONE)  GO TO 20
        IF (ABS(WI) .LT. ABS(RRJ)*BIG)  GO TO 20
          FAIL = .TRUE.
          RETURN
   20   CONTINUE
        W(I) = WI/RRJ
        IF (I.LT.K) CALL SAXPY(K-I,(-W(I)),RR(J+1),1,W(I+1),1)
        J = J + JDEL
        JDEL = JDEL - 1
   30 CONTINUE
C
C     ***************
C     NOW  (SOL) = (ZZ)*(W)
C     ***************
C
      CALL SCOPY(N,ZIP,0,SOL,1)
      DO 40 I=1,K
        CALL SAXPY(N,W(I),ZZ(1,I),1,SOL,1)
   40 CONTINUE
      FAIL = .FALSE.
      RETURN
      END
      SUBROUTINE ZDRGNV(N,K,ZZ,NZZR,RR,GV,SOL,FAIL)
C
      INTEGER K,N,NZZR
      LOGICAL FAIL
      REAL GV(1),RR(1),SOL(1),ZZ(NZZR,1)
C
C     ***************
C     PREPARED BY RICHARD BARTELS
C     THE UNIVERSITY OF WATERLOO
C     COMPUTER SCIENCE DEPARTMENT
C     LATEST UPDATE .... 30 NOVEMBER, 1979.
C
C     GIVEN THE FACTORIZATION
C
C          ZZ*DD*RR
C
C     OF SOME  N BY K  MATRIX
C
C       (1 .LE. K .LE. N)
C         (N .GE. 1),
C
C     WHERE
C
C          (ZZ-TRANSP)*(ZZ) = (DD-INV),
C          DD  IS DIAGONAL AND NONSINGULAR,
C     AND
C          RR  HAS ZEROS BELOW THE DIAGONAL,
C
C     AND GIVEN AN ARBITRARY VECTOR  GV  OF
C     APPROPRIATE DIMENSION, THIS ROUTINE FINDS THE
C     VECTOR  SOL  GIVEN BY
C
C          (SOL) = ((ZZ*DD*RR)-GEN.INV.)*(GV),
C
C     WHICH REPRESENTS THE LEAST SQUARES PROBLEM
C
C        (ZZ*DD*RR)*(SOL) = (GV).
C
C     THE ARRAY  DD  IS NOT NEEDED BY  ZDRGNV.
C
C     USE IS MADE OF ROUTINES FROM THE LIBRARY
C     OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS).
C
C     PARAMETERS...
C
C       NAME   TYPE   INPUT/    SUB-    DESCRIPTION
C                     OUTPUT/  SCRIPTS
C       -------------------------------------------
C       N      INT.      I              NUMBER OF ROWS
C
C       K      INT.     I/O             NUMBER OF COLUMNS
C
C       ZZ     REAL     I/O       2     SCALED ORTHOGONAL
C                                       MATRIX
C
C       NZZR   INT.      I              ROW DIMENSION OF ZZ
C
C       RR     REAL     I/O       1     RIGHT-TRIANGULAR
C                                       MATRIX IN COMPACT FORM.
C
C       GV     REAL      I        1     GIVEN VECTOR
C
C       SOL    REAL      O        1     SOLUTION
C
C       FAIL   LOG.      O              .TRUE. IF  N,K
C                                       ARE IMPROPER, OR IF
C                                       RR  IS SINGULAR
C       -------------------------------------------
C
C     THE  I-TH  SEGMENT OF THE ARRAY  RR  IS  N-I+2 SPACES
C     LONG AND CONTAINS  1  WORK SPACE FOLLOWED BY THE
C     K-I+1  ELEMENTS OF ROW  I  FOLLOWED BY  N-K
C     SCRATCH SPACES.
C     ***************
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C
C     BLAS  SDOT
C
C     BIG  IS THE LARGEST POSITIVE NUMBER
C     WHICH CAN BE REPRESENTED IN THE
C     PRECISION OF THE ARITHMETIC BEING USED.
C     +++++++++++++++
C
      INTEGER I,IX,J,JDEL
      REAL BIG,ONE,TDEN,TNUM
C
      REAL SDOT
C
      DATA BIG /1.0E+31/
      DATA ONE /1.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      IF (K .GE. 1 .AND. K .LE. N .AND. N .LE. NZZR)  GO TO 10
        FAIL = .TRUE.
        RETURN
   10 CONTINUE
C
C     ***************
C     FORM   (V) = (ZZ(1)-TRANSP)*(GV),   WHERE  ZZ(1)
C     IS THE MATRIX OF THE FIRST  K  COLUMNS OF  ZZ
C
C     V  CAN BE STORED IN THE ARRAY  SOL.
C     ***************
C
      DO 20 I=1,K
        SOL(I) = SDOT(N,ZZ(1,I),1,GV,1)
   20 CONTINUE
C
C     ***************
C     BACKSOLVE THE SYSTEM
C          (RR)*(SOL) = (V)
C     FOR THE VECTOR  SOL
C
C     NOTE THAT  SOL  AND  V
C     ARE STORED IN THE SAME ARRAY.
C     ***************
C
      J = (((N + 1)*(N + 2) - (N - K + 3)*(N - K + 2))/2) + 2
      JDEL = N - K + 3
      DO 40 IX=1,K
        I = K - IX + 1
        TDEN = RR(J)
        TNUM = SOL(I)
        IF (IX.GT.1) TNUM = TNUM-SDOT(IX-1,RR(J+1),1,SOL(I+1),1)
        IF (ABS(TDEN) .GE. ONE)  GO TO 30
        IF (ABS(TNUM) .LT. ABS(TDEN)*BIG)  GO TO 30
          FAIL = .TRUE.
          RETURN
   30   CONTINUE
        SOL(I) = TNUM/TDEN
        J = J - JDEL
        JDEL = JDEL + 1
   40 CONTINUE
      FAIL = .FALSE.
      RETURN
      END
      SUBROUTINE ZDRPOC(N,K,ZZ,NZZR,DD,GV,POC,FAIL)
C
      INTEGER K,N,NZZR
      LOGICAL FAIL
      REAL DD(1),POC(1),GV(1),ZZ(NZZR,1)
C
C     ***************
C     PREPARED BY RICHARD BARTELS
C     THE UNIVERSITY OF WATERLOO
C     COMPUTER SCIENCE DEPARTMENT
C     LATEST UPDATE .... 30 NOVEMBER, 1979.
C
C     ZZ IS AN  N BY N  (N .GE. 1)  SCALED
C     ORTHOGONAL MATRIX.  DD  CONTAINS THE
C     DIAGONAL ELEMENTS OF A DIAGONAL SCALING
C     MATRIX.  GV  IS A GIVEN VECTOR OF LENGTH  N.
C
C     WE HAVE
C
C          (ZZ-TRANSP.)*(ZZ) = (DD-INV.)
C
C     AND
C
C               ZZ*DD*RR = MAT
C
C     FOR SOME  N BY K  (0 .LE. K .LE. N)
C     MATRIX  RR  WITH ZEROS BELOW THE DIAGONAL
C     AND SOME GIVEN MATRIX  MAT.  (NIETHER  RR
C     NOR  MAT  ARE NEEDED BY  ZDRPOC.)
C
C     THEN
C
C    (PROJ(OC)) = (ZZ(2))*(DD(2))*(ZZ(2)-TRANSP.)
C
C     IS THE (ORTHOGONAL) PROJECTOR ON THE
C     COMPLEMENT OF THE RANGE SPACE OF  MAT,
C     WHERE  ZZ(2)  REPRESENTS THE LAST  N-K
C     COLUMNS OF  ZZ  AND  DD(2)  REPRESENTS THE
C     LOWER-RIGHT-HAND  N-K  ORDER SUBMATRIX OF  DD.
C
C     ZDRPOC  PRODUCES THE VECTOR
C
C               POC = (PROJ(OC))*GV .
C
C     USE IS MADE OF ROUTINES FROM THE LIBRARY
C     OF BASIC LINEAR ALGEBRA SUBROUTINES (BLAS).
C
C     PARAMETERS...
C
C                     INPUT/
C       NAME   TYPE   OUTPUT/   SUB-    DESCRIPTION
C                     SCRATCH  SCRIPTS
C       -------------------------------------------
C       N      INT.      I              ORDER OF  ZZ,DD
C
C       K      INT.     I/O             NUMBER OF COLUMNS
C                                       OF  ZZ  DEFINING
C                                       RANGE OF  MAT
C
C       ZZ     REAL     I/O       2     SCALED ORTHOGONAL
C                                       MATRIX
C
C       NZZR   INT.      I              ROW DIMENSION OF ZZ
C
C       DD     REAL     I/O       1     DIAGONAL SCALING
C                                       MATRIX (DIAGONAL
C                                       ELEMENTS ONLY)
C
C       GV     REAL      I        1     VECTOR TO BE PROJECTED
C
C       POC    REAL      O        1     PROJECTION
C
C       FAIL   LOG.      O              .TRUE.  IF  N,K
C                                       ARE IMPROPER
C
C       -------------------------------------------
C
C     ***************
C
C     +++++++++++++++
C     BLAS  SAXPY,SCOPY,SDOT
C     +++++++++++++++
C
      INTEGER I,KP1
      REAL WI,ZERO
C
      REAL SDOT
C
C     +++++++++++++++
C     THE FOLLOWING DECLARATIONS ARE NECESSARY
C     FOR PORTABILITY WHEN  SCOPY  IS USED, AS
C     IT IS BELOW, TO FILL ARRAYS WITH A SINGLE VALUE
C     (ZERO=ZIP  IN THIS CASE).
C     +++++++++++++++
C
      REAL ZIP(1)
      EQUIVALENCE (ZERO,ZIP(1))
C
      DATA ZERO /0.0E+00/
C
C     /////////////////  BEGIN PROGRAM  //////////////////
C
      KP1 = K + 1
      IF (K .GE. 0 .AND. K .LE. N
     *       .AND. N .GE. 1 .AND. N .LE. NZZR)  GO TO 5
        FAIL = .TRUE.
        RETURN
    5 CONTINUE
      IF (K .GT. 0)  GO TO 10
C
C     ***************
C     CASE 1 ... ZZ(2)=ZZ  (K=0)
C     ***************
C
        CALL SCOPY(N,GV,1,POC,1)
        FAIL = .FALSE.
        RETURN
   10 CONTINUE
      IF (K .LT. N)  GO TO 20
C
C     ***************
C     CASE 2 ... ZZ(2) IS VACUOUS  (K=N)
C     ***************
C
        CALL SCOPY(N,ZIP,0,POC,1)
        FAIL = .FALSE.
        RETURN
   20 CONTINUE
C
C     ***************
C     CASE 3 ... ZZ(2)  IS INTERMEDIATE
C     BETWEEN THE OTHER TWO CASES
C     (0 .LT. K .LT. N)
C     ***************
C
      CALL SCOPY(N,ZIP,0,POC,1)
      DO 30 I=KP1,N
        WI = SDOT(N,ZZ(1,I),1,GV,1)*DD(I)
        CALL SAXPY(N,WI,ZZ(1,I),1,POC,1)
   30 CONTINUE
      FAIL = .FALSE.
      RETURN
      END
C
C     ---------------------------------------------
C                  ** APPENDIX **
C
C     BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) --
C               SASUM,SAXPY,SCOPY,SDOT,
C               SROTM,SROTMG,SSCAL
C
C     THE  BLAS  ARE DESCRIBED IN
C
C              BASIC LINEAR ALGEBRA SUBPROGRAMS
C                    FOR FORTRAN USAGE
C
C                           BY
C              C. L. LAWSON, R. J. HANSON,
C              D. R. KINCAID, F. T. KROGH.
C              ACM TRANS. ON MATH SOFTWARE
C              VOL. 5, NO. 3, 308-325 (1979)
C
C              THE  BLAS  ARE TAKEN AS
C              PRIMITIVES.  THEIR CODE IS NOT
C              PART OF THE DEFINITION OF  CL1.
C              THEY ARE INCLUDED FOR REFERENCE
C              ONLY.  PLEASE REFER TO THE ABOVE
C              PUBLICATION FOR CURRENT INFORMATION.
C
C     IF THE  BLAS  ARE AVAILABLE INDEPENDANTLY
C     ON THE TARGET MACHINE, THIS APPENDIX SHOULD
C     BE REMOVED FROM THE CODE FOR  CL1.
C     ---------------------------------------------
C
      REAL FUNCTION SASUM(N,SX,INCX)
C
      INTEGER INCX,N
      REAL SX(1)
C
C     RETURNS SUM OF MAGNITUDES OF SINGLE PRECISION SX.
C
      INTEGER I,M,MP1,NS
      SASUM = 0.0E+00
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GOTO 20
C
C        CODE FOR INCREMENTS NOT EQUAL TO 1.
C
      NS = N*INCX
          DO 10 I=1,NS,INCX
          SASUM = SASUM + ABS(SX(I))
   10     CONTINUE
      RETURN
C
C        CODE FOR INCREMENTS EQUAL TO 1.
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6.
C
   20 M = N - (N/6)*6
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SASUM = SASUM + ABS(SX(I))
   30 CONTINUE
      IF( N .LT. 6 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,6
        SASUM = SASUM + ABS(SX(I)) + ABS(SX(I + 1)) + ABS(SX(I + 2))
     *  + ABS(SX(I + 3)) + ABS(SX(I + 4)) + ABS(SX(I + 5))
   50 CONTINUE
      RETURN
      END
      SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
C
      INTEGER INCX,INCY,N
      REAL SA,SX(1),SY(1)
C
C     OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY.
C
      INTEGER I,IX,IY,M,MP1,NS
C
      IF(N.LE.0.OR.SA.EQ.0.E0) RETURN
      IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
    5 CONTINUE
C
C        CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        SY(IY) = SY(IY) + SA*SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4.
C
   20 M = N - (N/4)*4
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SY(I) = SY(I) + SA*SX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        SY(I) = SY(I) + SA*SX(I)
        SY(I + 1) = SY(I + 1) + SA*SX(I + 1)
        SY(I + 2) = SY(I + 2) + SA*SX(I + 2)
        SY(I + 3) = SY(I + 3) + SA*SX(I + 3)
   50 CONTINUE
      RETURN
C
C        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
   60 CONTINUE
      NS = N*INCX
          DO 70 I=1,NS,INCX
          SY(I) = SA*SX(I) + SY(I)
   70     CONTINUE
      RETURN
      END
      SUBROUTINE  SCOPY(N,SX,INCX,SY,INCY)
C
      INTEGER INCX,INCY,N
      REAL SX(1),SY(1)
C
C     COPY SINGLE PRECISION SX TO SINGLE PRECISION SY.
C
      INTEGER I,IX,IY,M,MP1,NS
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60
    5 CONTINUE
C
C        CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        SY(IY) = SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7.
C
   20 M = N - (N/7)*7
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SY(I) = SX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        SY(I) = SX(I)
        SY(I + 1) = SX(I + 1)
        SY(I + 2) = SX(I + 2)
        SY(I + 3) = SX(I + 3)
        SY(I + 4) = SX(I + 4)
        SY(I + 5) = SX(I + 5)
        SY(I + 6) = SX(I + 6)
   50 CONTINUE
      RETURN
C
C        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
   60 CONTINUE
      NS = N*INCX
          DO 70 I=1,NS,INCX
          SY(I) = SX(I)
   70     CONTINUE
      RETURN
      END
      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C
      INTEGER INCX,INCY,N
      REAL SX(1),SY(1)
C
C     RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY.
C
      INTEGER I,IX,IY,M,MP1,NS
C
      SDOT = 0.0E+00
      IF(N.LE.0)RETURN
      IF(INCX.EQ.INCY) IF(INCX-1)5,20,60
    5 CONTINUE
C
C        CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        SDOT = SDOT + SX(IX)*SY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
   20 M = N - (N/5)*5
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SDOT = SDOT + SX(I)*SY(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        SDOT = SDOT + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
     *   SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4)
   50 CONTINUE
      RETURN
C
C        CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
   60 CONTINUE
      NS=N*INCX
      DO 70 I=1,NS,INCX
        SDOT = SDOT + SX(I)*SY(I)
   70   CONTINUE
      RETURN
      END
      SUBROUTINE SROTM (N,SX,INCX,SY,INCY,SPARAM)
C
      INTEGER INCX,INCY,N
      REAL SPARAM(5),SX(1),SY(1)
C
C     APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
C
C     (SX(1)     SX(N))
C     (      ...      )
C     (SY(1)     SY(N))
C
C     WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
C
C     SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
C
C       (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
C     H=(          )    (          )    (          )    (          )
C       (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
C
      INTEGER I,KX,KY,NSTEPS
      REAL SFLAG,SH11,SH12,SH21,SH22,TWO,W,Z,ZERO
C
      DATA ZERO,TWO /0.0E+00,2.0E+00/
C
      SFLAG=SPARAM(1)
      IF(N .LE. 0 .OR.(SFLAG+TWO.EQ.ZERO)) GO TO 140
          IF(.NOT.(INCX.EQ.INCY.AND. INCX .GT.0)) GO TO 70
C
               NSTEPS=N*INCX
               IF(SFLAG) 50,10,30
   10          CONTINUE
               SH12=SPARAM(4)
               SH21=SPARAM(3)
                    DO 20 I=1,NSTEPS,INCX
                    W=SX(I)
                    Z=SY(I)
                    SX(I)=W+Z*SH12
                    SY(I)=W*SH21+Z
   20               CONTINUE
               GO TO 140
   30          CONTINUE
               SH11=SPARAM(2)
               SH22=SPARAM(5)
                    DO 40 I=1,NSTEPS,INCX
                    W=SX(I)
                    Z=SY(I)
                    SX(I)=W*SH11+Z
                    SY(I)=-W+SH22*Z
   40               CONTINUE
               GO TO 140
   50          CONTINUE
               SH11=SPARAM(2)
               SH12=SPARAM(4)
               SH21=SPARAM(3)
               SH22=SPARAM(5)
                    DO 60 I=1,NSTEPS,INCX
                    W=SX(I)
                    Z=SY(I)
                    SX(I)=W*SH11+Z*SH12
                    SY(I)=W*SH21+Z*SH22
   60               CONTINUE
               GO TO 140
   70     CONTINUE
          KX=1
          KY=1
          IF(INCX .LT. 0) KX=1+(1-N)*INCX
          IF(INCY .LT. 0) KY=1+(1-N)*INCY
C
          IF(SFLAG)120,80,100
   80     CONTINUE
          SH12=SPARAM(4)
          SH21=SPARAM(3)
               DO 90 I=1,N
               W=SX(KX)
               Z=SY(KY)
               SX(KX)=W+Z*SH12
               SY(KY)=W*SH21+Z
               KX=KX+INCX
               KY=KY+INCY
   90          CONTINUE
          GO TO 140
  100     CONTINUE
          SH11=SPARAM(2)
          SH22=SPARAM(5)
               DO 110 I=1,N
               W=SX(KX)
               Z=SY(KY)
               SX(KX)=W*SH11+Z
               SY(KY)=-W+SH22*Z
               KX=KX+INCX
               KY=KY+INCY
  110          CONTINUE
          GO TO 140
  120     CONTINUE
          SH11=SPARAM(2)
          SH12=SPARAM(4)
          SH21=SPARAM(3)
          SH22=SPARAM(5)
               DO 130 I=1,N
               W=SX(KX)
               Z=SY(KY)
               SX(KX)=W*SH11+Z*SH12
               SY(KY)=W*SH21+Z*SH22
               KX=KX+INCX
               KY=KY+INCY
  130          CONTINUE
  140     CONTINUE
          RETURN
          END
      SUBROUTINE SROTMG (SD1,SD2,SX1,SY1,SPARAM)
C
      REAL SD1,SD2,SPARAM(5),SX1,SY1
C
C     CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H
C     WHICH ZEROS THE SECOND COMPONENT OF THE 2-VECTOR
C         (SQRT(SD1)*SX1,SQRT(SD2)*SY2)**T.
C     WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
C
C     SFLAG=-1.E0     SFLAG=0.E0        SFLAG=1.E0     SFLAG=-2.E0
C
C       (SH11  SH12)    (1.E0  SH12)    (SH11  1.E0)    (1.E0  0.E0)
C     H=(          )    (          )    (          )    (          )
C       (SH21  SH22),   (SH21  1.E0),   (-1.E0 SH22),   (0.E0  1.E0).
C
C     +++++++++++++++
C     SYSTEM ROUTINES  ABS
C     +++++++++++++++
C
      INTEGER IFLAG,IGO
      REAL GAM,GAMSQ,ONE,RGAM,RGAMSQ,SFLAG,SH11,SH12
      REAL SH21,SH22,SP1,SP2,SQ1,SQ2,STEMP,SU,TWO,ZERO
C
      DATA ZERO,ONE,TWO /0.0E+00,1.0E+00,2.0E+00/,IFLAG/1/
      DATA GAM/4096.0E+00/
      DATA GAMSQ /1.678E+07/
      DATA RGAM /2.441E-04/
      DATA RGAMSQ /5.960E-08/
C
      IF(.NOT. SD1 .LT. ZERO) GO TO 10
C       GO ZERO-H-D-AND-SX1..
          GO TO 60
   10 CONTINUE
C     CASE-SD1-NONNEGATIVE
      SP2=SD2*SY1
      IF(.NOT. SP2 .EQ. ZERO) GO TO 20
          SFLAG=-TWO
          GO TO 260
C     REGULAR-CASE..
   20 CONTINUE
      SP1=SD1*SX1
      SQ2=SP2*SY1
      SQ1=SP1*SX1
C
      IF(.NOT. ABS(SQ1) .GT. ABS(SQ2)) GO TO 40
          SH21=-SY1/SX1
          SH12=SP2/SP1
C
          SU=ONE-SH12*SH21
C
          IF(.NOT. SU .LE. ZERO) GO TO 30
C         GO ZERO-H-D-AND-SX1..
               GO TO 60
   30     CONTINUE
               SFLAG=ZERO
               SD1=SD1/SU
               SD2=SD2/SU
               SX1=SX1*SU
C         GO SCALE-CHECK..
               GO TO 100
   40 CONTINUE
          IF(.NOT. SQ2 .LT. ZERO) GO TO 50
C         GO ZERO-H-D-AND-SX1..
               GO TO 60
   50     CONTINUE
               SFLAG=ONE
               SH11=SP1/SP2
               SH22=SX1/SY1
               SU=ONE+SH11*SH22
               STEMP=SD2/SU
               SD2=SD1/SU
               SD1=STEMP
               SX1=SY1*SU
C         GO SCALE-CHECK
               GO TO 100
C     PROCEDURE..ZERO-H-D-AND-SX1..
   60 CONTINUE
          SFLAG=-ONE
          SH11=ZERO
          SH12=ZERO
          SH21=ZERO
          SH22=ZERO
C
          SD1=ZERO
          SD2=ZERO
          SX1=ZERO
C         RETURN..
          GO TO 220
C     PROCEDURE..FIX-H..
   70 CONTINUE
      IF(.NOT. SFLAG .GE. ZERO) GO TO 90
C
          IF(.NOT. SFLAG .EQ. ZERO) GO TO 80
          SH11=ONE
          SH22=ONE
          SFLAG=-ONE
          GO TO 90
   80     CONTINUE
          SH21=-ONE
          SH12=ONE
          SFLAG=-ONE
   90 CONTINUE
      GO TO IGO,(120,150,180,210)
C     PROCEDURE..SCALE-CHECK
  100 CONTINUE
  110     CONTINUE
          IF(.NOT. SD1 .LE. RGAMSQ) GO TO 130
               IF(.NOT. IFLAG.EQ.1) GO TO 105
C
C                   RECOMPUTE RESCALING PARAMETERS
C                   MORE ACCURATELY..
C
                    RGAM = ONE/GAM
                    GAMSQ = GAM**2
                    RGAMSQ = RGAM**2
                    IFLAG = 2
  105          CONTINUE
               IF(SD1 .EQ. ZERO) GO TO 160
               ASSIGN 120 TO IGO
C              FIX-H..
               GO TO 70
  120          CONTINUE
               SD1=SD1*GAMSQ
               SX1=SX1*RGAM
               SH11=SH11*RGAM
               SH12=SH12*RGAM
          GO TO 110
  130 CONTINUE
  140     CONTINUE
          IF(.NOT. SD1 .GE. GAMSQ) GO TO 160
               ASSIGN 150 TO IGO
C              FIX-H..
               GO TO 70
  150          CONTINUE
               SD1=SD1*RGAMSQ
               SX1=SX1*GAM
               SH11=SH11*GAM
               SH12=SH12*GAM
          GO TO 140
  160 CONTINUE
  170     CONTINUE
          IF(.NOT. ABS(SD2) .LE. RGAMSQ) GO TO 190
               IF(SD2 .EQ. ZERO) GO TO 220
               ASSIGN 180 TO IGO
C              FIX-H..
               GO TO 70
  180          CONTINUE
               SD2=SD2*GAMSQ
               SH21=SH21*RGAM
               SH22=SH22*RGAM
          GO TO 170
  190 CONTINUE
  200     CONTINUE
          IF(.NOT. ABS(SD2) .GE. GAMSQ) GO TO 220
               ASSIGN 210 TO IGO
C              FIX-H..
               GO TO 70
  210          CONTINUE
               SD2=SD2*RGAMSQ
               SH21=SH21*GAM
               SH22=SH22*GAM
          GO TO 200
  220 CONTINUE
          IF(SFLAG)250,230,240
  230     CONTINUE
               SPARAM(3)=SH21
               SPARAM(4)=SH12
               GO TO 260
  240     CONTINUE
               SPARAM(2)=SH11
               SPARAM(5)=SH22
               GO TO 260
  250     CONTINUE
               SPARAM(2)=SH11
               SPARAM(3)=SH21
               SPARAM(4)=SH12
               SPARAM(5)=SH22
  260 CONTINUE
          SPARAM(1)=SFLAG
          RETURN
      END
      SUBROUTINE  SSCAL(N,SA,SX,INCX)
C
      INTEGER INCX,N
      REAL SA,SX(1)
C
C     REPLACE SINGLE PRECISION SX BY SINGLE PRECISION SA*SX.
C
      INTEGER I,M,MP1,NS
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GOTO 20
C
C        CODE FOR INCREMENTS NOT EQUAL TO 1.
C
      NS = N*INCX
          DO 10 I = 1,NS,INCX
          SX(I) = SA*SX(I)
   10     CONTINUE
      RETURN
C
C        CODE FOR INCREMENTS EQUAL TO 1.
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
   20 M = N - (N/5)*5
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SX(I) = SA*SX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        SX(I) = SA*SX(I)
        SX(I + 1) = SA*SX(I + 1)
        SX(I + 2) = SA*SX(I + 2)
        SX(I + 3) = SA*SX(I + 3)
        SX(I + 4) = SA*SX(I + 4)
   50 CONTINUE
      RETURN
      END
ian@ian-HP-Stream-Notebook-PC-11:~/563dir$