1c program DRSBLAS4 2c>> 1994-10-19 DRSBLAS4 Krogh Changes to use M77CON 3c>> 1992-05-12 DRSBLAS4 CLL Improved consistency of various versions. 4c>> 1992-04-15 CLL 5c>> 1992-03-16 CLL 6c>> 1991-12-02 DRSBLAS4 CLL 7c>> 1991-07-31 DRSBLAS4 CLL 8c 9c Demonstrates the use of BLAS subroutines SROTMG, SROTM, SAXPY, 10c and SCOPY 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--S 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 real XTAB(MXY), YTAB(MXY), W(MC2), RG(MC1,MC2) 21 real C(MC), ZERO(1), X, Y, ESTSD, DIV, PARAM(5) 22 real ONE(1) 23 integer NXY, NC, NC1, NC2, IXY, J 24c 25 data XTAB / 0.0e0, .1e0, .2e0, .3e0, .4e0, .5e0, .6e0, .7e0, 26 * .8e0, .9e0, 1.0e0 / 27 data YTAB / 1.00e0, .91e0, .86e0, .82e0, .81e0, .82e0, .85e0, 28 * .89e0, .95e0, 1.02e0, 1.10e0 / 29 data NXY, NC / 11, 3 / 30 data ONE(1), ZERO(1) / 1.0e0, 0.0e0 / 31c ------------------------------------------------------------------ 32 print'(a/)',' DRSBLAS4.. Demo for SROTMG and SROTM.' 33 NC1 = NC + 1 34 NC2 = NC1 + 1 35 call SCOPY(MC1*MC2, ZERO, 0, RG, 1) 36 call SCOPY(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.0e0 44 W(2) = X 45 W(3) = EXP(-X) 46 W(4) = Y 47 W(5) = 1.0E0 48c 49c Accumulate W() into the triangular matrix [R:G] 50c using Givens rotations. 51c 52 do 10 J = 1, NC1 53 call SROTMG(RG(J,NC2), W(NC2), RG(J,J), W(J), PARAM) 54 if(J .lt. NC1) call SROTM(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 SCOPY(NC,RG(1,NC1),1,C,1) 62 do 30 J = NC, 1, -1 63 DIV = RG(J,J) 64 if (DIV .EQ. 0.0e0) then 65 print*,'ERROR:ZERO DIVISOR AT J =',J 66 stop 67 end if 68 C(J) = C(J) / DIV 69 call SAXPY(J-1,-C(J),RG(1,J),1,C,1) 70 30 continue 71c 72 call SVECP(C,NC,'0 Solution: C()=') 73 print*,' ' 74 ESTSD = abs(RG(NC1,NC1)) * sqrt(RG(NC1,NC2)) / sqrt(REAL(NXY-NC)) 75 print*,'Estimated STD. DEV. of data errors =',ESTSD 76 stop 77 end 78