| 1 |
aw0a |
1 |
************************************************************************ |
| 2 |
|
|
* |
| 3 |
|
|
* File mi65rmod fortran. |
| 4 |
|
|
* |
| 5 |
|
|
* m6bfgs m6bswp m6radd m6rcnd m6rdel |
| 6 |
|
|
* m6rmod m6rset m6rsol m6swap |
| 7 |
|
|
* |
| 8 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 9 |
|
|
|
| 10 |
|
|
subroutine m6bfgs( maxr, n, nr, r, g, g2, p, v, |
| 11 |
|
|
$ step, told, tolz, inform ) |
| 12 |
|
|
|
| 13 |
|
|
implicit double precision (a-h,o-z) |
| 14 |
|
|
double precision r(nr), g(n), g2(n), p(n), v(n) |
| 15 |
|
|
|
| 16 |
|
|
* ------------------------------------------------------------------ |
| 17 |
|
|
* m6bfgs applies the BFGS update to the upper-triangular matrix r, |
| 18 |
|
|
* which holds the cholesky factor of the quasi-newton approximation |
| 19 |
|
|
* of the (projected) hessian. |
| 20 |
|
|
* |
| 21 |
|
|
* r contains a triangle of size ncolr = min( n, maxr ). |
| 22 |
|
|
* If n .gt. maxr, r also contains a diagonal of size n - maxr. |
| 23 |
|
|
* |
| 24 |
|
|
* p holds the search direction. It is overwritten. |
| 25 |
|
|
* v must satisfy r(t)*v = g. It is overwritten. |
| 26 |
|
|
* |
| 27 |
|
|
* On exit, |
| 28 |
|
|
* inform = 0 if no update was performed, |
| 29 |
|
|
* = 1 if the update was successful, |
| 30 |
|
|
* = 2 if it was nearly singular. |
| 31 |
|
|
* ------------------------------------------------------------------ |
| 32 |
|
|
|
| 33 |
|
|
intrinsic abs, min, sqrt |
| 34 |
|
|
parameter ( one = 1.0d+0 ) |
| 35 |
|
|
|
| 36 |
|
|
inform = 0 |
| 37 |
|
|
ncolr = min( n, maxr ) |
| 38 |
|
|
gtp = ddot ( n, g , 1, p, 1 ) |
| 39 |
|
|
gtp2 = ddot ( n, g2, 1, p, 1 ) |
| 40 |
|
|
if (gtp2 .le. 0.91*gtp) return |
| 41 |
|
|
|
| 42 |
|
|
delta1 = one / sqrt( abs( gtp ) ) |
| 43 |
|
|
delta2 = one / sqrt( step*(gtp2 - gtp) ) |
| 44 |
|
|
|
| 45 |
|
|
* Normalize v and change its sign. |
| 46 |
|
|
|
| 47 |
|
|
call dscal ( n, (- delta1), v, 1 ) |
| 48 |
|
|
|
| 49 |
|
|
* Avoid cancellation error in forming the new vector p. |
| 50 |
|
|
|
| 51 |
|
|
if ( abs( delta1/delta2 - one ) .ge. 0.5d+0) then |
| 52 |
|
|
do 100 j = 1, n |
| 53 |
|
|
p(j) = delta2*( g2(j) - g(j) ) + delta1*g(j) |
| 54 |
|
|
100 continue |
| 55 |
|
|
else |
| 56 |
|
|
d = delta1 - delta2 |
| 57 |
|
|
do 120 j = 1, n |
| 58 |
|
|
p(j) = delta2*g2(j) + d*g(j) |
| 59 |
|
|
120 continue |
| 60 |
|
|
end if |
| 61 |
|
|
|
| 62 |
|
|
* Find the last nonzero in v. |
| 63 |
|
|
|
| 64 |
|
|
do 450 ii = 1, ncolr |
| 65 |
|
|
lastv = ncolr + 1 - ii |
| 66 |
|
|
if (abs( v(lastv) ) .gt. tolz) go to 500 |
| 67 |
|
|
450 continue |
| 68 |
|
|
lastv = 1 |
| 69 |
|
|
|
| 70 |
|
|
* Triangularize r + v * p(t). |
| 71 |
|
|
|
| 72 |
|
|
500 call m6rmod( ncolr, nr, r, v, p, lastv, told, tolz, inform ) |
| 73 |
|
|
|
| 74 |
|
|
* Deal with surplus diagonals of r. |
| 75 |
|
|
|
| 76 |
|
|
if (n .gt. maxr) then |
| 77 |
|
|
j1 = maxr + 1 |
| 78 |
|
|
l = maxr*j1/2 |
| 79 |
|
|
do 650 j = j1, n |
| 80 |
|
|
l = l + 1 |
| 81 |
|
|
r(l) = r(l) + v(j)*p(j) |
| 82 |
|
|
650 continue |
| 83 |
|
|
end if |
| 84 |
|
|
|
| 85 |
|
|
* end of m6bfgs |
| 86 |
|
|
end |
| 87 |
|
|
|
| 88 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 89 |
|
|
|
| 90 |
|
|
subroutine m6bswp( n, nr, r, v, w, lastv, told, tolz, inform ) |
| 91 |
|
|
|
| 92 |
|
|
implicit double precision (a-h,o-z) |
| 93 |
|
|
double precision r(nr), v(n), w(n) |
| 94 |
|
|
|
| 95 |
|
|
* ------------------------------------------------------------------ |
| 96 |
|
|
* m6bswp modifies the upper-triangular matrix r |
| 97 |
|
|
* to account for a basis exchange in which the lastv-th |
| 98 |
|
|
* superbasic variable becomes basic. r is changed to |
| 99 |
|
|
* r + v*w(t), which is triangularized by r1mod, |
| 100 |
|
|
* where v is the lastv-th column of r, and w is input. |
| 101 |
|
|
* ------------------------------------------------------------------ |
| 102 |
|
|
|
| 103 |
|
|
* Set v = lastv-th column of r and find its norm. |
| 104 |
|
|
|
| 105 |
|
|
l = (lastv - 1)*lastv/2 |
| 106 |
|
|
call dcopy ( lastv, r(l+1), 1, v, 1 ) |
| 107 |
|
|
vnorm = dasum ( lastv, v, 1 ) |
| 108 |
|
|
call m6rmod( n, nr, r, v, w, lastv, |
| 109 |
|
|
$ (vnorm*told), (vnorm*tolz), inform ) |
| 110 |
|
|
|
| 111 |
|
|
* end of m6bswp |
| 112 |
|
|
end |
| 113 |
|
|
|
| 114 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 115 |
|
|
|
| 116 |
|
|
subroutine m6radd( maxr, nr, ns, r ) |
| 117 |
|
|
|
| 118 |
|
|
implicit double precision (a-h,o-z) |
| 119 |
|
|
double precision r(nr) |
| 120 |
|
|
|
| 121 |
|
|
* ------------------------------------------------------------------ |
| 122 |
|
|
* m6radd adds column ns to the upper triangular matrix r. |
| 123 |
|
|
* In this version (Sep 1984) it is just a unit vector. |
| 124 |
|
|
* |
| 125 |
|
|
* Modified Jan 1983 to add a diagonal only, if ns .gt. maxr. |
| 126 |
|
|
* ------------------------------------------------------------------ |
| 127 |
|
|
|
| 128 |
|
|
parameter ( zero = 0.0d+0, one = 1.0d+0 ) |
| 129 |
|
|
|
| 130 |
|
|
if (ns .eq. 1) then |
| 131 |
|
|
r(1) = one |
| 132 |
|
|
else if (ns .le. maxr) then |
| 133 |
|
|
ldiag1 = (ns - 1)*ns/2 |
| 134 |
|
|
ldiag2 = ldiag1 + ns |
| 135 |
|
|
call dload ( ns, zero, r(ldiag1+1), 1 ) |
| 136 |
|
|
r(ldiag2) = one |
| 137 |
|
|
else |
| 138 |
|
|
ldiag2 = maxr*(maxr + 1)/2 + (ns - maxr) |
| 139 |
|
|
r(ldiag2) = one |
| 140 |
|
|
end if |
| 141 |
|
|
|
| 142 |
|
|
* end of m6radd |
| 143 |
|
|
end |
| 144 |
|
|
|
| 145 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 146 |
|
|
|
| 147 |
|
|
subroutine m6rcnd( maxr, nr, ns, r, dmax, dmin, cond ) |
| 148 |
|
|
|
| 149 |
|
|
implicit double precision (a-h,o-z) |
| 150 |
|
|
double precision r(nr) |
| 151 |
|
|
|
| 152 |
|
|
* ------------------------------------------------------------------ |
| 153 |
|
|
* m6rcnd finds the largest and smallest diagonals of the |
| 154 |
|
|
* upper triangular matrix r, and returns the square of their ratio. |
| 155 |
|
|
* This is a lower bound on the condition number of r(t)*r. |
| 156 |
|
|
* ------------------------------------------------------------------ |
| 157 |
|
|
|
| 158 |
|
|
intrinsic abs, max, min |
| 159 |
|
|
|
| 160 |
|
|
ncolr = min( maxr, ns ) |
| 161 |
|
|
dmax = abs( r(1) ) |
| 162 |
|
|
dmin = dmax |
| 163 |
|
|
|
| 164 |
|
|
if (ncolr .gt. 1) then |
| 165 |
|
|
l = 1 |
| 166 |
|
|
do 100 j = 2, ncolr |
| 167 |
|
|
l = l + j |
| 168 |
|
|
d = abs( r(l) ) |
| 169 |
|
|
dmax = max( dmax, d ) |
| 170 |
|
|
dmin = min( dmin, d ) |
| 171 |
|
|
100 continue |
| 172 |
|
|
end if |
| 173 |
|
|
|
| 174 |
|
|
cond = (dmax/dmin)**2 |
| 175 |
|
|
|
| 176 |
|
|
* end of m6rcnd |
| 177 |
|
|
end |
| 178 |
|
|
|
| 179 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 180 |
|
|
|
| 181 |
|
|
subroutine m6rdel( m, maxr, nr, ns, ms, |
| 182 |
|
|
$ kb, bbl, bbu, grd, r, rg, x, jq, rset ) |
| 183 |
|
|
|
| 184 |
|
|
implicit double precision (a-h,o-z) |
| 185 |
|
|
integer kb(ms) |
| 186 |
|
|
double precision bbl(ms), bbu(ms), grd(ms) |
| 187 |
|
|
double precision r(nr), rg(ns), x(ms) |
| 188 |
|
|
logical rset |
| 189 |
|
|
|
| 190 |
|
|
* ------------------------------------------------------------------ |
| 191 |
|
|
* m6rdel deletes the jq-th superbasic variable from r |
| 192 |
|
|
* and from various arrays kb, bbl, bbu, grd, rg, x. |
| 193 |
|
|
* ------------------------------------------------------------------ |
| 194 |
|
|
|
| 195 |
|
|
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 196 |
|
|
|
| 197 |
|
|
intrinsic abs, max, min, sqrt |
| 198 |
|
|
parameter ( zero = 0.0d+0 ) |
| 199 |
|
|
|
| 200 |
|
|
if (jq .eq. ns) return |
| 201 |
|
|
if (.not. rset) go to 500 |
| 202 |
|
|
|
| 203 |
|
|
* Delete the jq-th column of r and triangularize the |
| 204 |
|
|
* remaining columns using a partial forward sweep of rotations. |
| 205 |
|
|
|
| 206 |
|
|
ncolr = min( maxr, ns ) |
| 207 |
|
|
k = jq*(jq + 1)/2 |
| 208 |
|
|
|
| 209 |
|
|
do 50 i = jq + 1, ncolr |
| 210 |
|
|
k = k + i |
| 211 |
|
|
t1 = r(k-1) |
| 212 |
|
|
t2 = r(k) |
| 213 |
|
|
d = sqrt( t1*t1 + t2*t2 ) |
| 214 |
|
|
r(k-1) = d |
| 215 |
|
|
if (i .eq. ncolr) go to 20 |
| 216 |
|
|
sn = t2/d |
| 217 |
|
|
cs = t1/d |
| 218 |
|
|
if (abs( cs ) .le. eps0) cs = zero |
| 219 |
|
|
j1 = i + 1 |
| 220 |
|
|
k1 = k + i |
| 221 |
|
|
do 10 j = j1, ncolr |
| 222 |
|
|
t1 = r(k1-1) |
| 223 |
|
|
t2 = r(k1) |
| 224 |
|
|
r(k1-1) = cs*t1 + sn*t2 |
| 225 |
|
|
r(k1) = sn*t1 - cs*t2 |
| 226 |
|
|
k1 = k1 + j |
| 227 |
|
|
10 continue |
| 228 |
|
|
20 k1 = i - 1 |
| 229 |
|
|
j2 = k - i |
| 230 |
|
|
j1 = j2 - i + 2 |
| 231 |
|
|
do 30 j = j1, j2 |
| 232 |
|
|
r(j) = r(j+k1) |
| 233 |
|
|
30 continue |
| 234 |
|
|
50 continue |
| 235 |
|
|
|
| 236 |
|
|
* If necessary, clear out the last column of r. |
| 237 |
|
|
|
| 238 |
|
|
if (ns .gt. maxr) then |
| 239 |
|
|
lastr = maxr*(maxr + 1)/2 |
| 240 |
|
|
l = lastr - maxr + 1 |
| 241 |
|
|
if (jq .le. maxr) call dload ( maxr, zero, r(l), 1 ) |
| 242 |
|
|
|
| 243 |
|
|
* Shift surplus diagonals of r to the left. |
| 244 |
|
|
|
| 245 |
|
|
j1 = max( maxr, jq ) |
| 246 |
|
|
j2 = ns - 1 |
| 247 |
|
|
l = lastr + (j1 - maxr) |
| 248 |
|
|
do 320 j = j1, j2 |
| 249 |
|
|
r(l) = r(l+1) |
| 250 |
|
|
320 continue |
| 251 |
|
|
end if |
| 252 |
|
|
|
| 253 |
|
|
* Shift all the arrays one place to the left. |
| 254 |
|
|
|
| 255 |
|
|
500 j2 = ns - 1 |
| 256 |
|
|
do 550 j = jq, j2 |
| 257 |
|
|
rg(j) = rg(j+1) |
| 258 |
|
|
k = m + j |
| 259 |
|
|
kb(k) = kb(k+1) |
| 260 |
|
|
bbl(k) = bbl(k+1) |
| 261 |
|
|
bbu(k) = bbu(k+1) |
| 262 |
|
|
grd(k) = grd(k+1) |
| 263 |
|
|
x(k) = x(k+1) |
| 264 |
|
|
550 continue |
| 265 |
|
|
|
| 266 |
|
|
* end of m6rdel |
| 267 |
|
|
end |
| 268 |
|
|
|
| 269 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 270 |
|
|
|
| 271 |
|
|
subroutine m6rmod( n, nr, r, v, w, lastv, told, tolz, inform ) |
| 272 |
|
|
|
| 273 |
|
|
implicit double precision (a-h,o-z) |
| 274 |
|
|
double precision r(nr), v(n), w(n) |
| 275 |
|
|
|
| 276 |
|
|
* ------------------------------------------------------------------ |
| 277 |
|
|
* m6rmod modifies the upper-triangular matrix r so that |
| 278 |
|
|
* q*(r + v*w(t)) is upper triangular, where q is orthogonal, |
| 279 |
|
|
* v and w are vectors, and the new r overwrites the old. |
| 280 |
|
|
* |
| 281 |
|
|
* q is the product of two sweeps of plane rotations (not stored). |
| 282 |
|
|
* These affect the lastv-th row of r, which is temporarily held |
| 283 |
|
|
* in the vector v. Thus, v is overwritten. w is not altered. |
| 284 |
|
|
* |
| 285 |
|
|
* lastv points to the last nonzero of v. The value lastv = n |
| 286 |
|
|
* would always be ok, but sometimes it is known to be less |
| 287 |
|
|
* than n, so q reduces to two partial sweeps of |
| 288 |
|
|
* rotations. |
| 289 |
|
|
* |
| 290 |
|
|
* told is a tolerance on the lastv-th diagonal of r. |
| 291 |
|
|
* tolz is a tolerance for negligible elements in v. |
| 292 |
|
|
* |
| 293 |
|
|
* On exit, |
| 294 |
|
|
* inform = 1 if the diagonal of r is larger than told, |
| 295 |
|
|
* = 2 if not (in which case it is reset to told). |
| 296 |
|
|
* ------------------------------------------------------------------ |
| 297 |
|
|
|
| 298 |
|
|
intrinsic sqrt |
| 299 |
|
|
parameter ( zero = 0.0d+0 ) |
| 300 |
|
|
|
| 301 |
|
|
vlast = v(lastv) |
| 302 |
|
|
lm1 = lastv - 1 |
| 303 |
|
|
lastr = lastv*(lastv + 1)/2 |
| 304 |
|
|
|
| 305 |
|
|
* Copy the lastv-th row of r into the end of v. |
| 306 |
|
|
|
| 307 |
|
|
l = lastr |
| 308 |
|
|
do 100 j = lastv, n |
| 309 |
|
|
v(j) = r(l) |
| 310 |
|
|
l = l + j |
| 311 |
|
|
100 continue |
| 312 |
|
|
|
| 313 |
|
|
* Reduce v to a multiple of e(lastv) |
| 314 |
|
|
* using a partial backward sweep of rotations. |
| 315 |
|
|
* This fills in the lastv-th row of r (held in v). |
| 316 |
|
|
|
| 317 |
|
|
v2 = vlast**2 |
| 318 |
|
|
ldiag = lastr |
| 319 |
|
|
|
| 320 |
|
|
do 300 ii = 1, lm1 |
| 321 |
|
|
i = lastv - ii |
| 322 |
|
|
ldiag = ldiag - i - 1 |
| 323 |
|
|
s = v(i) |
| 324 |
|
|
v(i) = zero |
| 325 |
|
|
if (abs(s) .le. tolz) go to 300 |
| 326 |
|
|
v2 = s**2 + v2 |
| 327 |
|
|
root = sqrt(v2) |
| 328 |
|
|
cs = vlast/root |
| 329 |
|
|
sn = s/root |
| 330 |
|
|
vlast = root |
| 331 |
|
|
l = ldiag |
| 332 |
|
|
|
| 333 |
|
|
do 200 j = i, n |
| 334 |
|
|
s = v(j) |
| 335 |
|
|
t = r(l) |
| 336 |
|
|
v(j) = cs*s + sn*t |
| 337 |
|
|
r(l) = sn*s - cs*t |
| 338 |
|
|
l = l + j |
| 339 |
|
|
200 continue |
| 340 |
|
|
300 continue |
| 341 |
|
|
|
| 342 |
|
|
* Set v = v + vlast*w. |
| 343 |
|
|
|
| 344 |
|
|
call daxpy ( n, vlast, w, 1, v, 1 ) |
| 345 |
|
|
|
| 346 |
|
|
* Eliminate the front of the lastv-th row of r (held in v) |
| 347 |
|
|
* using a partial forward sweep of rotations. |
| 348 |
|
|
|
| 349 |
|
|
ldiag = 0 |
| 350 |
|
|
|
| 351 |
|
|
do 700 i = 1, lm1 |
| 352 |
|
|
ldiag = ldiag + i |
| 353 |
|
|
t = v(i) |
| 354 |
|
|
if (abs(t) .le. tolz) go to 700 |
| 355 |
|
|
s = r(ldiag) |
| 356 |
|
|
root = sqrt(s**2 + t**2) |
| 357 |
|
|
cs = s/root |
| 358 |
|
|
sn = t/root |
| 359 |
|
|
r(ldiag) = root |
| 360 |
|
|
l = ldiag + i |
| 361 |
|
|
|
| 362 |
|
|
do 600 j = i + 1, n |
| 363 |
|
|
s = r(l) |
| 364 |
|
|
t = v(j) |
| 365 |
|
|
r(l) = cs*s + sn*t |
| 366 |
|
|
v(j) = sn*s - cs*t |
| 367 |
|
|
l = l + j |
| 368 |
|
|
600 continue |
| 369 |
|
|
700 continue |
| 370 |
|
|
|
| 371 |
|
|
* Insert the new lastv-th row of r. |
| 372 |
|
|
|
| 373 |
|
|
l = lastr |
| 374 |
|
|
do 900 j = lastv, n |
| 375 |
|
|
r(l) = v(j) |
| 376 |
|
|
l = l + j |
| 377 |
|
|
900 continue |
| 378 |
|
|
|
| 379 |
|
|
* Test for (unlikely) singularity. |
| 380 |
|
|
|
| 381 |
|
|
inform = 1 |
| 382 |
|
|
if (abs( r(lastr) ) .le. told) then |
| 383 |
|
|
inform = 2 |
| 384 |
|
|
r(lastr) = told |
| 385 |
|
|
end if |
| 386 |
|
|
|
| 387 |
|
|
* end of m6rmod |
| 388 |
|
|
end |
| 389 |
|
|
|
| 390 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 391 |
|
|
|
| 392 |
|
|
subroutine m6rset( maxr, nr, ns, r, cond ) |
| 393 |
|
|
|
| 394 |
|
|
implicit double precision (a-h,o-z) |
| 395 |
|
|
double precision r(nr) |
| 396 |
|
|
|
| 397 |
|
|
* ------------------------------------------------------------------ |
| 398 |
|
|
* m6rset alters r, the upper-triangular factor of the |
| 399 |
|
|
* approximation to the reduced hessian. |
| 400 |
|
|
* |
| 401 |
|
|
* If r(1) = zero, r does not exist. |
| 402 |
|
|
* In this case, r is initialized to the identity matrix. |
| 403 |
|
|
* |
| 404 |
|
|
* Otherwise, r already exists and we attempt to make it better |
| 405 |
|
|
* conditioned by scaling its columns by the square roots of its |
| 406 |
|
|
* current diagonals. |
| 407 |
|
|
* ------------------------------------------------------------------ |
| 408 |
|
|
|
| 409 |
|
|
common /m1file/ iread,iprint,isumm |
| 410 |
|
|
|
| 411 |
|
|
intrinsic abs, max, min, sqrt |
| 412 |
|
|
parameter ( zero = 0.0d+0, one = 1.0d+0 ) |
| 413 |
|
|
|
| 414 |
|
|
cond = one |
| 415 |
|
|
ncolr = min( maxr, ns ) |
| 416 |
|
|
if (ncolr .eq. 0 ) return |
| 417 |
|
|
|
| 418 |
|
|
dmax = abs( r(1) ) |
| 419 |
|
|
if (dmax .eq. zero) then |
| 420 |
|
|
|
| 421 |
|
|
* Set r = the identity. |
| 422 |
|
|
|
| 423 |
|
|
l = ncolr*(ncolr + 1)/2 |
| 424 |
|
|
call dload ( l, zero, r, 1 ) |
| 425 |
|
|
ldiag = 0 |
| 426 |
|
|
|
| 427 |
|
|
do 100 k = 1, ncolr |
| 428 |
|
|
ldiag = ldiag + k |
| 429 |
|
|
r(ldiag) = one |
| 430 |
|
|
100 continue |
| 431 |
|
|
|
| 432 |
|
|
do 120 k = maxr + 1, ns |
| 433 |
|
|
ldiag = ldiag + 1 |
| 434 |
|
|
r(ldiag) = one |
| 435 |
|
|
120 continue |
| 436 |
|
|
else |
| 437 |
|
|
|
| 438 |
|
|
* Scale the columns of r. |
| 439 |
|
|
|
| 440 |
|
|
dmin = dmax |
| 441 |
|
|
ldiag = 0 |
| 442 |
|
|
|
| 443 |
|
|
do 250 k = 1, ncolr |
| 444 |
|
|
l = ldiag + 1 |
| 445 |
|
|
ldiag = ldiag + k |
| 446 |
|
|
diag = abs( r(ldiag) ) |
| 447 |
|
|
dmax = max( dmax, diag ) |
| 448 |
|
|
dmin = min( dmin, diag ) |
| 449 |
|
|
s = one / sqrt( diag ) |
| 450 |
|
|
do 240 i = l, ldiag |
| 451 |
|
|
r(i) = r(i)*s |
| 452 |
|
|
240 continue |
| 453 |
|
|
250 continue |
| 454 |
|
|
|
| 455 |
|
|
cond = dmax / dmin |
| 456 |
|
|
if (iprint .gt. 0) write(iprint, 2000) dmax, dmin, cond |
| 457 |
|
|
end if |
| 458 |
|
|
return |
| 459 |
|
|
|
| 460 |
|
|
2000 format(/ ' XXX Hessian modified. Diags =', 1p, 2e11.2, |
| 461 |
|
|
$ ', Cond R(t)*R =', e11.2) |
| 462 |
|
|
|
| 463 |
|
|
* end of m6rset |
| 464 |
|
|
end |
| 465 |
|
|
|
| 466 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 467 |
|
|
|
| 468 |
|
|
subroutine m6rsol( mode, maxr, nr, ns, r, y ) |
| 469 |
|
|
|
| 470 |
|
|
implicit double precision (a-h,o-z) |
| 471 |
|
|
double precision r(nr), y(ns) |
| 472 |
|
|
|
| 473 |
|
|
* ------------------------------------------------------------------ |
| 474 |
|
|
* m6rsol solves R*x = y or R(t)*x = y, where R is an |
| 475 |
|
|
* upper-triangular matrix stored by columns in r(nr). |
| 476 |
|
|
* The solution x overwrites y. |
| 477 |
|
|
* ------------------------------------------------------------------ |
| 478 |
|
|
|
| 479 |
|
|
intrinsic min |
| 480 |
|
|
|
| 481 |
|
|
ncolr = min( maxr, ns ) |
| 482 |
|
|
|
| 483 |
|
|
if (mode .eq. 1) then |
| 484 |
|
|
|
| 485 |
|
|
* mode = 1 -- solve R*y = y. |
| 486 |
|
|
|
| 487 |
|
|
l = ncolr*(ncolr + 1)/2 |
| 488 |
|
|
do 120 i = ncolr, 2, -1 |
| 489 |
|
|
t = y(i) / r(l) |
| 490 |
|
|
y(i) = t |
| 491 |
|
|
l = l - i |
| 492 |
|
|
call daxpy ( i-1, (-t), r(l+1), 1, y, 1 ) |
| 493 |
|
|
120 continue |
| 494 |
|
|
y(1) = y(1) / r(1) |
| 495 |
|
|
else |
| 496 |
|
|
|
| 497 |
|
|
* mode = 2 -- solve R(t)*y = y. |
| 498 |
|
|
|
| 499 |
|
|
y(1) = y(1) / r(1) |
| 500 |
|
|
l = 1 |
| 501 |
|
|
do 220 i = 2, ncolr |
| 502 |
|
|
t = ddot ( i-1, r(l+1), 1, y, 1 ) |
| 503 |
|
|
l = l + i |
| 504 |
|
|
y(i) = (y(i) - t) / r(l) |
| 505 |
|
|
220 continue |
| 506 |
|
|
end if |
| 507 |
|
|
|
| 508 |
|
|
* Deal with surplus diagonals of r. |
| 509 |
|
|
|
| 510 |
|
|
if (ns .gt. maxr) then |
| 511 |
|
|
l = maxr*i/2 |
| 512 |
|
|
do 320 j = maxr + 1, ns |
| 513 |
|
|
l = l + 1 |
| 514 |
|
|
y(j) = y(j) / r(l) |
| 515 |
|
|
320 continue |
| 516 |
|
|
end if |
| 517 |
|
|
|
| 518 |
|
|
* end of m6rsol |
| 519 |
|
|
end |
| 520 |
|
|
|
| 521 |
|
|
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 522 |
|
|
|
| 523 |
|
|
subroutine m6swap( m, maxr, nr, ns, ms, |
| 524 |
|
|
$ kb, bbl, bbu, grd, r, rg, x ) |
| 525 |
|
|
|
| 526 |
|
|
implicit double precision (a-h,o-z) |
| 527 |
|
|
integer kb(ms) |
| 528 |
|
|
double precision bbl(ms), bbu(ms), grd(ms) |
| 529 |
|
|
double precision r(nr), rg(ns), x(ms) |
| 530 |
|
|
|
| 531 |
|
|
* ------------------------------------------------------------------ |
| 532 |
|
|
* m6swap (superbasic swap) finds the largest reduced gradient in |
| 533 |
|
|
* the range rg(maxr+1), ..., rg(ns) and swaps it into position |
| 534 |
|
|
* maxr + 1 (so we know where it is). |
| 535 |
|
|
* ------------------------------------------------------------------ |
| 536 |
|
|
|
| 537 |
|
|
k1 = maxr + 1 |
| 538 |
|
|
if (ns .gt. k1) then |
| 539 |
|
|
nz2 = ns - maxr |
| 540 |
|
|
k2 = maxr + idamax( nz2, rg(k1), 1 ) |
| 541 |
|
|
if (k2 .gt. k1) then |
| 542 |
|
|
j = m + k1 |
| 543 |
|
|
k = m + k2 |
| 544 |
|
|
lastr = maxr*k1/2 |
| 545 |
|
|
ldiag1 = lastr + 1 |
| 546 |
|
|
ldiag2 = lastr + (k2 - maxr) |
| 547 |
|
|
|
| 548 |
|
|
rdiag1 = r(ldiag1) |
| 549 |
|
|
rg1 = rg(k1) |
| 550 |
|
|
j1 = kb(j) |
| 551 |
|
|
bl1 = bbl(j) |
| 552 |
|
|
bu1 = bbu(j) |
| 553 |
|
|
grd1 = grd(j) |
| 554 |
|
|
x1 = x(j) |
| 555 |
|
|
|
| 556 |
|
|
r(ldiag1) = r(ldiag2) |
| 557 |
|
|
rg(k1) = rg(k2) |
| 558 |
|
|
kb(j) = kb(k) |
| 559 |
|
|
bbl(j) = bbl(k) |
| 560 |
|
|
bbu(j) = bbu(k) |
| 561 |
|
|
grd(j) = grd(k) |
| 562 |
|
|
x(j) = x(k) |
| 563 |
|
|
|
| 564 |
|
|
r(ldiag2) = rdiag1 |
| 565 |
|
|
rg(k2) = rg1 |
| 566 |
|
|
kb(k) = j1 |
| 567 |
|
|
bbl(k) = bl1 |
| 568 |
|
|
bbu(k) = bu1 |
| 569 |
|
|
grd(k) = grd1 |
| 570 |
|
|
x(k) = x1 |
| 571 |
|
|
end if |
| 572 |
|
|
end if |
| 573 |
|
|
|
| 574 |
|
|
* end of m6swap |
| 575 |
|
|
end |
Parent Directory
| 
Revision Log