1c     program DRDBLAS2
2c>> 1996-06-18 DRDBLAS2 Krogh  Minor format change for C conversion.
3c>> 1994-10-19 DRDBLAS2 Krogh  Changes to use M77CON
4c>> 1991-11-27 DRDBLAS2 CLL
5c>> 1987-12-09 Lawson  Initial Code.
6c
7c     Demonstrates the use of BLAS subroutines DROTG, DROT, DAXPY,
8c     and DCOPY to implement an algorithm for solving a linear
9c     least squares problem using sequential accumulation of the
10c     data and Givens orthogonal transformations.
11c     YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X)
12c     -----------------------------------------------------------------
13c--D replaces "?": DR?BLAS2, ?ROTG, ?ROT, ?AXPY, ?COPY
14c     -----------------------------------------------------------------
15      integer MC, MC1, MXY
16      parameter ( MC=3, MC1=MC+1, MXY=11 )
17      integer IXY, J, NC, NC1, NXY
18      double precision X, XTAB(MXY), Y, YTAB(MXY), W(MC1)
19      double precision C, RG(MC1,MC1), S
20      double precision COEF(MC), DIV, ESTSD, ZERO(1)
21c
22      data XTAB / 0.0d0,  .1d0,  .2d0,    .3d0,   .4d0,  .5d0,
23     *             .6d0,  .7d0,  .8d0,    .9d0,  1.0d0 /
24      data YTAB / 1.00d0, .91d0, .86d0,   .82d0,  .81d0, .82d0,
25     *             .85d0, .89d0, .95d0,  1.02d0, 1.10d0 /
26      data NXY, NC / MXY, MC /
27      data ZERO(1) / 0.0d0 /
28c     -----------------------------------------------------------------
29      NC1 = NC + 1
30      call DCOPY(MC1*MC1, ZERO, 0, RG, 1)
31      do 20 IXY = 1, NXY
32        X = XTAB(IXY)
33        Y = YTAB(IXY)
34c                                      Build new row of [A:B] in W().
35        W(1) = 1.0d0
36        W(2) = X
37        W(3) = exp(-X)
38        W(4) = Y
39c                                      Process W() into [R:G].
40        do 10 J = 1, NC
41          call DROTG(RG(J,J),W(J),C,S)
42          call DROT(NC1-J,RG(J,J+1),MC1,W(J+1),1,C,S)
43   10   continue
44        call DROTG(RG(NC1,NC1),W(NC1),C,S)
45   20 continue
46c                                      Begin: Solve triangular system.
47        call DCOPY(NC,RG(1,NC1),1,COEF,1)
48        do 30 J = NC, 1, -1
49          DIV = RG(J,J)
50          if (DIV .eq. 0.0d0) then
51            print '(''ERROR:ZERO DIVISOR AT J ='', I2)', J
52            stop
53          end if
54          COEF(J) = COEF(J) / DIV
55          call DAXPY(J-1,-COEF(J),RG(1,J),1,COEF,1)
56   30   continue
57c                                        End: Solve triangular system.
58c
59      print'('' Solution: COEF() = '',3f8.3)',(COEF(J),J=1,NC)
60      ESTSD = abs(RG(NC1,NC1)) / sqrt(DBLE(NXY-NC))
61      print'(/'' Estimated Std. Dev. of data errors ='',f9.5)', ESTSD
62      stop
63      end
64