c program DRSBLAS4 c>> 1994-10-19 DRSBLAS4 Krogh Changes to use M77CON c>> 1992-05-12 DRSBLAS4 CLL Improved consistency of various versions. c>> 1992-04-15 CLL c>> 1992-03-16 CLL c>> 1991-12-02 DRSBLAS4 CLL c>> 1991-07-31 DRSBLAS4 CLL c c Demonstrates the use of BLAS subroutines SROTMG, SROTM, SAXPY, c and SCOPY to implement an algorithm for solving a linear c least squares problem using sequential accumulation of the c data and Modified (Fast) Givens orthogonal transformations. c YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X) c c ------------------------------------------------------------------ c--S replaces "?": DR?BLAS4, ?AXPY, ?ROTMG, ?ROTM, ?COPY, ?VECP c ------------------------------------------------------------------ integer MC, MC1, MC2, MXY parameter ( MC=3, MC1=MC+1, MC2=MC+2, MXY=11 ) real XTAB(MXY), YTAB(MXY), W(MC2), RG(MC1,MC2) real C(MC), ZERO(1), X, Y, ESTSD, DIV, PARAM(5) real ONE(1) integer NXY, NC, NC1, NC2, IXY, J 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 / 11, 3 / data ONE(1), ZERO(1) / 1.0e0, 0.0e0 / c ------------------------------------------------------------------ print'(a/)',' DRSBLAS4.. Demo for SROTMG and SROTM.' NC1 = NC + 1 NC2 = NC1 + 1 call SCOPY(MC1*MC2, ZERO, 0, RG, 1) call SCOPY(NC1, ONE, 0, RG(1,NC2), 1) do 20 IXY = 1, NXY X = XTAB(IXY) Y = YTAB(IXY) c c Set W() to be the next row of the least-squares problem. c W(1) = 1.0e0 W(2) = X W(3) = EXP(-X) W(4) = Y W(5) = 1.0E0 c c Accumulate W() into the triangular matrix [R:G] c using Givens rotations. c do 10 J = 1, NC1 call SROTMG(RG(J,NC2), W(NC2), RG(J,J), W(J), PARAM) if(J .lt. NC1) call SROTM(NC1-J,RG(J,J+1),MC1,W(J+1),1,PARAM) 10 continue c 20 continue c c Solve the triangular system. c call SCOPY(NC,RG(1,NC1),1,C,1) do 30 J = NC, 1, -1 DIV = RG(J,J) if (DIV .EQ. 0.0e0) then print*,'ERROR:ZERO DIVISOR AT J =',J stop end if C(J) = C(J) / DIV call SAXPY(J-1,-C(J),RG(1,J),1,C,1) 30 continue c call SVECP(C,NC,'0 Solution: C()=') print*,' ' ESTSD = abs(RG(NC1,NC1)) * sqrt(RG(NC1,NC2)) / sqrt(REAL(NXY-NC)) print*,'Estimated STD. DEV. of data errors =',ESTSD stop end