@@ -2802,7 +2802,7 @@ subroutine dlsei (w, mdw, me, ma, mg, n, prgopt, x, rnorme, &
2802
2802
integer :: i, imax, j, jp1, k, key, kranke, last, lchk, link, m, &
2803
2803
mapke1, mdeqc, mend, mep1, n1, n2, next, nlink, nopt, np1, &
2804
2804
ntimes
2805
- logical :: cov
2805
+ logical :: cov, done
2806
2806
character (len= 8 ) :: xern1, xern2, xern3, xern4
2807
2807
2808
2808
! Set the nominal tolerance used in the code for the equality
@@ -2817,11 +2817,11 @@ subroutine dlsei (w, mdw, me, ma, mg, n, prgopt, x, rnorme, &
2817
2817
write (xern3, ' (I8)' ) ma
2818
2818
write (xern4, ' (I8)' ) mg
2819
2819
write (* ,* ) ' ALL OF THE VARIABLES N, ME,' // &
2820
- ' MA, MG MUST BE >= 0$$ ENTERED ROUTINE WITH' // &
2821
- ' $$N = ' // xern1 // &
2822
- ' $$ ME = ' // xern2 // &
2823
- ' $$ MA = ' // xern3 // &
2824
- ' $$ MG = ' // xern4
2820
+ ' MA, MG MUST BE >= 0. ENTERED ROUTINE WITH: ' // &
2821
+ ' N = ' // trim ( adjustl ( xern1)) // &
2822
+ ' , ME = ' // trim ( adjustl ( xern2)) // &
2823
+ ' , MA = ' // trim ( adjustl ( xern3)) // &
2824
+ ' , MG = ' // trim ( adjustl ( xern4))
2825
2825
return
2826
2826
endif
2827
2827
@@ -2857,7 +2857,7 @@ subroutine dlsei (w, mdw, me, ma, mg, n, prgopt, x, rnorme, &
2857
2857
endif
2858
2858
2859
2859
if (mdw< m) then
2860
- write (* ,* ) ' MDW< ME+MA+MG IS AN ERROR'
2860
+ write (* ,* ) ' MDW < ME+MA+MG IS AN ERROR'
2861
2861
return
2862
2862
endif
2863
2863
@@ -2972,8 +2972,7 @@ subroutine dlsei (w, mdw, me, ma, mg, n, prgopt, x, rnorme, &
2972
2972
if (i/= imax) call dswap (np1, w(i,1 ), mdw, w(imax,1 ), mdw)
2973
2973
if (snmax> rnmax* tau** 2 ) then
2974
2974
! Eliminate elements I+1,...,N in row I.
2975
- call dh12 (1 , i, i+1 , n, w(i,1 ), mdw, ws(i), w(i+1 ,1 ), mdw, &
2976
- 1 , m- i)
2975
+ call dh12 (1 , i, i+1 , n, w(i,1 ), mdw, ws(i), w(i+1 ,1 ), mdw, 1 , m- i)
2977
2976
else
2978
2977
kranke = i - 1
2979
2978
exit
@@ -3033,6 +3032,7 @@ subroutine dlsei (w, mdw, me, ma, mg, n, prgopt, x, rnorme, &
3033
3032
x(kranke+1 ), rnorml, mode, ws(n2), ip(2 ))
3034
3033
3035
3034
! Test for consistency of equality constraints.
3035
+ done = .false.
3036
3036
3037
3037
if (me> 0 ) then
3038
3038
mdeqc = 0
@@ -3051,62 +3051,66 @@ subroutine dlsei (w, mdw, me, ma, mg, n, prgopt, x, rnorme, &
3051
3051
size = dasum(n,w(i,1 ),mdw)* xnorm + abs (w(i,np1))
3052
3052
if (w(i,np1)>tau* size) then
3053
3053
mode = mode + 2
3054
- go to 290
3054
+ done = .true.
3055
+ exit
3055
3056
endif
3056
3057
end do
3057
3058
endif
3058
3059
endif
3059
3060
3060
- ! Replace diagonal terms of lower trapezoidal matrix.
3061
+ if (.not. done) then
3062
+ ! Replace diagonal terms of lower trapezoidal matrix.
3061
3063
3062
- if (kranke> 0 ) then
3063
- call dcopy (kranke, ws(kranke+1 ), 1 , w, mdw+1 )
3064
- ! Reapply transformation to put solution in original coordinates.
3065
- do i = kranke,1 ,- 1
3066
- call dh12 (2 , i, i+1 , n, w(i,1 ), mdw, ws(i), x, 1 , 1 , 1 )
3067
- end do
3064
+ if (kranke> 0 ) then
3065
+ call dcopy (kranke, ws(kranke+1 ), 1 , w, mdw+1 )
3066
+ ! Reapply transformation to put solution in original coordinates.
3067
+ do i = kranke,1 ,- 1
3068
+ call dh12 (2 , i, i+1 , n, w(i,1 ), mdw, ws(i), x, 1 , 1 , 1 )
3069
+ end do
3068
3070
3069
- ! Compute covariance matrix of equality constrained problem.
3071
+ ! Compute covariance matrix of equality constrained problem.
3070
3072
3071
- if (cov) then
3072
- do j = min (kranke,n-1 ),1 ,- 1
3073
- rb = ws(j)* w(j,j)
3074
- if (rb/= 0.0_wp ) rb = 1.0_wp / rb
3075
- jp1 = j + 1
3076
- do i = jp1,n
3077
- w(i,j) = rb* ddot(n- j,w(i,jp1),mdw,w(j,jp1),mdw)
3078
- end do
3079
- gam = 0.5_wp * rb* ddot(n- j,w(jp1,j),1 ,w(j,jp1),mdw)
3080
- call daxpy (n- j, gam, w(j,jp1), mdw, w(jp1,j), 1 )
3081
- do i = jp1,n
3082
- do k = i,n
3083
- w(i,k) = w(i,k) + w(j,i)* w(k,j) + w(i,j)* w(j,k)
3084
- w(k,i) = w(i,k)
3073
+ if (cov) then
3074
+ do j = min (kranke,n-1 ),1 ,- 1
3075
+ rb = ws(j)* w(j,j)
3076
+ if (rb/= 0.0_wp ) rb = 1.0_wp / rb
3077
+ jp1 = j + 1
3078
+ do i = jp1,n
3079
+ w(i,j) = rb* ddot(n- j,w(i,jp1),mdw,w(j,jp1),mdw)
3085
3080
end do
3081
+ gam = 0.5_wp * rb* ddot(n- j,w(jp1,j),1 ,w(j,jp1),mdw)
3082
+ call daxpy (n- j, gam, w(j,jp1), mdw, w(jp1,j), 1 )
3083
+ do i = jp1,n
3084
+ do k = i,n
3085
+ w(i,k) = w(i,k) + w(j,i)* w(k,j) + w(i,j)* w(j,k)
3086
+ w(k,i) = w(i,k)
3087
+ end do
3088
+ end do
3089
+ uj = ws(j)
3090
+ vj = gam* uj
3091
+ w(j,j) = uj* vj + uj* vj
3092
+ do i = jp1,n
3093
+ w(j,i) = uj* w(i,j) + vj* w(j,i)
3094
+ end do
3095
+ call dcopy (n- j, w(j, jp1), mdw, w(jp1,j), 1 )
3086
3096
end do
3087
- uj = ws(j)
3088
- vj = gam* uj
3089
- w(j,j) = uj* vj + uj* vj
3090
- do i = jp1,n
3091
- w(j,i) = uj* w(i,j) + vj* w(j,i)
3092
- end do
3093
- call dcopy (n- j, w(j, jp1), mdw, w(jp1,j), 1 )
3094
- end do
3097
+ endif
3095
3098
endif
3096
- endif
3097
3099
3098
- ! Apply the scaling to the covariance matrix.
3100
+ ! Apply the scaling to the covariance matrix.
3099
3101
3100
- if (cov) then
3101
- do i = 1 ,n
3102
- call dscal (n, ws(i+ n1-1 ), w(i,1 ), mdw)
3103
- call dscal (n, ws(i+ n1-1 ), w(1 ,i), 1 )
3104
- end do
3105
- endif
3102
+ if (cov) then
3103
+ do i = 1 ,n
3104
+ call dscal (n, ws(i+ n1-1 ), w(i,1 ), mdw)
3105
+ call dscal (n, ws(i+ n1-1 ), w(1 ,i), 1 )
3106
+ end do
3107
+ endif
3108
+
3109
+ end if
3106
3110
3107
3111
! Rescale solution vector.
3108
3112
3109
- 290 if (mode<= 1 ) then
3113
+ if (mode<= 1 ) then
3110
3114
do j = 1 ,n
3111
3115
x(j) = x(j)* ws(n1+ j-1 )
3112
3116
end do
@@ -3891,7 +3895,7 @@ subroutine dwnlsm (w, mdw, mme, ma, n, l, prgopt, x, rnorm, mode, &
3891
3895
dope(2 ) = eanorm
3892
3896
dope(3 ) = tau
3893
3897
call dwnlit (w, mdw, m, n, l, ipivot, itype, h, scale, rnorm, &
3894
- idope, dope, done)
3898
+ idope, dope, done)
3895
3899
me = idope(1 )
3896
3900
krank = idope(2 )
3897
3901
niv = idope(3 )
0 commit comments