c     program DRSBLAS2
c>> 1996-06-18 DRSBLAS2 Krogh  Minor format change for C conversion.
c>> 1994-10-19 DRSBLAS2 Krogh  Changes to use M77CON
c>> 1991-11-27 DRSBLAS2 CLL
c>> 1987-12-09 Lawson  Initial Code.
c
c     Demonstrates the use of BLAS subroutines SROTG, SROT, SAXPY,
c     and SCOPY to implement an algorithm for solving a linear
c     least squares problem using sequential accumulation of the
c     data and Givens orthogonal transformations.
c     YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X)
c     -----------------------------------------------------------------
c--S replaces "?": DR?BLAS2, ?ROTG, ?ROT, ?AXPY, ?COPY
c     -----------------------------------------------------------------
      integer MC, MC1, MXY
      parameter ( MC=3, MC1=MC+1, MXY=11 )
      integer IXY, J, NC, NC1, NXY
      real             X, XTAB(MXY), Y, YTAB(MXY), W(MC1)
      real             C, RG(MC1,MC1), S
      real             COEF(MC), DIV, ESTSD, ZERO(1)
c
      data XTAB / 0.0e0,  .1e0,  .2e0,    .3e0,   .4e0,  .5e0,
     *             .6e0,  .7e0,  .8e0,    .9e0,  1.0e0 /
      data YTAB / 1.00e0, .91e0, .86e0,   .82e0,  .81e0, .82e0,
     *             .85e0, .89e0, .95e0,  1.02e0, 1.10e0 /
      data NXY, NC / MXY, MC /
      data ZERO(1) / 0.0e0 /
c     -----------------------------------------------------------------
      NC1 = NC + 1
      call SCOPY(MC1*MC1, ZERO, 0, RG, 1)
      do 20 IXY = 1, NXY
        X = XTAB(IXY)
        Y = YTAB(IXY)
c                                      Build new row of [A:B] in W().
        W(1) = 1.0e0
        W(2) = X
        W(3) = exp(-X)
        W(4) = Y
c                                      Process W() into [R:G].
        do 10 J = 1, NC
          call SROTG(RG(J,J),W(J),C,S)
          call SROT(NC1-J,RG(J,J+1),MC1,W(J+1),1,C,S)
   10   continue
        call SROTG(RG(NC1,NC1),W(NC1),C,S)
   20 continue
c                                      Begin: Solve triangular system.
        call SCOPY(NC,RG(1,NC1),1,COEF,1)
        do 30 J = NC, 1, -1
          DIV = RG(J,J)
          if (DIV .eq. 0.0e0) then
            print '(''ERROR:ZERO DIVISOR AT J ='', I2)', J
            stop
          end if
          COEF(J) = COEF(J) / DIV
          call SAXPY(J-1,-COEF(J),RG(1,J),1,COEF,1)
   30   continue
c                                        End: Solve triangular system.
c
      print'('' Solution: COEF() = '',3f8.3)',(COEF(J),J=1,NC)
      ESTSD = abs(RG(NC1,NC1)) / sqrt(REAL(NXY-NC))
      print'(/'' Estimated Std. Dev. of data errors ='',f9.5)', ESTSD
      stop
      end