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