To determine the least squares solution of a large overdetermined linear system.
Given the m by n matrix r (m >= n)
and the m-vector b,
this routine calculates the n-vector x such
that the euclidean norm of the residue (r*x-b)
is minimized. the subroutine accepts rows of
the matrix one by one so that the entire matrix
need not be stored at one time. this allows
large problems to be solved without peripheral
storage. the length of the rows is limited by
the amount of scratch storage which can be set
aside for use by the routine. there is no
restriction on the number of rows.
suprls is called once for
each row of the matrix. A final call returns
the solution vector and the euclidean norm
of the residual. This following sequence would
process the m by n matrix r and the right hand
side m-vector b
do i = 1,m
do j = 1,n
! here set rowi(j) to the (i,j) element of r
end do
! here set bi to the ith component of b.
call suprls(i,rowi,n,bi,a,nn,soln,err,ier)
end do
call suprls (0,rowi,n,bi,a,nn,soln,err,ier)
given the m by n matrix r (m>=n) and the
m-vector b, we wish to find an n-vector x
such that
e = l2-norm of r*x-b
is minimized. since the euclidean norm is
invariant under orthogonal transformation,
r and b may be premultiplied by any
orthogonal matrix without changing the
norm of the residual (r*x-b). r is reduced
to upper triangular form by premultiplying
r and b by a sequence of householder and
rotation matrices. when the reduction is
complete, the norm of the residual takes the
form
e = l2 norm(t*x-b(n))+l2 norm(b(m-n))
where t is an n by n upper triangular
matrix, b(n) is a vector of the first n
components of b, b(m-n) is a vector of
the remaining (m-n) components of b. e is
minimized by taking x to be the solution
of the system t*x=b(n). this triangular
system is therefore solved to give the
required least squares solution. the norm
of the residual is then the l2-norm of b(m-n).
at each phase of the reduction, as many rows
as space permits are entered into the scratch
area. householder transformations are then
used to zero out subdiagonal elements. space
is saved by eliminating storage for the
zero subdiagonal terms. if there is room
for only one new row, rotation rather than
householder matrices are used for greater
speed. when all m rows have been entered,
reduction is completed and the triangular
system solved.
This will depend upon the size and condition of the matrix. Near machine accuracy may be expected for well conditioned systems of moderate size. If ill conditioning is suspect, a version using pivoting may be necessary.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(splpak_type), | intent(inout) | :: | me | |||
| integer, | intent(in) | :: | i |
the index of the row being entered. ( |
||
| real(kind=wp), | intent(in) | :: | rowi(n) |
a vector which on the |
||
| integer, | intent(in) | :: | n |
the length of the rows of the matrix (i.e.,
the number of columns). |
||
| real(kind=wp), | intent(in) | :: | bi |
on the |
||
| real(kind=wp), | intent(inout) | :: | a(nn) |
a working array which must not be changed between the successive calls to suprls. |
||
| integer, | intent(in) | :: | nn |
length of scratch array |
||
| real(kind=wp), | intent(out) | :: | soln(n) |
the |
||
| real(kind=wp), | intent(out) | :: | err |
the euclidean norm of the residual is
returned in |
||
| integer, | intent(out) | :: | ier |
error parameter. fatal errors:
|
subroutine suprls(me,i,rowi,n,bi,a,nn,soln,err,ier) class(splpak_type),intent(inout) :: me integer,intent(in) :: i !! the index of the row being entered. (`i` is 1 !! for the first call, increases by 1 for each !! call, and is `m` when the final row is !! entered). After the final row has been !! entered, [[suprls]] is called with `i = 0` to !! complete the reduction and solution. integer,intent(in) :: nn !! length of scratch array `a`. `nn` must be at !! least `n*(n+5)/2+1`. For speed, `nn` should be !! as large as possible up to a maximum of !! `(n+1)*m`. integer,intent(in) :: n !! the length of the rows of the matrix (i.e., !! the number of columns). `n <= m`, where `m` is !! the number of rows. real(wp),intent(in) :: rowi(n) !! a vector which on the `i`th call contains the `n` !! components of the `i`th row of the matrix. the !! dimension of `rowi` in calling program must be !! at least `n`. real(wp),intent(in) :: bi !! on the `i`th call, `bi` contains the `i`th element !! of the right hand side vector `b`. real(wp),intent(inout) :: a(nn) !! a working array which must not be changed !! between the successive calls to [[suprls]]. real(wp),intent(out) :: soln(n) !! the `n`-components of the solution vector are !! returned in this array after the final call !! to [[suprls]]. real(wp),intent(out) :: err !! the euclidean norm of the residual is !! returned in `err` after the final call to !! [[suprls]]. integer,intent(out) :: ier !! error parameter. !! fatal errors: !! !! * 32 -- insufficient scratch storage provided, !! must have `nn >= n*(n+5)/2+1`. !! * 33 -- array has too few rows. must have !! `m >= n`. !! * 34 -- system is singular. !! * 35 -- values of `i` not in sequence. real(wp) :: s,temp,temp1,cn,sn integer :: j,ilj,ilnp,nreq,k,& idiag,i1,i2,ii,jp1,lmkm1,j1,jdel,idj,iijd,& i1jd,k11,k1m1,i11,np1mk,lmk,imov,iii,iiim,iim1,& ilk,npk,ilii,npii logical :: complete_reduction !! Routine entered with `I<=0` means complete !! the reduction and store the solution in `SOLN`. real(wp),parameter :: tol = 1.0e-18_wp !! small number tolerance ier = 0 complete_reduction = i <= 0 if (.not. complete_reduction) then if (i<=1) then ! Set up quantities on first call. me%iold = 0 me%np1 = n + 1 ! Compute how many rows can be input now. me%l = nn/me%np1 me%ilast = 0 me%il1 = 0 me%k = 0 me%k1 = 0 me%errsum = 0.0_wp nreq = ((n+5)*n+2)/2 ! Error exit if insufficient scratch storage provided. if (nn<nreq) then ier = 32 write(*,*) 'nn = ', nn write(*,*) 'nreq = ', nreq call cfaerr(ier, & ' suprls - insufficient scratch storage provided. '//& 'at least ((N+5)*N+2)/2 locations needed') return end if end if ! Error exit if (I-IOLD)/=1. if ((i-me%iold)/=1) then ier = 35 write(*,*) 'i =',i write(*,*) 'me%iold =',me%iold call cfaerr(ier,' suprls - values of I not in sequence') return end if ! Store the row in the scratch storage. me%iold = i do j = 1,n ilj = me%ilast + j a(ilj) = rowi(j) end do ilnp = me%ilast + me%np1 a(ilnp) = bi me%ilast = me%ilast + me%np1 me%isav = i if (i<me%l) return end if main : do if (.not. complete_reduction) then if (me%k/=0) then me%k1 = min(me%k,n) idiag = -me%np1 if (me%l-me%k==1) then ! Apply rotations to zero out the single new row. do j = 1,me%k1 idiag = idiag + (me%np1-j+2) i1 = me%il1 + j if (abs(a(i1))<=tol) then s = sqrt(a(idiag)*a(idiag)) else if (abs(a(idiag))<tol) then s = sqrt(a(i1)*a(i1)) else s = sqrt(a(idiag)*a(idiag)+a(i1)*a(i1)) end if if (s==0.0_wp) cycle temp = a(idiag) a(idiag) = s s = 1.0_wp/s cn = temp*s sn = a(i1)*s jp1 = j + 1 do j1 = jp1,me%np1 jdel = j1 - j idj = idiag + jdel temp = a(idj) i1jd = i1 + jdel a(idj) = cn*temp + sn*a(i1jd) a(i1jd) = -sn*temp + cn*a(i1jd) end do end do else ! Apply householder transformations to zero out new rows. do j = 1,me%k1 idiag = idiag + (me%np1-j+2) i1 = me%il1 + j i2 = i1 + me%np1* (me%l-me%k-1) s = a(idiag)*a(idiag) do ii = i1,i2,me%np1 s = s + a(ii)*a(ii) end do if (s==0.0_wp) cycle temp = a(idiag) a(idiag) = sqrt(s) if (temp>0.0_wp) a(idiag) = -a(idiag) temp = temp - a(idiag) temp1 = 1.0_wp / (temp*a(idiag)) jp1 = j + 1 do j1 = jp1,me%np1 jdel = j1 - j idj = idiag + jdel s = temp*a(idj) do ii = i1,i2,me%np1 iijd = ii + jdel s = s + a(ii)*a(iijd) end do s = s*temp1 a(idj) = a(idj) + s*temp do ii = i1,i2,me%np1 iijd = ii + jdel a(iijd) = a(iijd) + s*a(ii) end do end do end do end if if (me%k>=n) then lmkm1 = me%l - me%k ! Accumulate residual sum of squares. do ii = 1,lmkm1 ilnp = me%il1 + ii*me%np1 me%errsum = me%errsum + a(ilnp)*a(ilnp) end do if (i<=0) exit main me%k = me%l me%ilast = me%il1 ! Determine how many new rows may be input on next iteration. me%l = me%k + (nn-me%ilast)/me%np1 return end if end if k11 = me%k1 + 1 me%k1 = min(me%l,n) if (me%l-me%k/=1) then k1m1 = me%k1 - 1 if (me%l>n) k1m1 = n i1 = me%il1 + k11 - me%np1 - 1 ! Perform householder transformations to reduce rows to upper ! triangular form. do j = k11,k1m1 i1 = i1 + (me%np1+1) i2 = i1 + (me%l-j)*me%np1 s = 0.0_wp do ii = i1,i2,me%np1 s = s + a(ii)*a(ii) end do if (s==0.0_wp) cycle temp = a(i1) a(i1) = sqrt(s) if (temp>0.0_wp) a(i1) = -a(i1) temp = temp - a(i1) temp1 = 1.0_wp/ (temp*a(i1)) jp1 = j + 1 i11 = i1 + me%np1 do j1 = jp1,me%np1 jdel = j1 - j i1jd = i1 + jdel s = temp*a(i1jd) do ii = i11,i2,me%np1 iijd = ii + jdel s = s + a(ii)*a(iijd) end do s = s*temp1 i1jd = i1 + jdel a(i1jd) = a(i1jd) + s*temp do ii = i11,i2,me%np1 iijd = ii + jdel a(iijd) = a(iijd) + s*a(ii) end do end do end do if (me%l>n) then np1mk = me%np1 - me%k lmk = me%l - me%k ! Accumulate residual sum of squares. do ii = np1mk,lmk ilnp = me%il1 + ii*me%np1 me%errsum = me%errsum + a(ilnp)*a(ilnp) end do end if end if imov = 0 i1 = me%il1 + k11 - me%np1 - 1 ! Squeeze the unnecessary elements out of scratch storage to ! allow space for more rows. do ii = k11,me%k1 imov = imov + (ii-1) i1 = i1 + me%np1 + 1 i2 = i1 + me%np1 - ii do iii = i1,i2 iiim = iii - imov a(iiim) = a(iii) end do end do me%ilast = i2 - imov me%il1 = me%ilast if (i<=0) exit main me%k = me%l ! Determine how many new rows may be input on next iteration. me%l = me%k + (nn-me%ilast)/me%np1 return end if ! Complete reduction and store solution in SOLN. complete_reduction = .false. ! reset this flag me%l = me%isav ! Error exit if L less than N. if (me%l<n) then ier = 33 call cfaerr(ier,' suprls - array has too few rows.') return end if ! K/=ISAV means further reduction needed. if (me%k==me%isav) exit main end do main me%ilast = (me%np1* (me%np1+1))/2 - 1 if (a(me%ilast-1)==0.0_wp) then ! Error return if system is singular. ier = 34 call cfaerr(ier,' suprls - system is singular.') return end if ! Solve triangular system into ROWI. soln(n) = a(me%ilast)/a(me%ilast-1) do ii = 2,n iim1 = ii - 1 me%ilast = me%ilast - ii s = a(me%ilast) do k = 1,iim1 ilk = me%ilast - k npk = me%np1 - k s = s - a(ilk)*soln(npk) end do me%k = k ! JW : is this necessary ??? ilii = me%ilast - ii if (a(ilii)==0.0_wp) then ! Error return if system is singular. ier = 34 call cfaerr(ier,' suprls - system is singular.') return end if npii = me%np1 - ii soln(npii) = s/a(ilii) end do ! Store residual norm. err = sqrt(me%errsum) end subroutine suprls