diff options
Diffstat (limited to 'math/llsq/progs/prog2.f')
| -rw-r--r-- | math/llsq/progs/prog2.f | 125 |
1 files changed, 125 insertions, 0 deletions
diff --git a/math/llsq/progs/prog2.f b/math/llsq/progs/prog2.f new file mode 100644 index 00000000..5fb601cc --- /dev/null +++ b/math/llsq/progs/prog2.f @@ -0,0 +1,125 @@ +c prog2 +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c demonstrate algorithm hfti for solving least squares problems +c and algorithm cov for computing the associated unscaled +c covariance matrix. +c + dimension a(8,8),h(8),b(8),g(8) + real gen,anoise + integer ip(8) + double precision sm + data mda /8/ +c + do 180 noise=1,2 + anorm=500. + anoise=0. + tau=.5 + if (noise.eq.1) go to 10 + anoise=1.e-4 + tau=anorm*anoise*10. + 10 continue +c initialize the data generation function +c .. + dummy=gen(-1.) + write (6,230) + write (6,240) anoise,anorm,tau +c + do 180 mn1=1,6,5 + mn2=mn1+2 + do 180 m=mn1,mn2 + do 180 n=mn1,mn2 + write (6,250) m,n +c generate data +c .. + do 20 i=1,m + do 20 j=1,n + 20 a(i,j)=gen(anoise) + do 30 i=1,m + 30 b(i)=gen(anoise) +c +c ****** call hfti ****** +c + call hfti(a,mda,m,n,b,1,1,tau,krank,srsmsq,h,g,ip) +c +c + write (6,260) krank + write (6,200) (i,b(i),i=1,n) + write (6,190) srsmsq + if (krank.lt.n) go to 180 +c ****** algorithm cov bigins here ****** +c .. + do 40 j=1,n + 40 a(j,j)=1./a(j,j) + if (n.eq.1) go to 70 + nm1=n-1 + do 60 i=1,nm1 + ip1=i+1 + do 60 j=ip1,n + jm1=j-1 + sm=0.d0 + do 50 l=i,jm1 + 50 sm=sm+a(i,l)*dble(a(l,j)) + 60 a(i,j)=-sm*a(j,j) +c .. +c the upper triangle of a has been inverted +c upon itself. + 70 do 90 i=1,n + do 90 j=i,n + sm=0.d0 + do 80 l=j,n + 80 sm=sm+a(i,l)*dble(a(j,l)) + 90 a(i,j)=sm + if (n.lt.2) go to 160 + do 150 ii=2,n + i=n+1-ii + if (ip(i).eq.i) go to 150 + k=ip(i) + tmp=a(i,i) + a(i,i)=a(k,k) + a(k,k)=tmp + if (i.eq.1) go to 110 + do 100 l=2,i + tmp=a(l-1,i) + a(l-1,i)=a(l-1,k) + 100 a(l-1,k)=tmp + 110 ip1=i+1 + km1=k-1 + if (ip1.gt.km1) go to 130 + do 120 l=ip1,km1 + tmp=a(i,l) + a(i,l)=a(l,k) + 120 a(l,k)=tmp + 130 if (k.eq.n) go to 150 + kp1=k+1 + do 140 l=kp1,n + tmp=a(i,l) + a(i,l)=a(k,l) + 140 a(k,l)=tmp + 150 continue + 160 continue +c .. +c covariance has been computed and repermuted. +c the upper triangular part of the +c symmetric matrix (a**t*a)**(-1) has +c replaced the upper triangular part of +c the a array. + write (6,210) + do 170 i=1,n + 170 write (6,220) (i,j,a(i,j),j=i,n) + 180 continue + stop + 190 format (1h0,8x,17hresidual length =,e12.4) + 200 format (1h0,8x,34hestimated parameters, x=a**(+)*b,,22h computed + 1by 'hfti' //(9x,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8)) + 210 format (1h0,8x,31hcovariance matrix (unscaled) of,22h estimated pa + 1rameters.,19h computed by 'cov'./1x) + 220 format (9x,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8) + 230 format (52h1 prog2. this program demonstates the algorithms,16 + 1h hfti and cov.) + 240 format (1h0,54hthe relative noise level of the generated data will + 1 be,e16.4/33h0the matrix norm is approximately,e12.4/43h0the absol + 2ute pseudorank tolerance, tau, is,e12.4) + 250 format (1h0////9h0 m n/1x,2i4) + 260 format (1h0,8x,12hpseudorank =,i4) + end |