1c     program DRDBLAS4
2c>> 1994-10-19 DRDBLAS4 Krogh  Changes to use M77CON
3c>> 1992-05-12 DRDBLAS4 CLL Improved consistency of various versions.
4c>> 1992-04-15 CLL
5c>> 1992-03-16 CLL
6c>> 1991-12-02 DRDBLAS4 CLL
7c>> 1991-07-31 DRDBLAS4 CLL
8c
9c     Demonstrates the use of BLAS subroutines DROTMG, DROTM, DAXPY,
10c     and DCOPY to implement an algorithm for solving a linear
11c     least squares problem using sequential accumulation of the
12c     data and Modified (Fast) Givens orthogonal transformations.
13c     YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X)
14c
15c     ------------------------------------------------------------------
16c--D replaces "?": DR?BLAS4, ?AXPY, ?ROTMG, ?ROTM, ?COPY, ?VECP
17c     ------------------------------------------------------------------
18      integer MC, MC1, MC2, MXY
19      parameter ( MC=3, MC1=MC+1, MC2=MC+2, MXY=11 )
20      double precision XTAB(MXY), YTAB(MXY), W(MC2), RG(MC1,MC2)
21      double precision C(MC), ZERO(1), X, Y, ESTSD, DIV, PARAM(5)
22      double precision ONE(1)
23      integer NXY, NC, NC1, NC2, IXY, J
24c
25      data XTAB / 0.0d0, .1d0, .2d0, .3d0, .4d0, .5d0, .6d0, .7d0,
26     *           .8d0, .9d0, 1.0d0 /
27      data YTAB / 1.00d0, .91d0, .86d0, .82d0, .81d0, .82d0, .85d0,
28     * .89d0, .95d0, 1.02d0, 1.10d0 /
29      data NXY, NC / 11, 3 /
30      data ONE(1), ZERO(1) / 1.0d0, 0.0d0 /
31c     ------------------------------------------------------------------
32      print'(a/)',' DRDBLAS4..  Demo for DROTMG and DROTM.'
33      NC1 = NC + 1
34      NC2 = NC1 + 1
35      call DCOPY(MC1*MC2, ZERO, 0, RG, 1)
36      call DCOPY(NC1, ONE, 0, RG(1,NC2), 1)
37      do 20 IXY = 1, NXY
38        X = XTAB(IXY)
39        Y = YTAB(IXY)
40c
41c          Set W() to be the next row of the least-squares problem.
42c
43        W(1) = 1.0d0
44        W(2) = X
45        W(3) = EXP(-X)
46        W(4) = Y
47        W(5) = 1.0D0
48c
49c        Accumulate W() into the triangular matrix [R:G]
50c        using Givens rotations.
51c
52        do 10 J = 1, NC1
53          call DROTMG(RG(J,NC2), W(NC2), RG(J,J), W(J), PARAM)
54          if(J .lt. NC1) call DROTM(NC1-J,RG(J,J+1),MC1,W(J+1),1,PARAM)
55   10   continue
56c
57   20 continue
58c
59c               Solve the triangular system.
60c
61        call DCOPY(NC,RG(1,NC1),1,C,1)
62        do 30 J = NC, 1, -1
63          DIV = RG(J,J)
64          if (DIV .EQ. 0.0d0) then
65            print*,'ERROR:ZERO DIVISOR AT J =',J
66            stop
67          end if
68          C(J) = C(J) / DIV
69          call DAXPY(J-1,-C(J),RG(1,J),1,C,1)
70   30 continue
71c
72      call DVECP(C,NC,'0 Solution: C()=')
73      print*,' '
74      ESTSD = abs(RG(NC1,NC1)) * sqrt(RG(NC1,NC2)) / sqrt(DBLE(NXY-NC))
75      print*,'Estimated STD. DEV. of data errors =',ESTSD
76      stop
77      end
78