c program DRDBLAS2 c>> 1996-06-18 DRDBLAS2 Krogh Minor format change for C conversion. c>> 1994-10-19 DRDBLAS2 Krogh Changes to use M77CON c>> 1991-11-27 DRDBLAS2 CLL c>> 1987-12-09 Lawson Initial Code. c c Demonstrates the use of BLAS subroutines DROTG, DROT, DAXPY, c and DCOPY 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--D 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 double precision X, XTAB(MXY), Y, YTAB(MXY), W(MC1) double precision C, RG(MC1,MC1), S double precision COEF(MC), DIV, ESTSD, ZERO(1) c data XTAB / 0.0d0, .1d0, .2d0, .3d0, .4d0, .5d0, * .6d0, .7d0, .8d0, .9d0, 1.0d0 / data YTAB / 1.00d0, .91d0, .86d0, .82d0, .81d0, .82d0, * .85d0, .89d0, .95d0, 1.02d0, 1.10d0 / data NXY, NC / MXY, MC / data ZERO(1) / 0.0d0 / c ----------------------------------------------------------------- NC1 = NC + 1 call DCOPY(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.0d0 W(2) = X W(3) = exp(-X) W(4) = Y c Process W() into [R:G]. do 10 J = 1, NC call DROTG(RG(J,J),W(J),C,S) call DROT(NC1-J,RG(J,J+1),MC1,W(J+1),1,C,S) 10 continue call DROTG(RG(NC1,NC1),W(NC1),C,S) 20 continue c Begin: Solve triangular system. call DCOPY(NC,RG(1,NC1),1,COEF,1) do 30 J = NC, 1, -1 DIV = RG(J,J) if (DIV .eq. 0.0d0) then print '(''ERROR:ZERO DIVISOR AT J ='', I2)', J stop end if COEF(J) = COEF(J) / DIV call DAXPY(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(DBLE(NXY-NC)) print'(/'' Estimated Std. Dev. of data errors ='',f9.5)', ESTSD stop end