c program DRDBLAS4 c>> 1994-10-19 DRDBLAS4 Krogh Changes to use M77CON c>> 1992-05-12 DRDBLAS4 CLL Improved consistency of various versions. c>> 1992-04-15 CLL c>> 1992-03-16 CLL c>> 1991-12-02 DRDBLAS4 CLL c>> 1991-07-31 DRDBLAS4 CLL c c Demonstrates the use of BLAS subroutines DROTMG, DROTM, DAXPY, c and DCOPY 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--D 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 ) double precision XTAB(MXY), YTAB(MXY), W(MC2), RG(MC1,MC2) double precision C(MC), ZERO(1), X, Y, ESTSD, DIV, PARAM(5) double precision ONE(1) integer NXY, NC, NC1, NC2, IXY, J 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 / 11, 3 / data ONE(1), ZERO(1) / 1.0d0, 0.0d0 / c ------------------------------------------------------------------ print'(a/)',' DRDBLAS4.. Demo for DROTMG and DROTM.' NC1 = NC + 1 NC2 = NC1 + 1 call DCOPY(MC1*MC2, ZERO, 0, RG, 1) call DCOPY(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.0d0 W(2) = X W(3) = EXP(-X) W(4) = Y W(5) = 1.0D0 c c Accumulate W() into the triangular matrix [R:G] c using Givens rotations. c do 10 J = 1, NC1 call DROTMG(RG(J,NC2), W(NC2), RG(J,J), W(J), PARAM) if(J .lt. NC1) call DROTM(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 DCOPY(NC,RG(1,NC1),1,C,1) do 30 J = NC, 1, -1 DIV = RG(J,J) if (DIV .EQ. 0.0d0) then print*,'ERROR:ZERO DIVISOR AT J =',J stop end if C(J) = C(J) / DIV call DAXPY(J-1,-C(J),RG(1,J),1,C,1) 30 continue c call DVECP(C,NC,'0 Solution: C()=') print*,' ' ESTSD = abs(RG(NC1,NC1)) * sqrt(RG(NC1,NC2)) / sqrt(DBLE(NXY-NC)) print*,'Estimated STD. DEV. of data errors =',ESTSD stop end