1c program DRDBLAS2 2c>> 1996-06-18 DRDBLAS2 Krogh Minor format change for C conversion. 3c>> 1994-10-19 DRDBLAS2 Krogh Changes to use M77CON 4c>> 1991-11-27 DRDBLAS2 CLL 5c>> 1987-12-09 Lawson Initial Code. 6c 7c Demonstrates the use of BLAS subroutines DROTG, DROT, DAXPY, 8c and DCOPY to implement an algorithm for solving a linear 9c least squares problem using sequential accumulation of the 10c data and Givens orthogonal transformations. 11c YTAB() contains rounded values of -2 + 2*X + 3*Exp(-X) 12c ----------------------------------------------------------------- 13c--D replaces "?": DR?BLAS2, ?ROTG, ?ROT, ?AXPY, ?COPY 14c ----------------------------------------------------------------- 15 integer MC, MC1, MXY 16 parameter ( MC=3, MC1=MC+1, MXY=11 ) 17 integer IXY, J, NC, NC1, NXY 18 double precision X, XTAB(MXY), Y, YTAB(MXY), W(MC1) 19 double precision C, RG(MC1,MC1), S 20 double precision COEF(MC), DIV, ESTSD, ZERO(1) 21c 22 data XTAB / 0.0d0, .1d0, .2d0, .3d0, .4d0, .5d0, 23 * .6d0, .7d0, .8d0, .9d0, 1.0d0 / 24 data YTAB / 1.00d0, .91d0, .86d0, .82d0, .81d0, .82d0, 25 * .85d0, .89d0, .95d0, 1.02d0, 1.10d0 / 26 data NXY, NC / MXY, MC / 27 data ZERO(1) / 0.0d0 / 28c ----------------------------------------------------------------- 29 NC1 = NC + 1 30 call DCOPY(MC1*MC1, ZERO, 0, RG, 1) 31 do 20 IXY = 1, NXY 32 X = XTAB(IXY) 33 Y = YTAB(IXY) 34c Build new row of [A:B] in W(). 35 W(1) = 1.0d0 36 W(2) = X 37 W(3) = exp(-X) 38 W(4) = Y 39c Process W() into [R:G]. 40 do 10 J = 1, NC 41 call DROTG(RG(J,J),W(J),C,S) 42 call DROT(NC1-J,RG(J,J+1),MC1,W(J+1),1,C,S) 43 10 continue 44 call DROTG(RG(NC1,NC1),W(NC1),C,S) 45 20 continue 46c Begin: Solve triangular system. 47 call DCOPY(NC,RG(1,NC1),1,COEF,1) 48 do 30 J = NC, 1, -1 49 DIV = RG(J,J) 50 if (DIV .eq. 0.0d0) then 51 print '(''ERROR:ZERO DIVISOR AT J ='', I2)', J 52 stop 53 end if 54 COEF(J) = COEF(J) / DIV 55 call DAXPY(J-1,-COEF(J),RG(1,J),1,COEF,1) 56 30 continue 57c End: Solve triangular system. 58c 59 print'('' Solution: COEF() = '',3f8.3)',(COEF(J),J=1,NC) 60 ESTSD = abs(RG(NC1,NC1)) / sqrt(DBLE(NXY-NC)) 61 print'(/'' Estimated Std. Dev. of data errors ='',f9.5)', ESTSD 62 stop 63 end 64