1*----------------------------------------------------------------------* 2 subroutine optc_sbspman(lusbsp,lunew,facs,lunew2,nvec,nmaxvec, 3 & itask,ntask,iaddif,ndel_out, 4 & vec1,vec2) 5*----------------------------------------------------------------------* 6* 7* purpose: manage subspace of vectors on file lusbsp. at most nmaxvec 8* vectors will be kept on file. itask contains the tasks: 9* 10* 1: add n vectors (where n is contained in the succeeding 11* field of itask) 12* 2: delete n vectors (n: see above). the number of vectors 13* that need to be deleted anyway to keep the subspace 14* dimensions is subtracted from this 15* 3: restart (all vectors will be deleted) 16* 17* nvec is current size of subspace 18* 19* iaddif = 1: add sums lunew + lunew2 20* iaddif = 1: add diffs lunew - lunew2 21* 22*----------------------------------------------------------------------* 23 implicit none 24 25* constants 26 integer, parameter :: 27 & ntest = 00 28 29* input/output 30 integer, intent(in) :: 31 & lusbsp, lunew, lunew2, 32 & nmaxvec, ntask, 33 & itask(ntask), iaddif 34 integer, intent(inout) :: 35 & nvec 36 integer, intent(out) :: 37 & ndel_out 38 real(8), intent(in) :: 39 & facs(*) 40* scratch 41 real(8), intent(inout) :: 42 & vec1(*), vec2(*) 43* external function 44 integer, external :: 45 & iopen_nus 46 real(8), external :: 47 & inprdd 48* local 49 integer :: 50 & i, ierr, ndel, nadd, nnew, lblk, 51 & luscr, irec 52 real(8) :: 53 & xnrm 54 55 lblk = -1 56 57 if (ntest.ge.10) then 58 write(6,*) 'optc_sbspman:' 59 write(6,*) '==============' 60 write(6,*) ' nvec, nmaxvec : ', nvec, nmaxvec 61 write(6,*) ' ntask,itask :', ntask,'(',itask(1:ntask),')' 62 write(6,*) ' lusbsp, lunew: ',lusbsp, lunew 63 if (iabs(iaddif).eq.1) 64 & write(6,*) ' lnew2 :',lunew2 65 if (ntest.ge.100) then 66 write(6,*) 'Initial contents on subspace-file:' 67 write(6,*) ' record norm(vector)' 68 call rewino(lusbsp) 69 do irec = 1, nvec 70 xnrm = sqrt(inprdd(vec1,vec1,lusbsp,lusbsp,0,lblk)) 71 write(6,'(4x,i5,2x,e15.7)') irec, xnrm 72 end do 73 end if 74 end if 75 76 i = 0 77 ierr = 0 78 ndel = 0 79 nadd = 0 80 do while(i.lt.ntask) 81 i = i+1 82 if (itask(i).eq.1) then 83 i = i+1 84 if (i.gt.ntask) then 85 ierr = 1 86 exit 87 end if 88 nadd = nadd + itask(i) 89 else if (itask(i).eq.2) then 90 i = i+1 91 if (i.gt.ntask) then 92 ierr = 2 93 exit 94 end if 95 ndel = min(ndel+itask(i),nvec) 96 else if (itask(i).eq.3) then 97 ndel = nvec 98 else 99 ierr = 4 100 exit 101 end if 102 end do 103 104 if (ierr.ne.0) then 105 write(6,*) 'Error in optc_sbspman: ',ierr 106 stop 'error in optc_sbspman' 107 end if 108 109 nnew = nvec + nadd - ndel 110 if (nnew.gt.nmaxvec) then 111 ndel = ndel + nnew-nmaxvec 112 end if 113 114 if (ntest.ge.10) then 115 write(6,*) 'nadd, ndel: ',nadd, ndel 116 end if 117 118 if (ndel.gt.0.and.nvec-ndel.gt.0) then 119 luscr = iopen_nus('sbspmanscr') 120 call rewino(luscr) 121 call skpvcd(lusbsp,ndel,vec1,1,lblk) 122 do i = ndel+1,nvec 123 call copvcd(lusbsp,luscr,vec1,0,lblk) 124 end do 125 call rewino(luscr) 126 call rewino(lusbsp) 127 do i = ndel+1,nvec 128 call copvcd(luscr,lusbsp,vec1,0,lblk) 129 end do 130 call relunit(luscr,'delete') 131 else if (ndel.eq.nvec) then 132 call rewino(lusbsp) 133 else if (ndel.eq.0) then 134 call skpvcd(lusbsp,nvec,vec1,1,lblk) 135 else 136 write(6,*) 'Inconsistency in optc_sbspman' 137 stop 'Inconsistency in optc_sbspman' 138 end if 139 140 if (nadd.gt.0) then 141 call rewino(lunew) 142 if (iabs(iaddif).eq.1) call rewino(lunew2) 143 do i = 1, nadd 144 if (iaddif.eq.1) then 145 call vecsmd(vec1,vec2,facs(i),facs(i), 146 & lunew,lunew2,lusbsp,0,lblk) 147 else if (iaddif.eq.-1) then 148 call vecsmd(vec1,vec2,facs(i),-facs(i), 149 & lunew,lunew2,lusbsp,0,lblk) 150 else 151 call sclvcd(lunew,lusbsp,facs(i),vec1,0,lblk) 152 end if 153 end do 154 end if 155 156 ! return new subspace size ... 157 nvec = nvec - ndel + nadd 158 159 ! ... and how many vectors were deleted 160 ndel_out = ndel 161 162 if (ntest.ge.100) then 163 write(6,*) 'Final contents on subspace-file:' 164 write(6,*) ' record norm(vector)' 165 call rewino(lusbsp) 166 do irec = 1, nvec 167 xnrm = sqrt(inprdd(vec1,vec1,lusbsp,lusbsp,0,lblk)) 168 write(6,'(4x,i5,2x,e15.7)') irec, xnrm 169 end do 170 end if 171 172 return 173 174 end 175*----------------------------------------------------------------------* 176*----------------------------------------------------------------------* 177 subroutine optc_diagp(lugrvf,ludia,xmaxstp, 178 & luamp,xdamp,ldamp, 179 & vec1,vec2,nwfpar) 180*----------------------------------------------------------------------* 181 implicit none 182 183 integer, parameter :: 184 & ntest = 00 185 186 logical, intent(in) :: 187 & ldamp 188 integer, intent(in) :: 189 & lugrvf, ludia, luamp, 190 & nwfpar 191 real(8), intent(in) :: 192 & vec1(nwfpar), vec2(nwfpar), 193 & xmaxstp 194 real(8), intent(out) :: 195 & xdamp 196 integer :: 197 & inv, lblk 198 199 logical :: 200 & lconv 201 integer :: 202 & maxiter, iter 203 real(8) :: 204 & xnrm, xgg, xval, xder, xthr, xmx, xdamp_min, 205 & xdamp_old, xval_old, diamin, xstep, xhigh, xlow, 206 & xvhigh, xvlow 207 208 real(8), external :: 209 & inprdd, inprdd3, fdmnxd 210 211 inv=1 212 lblk = -1 213 diamin = fdmnxd(ludia,-1,vec1,1,lblk) 214 xdamp_min = max(-diamin,0d0) 215 if (ntest.ge.100) write(6,*) 'smallest diagonal element: ',diamin 216 217 call dmtvcd2(vec1,vec2,ludia,lugrvf,luamp,-1d0,0d0,1,inv,lblk) 218 219 xnrm = sqrt(inprdd(vec1,vec1,luamp,luamp,1,lblk)) 220 221 if (ntest.ge.100) write(6,*) 'norm of predicted step: ',xnrm 222 xdamp = 0d0 223 224 if ((xnrm.gt.xmaxstp.or.diamin.lt.0d0).and..not.ldamp) then 225 xdamp = 100d0 ! just used as a flag, another routine will 226 ! evaluate the actual damping 227 end if 228 229 xlow = xdamp_min 230 xvlow = +100 231 xhigh= 1d13 232 xvhigh= -100 233 234* step length OK? 235 if ((xnrm.gt.xmaxstp.or.diamin.lt.0d0).and.ldamp) then 236 ! else we have to search for the lowest positive value of 237 ! xdamp such that that condition is fulfilled 238 xgg = inprdd(vec1,vec1,lugrvf,lugrvf,1,lblk) 239 ! a guess, assuming sqrt(<1/H^2> ~ 0.4d0) 240 xdamp = max(sqrt(xgg)/xmaxstp - 0.4d0,0.1d0,-diamin+0.1d0) 241c xdamp = 0 242 maxiter = 100 243 xmx = 1.d0 244 xthr = 1d-12 245 do iter = 1, maxiter 246 ! get the value 247 xval=inprdd3(vec1,vec2,lugrvf,lugrvf,ludia,xdamp,-2d0,1,lblk) 248 & - xmaxstp*xmaxstp 249 250c? if (iter.gt.1.and.sign(1d0,xval).ne.sign(1d0,xval_old) 251c? & .and.abs(xval_old)/10d0.lt.abs(xval)) then 252c? xdamp = xdamp_old - 253c?c & xval_old*(xdamp_old-xdamp)/(xval_old-xval) 254c? & 0.5d0*(xdamp_old-xdamp) 255c? 256c? xval=inprdd3(vec1,vec2,lugrvf,lugrvf,ludia,xdamp,-2d0,1,lblk) 257c? & - xmaxstp*xmaxstp 258c? end if 259 260c if (iter.gt.1.and.abs(xval).gt.abs(xval_old)) then 261c xdamp = xdamp_old + (xdamp_old - xdamp) 262c xval=inprdd3(vec1,vec2,lugrvf,lugrvf,ludia,xdamp,-2d0,1,lblk) 263c & - xmaxstp*xmaxstp 264c end if 265 266 267 ! test convergence 268 lconv = abs(xval).lt.xthr .or. 269 & xval.lt.0d0 .and. abs(xdamp-xdamp_min).lt.xthr 270 if (lconv) exit 271 if (xval.gt.0d0) then 272 xlow = xdamp 273 xvlow = xval 274 else 275 xhigh = xdamp 276 xvhigh = xval 277 end if 278 279 if (xlow.gt.xhigh) stop 'AHA!' 280 281 ! get the derivative 282 xder = 283 & -2d0*inprdd3(vec1,vec2,lugrvf,lugrvf,ludia,xdamp,-3d0,1,lblk) 284 if (ntest.ge.100) then 285 write(6,'(x,a,i4,3(x,e15.6))') 286 & 'iter, damp, val, der ',iter, xdamp, xval, xder 287 end if 288 289 xdamp_old = xdamp 290 ! newton step: 291 xstep = - xval/xder 292 293 if (xdamp+xstep.ge.xhigh) then 294 xdamp = xdamp-xval*(xdamp-xhigh)/(xval-xvhigh) 295c xdamp = (xdamp+xhigh)/2d0 296 else if (xdamp+xstep.lt.xlow) then 297 if (xlow.eq.xdamp_min) xvlow=-xval 298 xdamp = xdamp-xval*(xdamp-xlow)/(xval-xvlow) 299c xdamp = (xdamp+xlow)/2d0 300 else 301 xdamp = xdamp + xstep 302 end if 303 304 xval_old = xval 305 306 end do 307 308 if (.not.lconv) then 309 write(6,*) 'problem with damping!!' 310 stop 'optc_diagp' 311 end if 312 313 if (ntest.ge.100) write(6,*) 'final damping: ',xdamp 314 315 call dmtvcd2(vec1,vec2,ludia,lugrvf,luamp,-1d0,xdamp,1,inv,lblk) 316 317 end if 318 319 return 320 end 321*----------------------------------------------------------------------* 322*----------------------------------------------------------------------* 323 subroutine optc_conjgrad(itype,ilin, 324 & lu_corstp, 325 & lu_ucrstp,lust_sbsp,nrec, 326 & lugrvf,lugrvfold,lusig, 327 & vec1,vec2,nwfpar,iprint) 328*----------------------------------------------------------------------* 329* 330* apply conjugate gradient correction to new search direction: 331* 332* |d> = |d_unc> - beta |d_last> 333* 334* |d_uncorr> is initially on lu_ucrstp 335* |d> will finally be on lu_corstp 336* 337* |d_last> is the last step (record nrec on file lust_sbsp) 338* 339* beta can be calculated according to the following approximations: 340* 341* 1) orthodox: 342* < g - g_last | d_unc > 343* beta = ------------------------ 344* < g - g_last | d_last> 345* 346* 2) Polak-Ribiere: 347* < g - g_last| d_unc > 348* beta = ------------------------- 349* < g_last | d_last > 350* 351* 3) Fletcher-Reeves: 352* < g | d_unc > 353* beta = ------------------------ 354* < g_last | d_last > 355* 356* if ilin==1: 357* < d_last A | d_unc > 358* beta = ------------------------ 359* < d_last A | d_last > 360* 361* where < d_last A | is on lusig 362* 363*----------------------------------------------------------------------* 364 365 implicit none 366 367* constants 368 integer, parameter :: 369 & ntest = 00 370 371* input/output 372 integer, intent(in) :: 373 & itype,ilin, 374 & lu_ucrstp,lu_corstp,lust_sbsp, 375 & lugrvf,lugrvfold,lusig, 376 & nrec,nwfpar,iprint 377 real(8), intent(inout) :: 378 & vec1(nwfpar), vec2(nwfpar) 379 380* local variables 381 integer :: 382 & lblk, luscr, iprintl 383 real(8) :: 384 & xnum, xdnm, beta 385 386*external functions 387 real(8), external :: 388 & inprdd 389 integer, external :: 390 & iopen_nus 391 392 lblk = -1 393 iprintl = max(iprint,ntest) 394 395 if (ntest.ge.10) then 396 write(6,*) 'optc_conjgrad:' 397 write(6,*) '==============' 398 write(6,*) ' itype, ilin : ', itype, ilin 399 write(6,*) ' lu_corstp, lu_ucrstp, lust_sbsp, nrec: ', 400 & lu_corstp, lu_ucrstp, lust_sbsp, nrec 401 write(6,*) ' lugrvf, lugrvfold, lusig: ', 402 & lugrvf, lugrvfold, lusig 403 end if 404 405c if (ilin.eq.0.and.(itype.eq.1.or.itype.eq.2)) then 406 if (itype.eq.1.or.itype.eq.2) then 407 luscr = iopen_nus('DLTGRD') 408 call vecsmd(vec1,vec2,1d0,-1d0,lugrvf,lugrvfold,luscr,1,lblk) 409 ! calc <dlt g|d_unc> 410 xnum = inprdd(vec1,vec2,luscr,lu_ucrstp,1,lblk) 411 if (itype.eq.1) then 412 call skpvcd(lust_sbsp,nrec-1,vec1,1,lblk) 413 call rewino(luscr) 414 ! calc <dlt g|d_last> 415 xdnm = inprdd(vec1,vec2,luscr,lust_sbsp,0,lblk) 416 end if 417 call relunit(luscr,'delete') 418 end if 419 420c if (ilin.eq.0.and.(itype.eq.2.or.itype.eq.3)) then 421 if (itype.eq.2.or.itype.eq.3) then 422 ! calc <g_last|d_last> 423 call skpvcd(lust_sbsp,nrec-1,vec1,1,lblk) 424 call rewino(lugrvfold) 425 xdnm = inprdd(vec1,vec2,lugrvfold,lust_sbsp,0,lblk) 426 if (itype.eq.3) then 427 ! calc < g | d_unc> 428 xnum = inprdd(vec1,vec2,lugrvf,lu_ucrstp,1,lblk) 429 end if 430 end if 431 432c if (ilin.eq.1) then 433c ! calc < d_last A | d_unc > 434c xnum = inprdd(vec1,vec2,lusig,lu_ucrstp,1,lblk) 435cc??? ! calc < d_unc | d_last> 436cc??? call skpvcd(lust_sbsp,nrec-1,vec1,1,lblk) 437cc??? call rewino(lu_ucrstp) 438cc??? xnum = xnum - inprdd(vec1,vec2,lu_ucrstp,lust_sbsp,0,lblk) 439c ! calc < d_last A | d_last> 440c call skpvcd(lust_sbsp,nrec-1,vec1,1,lblk) 441c call rewino(lusig) 442c xdnm = inprdd(vec1,vec2,lusig,lust_sbsp,0,lblk) 443c end if 444 445 if (ntest.ge.10) then 446 write(6,*) 'ilin = ', ilin 447 write(6,*) 'itype = ', itype 448 write(6,*) ' xnum = ', xnum 449 write(6,*) ' xdnm = ', xdnm 450 end if 451 452 beta = xnum / xdnm 453 454c?? 455 if (itype.eq.2.or.itype.eq.3) beta = -beta 456 457 if (iprintl.ge.1) then 458 write (6,'(x,a)') 'conjugate gradient correction:' 459 if (ilin.eq.1) 460 & write(6,'(4x,a,e10.4)') '(lin. equations): beta = ', beta 461 if (ilin.eq.0.and.itype.eq.1) 462 & write(6,'(4x,a,e10.4)') '(orthodox): beta = ', beta 463 if (ilin.eq.0.and.itype.eq.2) 464 & write(6,'(4x,a,e10.4)') '(Polak-Ribiere): beta = ', beta 465 if (ilin.eq.0.and.itype.eq.3) 466 & write(6,'(4x,a,e10.4)') '(Fletcher-Reeves): beta = ', beta 467 end if 468 469 ! obtain corrected step: 470 call rewino(lu_ucrstp) 471 call rewino(lu_corstp) 472 call skpvcd(lust_sbsp,nrec-1,vec1,1,lblk) 473 call vecsmd(vec1,vec2,1d0,-beta,lu_ucrstp,lust_sbsp, 474 & lu_corstp,0,lblk) 475 476 477 return 478 479 end 480*----------------------------------------------------------------------* 481*----------------------------------------------------------------------* 482 subroutine optc_linesearch(ilsrch,ilin,ivar,iprecnd,ipass, 483 & alpha,dalp,energy,de_pred,trrad,xnstp, 484 & vec1,vec2, 485 & lu_corstp,lu_ucrstp,lusig,ludia,iprint) 486*----------------------------------------------------------------------* 487* 488* purpose: obtain step length from an approximate second-order model 489* along the given search direction. 490* 491* One point model: 492* 493* The local model is 494* 495* E(alpha) = E_0 - alpha <g|d> + 1/2 alpha^2 <d|A|d> 496* 497* where the diagonal preconditioner C is used as approximation to A 498* 499* 500* < d | g > < d | g > 501* alpha = - --------------- ~ - --------------- 502* < d | A | d > < d | C | d > 503* 504* if ilin==1: exact < d A | on lusig 505* 506* for prec. steepest descend alpha then always is 1.0 507* 508* | d > is always on lu_corstp 509* | g > is found on lu_ucrstp 510* 511* 512* Two point model: (needs two passes) 513* 514* We predict an initial alpha as above, and set up the local function 515* (along our new step direction): 516* 517* E(alpha) = E_0 + alpha <g|d> + 1/2 alpha^2 h 518* 519* to find the optimum alpha as 520* 521* alpha = |g|/h 522* 523* where h can be estimated from our extra point E(alpha) as 524* 525* h = 2 (E(alpha) - E_0 - alpha <g|d>)/alpha^2 526* = 2 (E(alpha) - E_0 - dE(pred) ) / alpha^2 527* 528* dE(pred) being the energy estimate from the linear model 529* 530*----------------------------------------------------------------------* 531 implicit none 532 533* constants 534 integer, parameter :: 535 & ntest = 00 536 537* input/output 538 integer, intent(in) :: 539 & ilsrch,ilin,ivar,iprecnd,ipass, 540 & lu_corstp,lu_ucrstp,ludia,lusig, 541 & iprint 542 real(8), intent(in) :: 543 & energy, xnstp, trrad 544 real(8), intent(out) :: 545 & alpha, de_pred, dalp 546 real(8), intent(inout) :: 547 & vec1(*), vec2(*) 548 549* local save variables 550 integer, save :: 551 & last_ipass 552 real(8), save :: 553 & de_lin, de_qdr, e_old, alpha_old, alpha_max, xdg 554* local 555 integer :: 556 & lblk, inv, luscr, iprintl 557 real(8) :: 558 & xdad, xdad2, xh, de 559 560* external functions 561 real(8), external :: 562 & inprdd 563 integer, external :: 564 & iopen_nus 565 566 lblk = -1 567 568 iprintl = max(ntest,iprint) 569 570 if (iprintl.ge.1) write (6,*) 'line-search:' 571 572 if (ntest.ge.10) then 573 write(6,*) ' entered optc_linesearch with:' 574 write(6,*) ' ivar,iprecnd, ipass : ',ivar,iprecnd, ipass 575 write(6,*) ' energy,trrad,xnstp : ',energy,trrad,xnstp 576 write(6,*) ' lu_corstp,lu_ucrstp,ludia: ', 577 & lu_corstp,lu_ucrstp,ludia 578 end if 579 580 if (ipass.eq.1) then 581 582 last_ipass = 1 583 584 luscr = iopen_nus('DADSCR') 585 586 ! calc < d | g > or < d | C^-1 g >, respectively 587 xdg = inprdd(vec1,vec2,lu_corstp,lu_ucrstp,1,lblk) 588 589 inv = 0 590 ! calc < d | A | d > 591 call dmtvcd(vec1,vec2,ludia,lu_corstp,luscr,0d0,1,inv,lblk) 592 xdad = inprdd(vec1,vec2,lu_corstp,luscr,1,lblk) 593 alpha = 1d0 594 if (ilsrch.gt.0) 595 & alpha = - xdg / xdad ! interestingly, a factor of 2 is helpful 596 ! with conjugate gradient methods ... 597 ! get max alpha 598 599 if (ntest.ge.10) then 600 write(6,*) ' xdad, xdg: ', xdad, xdg 601 end if 602 603c if (alpha.lt.0d0) stop 'PANIC: alpha.lt.0d0!!!' 604 if (alpha.lt.0d0)write(6,*) '>>> WARNING: alpha.lt.0d0!!!' 605 606 if (ilsrch.eq.2) then 607 if (iprintl.ge.2) 608 & write(6,*) ' reduced trust radius for trial step (.75)' 609 alpha_max = 0.75*trrad/xnstp 610 else 611 alpha_max = trrad/xnstp 612 end if 613 if (abs(alpha).gt.alpha_max) then 614 if (iprintl.ge.2) then 615 write (6,*) ' alpha = ',alpha 616 write (6,*) ' norm step = ',alpha*xnstp 617 write (6,*) ' trust radius = ',trrad 618 end if 619 alpha = alpha_max*sign(1d0,alpha) 620 end if 621c TEST 622 alpha = 1d0 623 write(6,*) 'alpha set to 1d0' 624c 625 if (ivar.eq.1) then 626 alpha_old = alpha 627 e_old = energy 628 de_lin = alpha*xdg 629 de_qdr = alpha*xdg + .5d0 * alpha*alpha*xdad 630 de_pred = de_qdr 631 end if 632 633 if (iprintl.ge.1) then 634 write (6,*) ' alpha from one-point model: ',alpha 635 if (ivar.eq.1) then 636 write (6,*) ' predicted energy change: ',de_pred 637 end if 638 end if 639 640 call relunit(luscr,'delete') 641 642 else if (ipass.eq.2) then 643 if (ivar.eq.0) then 644 write (6,*) 645 & 'Error: 2 point line-search entered for non-variational model' 646 stop 'line-search' 647 end if 648 if (last_ipass.ne.1) then 649 write (6,*) 650 & 'Error: pass 2 without previous pass 1' 651 stop 'line-search' 652 end if 653 654 if (ntest.ge.10) then 655 write(6,*) 'energy, e_old, de_lin, alpha_old: ', 656 & energy, e_old, de_lin, alpha_old 657 end if 658 659 ! values have hopefully been set in previous pass 660 de = energy-e_old 661 if (abs(de).lt.1d-12) then 662 write(6,*) 'Energy difference too small : ', de 663 write(6,*) 'reverting to old alpha' 664 alpha = alpha_old 665 else 666 if (de.lt.0d0) then 667 alpha_max = min(alpha_max,5d0*alpha_old) 668 else 669 alpha_max = 0.5*alpha_old 670 end if 671 if (ntest.ge.10) then 672 write(6,*) 'alpha_max set to ',alpha_max 673 end if 674 xh = 2d0*(de - de_lin) / (alpha_old*alpha_old) 675 if (ntest.ge.10) then 676 write(6,*) ' xh, xdg: ', xh, xdg 677 end if 678 679 if (xh.le.0d0) then 680 alpha = alpha_max 681 if (ntest.ge.5) then 682 write(6,*)'h <= 0, so we step to the trust radius border' 683 write(6,*) 'alpha opt = ',alpha 684 end if 685 else 686 alpha = -xdg/xh 687c if (alpha.lt.0d0) stop 'PANIC: alpha.lt.0d0!!! (2)' 688 if (alpha.lt.0d0) 689 & write(6,*) '>>>WARNING: alpha.lt.0d0!!! (2)' 690 if (ntest.ge.5) then 691 write(6,*) 'alpha opt = ',alpha 692 write(6,*) 'norm of step = ',alpha*xnstp 693 write(6,*) 'trust rad = ',trrad 694 end if 695 end if 696 697 end if 698 699 alpha = sign(min(alpha_max,abs(alpha)),alpha) 700 if (ntest.ge.5) then 701 write(6,*) 'alpha taken = ',alpha 702 end if 703 704 de_pred = alpha * xdg + 0.5d0*xh*alpha*alpha 705 706 if (iprintl.ge.1) then 707 write (6,*) ' alpha from two-point model: ',alpha 708 write (6,*) ' predicted energy change: ',de_pred 709 end if 710 711 ! return also difference to old alpha 712 dalp = alpha - alpha_old 713 714 if (ntest.ge.10) then 715 write(6,*) 'returned alpha, alpha-alpha_old ',alpha,dalp 716 end if 717 718 end if 719 720 return 721 722 end 723*----------------------------------------------------------------------* 724*----------------------------------------------------------------------* 725 subroutine optc_diis(itype,thrsh,nvec_sbsp, 726 & navec,maxvec, 727 & nadd_in,navec_last,alpha_last, 728 & lu_step,lu_corstep,lu_sbsp,lu_sbsp2, 729 & luamp,lugrvf, 730 & bmat,scr,vec1,vec2,iprint) 731*----------------------------------------------------------------------* 732* 733* perform DIIS extrapolation of wavefunction vectors. 734* 735* itype = 1 take actual t_k+1-t_k as error vector for t_k+1 736* and interpolate in (t_k) basis 737* (*not* recommended!!) 738* 739* itype = 2 take perturbation steps delta t_k (= preconditioned 740* vector function) as error vectors (on lu_sbsp) 741* and interpolate in (t_k + delta t_k) basis 742* 743* itype = 3 take vector function as error vectors (on lu_sbsp) 744* and interpolate in (t_k + delta t_k) basis 745* 746* itype = 4 take vector function as error vectors (on lu_sbsp) 747* and interpolate (t_k) and (Omega_k) separately 748* 749*----------------------------------------------------------------------* 750 751 implicit none 752 753* constants 754 integer, parameter :: 755 & ntest = 00 756 757* input/output 758 integer, intent(in) :: 759 & itype, 760 & nvec_sbsp, ! number of vectors in subspace 761 & maxvec, ! max number of vectors allowed 762 & nadd_in, ! number of new vectors on lu_step 763 & lu_step, lu_corstep, lu_sbsp, 764 & lu_sbsp2, luamp, lugrvf, ! unit numbers 765 & iprint ! print level 766c integer, intent(out) :: 767c & ndel_out ! request to sbspman: delete these vectors 768 integer, intent(inout) :: 769 & navec, ! number active vectors in subspace 770 & ! (skip the first nvec-navec vectors) 771 & navec_last ! number active vectors in subspace in 772 ! previous call 773 real(8), intent(in) :: 774 & alpha_last, ! scaling factor for previous vector 775 & thrsh ! threshold for accepting DIIS step 776 real(8), intent(inout) :: 777 & bmat(*), scr(*), vec1(*), vec2(*) ! DIIS-matrix and scratch 778 779* local save 780c integer, save :: 781 782* local O(N) scratch 783 integer :: 784 & kpiv(navec+1) 785 real(8) :: 786 & scrvec(navec+1) 787 788* local 789 logical :: 790 & lincore, again 791 integer :: 792 & lblk, iprintl, navec_l, 793 & nskip, nold, ndel, ii, jj, iioff, iioff2, 794 & lu_sbsp_un, lu_curerr, lust, itype_l 795 real(8) :: 796 & cond, xcorsum 797 798 iprintl = max(iprint,ntest) 799 lblk = -1 800 lincore = .false. 801 802 if (iprintl.ge.5) then 803 write(6,*) 'DIIS-extrapolation' 804 write(6,*) '==================' 805 if (itype.eq.1) then 806 write(6,*) ' error vector: c^(k+1)-c^(k)' 807 else if (itype.eq.2) then 808 write(6,*) ' error vector: c^(k+1)_pert' 809 else if (itype.eq.3.or.itype.eq.4) then 810 write(6,*) ' error vector: grad_c^(k+1)' 811 end if 812 end if 813 814 if (ntest.ge.10) then 815 write(6,*) 'on entry in optc_diis:' 816 write(6,*) 'nvec_sbsp, navec, maxvec, nadd_in, navec_last: ', 817 & nvec_sbsp, navec, maxvec, nadd_in, navec_last 818 write(6,*) 'lu_step, lu_corstep, lu_sbsp: ', 819 & lu_step, lu_corstep, lu_sbsp 820 end if 821 822 if (itype.lt.1.or.itype.gt.4) then 823 write(6,*) 'DIIS: unknown itype = ',itype 824 stop 'DIIS arguments' 825 end if 826 827 828c if (nvec_sbsp+nadd_in.ne.navec) then 829c write(6,*) 'Panic 1: DIIS subspace dimensioning inconsistent!' 830c stop 'DIIS dimensions' 831c end if 832 if (navec.gt.maxvec) then 833 write(6,*) 'Panic 2: DIIS subspace dimensioning inconsistent!' 834 stop 'DIIS dimensions' 835 end if 836 837 nskip = nvec_sbsp+nadd_in-navec 838 ! vectors that are to be skipped on Bmat 839 ! (after deleting), normally 0 840 ndel = navec_last+nadd_in-navec 841 ! vectors that are on Bmat but are no more 842 ! needed 843 nold = navec_last-ndel ! vectors that are already on Bmat and are 844 ! still needed 845 846 navec_last = navec ! save for next iteration 847 848 if (nold.lt.0.or.nskip.lt.0) then 849 write(6,*) 'inconsistent dimensions in optc_diis!' 850 write(6,*) 'nold , nskip : ', nold, nskip 851 stop 'DIIS dimensions' 852 end if 853 854 if (ntest.ge.10) then 855 write(6,*) 'vectors to be deleted in B mat (ndel): ',ndel 856 write(6,*) 'old vectors on B matrix: (nold) ',nold 857 write(6,*) 'skipped vectors on B mat:(nskip)',nskip 858 end if 859 860* vectors deleted? move B matrix by the corresp. number of col.s and rows 861 if (ndel.gt.0) then 862c if (ndel.gt.nskip) then 863c write(6,*) 'Panic 2: DIIS subspace dimensioning inconsistent!' 864c stop 'DIIS dimensions' 865c end if 866 ! nold is the future dimension of the already initialized 867 ! DIIS matrix 868 ! we now move all entries according to 869 ! B(i,j) = B(i+ndel,j+ndel), where ij actually refer to a 870 ! upper triangular matrix 871 if (ntest.ge.100) then 872 write(6,*) 'DIIS matrix before deleting:' 873 call prtrlt(bmat,nold+ndel) 874 end if 875 876 do ii = 1, nold 877 iioff = ii*(ii-1)/2 878 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 879 do jj = 1, ii 880 bmat(iioff+jj) = bmat(iioff2+jj) 881 end do 882 end do 883 if (ntest.ge.100) then 884 write(6,*) 'DIIS matrix after deleting:' 885 call prtrlt(bmat,nold) 886 end if 887 888 end if 889 890* modify most recent vector (if alpha was != 1.0) 891 if (itype.eq.1.and.nold.gt.0.and.alpha_last.ne.1d0) then 892 iioff = nold*(nold-1)/2 893 do ii = 1, nold-1 894 bmat(iioff+ii) = alpha_last*bmat(iioff+ii) 895 end do 896 bmat(iioff+nold) = alpha_last*alpha_last*bmat(iioff+nold) 897 if (ntest.ge.100) then 898 write(6,*) 'DIIS matrix after modifying last vector:' 899 call prtrlt(bmat,nold) 900 end if 901 end if 902 903* add new vectors 904 ndel = 0 905 if (itype.eq.1.or.itype.eq.2) then 906 lu_curerr = lu_step 907 else 908 lu_curerr = lugrvf 909 end if 910 call optc_diis_bmat(bmat,nadd_in,nold,ndel,nskip, 911 & lu_sbsp,lu_curerr,vec1,vec2,lincore) 912 913 again = .true. 914 navec_l = navec 915 do while(again) 916* transfer to scratch array and augment 917 do ii = 1, navec_l 918 iioff = ii*(ii-1)/2 919 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 920 do jj = 1, ii 921 scr(iioff+jj) = bmat(iioff2+jj) 922 end do 923 end do 924 iioff = navec_l*(navec_l+1)/2 925 do jj = 1, navec_l 926 scr(iioff+jj) = -1d0 927 end do 928 scr(iioff+navec_l+1) = 0d0 929 930 if (ntest.ge.100) then 931 write(6,*) 'Augmented matrix passed to dspco:' 932 call prtrlt(scr,navec_l+1) 933 end if 934* solve DIIS problem to obtain new weights 935 ! first factorize the DIIS matrix and get its condition number 936 call dspco(scr,navec_l+1,kpiv,cond,scrvec) 937c call dppco(scr,navec+1,cond,scrvec,info) 938 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 939 & navec_l,cond 940 if (ntest.ge.100) then 941 write(6,*) 'Factorized matrix from dspco:' 942 call prtrlt(scr,navec_l+1) 943 write(6,*) 'Pivot array:' 944 call iwrtma(kpiv,1,navec_l+1,1,navec_l+1) 945 end if 946 947 again = .false. 948c if (itype.eq.1) again = cond.lt.1d-14 949 950 ! if everything is fine ... 951 if (.not.again) then 952 ! ... we set up the RHS vector and solve the DIIS equations 953 scrvec(1:navec_l) = 0d0 954 scrvec(navec_l+1) = -1d0 955 call dspsl(scr,navec_l+1,kpiv,scrvec) 956c call dppsl(scr,navec_l+1,kpiv,scrvec) 957 if (ntest.ge.10) then 958 write(6,*) 'result of dspsl:' 959 write(6,*) 'w:', scrvec(1:navec_l) 960 write(6,*) 'l:', scrvec(navec_l+1) 961 end if 962 else 963 ! ... else we have to reduce the number of active dimensions 964 nskip = nskip+1 965 ndel = ndel +1 966 navec_l = navec_l-1 967 end if 968* analyze solution and event. request deletion of vectors 969 if (.not.again) then 970 xcorsum = 0d0 971 do ii = 1, navec_l-1 972 xcorsum = xcorsum + abs(scrvec(ii)) 973 end do 974 if (xcorsum/abs(scrvec(navec_l)).gt.1.2d0) then 975 again = .true. 976 nskip = nskip + 1 977 ndel = ndel + 1 978 navec_l = navec_l -1 979 end if 980 end if 981 982 end do 983 984* get new step according to weigths 985 if (itype.eq.1.or.itype.eq.4) lu_sbsp_un = lu_sbsp 986 if (itype.eq.2.or. 987 & itype.eq.3) lu_sbsp_un = lu_sbsp2 988 if (itype.eq.1.or.itype.eq.2.or.itype.eq.3) lust = lu_step 989 if (itype.eq.4) lust = lugrvf 990 call optc_diis_nstp(itype,thrsh,scrvec,navec_l,nskip,lu_sbsp_un, 991 & lust,lu_corstep,luamp,vec1,vec2,lincore) 992 993 if (itype.eq.4) then 994 ! extrapolate here new internal step 995 itype_l = 5 996 lust = 0 997 call optc_diis_nstp(itype_l,thrsh,scrvec,navec_l,nskip,lu_sbsp2, 998 & lust,lu_step,luamp,vec1,vec2,lincore) 999 1000 end if 1001 1002 if (itype.eq.1) then 1003* reupdate DIIS B-matrix with actual new step 1004 call optc_diis_bmat(bmat,nadd_in,nold,ndel,nskip, 1005 & lu_sbsp,lu_corstep,vec1,vec2,lincore) 1006 end if 1007 1008* return the actual last DIIS dimension 1009 navec = navec_l 1010 1011 return 1012 end 1013*----------------------------------------------------------------------* 1014*----------------------------------------------------------------------* 1015 subroutine optc_diis_bmat(bmat,nadd,nold,nskipb,nskip, 1016 & lu_sbsp,lu_step, 1017 & vec1,vec2,lincore) 1018*----------------------------------------------------------------------* 1019* 1020* Update the DIIS B matrix ( B = <e_i|e_j> ) with nadd new vectors 1021* from lu_step. 1022* If lincore.eq..true., the vectors vec1 and vec2 can hold a complete 1023* vector from disc. 1024* 1025*----------------------------------------------------------------------* 1026 implicit none 1027 1028 integer, parameter :: 1029 & ntest = 00 1030 1031 logical, intent(in) :: 1032 & lincore 1033 integer, intent(in) :: 1034 & nadd,nold,lu_sbsp,lu_step,nskip,nskipb 1035 real(8), intent(inout) :: 1036 & bmat(*),vec1(*),vec2(*) 1037 1038 real(8) :: 1039 & bij, bii 1040 integer :: 1041 & ii, jj, iioff, iirec, lblk 1042 1043 real(8), external :: 1044 & inprdd 1045 1046 lblk = -1 1047 if (lincore) stop 'no lincore' 1048 if (nadd.gt.1.and..not.lincore) then 1049 write(6,*) 'optc_diis not completely read to add more than '// 1050 & 'vector at a time' 1051 stop 'DIIS: nadd.gt.1' 1052 end if 1053 1054 if (ntest.ge.10) then 1055 write(6,*) 'lu_sbsp, lu_step: ',lu_sbsp, lu_step 1056 write(6,*) 'updating B matrix with ', nadd,' vector(s)' 1057 write(6,*) 'old dimension was ',nold 1058 end if 1059 1060 do ii = nold + 1, nold+nadd 1061 iirec = ii-nold 1062 iioff = ii*(ii-1)/2 1063 call rewino(lu_sbsp) 1064 if (nskip.gt.0) call skpvcd(lu_sbsp,nskip,vec1,1,lblk) 1065 do jj = 1, nskipb 1066 bmat(iioff+jj) = 0d0 1067 end do 1068c if (lincore) call vec_from_disc(vec1,) 1069 do jj = nskipb + 1, nold 1070 call rewino(lu_step) 1071 if (iirec.gt.1) call skpvcd(lu_step,iirec-1,vec1,1,lblk) 1072 bij = inprdd(vec1,vec2,lu_sbsp,lu_step,0,lblk) 1073 bmat(iioff+jj) = bij 1074 end do 1075 ! well, here adding more than one vector becomes a bit tricky 1076 ! unless we can do that incore. we leave it as is for a while 1077 call rewino(lu_step) 1078 if (iirec.gt.1) call skpvcd(lu_step,iirec-1,vec1,1,lblk) 1079 bii = inprdd(vec1,vec2,lu_step,lu_step,0,lblk) 1080 bmat(iioff+nold+1) = bii 1081 end do 1082* 1083 if (ntest.ge.100) then 1084 write(6,*) 'DIIS matrix after adding new vector:' 1085 call prtrlt(bmat,nold+nadd) 1086 end if 1087 1088 return 1089 end 1090 1091*----------------------------------------------------------------------* 1092*----------------------------------------------------------------------* 1093 subroutine optc_diis_nstp(itype,thrsh,wvec,navec,nskip,lusbsp, 1094 & lu_step,lu_corstep,luamp, 1095 & vec1,vec2,lincore) 1096*----------------------------------------------------------------------* 1097* 1098* obtain new step from the DIIS-weigths on wvec, the old steps on 1099* lu_sbsp and the current perturbation estimate on lu_step 1100* 1101* if |diis_correction|/|perturbation_step| > thrsh, the DIIS is rejected 1102* 1103*----------------------------------------------------------------------* 1104 1105 implicit none 1106* constants 1107 integer, parameter :: 1108 & ntest = 00 1109 1110* input/output 1111 logical, intent(in) :: 1112 & lincore 1113 integer, intent(in) :: 1114 & itype,navec, nskip, 1115 & lusbsp,lu_step,lu_corstep,luamp 1116 real(8), intent(in) :: 1117 & thrsh, 1118 & wvec(navec), vec1(*), vec2(*) 1119 1120 integer :: 1121 & ii, jj, lblk, ndim, 1122 & luscr 1123 real(8) :: 1124 & swvec(1:navec), xfac, xnrm 1125 1126 integer, external :: 1127 & iopen_nus 1128 real(8), external :: 1129 & inprdd 1130 1131 if (lincore) stop 'no lincore in optc_diis_newstp' 1132 1133 lblk = -1 1134 luscr = iopen_nus('DIISCR') 1135 1136 if (itype.eq.1.or.itype.eq.5) then 1137 swvec(1:navec) = 0d0 1138 do jj = 1, navec-1 1139 do ii = jj,navec ! jj+1 -> jj ? 1140 swvec(jj) = swvec(jj) + wvec(ii) 1141 end do 1142 swvec(jj) = swvec(jj) - 1d0 1143 end do 1144 swvec(navec) = wvec(navec) 1145 if (itype.eq.5) then 1146 swvec(navec) = swvec(navec) - 1d0 1147 do ii = 1, navec 1148 swvec(ii) = swvec(ii+1) 1149 end do 1150 end if 1151 else if (itype.eq.2.or.itype.eq.3.or.itype.eq.4) then 1152 swvec(1:navec) = wvec(1:navec) 1153 else 1154 write(6,*) 'unexpected itype in optc_diis_nstep: ',itype 1155 stop 'optc_diis_nstep' 1156 end if 1157 1158 if (ntest.ge.10) then 1159 write(6,*) 1160 & 'optc_diis_nstep: weight vector for new step: (itype = ', 1161 & itype,')' 1162 ndim = navec 1163 if (itype.eq.5) ndim = navec-1 1164 call wrtmat(swvec,ndim,1,ndim,1) 1165 end if 1166 1167 call skpvcd(lusbsp,nskip,vec1,1,lblk) 1168 call rewino(lu_corstep) 1169 call rewino(luscr) 1170 1171 if (navec-1.gt.0) then 1172 call mvcsmd(lusbsp,swvec,luscr,lu_corstep, ! lu_corstep is only scratch 1173 & vec1,vec2,navec-1,0,lblk) 1174 if (lu_step.ne.0.and.itype.ne.5) then 1175 call vecsmd(vec1,vec2,1d0,swvec(navec),luscr,lu_step, 1176 & lu_corstep,1,lblk) 1177 else 1178 call copvcd(luscr,lu_corstep,vec1,1,lblk) 1179 end if 1180 if (itype.eq.2.or.itype.eq.3) then 1181* to be consistent, we add only (w_n - 1)times the current amplitudes 1182 call vecsmd(vec1,vec2,1d0,swvec(navec)-1d0, 1183 & lu_corstep,luamp, 1184 & luscr,1,lblk) 1185* too much copying, but for the moment it is OK 1186 call copvcd(luscr,lu_corstep,vec1,1,lblk) 1187 end if 1188 else 1189 if (lu_step.ne.0.and.itype.ne.5) 1190 & call sclvcd(lu_step,lu_corstep,swvec(navec),vec1,1,lblk) 1191 end if 1192 1193 ! get norm of DIIS correction: 1194 if (lu_step.ne.0.and.itype.ne.5.and.itype.ne.4) then 1195 call vecsmd(vec1,vec2,1d0,-1d0,lu_corstep,lu_step,luscr,1,lblk) 1196 xnrm = sqrt(inprdd(vec1,vec2,luscr,luscr,1,lblk)) 1197 if (ntest.ge.10) write(6,*) '|DIIS correction:| = ',xnrm 1198 if (xnrm.gt.thrsh) then 1199 if (ntest.ge.10) write(6,*) 'DIIS step not accepted!' 1200 call copvcd(lu_step,lu_corstep,vec1,1,lblk) 1201 end if 1202 end if 1203 1204 call relunit(luscr,'delete') 1205 1206 return 1207 end 1208*----------------------------------------------------------------------* 1209*----------------------------------------------------------------------* 1210 subroutine optc_sbspja(itype,thrsh,nstdim,ngvdim, 1211 & navec,maxvec, 1212 & nadd,navec_last, 1213 & lugrvf,ludia,trrad,lu_pertstp,lu_newstp, 1214 & xdamp,xdamp_last, 1215 & lust_sbsp,lugv_sbsp, 1216 & umat,scr,vec1,vec2,iprint) 1217*----------------------------------------------------------------------* 1218* 1219* subspace jacobian: A = A_0(1-P) + A_ex P 1220* 1221* if xdamp is 0d0, we try the direct inversion procedure; if that 1222* step is too large, we go over to the iterative solver to find the 1223* optimal damping (where we end up directly, if xdamp was != 0d0) 1224* 1225*----------------------------------------------------------------------* 1226 implicit none 1227 1228* constants 1229 integer, parameter :: 1230 & ntest = 00 1231 1232* input/output 1233 integer, intent(in) :: 1234 & itype, ! type (s.above) 1235 & nstdim,ngvdim, ! # records on sbsp files 1236 & maxvec, ! max dimension of sbsp 1237 & nadd, ! number of new vectors (usually 1) 1238 & lugrvf,ludia, ! current gradient, diagonal Hess/Jac. 1239 & lu_pertstp, ! step from diagonal Hess/Jac. 1240 & lu_newstp, ! new step (output) 1241 & lust_sbsp,lugv_sbsp, ! sbsp files 1242 & iprint ! print level 1243 integer, intent(inout) :: 1244 & navec, ! current dimension of sbsp 1245 & navec_last ! previous sbsp dimension 1246 real(8), intent(in) :: 1247 & thrsh, ! thresh for accepting low-rank correction 1248 & trrad, ! trust radius for step 1249 & xdamp_last ! previous damping 1250 real(8), intent(out) :: 1251 & xdamp ! current damping 1252 real(8), intent(inout) :: 1253 & umat(*), ! low rank Hessian/Jacobian 1254 & scr(*), vec1(*), vec2(*) ! scratch 1255 1256* local O(N) scratch 1257 integer :: 1258 & kpiv(navec) 1259 real(8) :: 1260 & scrvec(navec), scrvec2(navec) 1261 1262* local 1263 logical :: 1264 & lincore, again, accept, converged 1265 integer :: 1266 & lblk, iprintl, navec_l, isym, job, 1267 & nskipst, nskipgv, nnew, nold, nold2, ndel, 1268 & ii, jj, iioff, iioff2, irhsoff, 1269 & lu_sbsp_un, lu_curerr, 1270 & luvec, luAvec, nadd_l, nold_l, nskip_l, ludum, 1271 & iter 1272 real(8) :: 1273 & cond, xcorsum 1274 integer, external :: 1275 & iopen_nus 1276 1277 iprintl = max(iprint,ntest) 1278 lblk = -1 1279 lincore = .false. 1280 if (itype.eq.1) isym = 0 ! non-symmetric jacobian 1281 if (itype.eq.2) isym = 1 ! symmetric jacobian 1282 1283 if (iprintl.ge.5) then 1284 write(6,*) 'Subspace-Hessian/Jacobian' 1285 write(6,*) '=========================' 1286 if (itype.eq.1) then 1287 write(6,*) ' asymmetric update' 1288 else if (itype.eq.2) then 1289 write(6,*) ' symmetric update' 1290 end if 1291 end if 1292 1293 if (ntest.ge.10) then 1294 write(6,*) 'on entry in optc_sbspja:' 1295 write(6,*) 1296 & 'nstdim, ngvdim, navec, maxvec, nadd, navec_last: ', 1297 & nstdim, ngvdim, navec, maxvec, nadd, navec_last 1298 write(6,*) 1299 & 'lugrvf, lu_newstp, lu_pertstp, lust_sbsp, lugv_sbsp: ', 1300 & lugrvf, lu_newstp, lu_pertstp, lust_sbsp, lugv_sbsp 1301 end if 1302 1303 if (itype.lt.1.or.itype.gt.2) then 1304 write(6,*) 'SBSP hessian/jacobian: unknown itype = ',itype 1305 stop 'SBSPJA arguments' 1306 end if 1307 1308 if (navec.gt.maxvec) then 1309 write(6,*) 1310 & 'Panic 2: Hess./Jac. subspace dimensioning inconsistent!' 1311 stop 'SBSP dimensions' 1312 end if 1313 1314 nskipst = nstdim-navec 1315 nskipgv = ngvdim-navec 1316 ndel = navec_last+nadd-navec 1317 nold = navec_last-ndel 1318 1319 navec_last = navec ! save for next iteration 1320 1321 if (nold.lt.0.or.nskipst.lt.0.or.nskipgv.lt.0) then 1322 write(6,*) 'inconsistent dimensions in optc_sbspja!' 1323 write(6,*) 'nold , nskipst, nskipgv : ', nold, nskipst, nskipgv 1324 stop 'SBSPJA dimensions' 1325 end if 1326 1327 if (ntest.ge.10) then 1328 write(6,*) 'vectors to be deleted in U mat (ndel): ',ndel 1329 write(6,*) 'old vectors on U matrix: (nold) ',nold 1330 write(6,*) 'skipped vectors on step file:(nskipst)',nskipst 1331 write(6,*) 'skipped vectors on grad file:(nskipgv)',nskipgv 1332 end if 1333 1334* vectors deleted? move B matrix by the corresp. number of col.s and rows 1335 if (ndel.gt.0) then 1336 if (isym.eq.1) then ! symmetric U matrix 1337 ! nold is the future dimension of the already initialized 1338 ! low-rank matrix 1339 ! we now move all entries according to 1340 ! U(i,j) = U(i+ndel,j+ndel), where ij actually refer to a 1341 ! upper triangular matrix 1342 if (ntest.ge.100) then 1343 write(6,*) 'Low-rank matrix before deleting:' 1344 call prtrlt(umat,nold+ndel) 1345 end if 1346 1347 do ii = 1, nold 1348 iioff = ii*(ii-1)/2 1349 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 1350 do jj = 1, ii 1351 umat(iioff+jj) = umat(iioff2+jj) 1352 end do 1353 end do 1354 if (ntest.ge.100) then 1355 write(6,*) 'Low-rank matrix after deleting:' 1356 call prtrlt(umat,nold) 1357 end if 1358 else ! non-symmetric (and therefore full) matrix 1359 if (ntest.ge.100) then 1360 write(6,*) 'Low-rank matrix before deleting:' 1361 call wrtmat2(umat,nold+ndel,maxvec,nold+ndel,maxvec) 1362 end if 1363 1364 do ii = 1, nold 1365 iioff = (ii-1)*maxvec 1366 iioff2 = (ii+ndel-1)*maxvec + ndel 1367 do jj = 1, nold 1368 umat(iioff+jj) = umat(iioff2+jj) 1369 end do 1370 end do 1371 if (ntest.ge.100) then 1372 write(6,*) 'Low-rank matrix after deleting:' 1373 call wrtmat2(umat,nold,maxvec,nold,maxvec) 1374 end if 1375 1376 end if ! isym 1377 1378 end if ! ndel.gt.0 1379 1380* add new vectors 1381 ndel = 0 1382 nnew = nadd 1383 nold2 = nold 1384 call optc_sbspja_umat(itype, 1385 & umat,nnew,nold2,ndel,maxvec, 1386 & nskipst,nskipgv, 1387 & lust_sbsp,lugv_sbsp,ludia,xdamp, 1388 & vec1,vec2,lincore) 1389 1390* get projection of perturbation step in current subspace 1391 call optc_sbspja_prjpstp(scrvec2,navec,nskipst, 1392 & lu_pertstp,lust_sbsp, 1393 & vec1,vec2,lincore) 1394 1395 navec_l = navec 1396 1397* direct-inversion section: 1398* ========================= 1399 again = .true. 1400 do while(again) 1401* transfer to scratch array 1402 if (isym.eq.1) then 1403 do ii = 1, navec_l 1404 iioff = ii*(ii-1)/2 1405 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 1406 do jj = 1, ii 1407 scr(iioff+jj) = umat(iioff2+jj) 1408 end do 1409 end do 1410 if (ntest.ge.100) then 1411 write(6,*) 'Low-rank matrix passed to dspco:' 1412 call prtrlt(scr,navec_l) 1413 end if 1414 else 1415 do ii = 1, navec_l 1416 iioff = (ii-1)*navec_l 1417 iioff2 = (ii+ndel-1)*maxvec + ndel 1418 do jj = 1, navec_l 1419 scr(iioff+jj) = umat(iioff2+jj) 1420 end do 1421 end do 1422 if (ntest.ge.100) then 1423 write(6,*) 'Low-rank matrix passed to dgeco:' 1424 call wrtmat2(scr,navec_l,navec_l,navec_l,navec_l) 1425 end if 1426 end if 1427 1428* invert low-rank matrix: 1429* factorize and get condition 1430 if (isym.eq.1) then 1431 call dspco(scr,navec_l,kpiv,cond,scrvec) 1432 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 1433 & navec_l,cond 1434 if (ntest.ge.100) then 1435 write(6,*) 'Factorized matrix from dspco:' 1436 call prtrlt(scr,navec_l) 1437 write(6,*) 'Pivot array:' 1438 call iwrtma(kpiv,1,navec_l,1,navec_l) 1439 end if 1440 else 1441 call dgeco(scr,navec_l,navec_l,kpiv,cond,scrvec) 1442 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 1443 & navec_l,cond 1444 if (ntest.ge.100) then 1445 write(6,*) 'Factorized matrix from dgeco:' 1446 call wrtmat2(scr,navec_l,navec_l,navec_l,navec_l) 1447 write(6,*) 'Pivot array:' 1448 call iwrtma(kpiv,1,navec_l,1,navec_l) 1449 end if 1450 end if ! isym 1451 1452 again = cond.lt.1d-14 1453 1454 if (.not.again) then 1455 irhsoff = navec - navec_l + 1 1456 scrvec(1:navec_l) = scrvec2(irhsoff:irhsoff-1+navec_l) 1457* if OK, solve U x = rhs 1458 if (ntest.ge.100) then 1459 write(6,*) 'RHS: ',scrvec(1:navec_l) 1460 end if 1461 if (isym.eq.1) then 1462 call dspsl(scr,navec_l,kpiv,scrvec) 1463 else 1464 job = 0 1465 call dgesl(scr,navec_l,navec_l,kpiv,scrvec,job) 1466 end if 1467 if (ntest.ge.100) then 1468 write(6,*) 'Result for x:' 1469 call wrtmat(scrvec,navec_l,1,navec_l,1) 1470 end if 1471 end if 1472 1473 if (.not.again) then 1474* get low-rank contribution to new vector 1475 call optc_sbspja_lrstep(itype,thrsh,accept,scrvec,navec_l, 1476 & nskipst,nskipgv, 1477 & lust_sbsp,lugv_sbsp,ludia,xdamp, 1478 & lu_pertstp,lu_newstp, 1479 & vec1,vec2,lincore) 1480* and analyze it 1481 if (.not.accept.and.navec_l.gt.1) then 1482 again = .true. 1483 else if (.not.accept) then 1484 navec_l = 0 1485 end if 1486 end if 1487 1488 if (again) then 1489 nskipst = nskipst+1 1490 nskipgv = nskipgv+1 1491 ndel = ndel+1 1492 navec_l = navec_l-1 1493 end if 1494 1495 end do 1496 1497* return the actual last subspace dimension 1498 navec = navec_l 1499 1500 return 1501 1502 end 1503*----------------------------------------------------------------------* 1504*----------------------------------------------------------------------* 1505 subroutine optc_sbspja_new(itype,thrsh,nstdim,ngvdim, 1506 & navec,maxvec, 1507 & nadd,navec_last, 1508 & lugrvf,ludia,trrad,lu_pertstp,lu_newstp, 1509 & xdamp,xdamp_last, 1510 & lust_sbsp,lugv_sbsp, 1511 & umat,scr,vec1,vec2,iprint) 1512*----------------------------------------------------------------------* 1513* 1514* subspace jacobian: A(xdamp) = A_0(1-P) + A_ex P + xdamp*1 1515* 1516* xdamp is adjusted according to the trust radius on trrad 1517* 1518*----------------------------------------------------------------------* 1519 implicit none 1520 1521* constants 1522 integer, parameter :: 1523 & ntest = 00, mxd_iter = 100 1524 real(8), parameter :: 1525 & thrdamp = 1d-4, xinc = 1d-5, xfailfac = 1d5, 1526 & xmxdxdamp = 0.1d0 1527 1528* input/output 1529 integer, intent(in) :: 1530 & itype, ! type (s.above) 1531 & nstdim,ngvdim, ! # records on sbsp files 1532 & maxvec, ! max dimension of sbsp 1533 & nadd, ! number of new vectors (usually 1) 1534 & lugrvf,ludia, ! current gradient, diagonal Hess/Jac. 1535 & lu_pertstp, ! step from diagonal Hess/Jac. 1536 & lu_newstp, ! new step (output) 1537 & lust_sbsp,lugv_sbsp, ! sbsp files 1538 & iprint ! print level 1539 integer, intent(inout) :: 1540 & navec, ! current dimension of sbsp 1541 & navec_last ! previous sbsp dimension 1542 real(8), intent(in) :: 1543 & thrsh, ! thresh for accepting low-rank correction 1544 & trrad, ! trust radius for step 1545 & xdamp_last ! previous damping 1546 real(8), intent(out) :: 1547 & xdamp ! current damping 1548 real(8), intent(inout) :: 1549 & umat(*), ! low rank Hessian/Jacobian 1550 & scr(*), vec1(*), vec2(*) ! scratch 1551 1552* local O(N) scratch 1553 integer :: 1554 & kpiv(navec) 1555 real(8) :: 1556 & scrvec(navec), scrvec2(navec), xd(mxd_iter), xv(mxd_iter) 1557 1558* local 1559 logical :: 1560 & lincore, again, accept, converged 1561 integer :: 1562 & lblk, iprintl, navec_l, isym, job, 1563 & nskipst, nskipgv, nnew, nold, nold2, ndel, 1564 & ii, jj, iioff, iioff2, irhsoff, 1565 & lu_sbsp_un, lu_curerr, 1566 & luvec, luAvec, luvec_sbsp, luAvec_sbsp, lures, 1567 & nadd_l, nold_l, nskip_l, ludum, iter, maxiter, 1568 & kend, ksmat, ksred, kared, kevec, kscr, 1569 & itype_evp, itask(2), ntask, ndum, 1570 & ixd_iter,isubcnt 1571 real(8) :: 1572 & cond, xcorsum, xlam, xresnrm, fac, thrsh_l, xnrm, 1573 & xdiff, xgr, xhs, f0, f1, f2, a1, a2, xdis, dxdamp, xmxdxdmp2, 1574 & xval, xlow, xhigh, xvlow, xvhigh, xdamp_min 1575 integer, external :: 1576 & iopen_nus 1577 real(8), external :: 1578 & inprdd, fdmnxd 1579 1580 iprintl = max(iprint,ntest) 1581 lblk = -1 1582 lincore = .false. 1583 if (itype.eq.1) isym = 0 ! non-symmetric jacobian 1584 if (itype.eq.2) isym = 1 ! symmetric jacobian 1585 1586 if (iprintl.ge.5) then 1587 write(6,*) 'Subspace-Hessian/Jacobian' 1588 write(6,*) '=========================' 1589 if (itype.eq.1) then 1590 write(6,*) ' asymmetric update' 1591 else if (itype.eq.2) then 1592 write(6,*) ' symmetric update' 1593 end if 1594 end if 1595 1596 if (ntest.ge.10) then 1597 write(6,*) 'on entry in optc_sbspja:' 1598 write(6,*) 1599 & 'nstdim, ngvdim, navec, maxvec, nadd, navec_last: ', 1600 & nstdim, ngvdim, navec, maxvec, nadd, navec_last 1601 write(6,*) 1602 & 'lugrvf, lu_newstp, lu_pertstp, lust_sbsp, lugv_sbsp: ', 1603 & lugrvf, lu_newstp, lu_pertstp, lust_sbsp, lugv_sbsp 1604 end if 1605 1606 if (itype.lt.1.or.itype.gt.2) then 1607 write(6,*) 'SBSP hessian/jacobian: unknown itype = ',itype 1608 stop 'SBSPJA arguments' 1609 end if 1610 1611 if (navec.gt.maxvec) then 1612 write(6,*) 1613 & 'Panic 2: Hess./Jac. subspace dimensioning inconsistent!' 1614 stop 'SBSP dimensions' 1615 end if 1616 1617 nskipst = nstdim-navec 1618 nskipgv = ngvdim-navec 1619 ndel = navec_last+nadd-navec 1620 nold = navec_last-ndel 1621 1622 navec_last = navec ! save for next iteration 1623 1624 if (nold.lt.0.or.nskipst.lt.0.or.nskipgv.lt.0) then 1625 write(6,*) 'inconsistent dimensions in optc_sbspja!' 1626 write(6,*) 'nold , nskipst, nskipgv : ', nold, nskipst, nskipgv 1627 stop 'SBSPJA dimensions' 1628 end if 1629 1630 if (ntest.ge.10) then 1631 write(6,*) 'vectors to be deleted in U mat (ndel): ',ndel 1632 write(6,*) 'old vectors on U matrix: (nold) ',nold 1633 write(6,*) 'skipped vectors on step file:(nskipst)',nskipst 1634 write(6,*) 'skipped vectors on grad file:(nskipgv)',nskipgv 1635 end if 1636 1637* vectors deleted? move B matrix by the corresp. number of col.s and rows 1638 if (ndel.gt.0) then 1639 if (isym.eq.1) then ! symmetric U matrix 1640 ! nold is the future dimension of the already initialized 1641 ! low-rank matrix 1642 ! we now move all entries according to 1643 ! U(i,j) = U(i+ndel,j+ndel), where ij actually refer to a 1644 ! upper triangular matrix 1645 if (ntest.ge.100) then 1646 write(6,*) 'Low-rank matrix before deleting:' 1647 call prtrlt(umat,nold+ndel) 1648 end if 1649 1650 do ii = 1, nold 1651 iioff = ii*(ii-1)/2 1652 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 1653 do jj = 1, ii 1654 umat(iioff+jj) = umat(iioff2+jj) 1655 end do 1656 end do 1657 if (ntest.ge.100) then 1658 write(6,*) 'Low-rank matrix after deleting:' 1659 call prtrlt(umat,nold) 1660 end if 1661 else ! non-symmetric (and therefore full) matrix 1662 if (ntest.ge.100) then 1663 write(6,*) 'Low-rank matrix before deleting:' 1664 call wrtmat2(umat,nold+ndel,maxvec,nold+ndel,maxvec) 1665 end if 1666 1667 do ii = 1, nold 1668 iioff = (ii-1)*maxvec 1669 iioff2 = (ii+ndel-1)*maxvec + ndel 1670 do jj = 1, nold 1671 umat(iioff+jj) = umat(iioff2+jj) 1672 end do 1673 end do 1674 if (ntest.ge.100) then 1675 write(6,*) 'Low-rank matrix after deleting:' 1676 call wrtmat2(umat,nold,maxvec,nold,maxvec) 1677 end if 1678 1679 end if ! isym 1680 1681 end if ! ndel.gt.0 1682 1683 ixd_iter = 0 ! iteration counter for xdamp optimization 1684 isubcnt = 0 ! counter for optimizer steps 1685 converged = .false. 1686 1687 xdamp_min = max(-fdmnxd(ludia,-1,vec1,1,lblk),0d0) 1688 xlow = xdamp_min 1689 xhigh = 1d13 1690 xvlow = 100d0 1691 xvhigh = -100d0 1692 1693c test: if we damp, start from a somewhat larger value 1694CCC if (xdamp.gt.0d0) xdamp = xdamp+0.5d0 1695 1696 do while(.not.converged) 1697 1698* add new vectors 1699 ndel = 0 1700 if (xdamp.eq.0d0) then 1701c nnew = nadd 1702c nold2 = nold 1703c for the moment: rebuild everything, everytime 1704 nnew = nold+nadd 1705 nold2 = 0 1706 else 1707 nnew = nold+nadd 1708 nold2=0 1709 call optc_smat(umat,nnew,0,nold2, 1710 & lust_sbsp,nskipst,nold+nadd, 1711 & ludum,0,0, 1712 & vec1,vec2,lincore) 1713 end if 1714 1715 call optc_sbspja_umat(itype, 1716 & umat,nnew,nold2,ndel,maxvec, 1717 & nskipst,nskipgv, 1718 & lust_sbsp,lugv_sbsp,ludia,xdamp, 1719 & vec1,vec2,lincore) 1720 1721* for xdamp.ne.0d0, we produce the perturbation step here: 1722c if (xdamp.ne.0d0) then 1723 ! to be sure, rebuild the perturbation step always here: 1724 call dmtvcd2(vec1,vec2,ludia,lugrvf,lu_pertstp, 1725 & -1d0,xdamp,1,1,lblk) 1726c end if 1727 1728* get projection of perturbation step in current subspace 1729 call optc_sbspja_prjpstp(scrvec2,navec,nskipst, 1730 & lu_pertstp,lust_sbsp, 1731 & vec1,vec2,lincore) 1732 1733 navec_l = navec 1734 1735 again = .true. 1736 do while(again) 1737* transfer to scratch array 1738 if (isym.eq.1) then 1739 do ii = 1, navec_l 1740 iioff = ii*(ii-1)/2 1741 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 1742 do jj = 1, ii 1743 scr(iioff+jj) = umat(iioff2+jj) 1744 end do 1745 end do 1746 if (ntest.ge.100) then 1747 write(6,*) 'Low-rank matrix passed to dspco:' 1748 call prtrlt(scr,navec_l) 1749 end if 1750 else 1751 do ii = 1, navec_l 1752 iioff = (ii-1)*navec_l 1753 iioff2 = (ii+ndel-1)*maxvec + ndel 1754 do jj = 1, navec_l 1755 scr(iioff+jj) = umat(iioff2+jj) 1756 end do 1757 end do 1758 if (ntest.ge.100) then 1759 write(6,*) 'Low-rank matrix passed to dgeco:' 1760 call wrtmat2(scr,navec_l,navec_l,navec_l,navec_l) 1761 end if 1762 end if 1763 1764* invert low-rank matrix: 1765* factorize and get condition 1766 if (isym.eq.1) then 1767 call dspco(scr,navec_l,kpiv,cond,scrvec) 1768 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 1769 & navec_l,cond 1770 if (ntest.ge.100) then 1771 write(6,*) 'Factorized matrix from dspco:' 1772 call prtrlt(scr,navec_l) 1773 write(6,*) 'Pivot array:' 1774 call iwrtma(kpiv,1,navec_l,1,navec_l) 1775 end if 1776 else 1777 call dgeco(scr,navec_l,navec_l,kpiv,cond,scrvec) 1778 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 1779 & navec_l,cond 1780 if (ntest.ge.100) then 1781 write(6,*) 'Factorized matrix from dgeco:' 1782 call wrtmat2(scr,navec_l,navec_l,navec_l,navec_l) 1783 write(6,*) 'Pivot array:' 1784 call iwrtma(kpiv,1,navec_l,1,navec_l) 1785 end if 1786 end if ! isym 1787 1788c we do not care for the moment 1789c again = cond.lt.1d-14 1790 again = .false. 1791 1792 if (.not.again) then 1793 irhsoff = navec - navec_l + 1 1794 scrvec(1:navec_l) = scrvec2(irhsoff:irhsoff-1+navec_l) 1795* if OK, solve U x = rhs 1796 if (ntest.ge.100) then 1797 write(6,*) 'RHS: ',scrvec(1:navec_l) 1798 end if 1799 if (isym.eq.1) then 1800 call dspsl(scr,navec_l,kpiv,scrvec) 1801 else 1802 job = 0 1803 call dgesl(scr,navec_l,navec_l,kpiv,scrvec,job) 1804 end if 1805 if (ntest.ge.100) then 1806 write(6,*) 'Result for x:' 1807 call wrtmat(scrvec,navec_l,1,navec_l,1) 1808 end if 1809 end if 1810 1811 if (.not.again) then 1812* get low-rank contribution to new vector 1813c TEST unset thrsh 1814 thrsh_l = 1d10 1815 call optc_sbspja_lrstep(itype,thrsh_l,accept,scrvec,navec_l, 1816 & nskipst,nskipgv, 1817 & lust_sbsp,lugv_sbsp,ludia,xdamp, 1818 & lu_pertstp,lu_newstp, 1819 & vec1,vec2,lincore) 1820 1821 xnrm = sqrt(inprdd(vec1,vec1,lu_newstp,lu_newstp,1,lblk)) 1822 1823 xval = xnrm-trrad 1824 converged = xnrm.le.trrad.and.(xdamp-xdamp_min).le.thrdamp 1825 & .or.( xdamp.ne.0d0 .and. (abs(xnrm-trrad).lt.thrdamp)) 1826 1827 if (xval.gt.0d0) then 1828 xlow = xdamp 1829 xvlow = xval 1830 else 1831 xhigh = xdamp 1832 xvhigh = xval 1833 end if 1834 write(6,'(">>",4x,2i4,3(2x,e20.5))') 1835 & ixd_iter+1,isubcnt,xdamp,xnrm,xnrm-trrad 1836 1837* convergence control for xdamp: 1838* the step-length function should for any xdamp larger than the 1839* negative of the lowest eigenvalue of the approximate Jacobian 1840* be a monotonically decreasing konvex function. If our Quasi-Newton 1841* algorithm encounters problems, that can only mean that we are still 1842* not beyond the lowest eigenvalue. 1843 if (.not.converged) then 1844 if (xdamp.lt.0d0) then ! that should never ever happen!! 1845 write(6,*) 'What did you do? xdamp = ',xdamp 1846 stop 'optc_sbspja: xdamn' 1847 end if 1848 ixd_iter = ixd_iter+1 1849 if (ixd_iter.ge.mxd_iter) then 1850 write(6,*) 'unsolvable problems finding xdamp' 1851 stop 'optc_sbspja: xdamp' 1852 end if 1853 xd(ixd_iter) = xdamp 1854 xv(ixd_iter) = xnrm-trrad 1855 if (isubcnt.eq.0) then 1856 ! restart at a very new point 1857 ! make a small increment to get the gradient 1858 if (xv(ixd_iter).gt.0d0) then 1859 dxdamp = + xinc 1860 else 1861 dxdamp = - xinc 1862 end if 1863 if (ntest.ge.150) 1864 & write(6,*) 1865 & 'xdamp> initialized num. gradient with xinc = ', 1866 & xinc 1867 ! next time we see us in step two 1868 isubcnt = 2 1869 else if (isubcnt.eq.1) then 1870 ! this is step one of a numerical Newton optimization 1871 ! make a small increment to get the gradient 1872 if (xv(ixd_iter).gt.0d0) then 1873 dxdamp = + xinc 1874 else 1875 dxdamp = - xinc 1876 end if 1877 if (ntest.ge.150) 1878 & write(6,*) 1879 & 'xdamp> initialized num. gradient with xinc = ', 1880 & xinc 1881 ! next time we see us in step two 1882 isubcnt = 2 1883 else if (isubcnt.eq.2) then 1884 ! get gradient from finite difference 1885 if (ntest.ge.150) then 1886 write(6,*) 'xdamp> linear model ' 1887 write(6,*) ' points used:' 1888 do ii = 0, 1 1889 write(6,'(3x,i2,2(2x,e25.8))') 1890 & ii, xd(ixd_iter-ii), xv(ixd_iter-ii) 1891 end do 1892 end if 1893 xdiff = xd(ixd_iter)-xd(ixd_iter-1) 1894 xgr = (xv(ixd_iter)-xv(ixd_iter-1))/xdiff 1895 if (ntest.ge.150) 1896 & write(6,*) 'xdamp> current num. gradient: ',xgr 1897 ! gradient has to be negative 1898 if (xgr.lt.0d0) then 1899 if (ntest.ge.150) 1900 & write(6,*) 'xdamp> accepting step' 1901 ! make a Newton step 1902c dxdamp = - xv(ixd_iter-1)/xgr 1903 dxdamp = - xv(ixd_iter)/xgr 1904 isubcnt = 1 1905 if (ntest.eq.150) 1906 & write(6,*) 'xdamp> step = ',dxdamp 1907c switch to three point formula disabled 1908c if (dxdamp.lt.0.5d0) isubcnt = 3 1909 if (dxdamp.lt.0.5d0) isubcnt = 2 1910 else 1911 if (ntest.ge.150) 1912 & write(6,*) 'xdamp> search at a new place' 1913 ! retry with a larger xdamp 1914 dxdamp = + xfailfac*xinc 1915 ixd_iter = ixd_iter-1 1916 isubcnt = 0 1917 end if 1918 else if (isubcnt.eq.3) then 1919 if (ntest.ge.150) then 1920 write(6,*) 'xdamp> quadratic model' 1921 write(6,*) ' points used:' 1922 do ii = 0, 2 1923 write(6,'(3x,i2,2(2x,e25.8))') 1924 & ii, xd(ixd_iter-ii), xv(ixd_iter-ii) 1925 end do 1926 end if 1927 ! for convenience: 1928 f0 = xv(ixd_iter) 1929 f1 = xv(ixd_iter-1) 1930 f2 = xv(ixd_iter-2) 1931 a1 = xd(ixd_iter-1)-xd(ixd_iter) 1932 a2 = xd(ixd_iter-2)-xd(ixd_iter) 1933 ! the points should correspond to a monotonic 1934 ! decreasing function 1935 if (.not.(a1*(f0-f1).gt.0d0.and.a2*(f0-f2).gt.0d0)) then 1936 if (ntest.ge.150) 1937 & write(6,*) 'xdamp> search at a new place' 1938 ! search somewhere else 1939 dxdamp = + xfailfac*xinc 1940 xhigh = 1d13 1941 xvhigh = -100d0 1942 isubcnt = 0 1943 else 1944 ! get gradient and hessian from last three values 1945 ! gradient 1946 xgr = ((f1-f0)-(a1*a1)/(a2*a2)*(f2-f0))/(a1-(a1*a1)/a2) 1947 if (ntest.ge.150) 1948 & write(6,*) 'xdamp> num. gradient ',xgr 1949 ! hessian 1950 xhs = 2d0*((f1-f0)-a1/a2*(f2-f0))/(a1*a1-a1*a2) 1951 if (ntest.ge.150) 1952 & write(6,*) 'xdamp> num. hessian ',xhs 1953 ! discriminant 1954 xdis = (xgr*xgr)/(xhs*xhs) - 2d0*f0/xhs 1955 if (ntest.ge.150) 1956 & write(6,*) 'xdamp> discriminant ',xdis 1957 ! hessian positive? 1958 if (xhs.gt.0d0.and.xdis.gt.0d0) then 1959 if (ntest.ge.150) 1960 & write(6,*) 'xdamp> accepting solution ' 1961 ! take the lower solution 1962 dxdamp = -xgr/xhs - sqrt(xdis) 1963c dxdamp = -xgr/xhs + sqrt(xdis) 1964 ! we go on with the 3-point search 1965 isubcnt = 3 1966 else 1967 ! gradient negative (and to be trusted) 1968 if (xgr.lt.0d0.and.abs(a2).lt.xinc) then 1969 if (ntest.ge.150) 1970 & write(6,*) 'xdamp> follow gradient step' 1971 ! take only the gradient step 1972 dxdamp = - f0/xgr 1973 ! we go on with the 3-point search 1974 isubcnt = 3 1975 else 1976 ! xgr.gt.0 does not really mean something: 1977 ! we try to get a new gradient here 1978 if (ntest.ge.150) 1979 & write(6,*) 1980 & 'xdamp> prepare new gradient evaluation' 1981 if (f0.gt.0d0) then 1982 dxdamp = + xinc 1983 else 1984 dxdamp = - xinc 1985 end if 1986 isubcnt = 2 1987 end if 1988 end if ! hessian > 0? 1989 end if ! monotonic function? 1990 else 1991 write(6,*) 'unknown isubcnt = ',isubcnt 1992 stop 'optc_sbspja' 1993 end if 1994 1995c if (isubcnt.ne.0) then 1996c ! we use a simple step restriction 1997c xmxdxdmp2 = max(xmxdxdamp,xdamp-1d0/xdamp) 1998c dxdamp = sign(min(abs(dxdamp),xmxdxdmp2),dxdamp) 1999c end if 2000 ! unless isubcnt was reset to 0 2001 if (isubcnt.ne.0) then 2002 if (xdamp+dxdamp.ge.xhigh) then 2003 if (xval.eq.xvhigh) stop 'strange: xdamp (1)' 2004 xdamp = xdamp-xval*(xdamp-xhigh)/(xval-xvhigh) 2005 else if (xdamp+dxdamp.le.xlow) then 2006 if (xlow.eq.xdamp_min) xvlow=-xval 2007 if (xval.eq.xvlow) stop 'strange: xdamp (2)' 2008 xdamp = xdamp-xval*(xdamp-xlow)/(xval-xvlow) 2009 else 2010 xdamp = xdamp+dxdamp 2011 end if 2012 else 2013 xdamp = xdamp+dxdamp 2014 end if 2015 2016c if (isubcnt.eq.0) then 2017c xdamp = xdamp + dxdamp 2018c else 2019c ! we use a simple step restriction 2020c xmxdxdmp2 = max(xmxdxdamp,xdamp-1d0/xdamp) 2021c dxdamp = sign(min(abs(dxdamp),xmxdxdmp2),dxdamp) 2022c xdamp = max(0d0,xdamp + dxdamp) 2023c end if 2024 2025 if (ntest.ge.150) 2026 & write(6,*)'xdamp> next xdamp = ',xdamp 2027 2028c if (xdamp.eq.0d0) xdamp = 10d0 2029c xdamp = xdamp - 0.009d0 2030c if (xdamp.lt.-4d0) stop 'test' 2031 end if ! .not.converged 2032 2033 2034c* and analyze it 2035c if (.not.accept.and.navec_l.gt.1) then 2036c again = .true. 2037c else if (.not.accept) then 2038c navec_l = 0 2039c end if 2040 end if ! again 2041 2042 if (again) then 2043 nskipst = nskipst+1 2044 nskipgv = nskipgv+1 2045 ndel = ndel+1 2046 navec_l = navec_l-1 2047 end if 2048 2049 end do ! again 2050 2051 end do ! xdamp - optimization 2052 2053* return the actual last subspace dimension 2054 navec = navec_l 2055 2056* make an energy prediction (for variational cases) 2057c 2058c RHS : scrvec2(irhsoff:irhsoff-1+navec_l) 2059c c : scrvec 2060c 2061c call matvcb(work(khred),work(kscr1),work(kscr2), 2062c & imicdim,imicdim,0) 2063c de_pred = inprod(work(kgred),work(kscr1),imicdim) + 2064c & 0.5d0*inprod(work(kscr2),work(kscr1),imicdim) 2065c 2066 2067 return 2068 2069 end 2070*----------------------------------------------------------------------* 2071*----------------------------------------------------------------------* 2072 subroutine optc_sbspja_new_iter(itype,thrsh,nstdim,ngvdim, 2073 & navec,maxvec, 2074 & nadd,navec_last, 2075 & lugrvf,ludia,trrad,lu_pertstp,lu_newstp, 2076 & xdamp,xdamp_last, 2077 & lust_sbsp,lugv_sbsp, 2078 & umat,scr,vec1,vec2,iprint) 2079*----------------------------------------------------------------------* 2080* 2081* subspace jacobian: A = A_0(1-P) + A_ex P 2082* 2083* if xdamp is 0d0, we try the direct inversion procedure; if that 2084* step is too large, we go over to the iterative solver to find the 2085* optimal damping (where we end up directly, if xdamp was != 0d0) 2086* 2087*----------------------------------------------------------------------* 2088 implicit none 2089 2090* constants 2091 integer, parameter :: 2092 & ntest = 00 2093 2094* input/output 2095 integer, intent(in) :: 2096 & itype, ! type (s.above) 2097 & nstdim,ngvdim, ! # records on sbsp files 2098 & maxvec, ! max dimension of sbsp 2099 & nadd, ! number of new vectors (usually 1) 2100 & lugrvf,ludia, ! current gradient, diagonal Hess/Jac. 2101 & lu_pertstp, ! step from diagonal Hess/Jac. 2102 & lu_newstp, ! new step (output) 2103 & lust_sbsp,lugv_sbsp, ! sbsp files 2104 & iprint ! print level 2105 integer, intent(inout) :: 2106 & navec, ! current dimension of sbsp 2107 & navec_last ! previous sbsp dimension 2108 real(8), intent(in) :: 2109 & thrsh, ! thresh for accepting low-rank correction 2110 & trrad, ! trust radius for step 2111 & xdamp_last ! previous damping 2112 real(8), intent(out) :: 2113 & xdamp ! current damping 2114 real(8), intent(inout) :: 2115 & umat(*), ! low rank Hessian/Jacobian 2116 & scr(*), vec1(*), vec2(*) ! scratch 2117 2118* local O(N) scratch 2119 integer :: 2120 & kpiv(navec) 2121 real(8) :: 2122 & scrvec(navec), scrvec2(navec) 2123 2124* local 2125 logical :: 2126 & lincore, again, accept, converged, lin_dep 2127 integer :: 2128 & lblk, iprintl, navec_l, isym, job, 2129 & nskipst, nskipgv, nnew, nold, nold2, ndel, 2130 & ii, jj, iioff, iioff2, irhsoff, 2131 & lu_sbsp_un, lu_curerr, 2132 & luvec, luAvec, luvec_sbsp, luAvec_sbsp, lures, 2133 & nadd_l, nold_l, nskip_l, ludum, iter, maxiter, 2134 & kend, ksmat, ksred, kared, kevec, kscr, 2135 & itype_evp, itask(2), ntask, ndum 2136 real(8) :: 2137 & cond, xcorsum, xlam, xresnrm, fac, thrres, xnrm 2138 integer, external :: 2139 & iopen_nus 2140 real(8), external :: 2141 & inprdd 2142 2143 iprintl = max(iprint,ntest) 2144 lblk = -1 2145 lincore = .false. 2146 if (itype.eq.1) isym = 0 ! non-symmetric jacobian 2147 if (itype.eq.2) isym = 1 ! symmetric jacobian 2148 2149 if (iprintl.ge.5) then 2150 write(6,*) 'Subspace-Hessian/Jacobian' 2151 write(6,*) '=========================' 2152 if (itype.eq.1) then 2153 write(6,*) ' asymmetric update' 2154 else if (itype.eq.2) then 2155 write(6,*) ' symmetric update' 2156 end if 2157 end if 2158 2159 if (ntest.ge.10) then 2160 write(6,*) 'on entry in optc_sbspja:' 2161 write(6,*) 2162 & 'nstdim, ngvdim, navec, maxvec, nadd, navec_last: ', 2163 & nstdim, ngvdim, navec, maxvec, nadd, navec_last 2164 write(6,*) 2165 & 'lugrvf, lu_newstp, lu_pertstp, lust_sbsp, lugv_sbsp: ', 2166 & lugrvf, lu_newstp, lu_pertstp, lust_sbsp, lugv_sbsp 2167 end if 2168 2169 if (itype.lt.1.or.itype.gt.2) then 2170 write(6,*) 'SBSP hessian/jacobian: unknown itype = ',itype 2171 stop 'SBSPJA arguments' 2172 end if 2173 2174 if (navec.gt.maxvec) then 2175 write(6,*) 2176 & 'Panic 2: Hess./Jac. subspace dimensioning inconsistent!' 2177 stop 'SBSP dimensions' 2178 end if 2179 2180 nskipst = nstdim-navec 2181 nskipgv = ngvdim-navec 2182 ndel = navec_last+nadd-navec 2183 nold = navec_last-ndel 2184 2185 navec_last = navec ! save for next iteration 2186 2187 if (nold.lt.0.or.nskipst.lt.0.or.nskipgv.lt.0) then 2188 write(6,*) 'inconsistent dimensions in optc_sbspja!' 2189 write(6,*) 'nold , nskipst, nskipgv : ', nold, nskipst, nskipgv 2190 stop 'SBSPJA dimensions' 2191 end if 2192 2193 if (ntest.ge.10) then 2194 write(6,*) 'vectors to be deleted in U mat (ndel): ',ndel 2195 write(6,*) 'old vectors on U matrix: (nold) ',nold 2196 write(6,*) 'skipped vectors on step file:(nskipst)',nskipst 2197 write(6,*) 'skipped vectors on grad file:(nskipgv)',nskipgv 2198 end if 2199 2200* vectors deleted? move B matrix by the corresp. number of col.s and rows 2201 if (ndel.gt.0) then 2202 if (isym.eq.1) then ! symmetric U matrix 2203 ! nold is the future dimension of the already initialized 2204 ! low-rank matrix 2205 ! we now move all entries according to 2206 ! U(i,j) = U(i+ndel,j+ndel), where ij actually refer to a 2207 ! upper triangular matrix 2208 if (ntest.ge.100) then 2209 write(6,*) 'Low-rank matrix before deleting:' 2210 call prtrlt(umat,nold+ndel) 2211 end if 2212 2213 do ii = 1, nold 2214 iioff = ii*(ii-1)/2 2215 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 2216 do jj = 1, ii 2217 umat(iioff+jj) = umat(iioff2+jj) 2218 end do 2219 end do 2220 if (ntest.ge.100) then 2221 write(6,*) 'Low-rank matrix after deleting:' 2222 call prtrlt(umat,nold) 2223 end if 2224 else ! non-symmetric (and therefore full) matrix 2225 if (ntest.ge.100) then 2226 write(6,*) 'Low-rank matrix before deleting:' 2227 call wrtmat2(umat,nold+ndel,maxvec,nold+ndel,maxvec) 2228 end if 2229 2230 do ii = 1, nold 2231 iioff = (ii-1)*maxvec 2232 iioff2 = (ii+ndel-1)*maxvec + ndel 2233 do jj = 1, nold 2234 umat(iioff+jj) = umat(iioff2+jj) 2235 end do 2236 end do 2237 if (ntest.ge.100) then 2238 write(6,*) 'Low-rank matrix after deleting:' 2239 call wrtmat2(umat,nold,maxvec,nold,maxvec) 2240 end if 2241 2242 end if ! isym 2243 2244 end if ! ndel.gt.0 2245 2246* add new vectors 2247 ndel = 0 2248 nnew = nadd 2249 nold2 = nold 2250 call optc_sbspja_umat(itype, 2251 & umat,nnew,nold2,ndel,maxvec, 2252 & nskipst,nskipgv, 2253 & lust_sbsp,lugv_sbsp,ludia,xdamp, 2254 & vec1,vec2,lincore) 2255 2256* get projection of perturbation step in current subspace 2257 call optc_sbspja_prjpstp(scrvec2,navec,nskipst, 2258 & lu_pertstp,lust_sbsp, 2259 & vec1,vec2,lincore) 2260 2261 navec_l = navec 2262 2263 if (xdamp.eq.0d0) then 2264 2265* direct-inversion section: 2266* ========================= 2267 again = .true. 2268 do while(again) 2269* transfer to scratch array 2270 if (isym.eq.1) then 2271 do ii = 1, navec_l 2272 iioff = ii*(ii-1)/2 2273 iioff2 = (ii+ndel)*(ii+ndel-1)/2 + ndel 2274 do jj = 1, ii 2275 scr(iioff+jj) = umat(iioff2+jj) 2276 end do 2277 end do 2278 if (ntest.ge.100) then 2279 write(6,*) 'Low-rank matrix passed to dspco:' 2280 call prtrlt(scr,navec_l) 2281 end if 2282 else 2283 do ii = 1, navec_l 2284 iioff = (ii-1)*navec_l 2285 iioff2 = (ii+ndel-1)*maxvec + ndel 2286 do jj = 1, navec_l 2287 scr(iioff+jj) = umat(iioff2+jj) 2288 end do 2289 end do 2290 if (ntest.ge.100) then 2291 write(6,*) 'Low-rank matrix passed to dgeco:' 2292 call wrtmat2(scr,navec_l,navec_l,navec_l,navec_l) 2293 end if 2294 end if 2295 2296* invert low-rank matrix: 2297* factorize and get condition 2298 if (isym.eq.1) then 2299 call dspco(scr,navec_l,kpiv,cond,scrvec) 2300 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 2301 & navec_l,cond 2302 if (ntest.ge.100) then 2303 write(6,*) 'Factorized matrix from dspco:' 2304 call prtrlt(scr,navec_l) 2305 write(6,*) 'Pivot array:' 2306 call iwrtma(kpiv,1,navec_l,1,navec_l) 2307 end if 2308 else 2309 call dgeco(scr,navec_l,navec_l,kpiv,cond,scrvec) 2310 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 2311 & navec_l,cond 2312 if (ntest.ge.100) then 2313 write(6,*) 'Factorized matrix from dgeco:' 2314 call wrtmat2(scr,navec_l,navec_l,navec_l,navec_l) 2315 write(6,*) 'Pivot array:' 2316 call iwrtma(kpiv,1,navec_l,1,navec_l) 2317 end if 2318 end if ! isym 2319 2320 again = cond.lt.1d-14 2321 2322 if (.not.again) then 2323 irhsoff = navec - navec_l + 1 2324 scrvec(1:navec_l) = scrvec2(irhsoff:irhsoff-1+navec_l) 2325* if OK, solve U x = rhs 2326 if (ntest.ge.100) then 2327 write(6,*) 'RHS: ',scrvec(1:navec_l) 2328 end if 2329 if (isym.eq.1) then 2330 call dspsl(scr,navec_l,kpiv,scrvec) 2331 else 2332 job = 0 2333 call dgesl(scr,navec_l,navec_l,kpiv,scrvec,job) 2334 end if 2335 if (ntest.ge.100) then 2336 write(6,*) 'Result for x:' 2337 call wrtmat(scrvec,navec_l,1,navec_l,1) 2338 end if 2339 end if 2340 2341 if (.not.again) then 2342* get low-rank contribution to new vector 2343 call optc_sbspja_lrstep(itype,thrsh,accept,scrvec,navec_l, 2344 & nskipst,nskipgv, 2345 & lust_sbsp,lugv_sbsp,ludia,xdamp, 2346 & lu_pertstp,lu_newstp, 2347 & vec1,vec2,lincore) 2348 2349 xnrm = sqrt(inprdd(vec1,vec1,lu_newstp,lu_newstp,1,lblk)) 2350 accept = xnrm.le.trrad 2351 2352c* and analyze it 2353c if (.not.accept.and.navec_l.gt.1) then 2354c again = .true. 2355c else if (.not.accept) then 2356c navec_l = 0 2357c end if 2358 end if 2359 2360 if (again) then 2361 nskipst = nskipst+1 2362 nskipgv = nskipgv+1 2363 ndel = ndel+1 2364 navec_l = navec_l-1 2365 end if 2366 2367 end do 2368 2369 ! abuse xdamp as flag 2370 if (.not.accept) xdamp = 100d0 2371 2372* end of direct-inversion section 2373 end if 2374 2375* here starts the iterative inversion section with damping: 2376* ========================================================= 2377 if (xdamp.ne.0d0) then 2378 2379 kend = 1 2380 2381 ksmat = kend 2382 kend = kend + navec_l*(navec_l+2)/2 2383 2384 kscr = kend 2385 kend = kend + 27*(navec_l+2)**2 2386 2387 ksred = kend 2388 kend = kend + 27*(navec_l+2)**2 2389 2390 kared = kend 2391 kend = kend + 27*(navec_l+2)**2 2392 2393 kevec = kend 2394 kend = kend + 27*(navec_l+2)**2 2395 2396 luvec = iopen_nus('SSJA_VEC') 2397 luAvec = iopen_nus('SSJAAVEC') 2398 luvec_sbsp = iopen_nus('SSJA_V_SP') 2399 luAvec_sbsp = iopen_nus('SSJAAV_SP') 2400 lures = iopen_nus('SSJARVEC') 2401 2402 ! set up the S matrix 2403 nadd_l = navec_l 2404 nold_l = 0 2405 nskip_l = 0 2406 call optc_smat(scr(ksmat),navec_l,nold_l,nskip_l, 2407 & lust_sbsp,nskipst,nstdim, 2408 & ludum,0,0, 2409 & vec1,vec2,lincore) 2410 2411 ! decompose it 2412 call dspco(scr(ksmat),navec_l,kpiv,cond,scrvec) 2413 if (ntest.ge.10) write(6,*)'navec_l,condition number: ', 2414 & navec_l,cond 2415 if (ntest.ge.100) then 2416 write(6,*) 'Factorized matrix from dspco:' 2417 call prtrlt(scr(ksmat),navec_l) 2418 write(6,*) 'Pivot array:' 2419 call iwrtma(kpiv,1,navec_l,1,navec_l) 2420 end if 2421 ! delete vectors in subspace if S becomes singular 2422 ! to be done later ... 2423 2424 ! microiterations: 2425 converged = .false. 2426 iter = 0 2427 maxiter = 3*(navec_l+1) 2428 thrres = 1d-6 2429 2430 call sclvcd(lu_pertstp,luvec,-1d0,vec1,1,lblk) 2431 2432 do while(.not.converged) 2433 iter = iter+1 2434 ! orthonormalize trial vector 2435 call optc_orthvec(scr(kscr),scr(ksred), 2436 & iter-1,1,lin_dep, 2437 & luvec_sbsp,luvec, 2438 & vec1,vec2) 2439 2440 ! get Mv product 2441 call optc_sbspja_mvp(itype, 2442 & luvec,luAvec, 2443 & scr(ksmat),kpiv,navec_l, 2444 & lust_sbsp,nskipst,lugv_sbsp,nskipgv, 2445 & ludia, 2446 & vec1,vec2,lincore) 2447 2448 ! solve reduced EVP 2449 itype_evp = 0 ! trust-radius method 2450 !itype_evp = 1 ! augmented hessian method 2451 call optc_sbspja_redevp(itype_evp, 2452 & scr(kared),scr(ksred),scr(kscr),scr(kevec),trrad, 2453 & xlam,lures,xresnrm, 2454 & iter, 2455 & lugrvf, 2456 & luvec,luAvec,luvec_sbsp,luAvec_sbsp, 2457 & vec1,vec2,lincore) 2458 2459 ! test convergence 2460 if (abs(xresnrm).lt.thrres) converged = .true. 2461 2462 write(6,'(">>",i4,2(2x,e12.2),2x,l)') 2463 & iter,xlam,xresnrm,converged 2464 2465 if (iter.eq.maxiter) exit 2466 2467 if (.not.converged) then 2468 ! put luvec on subspace-file 2469 fac = 1d0 2470 itask(1) = 1 2471 itask(2) = 1 2472 ntask = 2 2473 call optc_sbspman(luvec_sbsp,luvec,fac,ludum,iter-1,maxiter, 2474 & itask,ntask,0,ndum,vec1,vec2) 2475 2476 ! put luAvec on subspace file 2477 call optc_sbspman(luAvec_sbsp,luAvec,fac, 2478 & ludum,iter-1,maxiter, 2479 & itask,ntask,0,ndum,vec1,vec2) 2480 2481 ! precondition new direction: 2482 call dmtvcd2(vec1,vec2,ludia,lures,luvec,1d0,0d0,1,1,lblk) 2483 2484 end if 2485 2486 end do 2487 2488 if (.not.converged) then 2489 write(6,*) 'WARNING: Microiterations not converged' 2490 write(6,*) ' final residual : ',xresnrm 2491 end if 2492 2493 ! assemble new step from first vector in scr(kevec) and 2494 ! the subspace vectors in luvec_sbsp and the last 2495 ! trial vector (which we didn't put on the subspace file) 2496 if (iter.gt.1) then 2497 ! luAvec and luAvec_sbsp used as scratch 2498 call mvcsmd(luvec_sbsp,scr(kevec),luAvec,luAvec_sbsp, 2499 & vec1,vec2,iter-1,1,lblk) 2500c TEST 2501c call vecsmd(vec1,vec2,-1d0,-scr(kevec-1+iter), 2502 call vecsmd(vec1,vec2,1d0,scr(kevec-1+iter), 2503 & luAvec,luvec,lu_newstp,1,lblk) 2504 else 2505c TEST 2506c call sclvcd(luvec,lu_newstp,-scr(kevec),vec1,1,lblk) 2507 call sclvcd(luvec,lu_newstp,scr(kevec),vec1,1,lblk) 2508 end if 2509 2510 call relunit(lures,'delete') 2511 call relunit(luvec,'delete') 2512 call relunit(luAvec,'delete') 2513 call relunit(luvec_sbsp,'delete') 2514 call relunit(luAvec_sbsp,'delete') 2515 2516 end if 2517 2518* return the actual last subspace dimension 2519 navec = navec_l 2520 2521 return 2522 2523 end 2524*----------------------------------------------------------------------* 2525*----------------------------------------------------------------------* 2526 subroutine optc_orthvec(scr1,scr2,ndim1,ndim2,lin_dep, 2527 & luvec1,luvec2, 2528 & vec1,vec2) 2529*----------------------------------------------------------------------* 2530* 2531* luvec1 contains a set of ndim1 orthonormal vectors, 2532* luvec2 a set of ndim2 non-orthonormal vectors. On output luvec2 2533* is orthonormalized against luvec1 and among itself. 2534* 2535*----------------------------------------------------------------------* 2536 implicit none 2537 2538 integer, parameter :: 2539 & ntest = 00 2540 2541 integer, intent(in) :: 2542 & ndim1, ndim2, 2543 & luvec1, luvec2 2544 logical, intent(out) :: 2545 & lin_dep(ndim2) 2546 real(8), intent(inout) :: 2547 & scr1(ndim1+ndim2,ndim1+ndim2), 2548 & scr2(ndim1+ndim2,ndim1+ndim2), 2549 & vec1(*), vec2(*) 2550 2551 integer :: 2552 & lblk, ii, jj, ndim, 2553 & lusc1, lusc2, lusc3, lusc4 2554 real(8) :: 2555 & xs, xnorm 2556 2557 real(8) :: 2558 & scrvec(ndim1+ndim2) 2559 2560 integer, external :: 2561 & iopen_nus 2562 real(8), external :: 2563 & inprdd, inprod 2564 2565 lblk = -1 2566 2567 lusc1 = iopen_nus('OPTC_ORTHSC1') 2568 lusc2 = iopen_nus('OPTC_ORTHSC2') 2569 lusc3 = iopen_nus('OPTC_ORTHSC3') 2570 lusc4 = iopen_nus('OPTC_ORTHSC4') 2571 2572 ! set up metric on s 2573 ndim = ndim1+ndim2 2574 scr1(1:ndim,1:ndim) = 0d0 2575 do ii = 1, ndim1 2576 scr1(ii,ii) = 1d0 2577 end do 2578 2579 call rewino(luvec2) 2580 do jj = ndim1+1, ndim 2581 call rewino(lusc1) 2582 call copvcd(luvec2,lusc1,vec1,0,lblk) 2583 call rewino(luvec1) 2584 do ii = 1, ndim1 2585 call rewino(lusc1) 2586 xs = inprdd(vec1,vec2,lusc1,luvec1,0,lblk) 2587 scr1(ii,jj) = xs 2588 scr1(jj,ii) = xs 2589 end do 2590 end do 2591 2592 ! get first element 2593 scr1(ndim1+1,ndim1+1) = inprdd(vec1,vec1,luvec2,luvec2,1,lblk) 2594 do ii = ndim1+2, ndim 2595 ! we are now at element ii, copy to scratch 2596 call copvcd(luvec2,lusc1,vec1,0,lblk) 2597 ! and rewind 2598 call rewino(luvec2) 2599 do jj = ndim1+2, ii 2600 call rewino(lusc1) 2601 xs = inprdd(vec1,vec2,luvec2,lusc1,0,lblk) 2602 scr1(ii,jj) = xs 2603 scr1(jj,ii) = xs 2604 end do 2605 end do 2606 2607 if (ntest.ge.100) then 2608 write(6,*) 'overlap matrix:' 2609 call wrtmat2(scr1,ndim,ndim,ndim,ndim) 2610 end if 2611 2612 ! call modified Gram-Schmidt orthonormalization 2613 call mgs3(scr2,scr1,ndim,scrvec) 2614 2615 if (ntest.ge.100) then 2616 write(6,*) 'orthonormalization matrix:' 2617 call wrtmat2(scr2,ndim,ndim,ndim,ndim) 2618 end if 2619 2620 do ii = ndim1+1, ndim 2621 xnorm = inprod(scr2(1,ii),scr2(1,ii),ndim) 2622 if (xnorm.lt.epsilon(1d0)) then 2623 if (ntest.ge.10) write(6,*) 2624 & ' linear dependency detected for vector ',ii 2625 lin_dep(ii-ndim1) = .true. 2626 ! leave this vector unmodified 2627 scr2(ii,ii) = 1d0 2628 else 2629 lin_dep(ii-ndim1) = .false. 2630 end if 2631 end do 2632 2633 ! update luvec2 2634 call rewino(luvec2) 2635 call rewino(lusc2) 2636 do ii = 1, ndim2 2637 call copvcd(luvec2,lusc2,vec1,0,lblk) 2638 end do 2639 2640 call rewino(luvec2) 2641 do ii = ndim1+1, ndim 2642 ! assemble new vector ii for luvec2 2643 ! contributions from luvec1 2644 if (ndim1.gt.0) 2645 & call mvcsmd(luvec1,scr2(1,ii),lusc1,lusc3, 2646 & vec1,vec2,ndim1,1,lblk) 2647 2648 ! contributions from luvec2 (now on lusc2) 2649 if (ndim2.gt.0) 2650 & call mvcsmd(lusc2,scr2(ndim1+1,ii),lusc4,lusc3, 2651 & vec1,vec2,ndim2,1,lblk) 2652 2653 if (ndim1.gt.0.and.ndim2.gt.0) then 2654 call rewino(lusc1) 2655 call rewino(lusc4) 2656 call vecsmd(vec1,vec2,1d0,1d0,lusc1,lusc4,luvec2,0,lblk) 2657 else if (ndim1.gt.0) then 2658 call rewino(lusc1) 2659 call copvcd(lusc1,luvec2,vec1,0,lblk) 2660 else if (ndim2.gt.0) then 2661 call rewino(lusc4) 2662 call copvcd(lusc4,luvec2,vec1,0,lblk) 2663 end if 2664 2665 end do 2666 2667 call relunit(lusc1,'delete') 2668 call relunit(lusc2,'delete') 2669 call relunit(lusc3,'delete') 2670 call relunit(lusc4,'delete') 2671 2672 return 2673 2674 end 2675*----------------------------------------------------------------------* 2676*----------------------------------------------------------------------* 2677 subroutine optc_sbspja_redevp(itype, 2678 & amat,smat,scr,xvec,xs, 2679 & xlam,lures,xresnrm, 2680 & ndim, 2681 & lugvec, 2682 & luvec,luAvec,luvec_sbsp,luAvec_sbsp, 2683 & vec1,vec2,lincore) 2684*----------------------------------------------------------------------* 2685* 2686* solve 2687* 2688* A_aug |x_aug> = 2689* 2690* / \ / \ / \ 2691* | <1|A|1> .. <1|A|n> <1|g> | | x1 | | x1 | 2692* | . . . | | . | | . | 2693* | . . . | | . | = lambda S | . | 2694* | <n|A|1> .. <n|A|n> <n|g> | | xn | | xn | 2695* | <g|1> .. <g|n> 0 | | 1 | | 1 | 2696* \ / \ / \ / 2697* 2698* (i.e. |x_aug> will be normalized such that its last coeff. is +1) 2699* 2700* where for itype==0 (trust radius method) 2701* 2702* / \ 2703* | 1 0 . 0 0 | 2704* | 0 . . | 2705* S = | . . 0 . | 2706* | 0 . 0 1 0 | 2707* | 0 ... 0 -s**2 | 2708* \ / 2709* 2710* or for itype==1 (rational function method) 2711* 2712* / \ 2713* | s 0 . 0 0 | 2714* | 0 . . | 2715* S = | . . 0 . | 2716* | 0 . 0 s 0 | 2717* | 0 ... 0 1 | 2718* \ / 2719* 2720* on input: on output: 2721* amat A_aug in (n-1)-dim subspc. A_aug in n-dim subspace 2722* smat scratch holding ---> S in n-dim subspace 2723* xvec eigenvectors of reduced EVP 2724* xs scalar defining the metric 2725* xlam lowest eigenvalue 2726* xresnrm residual norm in full space 2727* lures file: residual vector 2728* ndim current dimension of 2729* subspace 2730* lugvec file: g in full space 2731* luvec file: last trial vector |n> 2732* luAvec file: last A|n> 2733* luvec_sbsp file: (n-1)-dim subspace {|n>} 2734* luAvec_sbsp file: (n-1)-dim subspace {A|n>} 2735* vec1,vec2 scratch (size: at least largest block on vector file) 2736* lincore .true. if vec1/vec2 can hold one complete vector 2737*----------------------------------------------------------------------* 2738 implicit none 2739 2740 integer, parameter :: 2741 & ntest = 00 2742 2743* interface: 2744 logical, intent(in) :: 2745 & lincore 2746 integer, intent(in) :: 2747 & itype, ndim, 2748 & lures, 2749 & lugvec,luvec,luAvec,luvec_sbsp,luAvec_sbsp 2750 real(8), intent(in) :: 2751 & xs 2752 real(8), intent(out) :: 2753 & xlam, xresnrm 2754 real(8), intent(inout) :: 2755 & amat(*), smat(ndim+1,*), scr(*), 2756 & xvec(ndim+1,*), vec1(*), vec2(*) 2757 2758* local 2759 integer :: 2760 & lblk, ndimold, ndimp1, ndimp1old, matz, ierr, 2761 & iad, iadold, iminev, 2762 & ii, jj, iiold, jjold, 2763 & lusc1, lusc2 2764 2765 real(8) :: 2766 & xx, xminev 2767 2768* local order(N)-scratch, allocated on the stack: 2769 real(8) :: 2770 & evr(ndim+1),evi(ndim+1),evd(ndim+1) 2771 2772* external 2773 integer, external :: 2774 & iopen_nus 2775 real(8), external :: 2776 & inprdd, inprod 2777 2778c TEST 2779 real(8) xtest 2780 2781 lblk = -1 2782 lusc1 = iopen_nus('OPTCEVPSC1') 2783 lusc2 = iopen_nus('OPTCEVPSC2') 2784 2785 ndimp1 = ndim+1 2786 ! set up projected metric matrix 2787 smat(1:ndim+1,1:ndim+1) = 0d0 2788 if (itype.eq.0) then 2789 xx = 1d0 2790 else 2791 xx = 1d0/xs 2792 end if 2793 do ii = 1, ndim 2794 smat(ii,ii) = xx 2795 end do 2796 if (itype.eq.0) then 2797 smat(ndim+1,ndim+1) = 1d0!-xs*xs 2798 else 2799 smat(ndim+1,ndim+1) = 1d0 2800 end if 2801 2802 if (ntest.ge.100) then 2803 write(6,*) 'The S matrix: ' 2804 call wrtmat(smat,ndimp1,ndimp1,ndimp1,ndimp1) 2805 end if 2806 2807 ! update A matrix 2808 ! the ndim+1,ndim+1 element is zero 2809 amat((ndimp1-1)*ndimp1+ndimp1) = 0d0 2810 ! first rearrange from (ndim+1)-1 to (ndim+1) 2811 ! a) the g column 2812 ndimold = ndim-1 2813 ndimp1old = ndimp1-1 2814 if (ntest.ge.100.and.ndimp1old.gt.1) then 2815 write(6,*) 'The previous augmented A matrix: ' 2816 call wrtmat(amat,ndimp1old,ndimp1old,ndimp1old,ndimp1old) 2817 end if 2818 jjold = ndimp1old 2819 jj = ndimp1 2820 do iiold = ndimold, 1, -1 2821 iadold = ((jjold-1)*ndimp1old)+iiold 2822 iad = ((jj-1)*ndimp1)+iiold 2823 amat(iad) = amat(iadold) 2824 end do 2825 2826 do jjold = ndimold, 1, -1 2827 ! b) one element of the g row 2828 iadold = ((jjold-1)*ndimp1old)+ndimp1old 2829 iad = ((jjold-1)*ndimp1old)+ndimp1 2830 amat(iad) = amat(iadold) 2831 2832 ! c) the elements of A 2833 do iiold = ndimold, 1, -1 2834 iadold = ((jjold-1)*ndimp1old)+iiold 2835 iad = ((jjold-1)*ndimp1)+iiold 2836 amat(iad) = amat(iadold) 2837 end do 2838 2839 end do 2840 2841 if (ntest.ge.100) then 2842 write(6,*) 'The previous augm. A matrix after rearranging: ' 2843 call wrtmat(amat,ndimp1,ndimp1,ndimp1,ndimp1) 2844 end if 2845 2846 ! update the new row in A_aug 2847 call rewino(luvec_sbsp) 2848 do ii = 1, ndim-1 2849 call rewino(luAvec) 2850 amat(((ndim-1)*ndimp1)+ii) 2851 & = inprdd(vec1,vec2,luvec_sbsp,luAvec,0,lblk) 2852 end do 2853 2854 ! update the new row in A_aug 2855 call rewino(luAvec_sbsp) 2856 do jj = 1, ndim-1 2857 call rewino(luvec) 2858 amat(((jj-1)*ndimp1)+ndim) 2859 & = inprdd(vec1,vec2,luvec,luAvec_sbsp,0,lblk) 2860 end do 2861 2862 ! the (ndim,ndim)-element: 2863 amat((ndim-1)*ndimp1+ndim) = inprdd(vec1,vec2,luvec,luAvec,1,lblk) 2864 2865 ! update with <g|x_n> 2866 xx = inprdd(vec1,vec2,lugvec,luvec,1,lblk) 2867c TEST 2868 amat((ndimp1-1)*ndimp1+ndim) = 0d0!xx 2869 amat((ndim-1)*ndimp1+ndimp1) = 0d0!xx 2870 2871 if (ntest.ge.100) then 2872 write(6,*) 'The new augmented A matrix: ' 2873 call wrtmat(amat,ndimp1,ndimp1,ndimp1,ndimp1) 2874 end if 2875 2876 ! solve reduced EVP 2877 matz = 1 2878 ! transfer amat to scratch array (as rgg destroys it) 2879 scr(1:ndimp1**2) = amat(1:ndimp1**2) 2880 ! call the general EVP-solver from eispack: 2881c TEST 2882 call rgg(ndimp1,ndim,scr,smat,evr,evi,evd, 2883c call rgg(ndimp1,ndimp1,scr,smat,evr,evi,evd, 2884 & matz,xvec,ierr) 2885 if (ierr.ne.0) then 2886 write(6,*) 'error code from rgg: ',ierr 2887 stop 'optc_sbspja_redevp' 2888 end if 2889 2890 xminev = 1d100 2891 iminev = 0 2892c TEST 2893 do ii = 1, ndim 2894c do ii = 1, ndimp1 2895 evr(ii) = evr(ii)/evd(ii) 2896 if (evr(ii).lt.xminev) then 2897 xminev = evr(ii) 2898 iminev = ii 2899 end if 2900 evi(ii) = evi(ii)/evd(ii) 2901 end do 2902 if (iminev.eq.0.or.iminev.gt.ndimp1) then 2903 write(6,*) 'internal error' 2904 stop 'optc_sbspja_redevp' 2905 end if 2906 2907 if (ntest.ge.100) then 2908 write(6,*) 'The eigenvalues from RGG:' 2909 do ii = 1, ndimp1 2910 write(6,'(x,i4,2f12.6)') ii,evr(ii),evi(ii) 2911 end do 2912 write(6,*) 'The eigenvector with lowest real component: ', 2913 & evr(iminev) 2914 call wrtmat(xvec(1,iminev),ndimp1,1,ndimp1,1) 2915 if (evi(1).ne.0d0) then 2916 write(6,*) 'imaginary component:' 2917 call wrtmat(xvec(2,iminev),ndimp1,1,ndimp1,1) 2918 end if 2919 end if 2920 2921 xlam = evr(iminev) 2922 if (evi(iminev).ne.0d0) then 2923 write(6,*) 'Huah! Imaginary lowest eigenvalue encountered!' 2924 write(6,*) evr(iminev),' +/- i * ',evi(iminev) 2925 write(6,*) 'Don''t know what to do and give up ...' 2926 stop 'optc_sbspja_redevp' 2927 end if 2928 ! renormalize 2929c TEST 2930 xx = inprod(xvec(1,iminev),xvec(1,iminev),ndim) 2931c xx = xvec(ndimp1,iminev) 2932 do ii = 1, ndimp1 2933 xvec(ii,iminev) = xvec(ii,iminev)/xx 2934 end do 2935 2936 if (ntest.ge.100) then 2937 write(6,*) 'The lowest eigenvector (renormalized):' 2938 call wrtmat(xvec(1,iminev),ndimp1,1,ndimp1,1) 2939 end if 2940 2941 ! get residual in full space 2942 if (itype.eq.0) then 2943 xx = 1d0 2944 else 2945 xx = 1d0/xs 2946 end if 2947 2948 if (ndim.gt.1) then 2949 ! (A - lambda S) x ... 2950 call mvcsmd(luAvec_sbsp,xvec(1,iminev),lusc1,lusc2, 2951 & vec1,vec2,ndim-1,1,lblk) 2952 call vecsmd(vec1,vec2,1d0,xvec(ndim,iminev), 2953 & lusc1,luAvec,lusc2,1,lblk) 2954 ! lures used as scratch 2955 call mvcsmd(luvec_sbsp,xvec(1,iminev),lusc1,lures, 2956 & vec1,vec2,ndim-1,1,lblk) 2957 2958 call vecsmd(vec1,vec2,-xlam*xx,1d0, 2959 & lusc1,lusc2,lures,1,lblk) 2960 call vecsmd(vec1,vec2,1d0,-xlam*xx*xvec(ndim,iminev), 2961 & lures,luvec,lusc2,1,lblk) 2962 else 2963 ! (A - lambda S) x ... 2964 call vecsmd(vec1,vec2, 2965 & xvec(ndim,iminev),-xlam*xx*xvec(ndim,iminev), 2966 & luAvec,luvec,lusc2,1,lblk) 2967 end if 2968 2969 call copvcd(lusc2,lures,vec1,1,lblk) 2970 2971 xresnrm = sqrt(inprdd(vec1,vec1,lures,lures,1,lblk)) 2972 2973 call relunit(lusc1,'delete') 2974 call relunit(lusc2,'delete') 2975 2976 return 2977 end 2978*----------------------------------------------------------------------* 2979*----------------------------------------------------------------------* 2980 subroutine optc_sbspja_mvp(itype, 2981 & luvec,luAvec, 2982 & smat,kpiv,ndim, 2983 & lust_sbsp,nskipst,lujt_sbsp,nskipjt, 2984 & ludia, 2985 & vec1,vec2,lincore) 2986*----------------------------------------------------------------------* 2987* 2988*----------------------------------------------------------------------* 2989 implicit none 2990 2991 integer, parameter :: 2992 & ntest = 00 2993 2994 logical, intent(in) :: 2995 & lincore 2996 2997 integer, intent(in) :: 2998 & itype, luvec, luAvec, 2999 & ndim, lust_sbsp, lujt_sbsp, 3000 & nskipst, nskipjt, ludia 3001 3002 integer, intent(in) :: 3003 & kpiv(ndim) 3004 3005 real(8), intent(in) :: 3006 & smat(*) 3007 3008 real(8), intent(inout) :: 3009 & vec1(*), vec2(*) 3010 3011 integer :: 3012 & lusc1, lusc2, lblk, 3013 & ii, jj 3014 3015 real(8) :: 3016 & v(ndim) 3017 3018 real(8), external :: 3019 & inprdd 3020 integer, external :: 3021 & iopen_nus 3022 3023 lblk = -1 3024 3025 lusc1 = iopen_nus('MATVSCR1') 3026 lusc2 = iopen_nus('MATVSCR2') 3027 3028 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3029 ! v(j) = <x_j|v> 3030 do jj = 1, ndim 3031 call rewino(luvec) 3032 v(jj) = inprdd(vec1,vec2,luvec,lust_sbsp,0,lblk) 3033 end do 3034 3035 if (ntest.ge.100) then 3036 write(6,*) 'v:' 3037 call wrtmat(v,ndim,1,ndim,1) 3038 write(6,*) 'smat:' 3039 call prtrlt(smat,ndim) 3040 end if 3041 3042 ! S_{ij}^{-1}v_j 3043 call dspsl(smat,ndim,kpiv,v) 3044 3045 if (ntest.ge.100) then 3046 write(6,*) 'vbar:' 3047 call wrtmat(v,ndim,1,ndim,1) 3048 end if 3049 3050 call skpvcd(lujt_sbsp,nskipjt,vec1,1,lblk) 3051 call rewino(luAvec) 3052 call rewino(lusc2) 3053 ! v(i)*|Jx_i> --> luAvec 3054 call mvcsmd(lujt_sbsp,v,luAvec,lusc2,vec1,vec2, 3055 & ndim,0,lblk) 3056 3057 if (ntest.ge.1000) then 3058 write(6,*) 'v(i)|Jx_i>:' 3059 call wrtvcd(vec1,luAvec,1,lblk) 3060 end if 3061 3062 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3063 do ii = 1, ndim 3064 ! D|x_i> --> lusc2 3065 call rewino(ludia) 3066 call rewino(lusc2) 3067 call dmtvcd2(vec1,vec2,ludia,lust_sbsp,lusc2,1d0,0d0,0,0,lblk) 3068 ! add -v(i)D|x_i> to luAvec 3069 call vecsmd(vec1,vec2,1d0,-v(ii),luAvec,lusc2,lusc1,1,lblk) 3070 call copvcd(lusc1,luAvec,vec1,1,lblk) 3071 end do 3072 3073 if (ntest.ge.1000) then 3074 write(6,*) 'v(i)(|Jx_i>-D|x_i>):' 3075 call wrtvcd(vec1,luAvec,1,lblk) 3076 end if 3077 3078 ! D|v> 3079 call dmtvcd2(vec1,vec2,ludia,luvec,lusc2,1d0,0d0,1,0,lblk) 3080 ! add D|v> to luAvec 3081 call vecsmd(vec1,vec2,1d0,1d0,luAvec,lusc2,lusc1,1,lblk) 3082 call copvcd(lusc1,luAvec,vec1,1,lblk) 3083 3084 if (ntest.ge.1000) then 3085 write(6,*) 'D|v> + v(i)(|Jx_i>-D|x_i>):' 3086 call wrtvcd(vec1,luAvec,1,lblk) 3087 end if 3088 3089 call relunit(lusc1,'delete') 3090 call relunit(lusc2,'delete') 3091 3092 return 3093 3094 end 3095*----------------------------------------------------------------------* 3096*----------------------------------------------------------------------* 3097 subroutine optc_sbspja_umat(itype, 3098 & umat,nnew,nold,nskip_u,maxvec, 3099 & nskipst,nskipjt, 3100 & lust_sbsp,lujt_sbsp,ludia,xdamp, 3101 & vec1,vec2,lincore) 3102*----------------------------------------------------------------------* 3103* 3104* update the low-rank hessian/jacobian 3105* 3106* itype = 1: U_ij = <dt^(i)|(D+damp)^-1|Adt^(j)> 3107* 3108* itype = 2: dto. but symmetrized 3109* 3110* |dt^(i)> is on lust_sbsp 3111* |Adt^(i)> is on lujt_sbsp (usually gradient differences) 3112* D is on ludia 3113* 3114* if (xdamp.ne.0d0) the complete matrix needs to be rebuild. we 3115* expect the overlap matrix in packed triangular form on input. 3116* 3117*----------------------------------------------------------------------* 3118 implicit none 3119 3120* constants 3121 integer, parameter :: 3122 & ntest = 00 3123 logical, parameter :: 3124 & ldo_extra = .true. 3125 3126* input/output 3127 logical, intent(in) :: 3128 & lincore 3129 integer, intent(in) :: 3130 & itype, 3131 & nnew, nold, nskip_u, maxvec, 3132 & nskipst, nskipjt, 3133 & lust_sbsp, lujt_sbsp, ludia 3134 real(8), intent(in) :: 3135 & xdamp 3136 real(8), intent(inout) :: 3137 & umat(*), 3138 & vec1(*), vec2(*) 3139 3140* local 3141 integer :: 3142 & ii, jj, ijdx, isym, lblk, luscr 3143 real(8) :: 3144 & uij, fac 3145 3146* external 3147 integer, external :: 3148 & iopen_nus 3149 real(8), external :: 3150 & inprdd, inprdd3 3151 3152 lblk = -1 3153 3154 if (itype.lt.1.or.itype.gt.2) then 3155 write(6,*) 'optc_sbspja_umat: not prepared for itype = ',itype 3156 stop 'optc_sbspja_umat' 3157 end if 3158 3159 luscr = iopen_nus('UMATSCR') 3160 3161 if (itype.eq.1) isym = 0 3162 if (itype.eq.2) isym = 1 3163 3164 if (itype.eq.1) fac = 1d0 3165 if (itype.eq.2) fac = .5d0 3166 3167 if (xdamp.ne.0d0.and.ldo_extra) then 3168 if (nold.ne.0) then 3169 write(6,*) 'error: nold must be 0 if xdamp.ne.0d0!' 3170 stop 'optc_sbspja_umat' 3171 end if 3172 if (ntest.ge.100) then 3173 write(6,*) 'on entry: S packed' 3174 call prtrlt(umat,nnew) 3175 end if 3176 if (isym.eq.0) then 3177 call uptripak(umat,umat,2,nnew,maxvec) 3178 if (ntest.ge.100) then 3179 write(6,*) 'S unpacked' 3180 call wrtmat(umat,nnew,nnew,maxvec,maxvec) 3181 end if 3182 end if 3183 else 3184 3185 ! zero the new elements: 3186 if (isym.eq.0) then 3187 do jj = nskip_u+nold+1, nskip_u+nold+nnew 3188 do ii = nskip_u+1, nskip_u+nold+nnew 3189 ijdx = (jj-1)*maxvec + ii 3190 umat(ijdx) = 0d0 3191 end do 3192 end do 3193 do ii = nskip_u+nold+1, nskip_u+nold+nnew 3194 do jj = nskip_u+1, nskip_u+nold 3195 ijdx = (jj-1)*maxvec + ii 3196 umat(ijdx) = 0d0 3197 end do 3198 end do 3199 else 3200 do jj = nskip_u+nold+1, nskip_u+nold+nnew 3201 do ii = nskip_u+1, jj 3202 ijdx = jj*(jj-1)/2 + ii 3203 umat(ijdx) = 0d0 3204 end do 3205 end do 3206 end if 3207 3208 end if ! if (xdamp.ne.0d0) 3209 3210c a) update columns jj 3211c first let�s do the damping terms (to see if they are significant) 3212c u(i,j) -= <dt^i|(A_0+xdamp)^-1 A_0 |dt^j> 3213 if (xdamp.ne.0d0.and.ldo_extra) then 3214 do jj = nskip_u+nold+1, nskip_u+nold+nnew 3215 call rewino(ludia) 3216 call rewino(luscr) 3217 ! well, not ideal ... : 3218 call skpvcd(lust_sbsp,nskipst-1+jj,vec1,1,lblk) 3219 call dmtvcd(vec1,vec2,ludia,lust_sbsp,luscr,0d0,0,0,lblk) 3220 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3221 do ii = nskip_u+1, jj 3222 call rewino(luscr) 3223 call rewino(ludia) 3224 uij = inprdd3(vec1,vec2,lust_sbsp,luscr,ludia, 3225 & xdamp,-1d0,0,lblk) 3226 if (isym.eq.0) then 3227 ijdx = (jj-1)*maxvec + ii 3228 umat(ijdx) = umat(ijdx) - uij 3229 if (ii.ne.jj) then 3230 ijdx = (ii-1)*maxvec + jj 3231 umat(ijdx) = umat(ijdx) - uij 3232 end if 3233 else 3234 ijdx = jj*(jj-1)/2 + ii 3235 umat(ijdx) = umat(ijdx) - uij 3236 end if 3237 end do 3238 end do 3239 end if 3240 3241 if (ntest.ge.100.and.xdamp.ne.0d0.and.ldo_extra) then 3242 write(6,*) 'S_ij-<i|A_0(l)^-1 A_0|j> contrib. to U matrix: ' 3243 if (isym.eq.1) then 3244 call prtrlt(umat,nold+nnew) 3245 else 3246 call wrtmat(umat,nold+nnew,nold+nnew,maxvec,maxvec) 3247 end if 3248 end if 3249 3250 call skpvcd(lujt_sbsp,nskipjt+nold,vec1,1,lblk) 3251c and now the important contribution: 3252c u(i,j) += <dt^i|A_0^-1|Adt^j> 3253 do jj = nskip_u+nold+1, nskip_u+nold+nnew 3254 call rewino(ludia) 3255 call rewino(luscr) 3256 call dmtvcd(vec1,vec2,ludia,lujt_sbsp,luscr,xdamp,0,1,lblk) 3257 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3258 do ii = nskip_u+1, nskip_u+nold+nnew 3259 call rewino(luscr) 3260 uij = inprdd(vec1,vec2,lust_sbsp,luscr,0,lblk) 3261 if (isym.eq.0) then 3262 ijdx = (jj-1)*maxvec + ii 3263 else 3264 if (ii.le.jj) then 3265 ijdx = jj*(jj-1)/2 + ii 3266 else 3267 ijdx = ii*(ii-1)/2 + jj 3268 end if 3269 end if 3270 umat(ijdx) = umat(ijdx) + fac*uij 3271 if (isym.eq.1.and.ii.eq.jj) umat(ijdx) = umat(ijdx) + fac*uij 3272 end do 3273 end do 3274* b) update rows ii 3275* (as we assume nold.eq.0 for xdamp.ne.0d0, here nothing is 3276* to be done for this case) 3277 if (nold.gt.0) then 3278 call skpvcd(lust_sbsp,nskipst+nold,vec1,1,lblk) 3279 do ii = nskip_u+nold+1, nskip_u+nold+nnew 3280 call rewino(ludia) 3281 call rewino(luscr) 3282 call dmtvcd(vec1,vec2,ludia,lust_sbsp,luscr,xdamp,0,1,lblk) 3283 call skpvcd(lujt_sbsp,nskipjt,vec1,1,lblk) 3284 do jj = nskip_u+1, nskip_u+nold ! do not doubly visit the elements 3285 ! we already updated in a) 3286 call rewino(luscr) 3287 uij = inprdd(vec1,vec2,lujt_sbsp,luscr,0,lblk) 3288 if (isym.eq.0) then 3289 ijdx = (jj-1)*maxvec + ii 3290 else 3291 if (ii.le.jj) then 3292 ijdx = jj*(jj-1)/2 + ii 3293 else 3294 ijdx = ii*(ii-1)/2 + jj 3295 end if 3296 end if 3297 umat(ijdx) = umat(ijdx) + fac*uij 3298 end do 3299 end do 3300 3301 end if ! nold.gt.0 3302 3303 if (ntest.ge.100) then 3304 write(6,*) 'Updated U matrix: ',isym 3305 if (isym.eq.1) then 3306 call prtrlt(umat,nold+nnew) 3307 else 3308 call wrtmat(umat,nold+nnew,nold+nnew,maxvec,maxvec) 3309 end if 3310 end if 3311 call relunit(luscr,'delete') 3312 3313 return 3314 end 3315*----------------------------------------------------------------------* 3316*----------------------------------------------------------------------* 3317 subroutine optc_sbspja_prjpstp(vecout,navec,nskipst, 3318 & lu_pertstp,lust_sbsp, 3319 & vec1,vec2,lincore) 3320*----------------------------------------------------------------------* 3321* 3322* get projection of vector on lu_pertstp in subspace lust_sbsp 3323* 3324*----------------------------------------------------------------------* 3325 implicit none 3326 3327 integer, parameter :: 3328 & ntest = 00 3329 3330 logical, intent(in) :: 3331 & lincore 3332 integer, intent(in) :: 3333 & navec, nskipst, 3334 & lu_pertstp,lust_sbsp 3335 real(8), intent(out) :: 3336 & vecout(navec) 3337 real(8), intent(inout) :: 3338 & vec1(*), vec2(*) 3339 3340 integer :: 3341 & ii, lblk 3342 3343 real(8), external :: 3344 & inprdd 3345 3346 lblk = -1 3347 3348 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3349 do ii = 1, navec 3350 call rewino(lu_pertstp) 3351 vecout(ii) = inprdd(vec1,vec2,lust_sbsp,lu_pertstp,0,lblk) 3352 end do 3353 3354 if (ntest.ge.100) then 3355 write(6,*) 'projected RHS:' 3356 call wrtmat(vecout,navec,1,navec,1) 3357 end if 3358 3359 return 3360 end 3361*----------------------------------------------------------------------* 3362*----------------------------------------------------------------------* 3363 subroutine optc_sbspja_lrstep(itype,thrsh,accept,vecin,navec, 3364 & nskipst,nskipjt, 3365 & lust_sbsp,lujt_sbsp,ludia,xdamp, 3366 & lu_pertstp,lu_newstp, 3367 & vec1,vec2,lincore) 3368*----------------------------------------------------------------------* 3369 implicit none 3370 3371 integer, parameter :: 3372 & ntest = 00 3373 3374 logical, intent(in) :: 3375 & lincore 3376 logical, intent(out) :: 3377 & accept 3378 integer, intent(in) :: 3379 & itype, navec, nskipst, nskipjt, 3380 & lust_sbsp, lujt_sbsp, ludia, 3381 & lu_pertstp, lu_newstp 3382 real(8), intent(in) :: 3383 & vecin(navec), xdamp, thrsh 3384 real(8), intent(inout) :: 3385 & vec1(*), vec2(*) 3386 3387 integer :: 3388 & lblk, luscr, luscr2 3389 real(8) :: 3390 & xnrm 3391 3392 integer, external :: 3393 & iopen_nus 3394 real(8), external :: 3395 & inprdd 3396 3397 lblk = -1 3398 3399 if (ntest.ge.100) then 3400 write(6,*) 'assembling low-rank step' 3401 write(6,*) '------------------------' 3402 write(6,*) ' navec, nskipst, nskipjt: ',navec, nskipst, nskipjt 3403 write(6,*) ' lust_sbsp, lujt_sbsp, ludia: ', 3404 & lust_sbsp, lujt_sbsp, ludia 3405 write(6,*) ' lu_pertstp, lu_newstp: ',lu_pertstp, lu_newstp 3406 end if 3407 3408 luscr = iopen_nus('LRCSCR1') 3409 luscr2 = iopen_nus('LRCSCR2') 3410 3411 call rewino(luscr) 3412 call rewino(luscr2) 3413 call skpvcd(lujt_sbsp,nskipjt,vec1,1,lblk) 3414 ! Omg = vec(i) * Omg(i) --> luscr 3415 call mvcsmd(lujt_sbsp,vecin,luscr,luscr2,vec1,vec2, 3416 & navec,0,lblk) 3417 ! D^-1 Omg --> luscr2 3418 call dmtvcd(vec1,vec2,ludia,luscr,luscr2,xdamp,1,1,lblk) 3419 3420 call rewino(luscr) 3421 call rewino(luscr2) 3422 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3423 ! dt = vec(i) * dt(i) --> lu_newstp (only intermediate scratch) 3424 call mvcsmd(lust_sbsp,vecin,lu_newstp,luscr,vec1,vec2, 3425 & navec,0,lblk) 3426* in the case of damping, we have to modify the result with A_0/(A_0+xdamp) 3427 if (xdamp.ne.0d0) then 3428 call dmtvcd(vec1,vec2,ludia,lu_newstp,luscr,xdamp,1,1,lblk) 3429 call dmtvcd(vec1,vec2,ludia,luscr,lu_newstp,0d0,1,0,lblk) 3430 end if 3431 ! assemble total low-rank correction on luscr 3432 call vecsmd(vec1,vec2,1d0,-1d0,lu_newstp,luscr2,luscr,1,lblk) 3433 3434 xnrm = sqrt(inprdd(vec1,vec1,luscr,luscr,1,lblk)) 3435 3436 if (ntest.ge.10) write(6,*) '|low-rank correction:| =',xnrm 3437 if (xnrm.gt.thrsh) then 3438 if (ntest.ge.10) write(6,*) 'low-rank correction not accepted!' 3439 call copvcd(lu_pertstp,lu_newstp,vec1,1,lblk) 3440 accept = .false. 3441 else 3442 call vecsmd(vec1,vec2,1d0,1d0,lu_pertstp,luscr,lu_newstp,1,lblk) 3443 accept = .true. 3444 end if 3445 3446 call relunit(luscr,'delete') 3447 call relunit(luscr2,'delete') 3448 3449 return 3450 end 3451*----------------------------------------------------------------------* 3452*----------------------------------------------------------------------* 3453 subroutine optc_smat(smat,ndim,nold,nskips, 3454 & lu_sbsp1,nskip1,nvec1, 3455 & lu_sbsp2,nskip2,nvec2, 3456 & vec1,vec2,lincore) 3457*----------------------------------------------------------------------* 3458* update or build S matrix = <i|j> as upper triangular matrix 3459* ndim the new S dimension 3460* nold vectors already reside in the smat array 3461* nskips vectors for that the new <i|new> contributions are not 3462* needed 3463* lu_sbsp1 the first file containing (usually) the old vectors 3464* of the previously initialized subspace 3465* nskip1 number of vectors to be skipped on that file 3466* nvec1 total (!) number of vectors on that file 3467* lu_sbsp2 the second file containing (usually) the new vector(s) 3468* to be added to the subspace 3469* nskip2 same as nskip1 but for second file 3470* nvec2 dto. 3471*----------------------------------------------------------------------* 3472 implicit none 3473 3474 integer, parameter :: 3475 & ntest = 00 3476 3477 logical, intent(in) :: 3478 & lincore 3479 integer, intent(in) :: 3480 & ndim,nold,nskips, 3481 & lu_sbsp1,lu_sbsp2,nskip1,nskip2,nvec1,nvec2 3482 real(8), intent(inout) :: 3483 & smat(*),vec1(*),vec2(*) 3484 3485 integer :: 3486 & ii, jj, idxii, lblk, lusc1, 3487 & nadd, nadd1, nadd2, nvecs 3488 3489 real(8), external :: 3490 & inprdd 3491 integer, external :: 3492 & iopen_nus 3493 3494 if (ntest.ge.100) then 3495 write(6,*) 'Welcome to optc_smat:' 3496 write(6,*) '=====================' 3497 write(6,*) ' ndim,nold,nskips: ',ndim,nold,nskips 3498 write(6,*) ' lu_sbsp1,nskip1,nvec1: ',lu_sbsp1,nskip1,nvec1 3499 write(6,*) ' lu_sbsp2,nskip2,nvec2: ',lu_sbsp2,nskip2,nvec2 3500 write(6,*) ' infos on units:' 3501 if (nvec1.gt.0) then 3502 print *,'lu_sbsp1:' 3503 call unit_info(lu_sbsp1) 3504 end if 3505 if (nvec2.gt.0) then 3506 print *,'lu_sbsp2:' 3507 call unit_info(lu_sbsp2) 3508 end if 3509 end if 3510 3511 lblk = -1 3512 ! consistency check: 3513 nvecs = nvec1+nvec2 - (nskip1+nskip2) 3514 if (nvecs.ne.ndim) then 3515 write(6,*) 'consistency error in optc_smat: ' 3516 write(6,*) ' nvec1, nvec2, sum1 : ',nvec1, nvec2, nvec1+nvec2 3517 write(6,*) ' nskip1, nskip2, sum2 : ', 3518 & nskip1, nskip2, nskip1+nskip2 3519 write(6,*) ' sum1 - sum2, ndim : ',nvecs,ndim 3520 stop 'optc_smat' 3521 end if 3522 3523 lusc1 = iopen_nus('SMATSCR') 3524 3525 if (ntest.ge.100) then 3526 write(6,*) 'S on input:' 3527 call prtrlt(smat,ndim) 3528 end if 3529 3530 ! number of rows to be added 3531 nadd = ndim - nold 3532 nadd2 = nvec2-nskip2 3533 nadd1 = nadd-nadd2 3534 3535 ! contributions within lu_sbsp1 3536 do ii = nold+1, nold+nadd1 3537 ! goto the appropriate vector on lu_sbsp1 3538 call skpvcd(lu_sbsp1,nskip1+ii-1,vec1,1,lblk) 3539 ! get a copy 3540 call rewino(lusc1) 3541 call copvcd(lu_sbsp1,lusc1,vec1,0,lblk) 3542 ! rewind 3543 call skpvcd(lu_sbsp1,nskip1,vec1,1,lblk) 3544 idxii = (ii-1)*ii/2 3545 do jj = 1+nskips, ii 3546 call rewino(lusc1) 3547 smat(idxii+jj) = inprdd(vec1,vec2,lusc1,lu_sbsp1,0,lblk) 3548 end do 3549 end do 3550 3551 if (ntest.ge.100) then 3552 write(6,*) 'S after 1st block:' 3553 call prtrlt(smat,ndim) 3554 end if 3555 3556 ! contributions from lu_sbsp1/lu_sbsp2 3557 if (nadd2.gt.0) 3558 & call skpvcd(lu_sbsp2,nskip2,vec1,1,lblk) 3559 do ii = nold+nadd1+1, nold+nadd1+nadd2 3560 call rewino(lusc1) 3561 call copvcd(lu_sbsp2,lusc1,vec1,0,lblk) 3562 call skpvcd(lu_sbsp1,nskip1,vec1,1,lblk) 3563 idxii = (ii-1)*ii/2 3564 do jj = 1+nskips, nold+nadd1 3565 call rewino(lusc1) 3566 smat(idxii+jj) = inprdd(vec1,vec2,lusc1,lu_sbsp1,0,lblk) 3567 end do 3568 end do 3569 3570 if (ntest.ge.100) then 3571 write(6,*) 'S after 2nd block:' 3572 call prtrlt(smat,ndim) 3573 end if 3574 3575 ! contributions within lu_sbsp2 3576 do ii = nold+nadd1+1, nold+nadd1+nadd2 3577 ! go to the appropriate vector on lu_sbsp2 3578 call skpvcd(lu_sbsp2,nskip1+ii-1,vec1,1,lblk) 3579 ! get a copy 3580 call rewino(lusc1) 3581 call copvcd(lu_sbsp2,lusc1,vec1,0,lblk) 3582 ! rewind 3583 call skpvcd(lu_sbsp2,nskip1,vec1,1,lblk) 3584 idxii = (ii-1)*ii/2 3585 do jj = nold+nadd1+1, ii 3586 call rewino(lusc1) 3587 smat(idxii+jj) = inprdd(vec1,vec2,lusc1,lu_sbsp2,0,lblk) 3588 end do 3589 end do 3590 3591 if (ntest.ge.100) then 3592 write(6,*) 'S after 3rd block:' 3593 call prtrlt(smat,ndim) 3594 end if 3595 3596 call relunit(lusc1,'delete') 3597 3598 return 3599 3600 end 3601*----------------------------------------------------------------------* 3602*----------------------------------------------------------------------* 3603 subroutine optc_updtja(itype,nrank,thrsh, 3604 & nstdim,nhgdim, 3605 & navec,maxvec, 3606 & nadd,navec_last, 3607 & lugrvf,lugrvf_last, 3608 & ludia,trrad,lu_pertstp,lu_newstp, 3609 & luhg_last,luhgam_new, 3610 & xdamp,xdamp_last, 3611 & lust_sbsp,luhg_sbsp, 3612 & hkern,vec1,vec2,iprint) 3613*----------------------------------------------------------------------* 3614* 3615* calculate the matrix-vector product of the current gradient passed 3616* in with an updated inverse Hessian matrix H^-1(k): 3617* 3618* |w> = H^-1(k)|v> = H_0^-1|v> + sum(i=2,k) E_i |v> 3619* 3620* where the matrices E_i are obtained according to the update 3621* formulae of 3622* 3623* 11-14: Broyden-Family: 3624* 11 Broyden-Fletcher-Goldfarb-Shanno 3625* 12 Davidon-Fletcher-Powell 3626* 13 Broyden Rank 1 3627* 14 Hoshino 3628* 3629* 15 Broyden assymetric rank 1 3630* 3631* |v> comes in on lugrvf 3632* |w> goes into the wide world through unit lu_newstp 3633* 3634* |delta_i> previous steps on lust_sbsp 3635* |H gamma_i> previous grad. diffs. times prev. H on luhg_sbsp 3636* for Powell: |gamma - H^{-1}delta>, instead 3637* H_0 = H_1 on ludia 3638* H_0^-1|g> is assumed to be on lu_pertstp 3639* H_(k-1)^-1|gamma_(k-1)> will leave on luhgnew 3640* for Powell: |gamma - H_(k-1)^{-1}delta> 3641* 3642* H_(k-1)^-1|g_(k-1)> comes from previous iteration thru luhg_last 3643* H_(k)^-1|g_(k)> goes to next iteration thru luhg_last 3644* ( not used for Powell ) 3645* 3646*----------------------------------------------------------------------* 3647 implicit none 3648 3649* constants 3650 integer, parameter :: 3651 & ntest = 00 3652 3653* input/output 3654 integer, intent(in) :: 3655 & itype, ! type (s.above) 3656 & nrank, ! rank of update 3657 & nstdim,nhgdim, ! # records on sbsp files 3658 & maxvec, ! max dimension of sbsp 3659 & nadd, ! number of new vectors (usually 1) 3660 & lugrvf,ludia, ! current gradient, diagonal Hess/Jac. 3661 & lugrvf_last, ! previous gradient 3662 & lu_pertstp, ! step from diagonal Hess/Jac. 3663 & lu_newstp,luhgam_new, ! new step (output) 3664 & luhg_last, ! H_(k-1)^-1|g_(k-1)> from last iteration 3665 & lust_sbsp, ! sbsp files 3666 & luhg_sbsp, 3667 & iprint ! print level 3668 integer, intent(inout) :: 3669 & navec, ! current dimension of sbsp 3670 & navec_last ! previous sbsp dimension 3671 real(8), intent(in) :: 3672 & thrsh, ! thresh for accepting low-rank correction 3673 & trrad, ! trust radius for step 3674 & xdamp_last ! previous damping 3675 real(8), intent(inout) :: 3676 & xdamp ! current damping 3677 real(8), intent(inout) :: 3678 & hkern(nrank*(nrank+1)/2,navec), 3679 & ! low-rank kernels of previous updates 3680 & vec1(*), vec2(*) ! scratch 3681 3682* local O(N) scratch 3683c integer :: 3684c & kpiv(navec) 3685 real(8) :: 3686 & v1(navec), v2(navec), 3687 & hv1(navec), hv2(navec) 3688 3689* local 3690c logical :: 3691c & lincore, again, accept 3692 integer :: 3693 & ii, iprintl, ndel1, ndel2, nn, 3694 & nskipst, nskiphg, lblk, 3695 & luscr1, luscr2, luscr3, ludlcnt, lupnt1, lupnt2 3696 real(8) :: 3697 & x4, x2, x3, x1, phi, 3698 & f1, f2, fac 3699 3700 integer, external :: 3701 & iopen_nus 3702 real(8), external :: 3703 & inprdd 3704 3705 luscr1 = iopen_nus('RUHSCR1') 3706 luscr2 = iopen_nus('RUHSCR2') 3707 luscr3 = iopen_nus('RUHSCR3') 3708 3709 nskipst = nstdim - navec 3710 nskiphg = nhgdim - (navec-1) 3711 3712 iprintl = max(iprint,ntest) 3713 lblk = -1 3714 3715 if (itype.lt.11.or.itype.gt.15) then 3716 write(6,*) 'unknown update method!' 3717 stop 'optc_updtja' 3718 end if 3719 3720 if ((itype.ge.11.and.itype.le.14.and.nrank.ne.2).or. 3721 & (itype.eq.15 .and.nrank.ne.1)) then 3722 write(6,*) 'illegal combination of type and rank: ',itype,nrank 3723 stop 'optc_updtja' 3724 end if 3725 3726 if (iprintl.ge.5) then 3727 write(6,*) 'Updated-Hessian/Jacobian' 3728 write(6,*) '========================' 3729 if (itype.eq.11) then 3730 write(6,*) ' BFGS update' 3731 else if (itype.eq.12) then 3732 write(6,*) ' DFP update' 3733 else if (itype.eq.13) then 3734 write(6,*) ' Broyden rank 1 update' 3735 else if (itype.eq.14) then 3736 write(6,*) ' Hoshino update' 3737 else if (itype.eq.15) then 3738 write(6,*) ' Broyden asymmetric update' 3739 end if 3740 end if 3741 3742 if (ntest.ge.10) then 3743 write(6,*) 'on entry in optc_sbspja:' 3744 write(6,*) 3745 & 'nstdim, nhgdim, navec, maxvec, nadd, navec_last: ', 3746 & nstdim, nhgdim, navec, maxvec, nadd, navec_last 3747 write(6,*) 3748 & 'lugrvf, lugrvf_last, lu_newstp, lu_pertstp: ', 3749 & lugrvf, lugrvf_last, lu_newstp, lu_pertstp 3750 write(6,*) 3751 & 'luhgam_new, luhg_last: ',luhgam_new, luhg_last 3752 write(6,*) 3753 & 'lust_sbsp, luhg_sbsp: ', lust_sbsp, luhg_sbsp 3754 end if 3755 3756 if (ntest.ge.10) then 3757 write(6,*) 'nskipst, nskiphg: ', 3758 & nskipst, nskiphg 3759 end if 3760 3761 if (ntest.ge.100) then 3762 write(6,*) 'previous low-rank matrices:' 3763 do ii = 1, navec-1 3764 write(6,*) '(',ii,')' 3765 call prtrlt(hkern(1,ii),nrank) 3766 end do 3767 end if 3768 3769 ndel1 = 0 3770 ndel2 = 0 3771 if (navec.eq.maxvec) then 3772 ! be prepared to remove the first vector from the subspace 3773 ndel1 = 1 ! for Hinv_{k}g_{k} contributions 3774 ndel2 = 1 ! for moving h^{k} around 3775c if (itype.eq.15) ndel1 = 0 ! not necessary 3776 if (ndel1.ne.0) ludlcnt = iopen_nus('RUHDLCNT') 3777 end if 3778 3779* 1) calculate 3780* 3781* Hinv_{k-1} gamma_{k-1} = Hinv_{k-1} g_{k} - Hinv_{k-1} g_{k-1} 3782* 3783* where 3784* 3785* Hinv_{k-1} g_(k) = Hinv_{1}g_k 3786* 3787* + sum_{i=2,k-1} (delta_{i-1}, Hgamma_{i-1}) x 3788* 3789* / h11^(i) h12^(i) \ / <delta_{i-1}|g_(k)>\ 3790* | | | | 3791* \ h12^(i) h22^(i) / \<Hgamma_{i-1}|g_(k)>/ 3792* 3793* the low-rank kernel h^(i) is on hkern(1..3,i) 3794* -Hinv_{1}g_k is on lu_pertstp 3795* -Hinv_{k-1}g_{k-1} is on luhg_last 3796* 3797 ! get this file into position 3798 if (navec.gt.0) then 3799 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3800 end if 3801 3802 if (navec.gt.1.and.itype.ge.11.and.itype.le.14) then 3803 3804 ! <delta_{i}|g_(k)> 3805 ! this has been done outside the 'if': 3806 ! call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3807 do ii = 1, navec-1 3808 call rewino(lugrvf) 3809 v1(ii) = inprdd(vec1,vec2,lust_sbsp,lugrvf,0,lblk) 3810 end do 3811 3812 ! <Hinv gamma_{i}|g_(k)> 3813 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 3814 do ii = 1, navec-1 3815 call rewino(lugrvf) 3816 v2(ii) = inprdd(vec1,vec2,luhg_sbsp,lugrvf,0,lblk) 3817 end do 3818 3819 if (ntest.ge.100) then 3820 write(6,*) '<delta_{i}|g_(k)>:' 3821 call wrtmat(v1,navec-1,1,navec-1,1) 3822 write(6,*) '<Hgamma_{i}|g_(k)>:' 3823 call wrtmat(v2,navec-1,1,navec-1,1) 3824 end if 3825 3826 do ii = 1, navec-1 3827 hv1(ii) = hkern(1,ii)*v1(ii) + hkern(2,ii)*v2(ii) 3828 hv2(ii) = hkern(2,ii)*v1(ii) + hkern(3,ii)*v2(ii) 3829 end do 3830 3831 if (ntest.ge.100) then 3832 write(6,*) 'hv1:' 3833 call wrtmat(hv1,navec-1,1,navec-1,1) 3834 write(6,*) 'hv2:' 3835 call wrtmat(hv2,navec-1,1,navec-1,1) 3836 end if 3837 3838 ! bring subspace-files into position: 3839 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 3840 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 3841 if (ndel1.gt.0) then 3842 ! store contributions from 'to-be-deleted' vectors on 3843 ! separate file ludlcnt 3844 call mvcsmd(lust_sbsp,hv1(1),luscr1,luscr2, 3845 & vec1,vec2,ndel1,0,lblk) 3846 call mvcsmd(luhg_sbsp,hv2(1),luscr2,luscr3, 3847 & vec1,vec2,ndel1,0,lblk) 3848 call vecsmd(vec1,vec2,1d0,1d0,luscr1,luscr2,ludlcnt,1,lblk) 3849 end if 3850 3851 ! assemble delta contributions 3852 call rewino(luscr1) 3853 call rewino(luscr2) 3854 call mvcsmd(lust_sbsp,hv1(1+ndel1),luscr1,luscr2, 3855 & vec1,vec2,navec-1-ndel1,0,lblk) 3856c call mvcsmd(lust_sbsp,hv1,luscr1,luscr2, 3857c & vec1,vec2,navec-1,0,lblk) 3858 ! add to Hinv_1g_k on lu_pertstp --> luscr3 3859 call vecsmd(vec1,vec2,1d0,-1d0,lu_pertstp,luscr1,luscr3,1,lblk) 3860 3861 ! assemble Hinv gamma contributions 3862 call rewino(luscr1) 3863 call rewino(luscr2) 3864 call mvcsmd(luhg_sbsp,hv2(1+ndel1),luscr1,luscr2, 3865 & vec1,vec2,navec-1-ndel1,0,lblk) 3866c call mvcsmd(luhg_sbsp,hv2,luscr1,luscr2, 3867c & vec1,vec2,navec-1,0,lblk) 3868 ! add to luscr3 --> lu_newstp (used as intermediate scratch) 3869 if (ndel1.eq.0) then 3870 call vecsmd(vec1,vec2,1d0,-1d0,luscr3,luscr1,lu_newstp,1,lblk) 3871 else 3872 ! add contribution from 'to-be-deleted' vectors here 3873 call vecsmd(vec1,vec2,1d0,1d0,luscr1,ludlcnt,luscr2,1,lblk) 3874 call vecsmd(vec1,vec2,1d0,-1d0,luscr3,luscr2,lu_newstp,1,lblk) 3875 end if 3876 3877c if (itype.ge.11.and.itype.le.14) then 3878c ! BFGS and DFP: 3879c ! so far, we have built -Hinv_{k-1}g_k on lu_newstp 3880c ! combine with -Hinv_{k-1}g_{k-1} on luhg_last to get 3881c ! Hinv_{k-1}gamma_{k-1} --> luhgam_new 3882c call vecsmd(vec1,vec2,-1d0,1d0, 3883c & lu_newstp,luhg_last,luhgam_new,1,lblk) 3884c end if 3885 3886 else if (navec.gt.1.and.itype.eq.15) then 3887* 3888* assymetric Broyden: we have to recursively update Hinv g 3889* 3890 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 3891 ! point to lu_newstp as starting point 3892 lupnt1 = luscr2 3893 lupnt2 = lu_pertstp 3894 ! initial prefactor of -1 as lu_newstp carries -|Hinv0 g> 3895 fac = -1d0 3896 do ii = 1, navec-1 3897 ! <delta^{ii}|Hinv{ii}g> 3898 call rewino(lupnt2) 3899 x1 = fac*inprdd(vec1,vec2,lust_sbsp,lupnt2,0,lblk) 3900 fac = -1d0 3901 ! multiply with rank-1 kernel 3902 x2 = hkern(1,ii)*x1 3903 3904 if (ntest.ge.100) then 3905 write(6,*) ii,'/',navec-1, ' x1 = ',x1,' x2 = ',x2 3906 end if 3907 3908 if (ii.eq.navec-1) lupnt1 = lu_newstp 3909 3910 call rewino(lupnt1) 3911 call rewino(lupnt2) 3912 call vecsmd(vec1,vec2,1d0,-x2,lupnt2,luhg_sbsp,lupnt1,0,lblk) 3913 if (mod(ii,2).eq.1) then 3914 lupnt1 = luscr1 3915 lupnt2 = luscr2 3916 else 3917 lupnt1 = luscr2 3918 lupnt2 = luscr1 3919 end if 3920 3921 end do 3922 ! now we should have -Hinv^{k-1}g^{k} on lu_newstp (as before) 3923c ! combine with old -Hinv^{k-1}g^{k-1} 3924c call vecsmd(vec1,vec2,-1d0,1d0, 3925c & lu_newstp,luhg_last,luhgam_new,1,lblk) 3926 3927 else ! navec.gt.1 3928 3929 ! first two iterations: just copy 3930 call copvcd(lu_pertstp,lu_newstp,vec1,1,lblk) 3931 end if 3932 3933 if (navec.gt.0) then 3934 ! so far, we have built -Hinv_{k-1}g_k on lu_newstp 3935 ! combine with -Hinv_{k-1}g_{k-1} on luhg_last to get 3936 ! Hinv_{k-1}gamma_{k-1} --> luhgam_new 3937 call vecsmd(vec1,vec2,-1d0,1d0, 3938 & lu_newstp,luhg_last,luhgam_new,1,lblk) 3939 3940 else 3941 ! else Hinv_{k-1}gamma{k-1} is identical with Hinv_1 g_{k} 3942 call sclvcd(lu_pertstp,luhgam_new,-1d0,vec1,1,lblk) 3943 end if 3944* 3945* what is expected at this position: 3946* -Hinv^{k-1}g^{k} on lu_newstp 3947* lust_sbsp positioned on delta^{k-1} 3948* 3949 if (navec.gt.0) then 3950 ! we take advantage of the file lust_sbsp being 3951 ! the correct postion for this 3952 ! delta^{k-1} --> luscr2 3953 call rewino(luscr2) 3954 call copvcd(lust_sbsp,luscr2,vec1,0,lblk) 3955 3956 end if 3957 3958* 3959* 2) calculate new h^{k} 3960* 3961* we expect that delta^(k-1) is still save(d) on luscr2 3962* 3963 if (navec.gt.0.and.(itype.ge.11.and.itype.le.14)) 3964 & then 3965 3966 ! calculate the current gamma from current and previous gradient 3967 ! --> luscr3 3968 call vecsmd(vec1,vec2,1d0,-1d0,lugrvf,lugrvf_last,luscr3,1,lblk) 3969 3970 ! <gamma|Hgamma> 3971 x1 = inprdd(vec1,vec2,luscr3,luhgam_new,1,lblk) 3972 ! <delta^(k-1)|gamma^{k-1}> 3973 x2 = inprdd(vec1,vec2,luscr2,luscr3,1,lblk) 3974 ! <delta^(k-1)|g^{k}> 3975 x3 = inprdd(vec1,vec2,luscr2,lugrvf,1,lblk) 3976 ! <Hgamma|g^{k}> 3977 x4 = inprdd(vec1,vec2,luhgam_new,lugrvf,1,lblk) 3978 3979 if (itype.eq.11) then 3980 phi = 1d0 ! BFGS 3981 else if (itype.eq.12) then 3982 phi = 0d0 ! DFP 3983 else if (itype.eq.13) then 3984 phi = 1d0/(1d0 - x1/x2) ! Broyden rank 1 3985 else if (itype.eq.14) then 3986 phi = 1d0/(1d0 + x1/x2) ! Hoshino 3987 end if 3988 3989 hkern(1,navec) = (1d0+phi*x1/x2)/x2 3990 hkern(2,navec) = -phi/x2 3991 hkern(3,navec) = (phi-1d0)/x1 3992 3993 if (ntest.ge.100) then 3994 write(6,*) '<gamma^{k-1}|Hgamma^{k-1}> = ',x1 3995 write(6,*) ' <delta^(k-1)|gamma^{k-1}> = ',x2 3996 write(6,*) ' <delta^(k-1)|g^{k}> = ',x3 3997 write(6,*) ' <Hgamma^{k-1}|g^{k}> = ',x4 3998 end if 3999 else if (navec.gt.0.and.itype.eq.15) then 4000 ! <delta^{k-1}|Hgamma^{k-1}> 4001 x1 = inprdd(vec1,vec2,luscr2,luhgam_new,1,lblk) 4002 ! |delta^{k-1}> - |Hgamma^{k-1}> --> luhgam_new 4003 call vecsmd(vec1,vec2,1d0,-1d0, 4004 & luscr2,luhgam_new,luscr3,1,lblk) 4005 call copvcd(luscr3,luhgam_new,vec1,1,lblk) 4006 4007 ! <delta^{k-1}|Hg^{k}> 4008 x3 = -inprdd(vec1,vec2,luscr2,lu_newstp,1,lblk) 4009 4010 if (ntest.ge.100) then 4011 print *,'<delta^{k-1}|Hgamma^{k-1}> = ',x1 4012 print *,'<delta^{k-1}|Hg^{k}> = ',x3 4013 end if 4014 4015 hkern(1,navec) = 1d0/x1 4016 end if 4017 4018 if (navec.gt.0.and.ntest.ge.100) then 4019 write(6,*) 'new rank n kernel (rank = ',nrank,'):' 4020 call prtrlt(hkern(1,navec),nrank) 4021 end if 4022* 4023* 3) complete Hinv^{k} g^{k} 4024* 4025 if (navec.gt.0.and.nrank.eq.2) then 4026 f1 = hkern(1,navec) * x3 + hkern(2,navec) * x4 4027 f2 = hkern(2,navec) * x3 + hkern(3,navec) * x4 4028 4029 if (ntest.ge.100) then 4030 write(6,*) 'f1 = ', f1 4031 write(6,*) 'f2 = ', f2 4032 end if 4033 4034 ! add to what we already have on lu_newstp 4035 call vecsmd(vec1,vec2,1d0,-f1,lu_newstp,luscr2,luscr3,1,lblk) 4036 call vecsmd(vec1,vec2,1d0,-f2, 4037 & luscr3,luhgam_new,lu_newstp,1,lblk) 4038 else if (navec.gt.0.and.nrank.eq.1) then 4039 f1 = hkern(1,navec)*x3 4040 4041 if (ntest.ge.100) then 4042 write(6,*) 'f1 = ', f1 4043 end if 4044 4045 call vecsmd(vec1,vec2,1d0,-f1, 4046 & lu_newstp,luhgam_new,luscr3,1,lblk) 4047 call copvcd(luscr3,lu_newstp,vec1,1,lblk) 4048 4049 end if ! navec.gt.0 4050 4051 ! save Hinv_k|g_k> for next iteration (only BFGS, DFP) 4052 if (itype.ge.11.and.itype.le.15) then 4053 if (ndel1.eq.0) then 4054 call copvcd(lu_newstp,luhg_last,vec1,1,lblk) 4055 else 4056 ! remove contributions from 'to-be-deleted' vectors 4057 ! such that in the next iteration we really calculate 4058 ! H^{k-1}gamma{k-1} = H^{k-1}g^k - H^{k-1}g^{k-1} 4059 ! as this ------------^^^ won't have these vectors anymore 4060 ! yes, the signs are ok: 4061 call vecsmd(vec1,vec2,1d0,1d0, 4062 & lu_newstp,ludlcnt,luhg_last,1,lblk) 4063 end if 4064 end if 4065 4066 ! if necessary, move h matrices down by ndel2 4067 if (ndel2.gt.0) then 4068 nn = nrank*(nrank+1)/2 4069 do ii = ndel2+1, navec 4070 hkern(1:nn,ii-ndel2) = hkern(1:nn,ii) 4071 end do 4072 end if 4073 4074 ! remove scratch files 4075 call relunit(luscr1,'delete') 4076 call relunit(luscr2,'delete') 4077 call relunit(luscr3,'delete') 4078 if (ndel1.gt.0) call relunit(ludlcnt,'delete') 4079 4080 return 4081 4082 end 4083*----------------------------------------------------------------------* 4084* the version with wrong Powell: 4085*----------------------------------------------------------------------* 4086 subroutine optc_updtja_old(itype,nrank,thrsh, 4087 & nstdim,nhgdim, 4088 & navec,maxvec, 4089 & nadd,navec_last, 4090 & lugrvf,lugrvf_last, 4091 & ludia,trrad,lu_pertstp,lu_newstp, 4092 & luhg_last,luhgam_new, 4093 & xdamp,xdamp_last, 4094 & lust_sbsp,luhg_sbsp, 4095 & hkern,vec1,vec2,iprint) 4096*----------------------------------------------------------------------* 4097* 4098* calculate the matrix-vector product of the current gradient passed 4099* in with an updated inverse Hessian matrix H^-1(k): 4100* 4101* |w> = H^-1(k)|v> = H_0^-1|v> + sum(i=2,k) E_i |v> 4102* 4103* where the matrices E_i are obtained according to the update 4104* formulae of 4105* 4106* 11 Broyden-Fletcher-Goldfarb-Shanno 4107* 12 Davidon-Fletcher-Powell 4108* 13 Powell symmetric Broyden 4109* 4110* |v> comes in on lugrvf 4111* |w> goes into the wide world through unit lu_newstp 4112* 4113* |delta_i> previous steps on lust_sbsp 4114* |H gamma_i> previous grad. diffs. times prev. H on luhg_sbsp 4115* for Powell: |gamma - H^{-1}delta>, instead 4116* H_0 = H_1 on ludia 4117* H_0^-1|g> is assumed to be on lu_pertstp 4118* H_(k-1)^-1|gamma_(k-1)> will leave on luhgnew 4119* for Powell: |gamma - H_(k-1)^{-1}delta> 4120* 4121* H_(k-1)^-1|g_(k-1)> comes from previous iteration thru luhg_last 4122* H_(k)^-1|g_(k)> goes to next iteration thru luhg_last 4123* ( not used for Powell ) 4124* 4125*----------------------------------------------------------------------* 4126 implicit none 4127 4128* constants 4129 integer, parameter :: 4130 & ntest = 00 4131 4132* input/output 4133 integer, intent(in) :: 4134 & itype, ! type (s.above) 4135 & nrank, ! rank of update 4136 & nstdim,nhgdim, ! # records on sbsp files 4137 & maxvec, ! max dimension of sbsp 4138 & nadd, ! number of new vectors (usually 1) 4139 & lugrvf,ludia, ! current gradient, diagonal Hess/Jac. 4140 & lugrvf_last, ! previous gradient 4141 & lu_pertstp, ! step from diagonal Hess/Jac. 4142 & lu_newstp,luhgam_new, ! new step (output) 4143 & luhg_last, ! H_(k-1)^-1|g_(k-1)> from last iteration 4144 & lust_sbsp, ! sbsp files 4145 & luhg_sbsp, 4146 & iprint ! print level 4147 integer, intent(inout) :: 4148 & navec, ! current dimension of sbsp 4149 & navec_last ! previous sbsp dimension 4150 real(8), intent(in) :: 4151 & thrsh, ! thresh for accepting low-rank correction 4152 & trrad, ! trust radius for step 4153 & xdamp_last ! previous damping 4154 real(8), intent(inout) :: 4155 & xdamp ! current damping 4156 real(8), intent(inout) :: 4157 & hkern(nrank*(nrank+1)/2,navec), 4158 & ! low-rank kernels of previous updates 4159 & vec1(*), vec2(*) ! scratch 4160 4161* local O(N) scratch 4162c integer :: 4163c & kpiv(navec) 4164 real(8) :: 4165 & v1(navec), v2(navec), 4166 & hv1(navec), hv2(navec) 4167 4168* local 4169c logical :: 4170c & lincore, again, accept 4171 integer :: 4172 & ii, iprintl, ndel1, ndel2, nn, 4173 & nskipst, nskiphg, lblk, 4174 & luscr1, luscr2, luscr3, ludlcnt 4175 real(8) :: 4176 & x4, x2, x3, x1, phi, 4177 & f1, f2 4178 4179 integer, external :: 4180 & iopen_nus 4181 real(8), external :: 4182 & inprdd 4183 4184 luscr1 = iopen_nus('RUHSCR1') 4185 luscr2 = iopen_nus('RUHSCR2') 4186 luscr3 = iopen_nus('RUHSCR3') 4187 4188 nskipst = nstdim - navec 4189 nskiphg = nhgdim - (navec-1) 4190 4191 iprintl = max(iprint,ntest) 4192 lblk = -1 4193 4194 if (iprintl.ge.5) then 4195 write(6,*) 'Updated-Hessian/Jacobian' 4196 write(6,*) '========================' 4197 if (itype.eq.11) then 4198 write(6,*) ' BFGS update' 4199 else if (itype.eq.12) then 4200 write(6,*) ' DFP update' 4201 else if (itype.eq.13) then 4202 write(6,*) ' Broyden rank 1 update' 4203 else if (itype.eq.14) then 4204 write(6,*) ' Hoshino update' 4205 else if (itype.eq.15) then 4206 write(6,*) ' PSB update' 4207 end if 4208 end if 4209 4210 if (ntest.ge.10) then 4211 write(6,*) 'on entry in optc_sbspja:' 4212 write(6,*) 4213 & 'nstdim, nhgdim, navec, maxvec, nadd, navec_last: ', 4214 & nstdim, nhgdim, navec, maxvec, nadd, navec_last 4215 write(6,*) 4216 & 'lugrvf, lugrvf_last, lu_newstp, lu_pertstp: ', 4217 & lugrvf, lugrvf_last, lu_newstp, lu_pertstp 4218 write(6,*) 4219 & 'luhgam_new, luhg_last: ',luhgam_new, luhg_last 4220 write(6,*) 4221 & 'lust_sbsp, luhg_sbsp: ', lust_sbsp, luhg_sbsp 4222 end if 4223 4224 if (ntest.ge.10) then 4225 write(6,*) 'nskipst, nskiphg: ', 4226 & nskipst, nskiphg 4227 end if 4228 4229 if (ntest.ge.100) then 4230 write(6,*) 'previous low-rank matrices:' 4231 do ii = 1, navec-1 4232 write(6,*) '(',ii,')' 4233 call prtrlt(hkern(1,ii),nrank) 4234 end do 4235 end if 4236 4237 ndel1 = 0 4238 ndel2 = 0 4239 if (navec.eq.maxvec) then 4240 ! be prepared to remove the first vector from the subspace 4241 ndel1 = 1 ! for Hinv_{k}g_{k} contributions 4242 ndel2 = 1 ! for moving h^{k} around 4243 if (itype.eq.15) ndel1 = 0 ! not necessary 4244 if (ndel1.ne.0) ludlcnt = iopen_nus('RUHDLCNT') 4245 end if 4246 4247* 1) calculate 4248* 4249* Hinv_{k-1} gamma_{k-1} = Hinv_{k-1} g_{k} - Hinv_{k-1} g_{k-1} 4250* 4251* where 4252* 4253* Hinv_{k-1} g_(k) = Hinv_{1}g_k 4254* 4255* + sum_{i=2,k-1} (delta_{i-1}, Hgamma_{i-1}) x 4256* 4257* / h11^(i) h12^(i) \ / <delta_{i-1}|g_(k)>\ 4258* | | | | 4259* \ h12^(i) h22^(i) / \<Hgamma_{i-1}|g_(k)>/ 4260* 4261* the low-rank kernel h^(i) is on hkern(1..3,i) 4262* -Hinv_{1}g_k is on lu_pertstp 4263* -Hinv_{k-1}g_{k-1} is on luhg_last 4264* 4265 ! get this file into position 4266 if (navec.gt.0) then 4267 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 4268 end if 4269 4270 if (navec.gt.1) then 4271 4272 ! <delta_{i}|g_(k)> 4273 ! this has been done outside the 'if': 4274 ! call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 4275 do ii = 1, navec-1 4276 call rewino(lugrvf) 4277 v1(ii) = inprdd(vec1,vec2,lust_sbsp,lugrvf,0,lblk) 4278 end do 4279 4280 ! <Hinv gamma_{i}|g_(k)> 4281 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 4282 do ii = 1, navec-1 4283 call rewino(lugrvf) 4284 v2(ii) = inprdd(vec1,vec2,luhg_sbsp,lugrvf,0,lblk) 4285 end do 4286 4287 if (ntest.ge.100) then 4288 write(6,*) '<delta_{i}|g_(k)>:' 4289 call wrtmat(v1,navec-1,1,navec-1,1) 4290 write(6,*) '<Hgamma_{i}|g_(k)>:' 4291 call wrtmat(v2,navec-1,1,navec-1,1) 4292 end if 4293 4294 do ii = 1, navec-1 4295 hv1(ii) = hkern(1,ii)*v1(ii) + hkern(2,ii)*v2(ii) 4296 hv2(ii) = hkern(2,ii)*v1(ii) + hkern(3,ii)*v2(ii) 4297 end do 4298 4299 if (ntest.ge.100) then 4300 write(6,*) 'hv1:' 4301 call wrtmat(hv1,navec-1,1,navec-1,1) 4302 write(6,*) 'hv2:' 4303 call wrtmat(hv2,navec-1,1,navec-1,1) 4304 end if 4305 4306 ! bring subspace-files into position: 4307 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 4308 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 4309 if (ndel1.gt.0) then 4310 ! store contributions from 'to-be-deleted' vectors on 4311 ! separate file ludlcnt 4312 call mvcsmd(lust_sbsp,hv1(1),luscr1,luscr2, 4313 & vec1,vec2,ndel1,0,lblk) 4314 call mvcsmd(luhg_sbsp,hv2(1),luscr2,luscr3, 4315 & vec1,vec2,ndel1,0,lblk) 4316 call vecsmd(vec1,vec2,1d0,1d0,luscr1,luscr2,ludlcnt,1,lblk) 4317 end if 4318 4319 ! assemble delta contributions 4320 call rewino(luscr1) 4321 call rewino(luscr2) 4322 call mvcsmd(lust_sbsp,hv1(1+ndel1),luscr1,luscr2, 4323 & vec1,vec2,navec-1-ndel1,0,lblk) 4324c call mvcsmd(lust_sbsp,hv1,luscr1,luscr2, 4325c & vec1,vec2,navec-1,0,lblk) 4326 ! add to Hinv_1g_k on lu_pertstp --> luscr3 4327 call vecsmd(vec1,vec2,1d0,-1d0,lu_pertstp,luscr1,luscr3,1,lblk) 4328 4329 ! assemble Hinv gamma contributions 4330 call rewino(luscr1) 4331 call rewino(luscr2) 4332 call mvcsmd(luhg_sbsp,hv2(1+ndel1),luscr1,luscr2, 4333 & vec1,vec2,navec-1-ndel1,0,lblk) 4334c call mvcsmd(luhg_sbsp,hv2,luscr1,luscr2, 4335c & vec1,vec2,navec-1,0,lblk) 4336 ! add to luscr3 --> lu_newstp (used as intermediate scratch) 4337 if (ndel1.eq.0) then 4338 call vecsmd(vec1,vec2,1d0,-1d0,luscr3,luscr1,lu_newstp,1,lblk) 4339 else 4340 ! add contribution from 'to-be-deleted' vectors here 4341 call vecsmd(vec1,vec2,1d0,1d0,luscr1,ludlcnt,luscr2,1,lblk) 4342 call vecsmd(vec1,vec2,1d0,-1d0,luscr3,luscr2,lu_newstp,1,lblk) 4343 end if 4344 4345 if (itype.ge.11.and.itype.le.14) then 4346 ! BFGS and DFP: 4347 ! so far, we have built -Hinv_{k-1}g_k on lu_newstp 4348 ! combine with -Hinv_{k-1}g_{k-1} on luhg_last to get 4349 ! Hinv_{k-1}gamma_{k-1} --> luhgam_new 4350 call vecsmd(vec1,vec2,-1d0,1d0, 4351 & lu_newstp,luhg_last,luhgam_new,1,lblk) 4352 end if 4353 4354 else ! navec.gt.1 4355 4356 ! else Hinv_{k-1}gamma{k-1} is identical with Hinv_1 g_{k} 4357 if (itype.ge.11.and.itype.le.14) 4358 & call sclvcd(lu_pertstp,luhgam_new,-1d0,vec1,1,lblk) 4359 call copvcd(lu_pertstp,lu_newstp,vec1,1,lblk) 4360 4361 end if 4362* 4363* what is expected at this position: 4364* -Hinv^{k-1}g^{k} on lu_newstp 4365* lust_sbsp positioned on delta^{k-1} 4366* 4367 if (navec.gt.0) then 4368 ! we take advantage of the file lust_sbsp being 4369 ! the correct postion for this 4370 ! delta^{k-1} --> luscr2 4371 call rewino(luscr2) 4372 call copvcd(lust_sbsp,luscr2,vec1,0,lblk) 4373 4374 end if 4375* 4376* 1) a. for Powell-Update: 4377* 4378* calculate H^{k-1} |delta^{k-1}> 4379* 4380 if (navec.gt.1.and.itype.eq.15) then 4381 4382 ! <delta_{i}|delta_(k-1)> 4383 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 4384 do ii = 1, navec-1 4385 call rewino(luscr2) 4386 v1(ii) = inprdd(vec1,vec2,lust_sbsp,luscr2,0,lblk) 4387 end do 4388 4389 ! <tau_{i}|delta_(k-1)> 4390 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 4391 do ii = 1, navec-1 4392 call rewino(luscr2) 4393 v2(ii) = inprdd(vec1,vec2,luhg_sbsp,luscr2,0,lblk) 4394 end do 4395 4396 if (ntest.ge.100) then 4397 write(6,*) '<delta_{i}|delta_(k-1)>:' 4398 call wrtmat(v1,navec-1,1,navec-1,1) 4399 write(6,*) '<tau_{i}|delta_(k-1)>:' 4400 call wrtmat(v2,navec-1,1,navec-1,1) 4401 end if 4402 4403 do ii = 1, navec-1 4404 hv1(ii) = hkern(1,ii)*v1(ii) + hkern(2,ii)*v2(ii) 4405 hv2(ii) = hkern(2,ii)*v1(ii) + hkern(3,ii)*v2(ii) 4406 end do 4407 4408 if (ntest.ge.100) then 4409 write(6,*) 'hv1:' 4410 call wrtmat(hv1,navec-1,1,navec-1,1) 4411 write(6,*) 'hv2:' 4412 call wrtmat(hv2,navec-1,1,navec-1,1) 4413 end if 4414 4415 ! bring subspace-files into position again: 4416 call skpvcd(lust_sbsp,nskipst,vec1,1,lblk) 4417 call skpvcd(luhg_sbsp,nskiphg,vec1,1,lblk) 4418 4419 ! assemble delta contributions 4420 call rewino(luscr1) 4421 call rewino(luscr2) 4422 call mvcsmd(lust_sbsp,hv1,luscr1,luscr3, 4423 & vec1,vec2,navec-1,0,lblk) 4424 4425 ! assemble gamma - Hinv delta contributions 4426 call rewino(luscr1) 4427 call rewino(luscr2) 4428 ! luhgam_new used as scratch: 4429 call mvcsmd(luhg_sbsp,hv2,luscr3,luhgam_new, 4430 & vec1,vec2,navec-1,0,lblk) 4431 ! add luscr1 and luscr3 --> luhgam_new 4432 call vecsmd(vec1,vec2,1d0,1d0,luscr1,luscr3,luhgam_new,1,lblk) 4433 4434 ! luhgam_new is still missing the contributions 4435 ! H_0^-1 delta^{k-1} and gamma^{k-1} 4436 ! coming soon ... 4437 4438 end if 4439 4440* 4441* 2) calculate new h^{k} 4442* 4443* we expect that delta^(k-1) is still save(d) on luscr2 4444* 4445 if (navec.gt.0.and.(itype.ge.11.and.itype.le.14)) 4446 & then 4447 4448 ! calculate the current gamma from current and previous gradient 4449 ! --> luscr3 4450 call vecsmd(vec1,vec2,1d0,-1d0,lugrvf,lugrvf_last,luscr3,1,lblk) 4451 4452 ! <gamma|Hgamma> 4453 x1 = inprdd(vec1,vec2,luscr3,luhgam_new,1,lblk) 4454 ! <delta^(k-1)|gamma^{k-1}> 4455 x2 = inprdd(vec1,vec2,luscr2,luscr3,1,lblk) 4456 ! <delta^(k-1)|g^{k}> 4457 x3 = inprdd(vec1,vec2,luscr2,lugrvf,1,lblk) 4458 ! <Hgamma|g^{k}> 4459 x4 = inprdd(vec1,vec2,luhgam_new,lugrvf,1,lblk) 4460 4461 if (itype.eq.11) then 4462 phi = 1d0 ! BFGS 4463 else if (itype.eq.12) then 4464 phi = 0d0 ! DFP 4465 else if (itype.eq.13) then 4466 phi = 1d0/(1d0 - x1/x2) ! Broyden rank 1 4467 else if (itype.eq.14) then 4468 phi = 1d0/(1d0 + x1/x2) ! Hoshino 4469 end if 4470 4471 hkern(1,navec) = (1d0+phi*x1/x2)/x2 4472 hkern(2,navec) = -phi/x2 4473 hkern(3,navec) = (phi-1d0)/x1 4474 4475 if (ntest.ge.100) then 4476 write(6,*) '<gamma^{k-1}|Hgamma^{k-1}> = ',x1 4477 write(6,*) ' <delta^(k-1)|gamma^{k-1}> = ',x2 4478 write(6,*) ' <delta^(k-1)|g^{k}> = ',x3 4479 write(6,*) ' <Hgamma^{k-1}|g^{k}> = ',x4 4480 end if 4481 else if (navec.gt.0.and.itype.eq.15) then 4482 4483 if (navec.gt.1) then 4484 ! Hinv^{1} delta^{k-1} on luscr3 4485 call dmtvcd(vec1,vec2,ludia,luscr2,luscr3,0d0,1,1,lblk) 4486 ! add to luhgam_new --> luscr1 4487 call vecsmd(vec1,vec2,1d0,1d0,luhgam_new,luscr3,luscr1,1,lblk) 4488 else 4489 ! no previous contributions? 4490 ! Hinv^{1} delta^{k-1} directly on luscr1 4491 call dmtvcd(vec1,vec2,ludia,luscr2,luscr1,0d0,1,1,lblk) 4492 end if 4493 4494 ! gamma^{k-1} on luscr3 4495 call vecsmd(vec1,vec2,1d0,-1d0,lugrvf,lugrvf_last,luscr3,1,lblk) 4496 4497 ! subtract luscr1 --> final result on luhgam_new 4498 call vecsmd(vec1,vec2,-1d0,1d0,luscr1,luscr3,luhgam_new,1,lblk) 4499 4500 ! <delta^{k-1}|gamma^{k-1}-Hinv^{k-1}delta^{k-1}> 4501 ! = <delta^{k-1}|tau^{k-1}> 4502 x1 = inprdd(vec1,vec2,luscr2,luhgam_new,1,lblk) 4503 ! <delta^{k-1}|delta^{k-1}> 4504 x2 = inprdd(vec1,vec1,luscr2,luscr2,1,lblk) 4505 ! <delta^{k-1}|g^k> 4506 x3 = inprdd(vec1,vec2,luscr2,lugrvf,1,lblk) 4507 ! <tau^{k-1}|g^k> 4508 x4 = inprdd(vec1,vec2,luhgam_new,lugrvf,1,lblk) 4509 4510c hkern(1,navec) = -x1/(x2*x2) 4511c hkern(2,navec) = 1d0/x2 4512c hkern(3,navec) = 0d0 4513c TEST 4514 hkern(1,navec) = -x1 4515 hkern(2,navec) = x2 4516 hkern(3,navec) = 0d0 4517 4518 if (ntest.ge.100) then 4519 write(6,*) ' <delta^{k-1}|tau^{k-1}> = ',x1 4520 write(6,*) '<delta^(k-1)|delta^{k-1}> = ',x2 4521 write(6,*) ' <delta^(k-1)|g^{k}> = ',x3 4522 write(6,*) ' <tau^{k-1}|g^{k}> = ',x4 4523 end if 4524 4525 end if 4526 4527 if (navec.gt.0.and.ntest.ge.100) then 4528 write(6,*) 'new rank n kernel (rank = ',nrank,'):' 4529 call prtrlt(hkern(1,navec),nrank) 4530 end if 4531* 4532* 3) complete Hinv^{k} g^{k} 4533* 4534 if (navec.gt.0) then 4535 f1 = hkern(1,navec) * x3 + hkern(2,navec) * x4 4536 f2 = hkern(2,navec) * x3 + hkern(3,navec) * x4 4537 4538 if (ntest.ge.100) then 4539 write(6,*) 'f1 = ', f1 4540 write(6,*) 'f2 = ', f2 4541 end if 4542 4543 ! add to what we already have on lu_newstp 4544 call vecsmd(vec1,vec2,1d0,-f1,lu_newstp,luscr2,luscr3,1,lblk) 4545 call vecsmd(vec1,vec2,1d0,-f2, 4546 & luscr3,luhgam_new,lu_newstp,1,lblk) 4547 end if ! navec.gt.0 4548 4549 ! save Hinv_k|g_k> for next iteration (only BFGS, DFP) 4550 if (itype.ge.11.and.itype.le.14) then 4551 if (ndel1.eq.0) then 4552 call copvcd(lu_newstp,luhg_last,vec1,1,lblk) 4553 else 4554 ! remove contributions from 'to-be-deleted' vectors 4555 ! such that in the next iteration we really calculate 4556 ! H^{k-1}gamma{k-1} = H^{k-1}g^k - H^{k-1}g^{k-1} 4557 ! as this ------------^^^ won't have these vectors anymore 4558 ! yes, the signs are ok: 4559 call vecsmd(vec1,vec2,1d0,1d0, 4560 & lu_newstp,ludlcnt,luhg_last,1,lblk) 4561 end if 4562 end if 4563 4564 ! if necessary, move h matrices down by ndel2 4565 if (ndel2.gt.0) then 4566 nn = nrank*(nrank+1)/2 4567 do ii = ndel2+1, navec 4568 hkern(1:nn,ii-ndel2) = hkern(1:nn,ii) 4569 end do 4570 end if 4571 4572 ! remove scratch files 4573 call relunit(luscr1,'delete') 4574 call relunit(luscr2,'delete') 4575 call relunit(luscr3,'delete') 4576 if (ndel1.gt.0) call relunit(ludlcnt,'delete') 4577 4578 return 4579 4580 end 4581*----------------------------------------------------------------------* 4582*----------------------------------------------------------------------* 4583 subroutine cmbamp(imode,lu1,lu2,lu3,lu_comb, 4584 & vec,namp1,namp2,namp3) 4585*----------------------------------------------------------------------* 4586* 4587* imode = 11 : combine vectors on lu1, lu2, lu3 into one single vector 4588* imode = 01 : reverse procedure 4589* 4590*----------------------------------------------------------------------* 4591 4592 implicit none 4593 4594 integer, parameter :: 4595 & ntest = 00 4596 4597 integer, intent(in) :: 4598 & imode, 4599 & lu1, lu2, lu3, lu_comb, 4600 & namp1, namp2, namp3 4601 real(8), intent(inout) :: 4602 & vec(*) 4603 4604 logical :: 4605 & l_end 4606 integer :: 4607 & ipass, 4608 & lblk, lu, iamzero, iampacked, namp, namp_r, namp_sum!, imone 4609 real(8), external :: 4610 & inprod 4611 4612 lblk = -1 4613 4614 if (ntest.ge.10) then 4615 write(6,*) 'cmbamp at work!' 4616 write(6,*) '===============' 4617 write(6,*) 'imode = ',imode 4618 if (imode.lt.12)write(6,*) 'lu1, lu2, lu3 = ', lu1, lu2, lu3 4619 write(6,*) 'lu_comb = ',lu_comb 4620 end if 4621 4622 if (imode.ge.11) then 4623* lu1, lu2, lu3 --> lu_comb 4624 4625 call rewino(lu_comb) 4626 do ipass = 1, 3 4627 if (ipass.eq.1) lu = lu1 4628 if (ipass.eq.2) lu = lu2 4629 if (ipass.eq.3) lu = lu3 4630 if (ipass.eq.1) namp = namp1 4631 if (ipass.eq.2) namp = namp2 4632 if (ipass.eq.3) namp = namp3 4633 if (ntest.ge.50) then 4634 write(6,*) 'ipass, lu, namp: ', ipass, lu, namp 4635 end if 4636 if (lu.gt.0) then 4637 call rewino(lu) 4638 namp_sum = 0 4639 l_end = .false. 4640 do while(.not.l_end) 4641 call ifrmds(namp_r,1,lblk,lu) 4642 if (namp_r.eq.-1) then 4643 if (ntest.ge.50) write(6,*) 'end of vector on lu ',lu 4644 l_end = .true. 4645 else 4646 namp_sum = namp_sum + namp_r 4647 if (ntest.ge.50) 4648 & write(6,*) 'transfer ',namp_r, 4649 & ' words from lu ',lu,' to lu_comb ',lu_comb 4650 call frmdsc(vec,namp_r,lblk,lu,iamzero,iampacked) 4651 if (ntest.ge.100) 4652 & write(6,*) ' norm of that block: ', 4653 & sqrt(inprod(vec,vec,namp_r)) 4654 call itods(namp_r,1,lblk,lu_comb) 4655 call todsc(vec,namp_r,lblk,lu_comb) 4656 end if 4657 end do 4658 if (namp_sum.ne.namp) then 4659 write(6,*) 'WARNING from cmbamp:' 4660 write(6,'(x,a,i2,2(a,i10))') ' File ',ipass, 4661 & ': read ',namp_sum,' expected: ',namp 4662 end if 4663 end if 4664 end do ! ipass 4665 ! mark end of record 4666 call itods(-1,1,lblk,lu_comb) 4667 4668 else 4669* lu_comb --> lu1, lu2 4670 4671 call rewino(lu_comb) 4672 pass_loop: do ipass = 1, 3 4673 if (ipass.eq.1) lu = lu1 4674 if (ipass.eq.2) lu = lu2 4675 if (ipass.eq.3) lu = lu3 4676 if (ipass.eq.1) namp = namp1 4677 if (ipass.eq.2) namp = namp2 4678 if (ipass.eq.3) namp = namp3 4679 4680 if (lu.gt.0) then 4681 call rewino(lu) 4682 namp_sum = 0 4683 l_end = .false. 4684 do while(.not.l_end) 4685 call ifrmds(namp_r,1,lblk,lu_comb) 4686 if (namp_r.eq.-1) then 4687 write(6,*) 4688 & 'WARNING from cmbamp: Unexpected end of vector' 4689 write(6,'(x,a,i2,2(a,i10))') ' File ',ipass, 4690 & ': read ',namp_sum,' expected: ',namp 4691 exit pass_loop 4692 else if (namp_r+namp_sum.gt.namp) then 4693 write(6,*) 4694 & 'WARNING from cmbamp: Too long new record ',namp_r 4695 write(6,'(x,a,i2,2(a,i10))') ' File ',ipass, 4696 & ': read ',namp_r,' expected: ',namp 4697 ! well, this is fatal 4698 stop 'cmbamp' 4699 l_end = .true. 4700 else 4701 namp_sum = namp_sum + namp_r 4702 if (ntest.ge.50) 4703 & write(6,*) 'transfer ',namp_r, 4704 & ' words from lu_comb ',lu_comb,' to lu ',lu 4705 call frmdsc(vec,namp_r,lblk,lu_comb,iamzero,iampacked) 4706 if (ntest.ge.100) 4707 & write(6,*) ' norm of that block: ', 4708 & sqrt(inprod(vec,vec,namp_r)) 4709 call itods(namp_r,1,lblk,lu) 4710 call todsc(vec,namp_r,lblk,lu) 4711 if (namp_sum.eq.namp) then 4712 if (ntest.ge.50) 4713 & write(6,*) 'vector on lu ',lu,' complete' 4714 l_end = .true. 4715 call itods(-1,1,lblk,lu) 4716 end if 4717 end if 4718 end do ! while (.not.l_end) 4719 end if ! lu.gt.0 4720 end do pass_loop 4721 ! skip the -1 4722c call ifrmds(imone,1,lblk,lu_comb) 4723 4724 end if 4725 4726 return 4727 end 4728*----------------------------------------------------------------------* 4729* follow: routines for 2nd-order optimization 4730*----------------------------------------------------------------------* 4731 subroutine optc_redh(isymm,hred,ndim,ndim_o, 4732 & luvec,lumv,luvec_sbsp,lumv_sbsp, 4733 & vec1,vec2) 4734*----------------------------------------------------------------------* 4735* update/set up reduced H matrix 4736* 4737* isymm != 0 : symmetrize (for more numerical stability) 4738* 4739*----------------------------------------------------------------------* 4740 implicit none 4741 4742 integer, parameter :: 4743 & ntest = 00 4744 4745 integer, intent(in) :: 4746 & isymm, 4747 & ndim, ndim_o, 4748 & luvec, lumv, luvec_sbsp, lumv_sbsp 4749 4750 real(8), intent(inout) :: 4751 & hred(ndim*ndim), vec1(*), vec2(*) 4752 4753 real(8) :: 4754 & xel, fac 4755 integer :: 4756 & ii, jj, 4757 & lblk, luscr 4758 4759 integer, external :: 4760 & iopen_nus 4761 real(8), external :: 4762 & inprdd 4763 4764 lblk = -1 4765 4766 if (ntest.ge.10) then 4767 write(6,*) '-----------' 4768 write(6,*) ' optc_redh' 4769 write(6,*) '-----------' 4770 write(6,*) ' ndim, ndim_o : ',ndim, ndim_o 4771 write(6,*) ' luvec,lumv,luvec_sbsp,lumv_sbsp: ', 4772 & luvec,lumv,luvec_sbsp,lumv_sbsp 4773 end if 4774 4775 if (ndim_o.gt.0) then 4776 if (ntest.ge.100) then 4777 write(6,*) 'H(red) on input:' 4778 call wrtmat2(hred,ndim_o,ndim_o,ndim_o,ndim_o) 4779 end if 4780 4781 do ii = ndim_o, 1, -1 4782 do jj = ndim_o, 1, -1 4783 hred((ii-1)*ndim+jj) = hred((ii-1)*ndim_o+jj) 4784 end do 4785 end do 4786 4787 if (ntest.ge.100) then 4788 write(6,*) 'H(red) on input (resorted):' 4789 call wrtmat2(hred,ndim,ndim,ndim,ndim) 4790 end if 4791 4792 end if 4793 4794 if (ndim-ndim_o.gt.1) then 4795 stop 'i shall not enter this section!' 4796 luscr = iopen_nus('HRED_SCR') 4797 4798 call skpvcd(luvec_sbsp,ndim_o,vec1,1,lblk) 4799 do ii = ndim_o+1, ndim-1 4800 call copvcd(luvec_sbsp,luscr,vec1,0,lblk) 4801 call rewino(lumv_sbsp) 4802 do jj = 1, ndim-1 4803 call rewino(luscr) 4804 xel = inprdd(vec1,vec2,luscr,lumv_sbsp,0,lblk) 4805c hred((ii-1)*ndim+jj) = xel 4806 hred((jj-1)*ndim+ii) = xel 4807 end do 4808 end do 4809 4810 call skpvcd(lumv_sbsp,ndim_o,vec1,1,lblk) 4811 do jj = ndim_o+1, ndim-1 4812 call copvcd(lumv_sbsp,luscr,vec1,0,lblk) 4813 call rewino(luvec_sbsp) 4814 do ii = 1, ndim_o 4815 call rewino(luscr) 4816 xel = inprdd(vec1,vec2,luscr,luvec_sbsp,0,lblk) 4817c hred((ii-1)*ndim+jj) = xel 4818 hred((jj-1)*ndim+ii) = xel 4819 end do 4820 end do 4821 4822 call relunit(luscr,'delete') 4823 4824 end if 4825 4826 if (ndim.gt.1) then 4827 fac = 1d0 4828 if (isymm.ne.0) fac = 0.5d0 4829 call rewino(lumv_sbsp) 4830 do jj = 1, ndim-1 4831 call rewino(luvec) 4832 xel = inprdd(vec1,vec2,luvec,lumv_sbsp,0,lblk) 4833c hred((ndim-1)*ndim+jj) = xel 4834 hred((jj-1)*ndim+ndim) = fac*xel 4835 if (isymm.ne.0) 4836 & hred((ndim-1)*ndim+jj) = fac*xel 4837 end do 4838 if (isymm.eq.0) hred((ndim-1)*ndim+1:(ndim-1)*ndim+ndim) = 0d0 4839 call rewino(luvec_sbsp) 4840 do ii = 1, ndim-1 4841 call rewino(lumv) 4842 xel = inprdd(vec1,vec2,lumv,luvec_sbsp,0,lblk) 4843 hred((ndim-1)*ndim+ii) = hred((ndim-1)*ndim+ii) + fac*xel 4844 if (isymm.ne.0) 4845 & hred((ii-1)*ndim+ndim) = hred((ii-1)*ndim+ndim) +fac*xel 4846 end do 4847 end if 4848 4849 xel = inprdd(vec1,vec2,lumv,luvec,1,lblk) 4850 hred(ndim*ndim) = xel 4851 4852 if (ntest.ge.100) then 4853 write(6,*) 'Final H(red):' 4854 call wrtmat2(hred,ndim,ndim,ndim,ndim) 4855 end if 4856 4857 return 4858 end 4859*----------------------------------------------------------------------* 4860 subroutine optc_trnewton(imode,iret,isymmet, 4861 & hred,gred,cred,scr,scr2,ndim,xlamb,gamma,trrad,de_pred) 4862*----------------------------------------------------------------------* 4863* solve reduced trust-radius newton equations 4864* 4865* imode = 1: try bare newton step first 4866* the subspace Hessian will be tested for negative 4867* eigenvalues; if so, iret!=0 will signal that this 4868* does not work 4869* imode = 2: solve augmented Hessian equations with adapted damping 4870* (reminiscent of the NEO, see "the book") If a too 4871* low gamma is found (which is, together with the lowest 4872* eigenvalue of the augmented Hessian, nu, related to 4873* the usual damping parameter lambda = nu*gamma**2) 4874* the routine will use the lowest possible gamma and 4875* suggest to take the newton step by setting iret!=0 4876* 4877*----------------------------------------------------------------------* 4878 4879 implicit none 4880 4881 integer, parameter :: 4882 & ntest = 00, maxit = 400 4883 real(8), parameter :: 4884 & thrsh = 1d-6, xmxstp = 1.0d0, gamma_min = 1d-1 4885 4886 integer, intent(in) :: 4887 & ndim, imode 4888 real(8), intent(in) :: 4889 & hred(ndim*ndim),gred(ndim), 4890 & trrad 4891 4892 integer, intent(inout) :: 4893 & isymmet 4894 integer, intent(out) :: 4895 & iret 4896 real(8), intent(out) :: 4897 & de_pred,xlamb, cred(ndim) 4898 real(8), intent(inout) :: 4899 & gamma, 4900 & scr((ndim+1)*(ndim+1)),scr2((ndim+1)*(ndim+1)) 4901 4902 logical :: 4903 & opt, conv 4904 integer :: 4905 & iter, ii, jj, idx_min, ierr, irg 4906 real(8) :: 4907 & eig_min, dlt, fac, xel, 4908 & xnrm, xnrm_prev, xdum, gamma_prev, gamma_new, 4909 & gamma_high, gamma_low 4910 ! local O(N) scratch 4911 integer :: 4912 & iscr(ndim+1) 4913 real(8) :: 4914 & xscr(ndim+1), eigr(ndim+1), eigi(ndim+1) 4915 4916 integer, external :: 4917 & iopen_nus 4918 real(8), external :: 4919 & inprdd, inprod 4920 4921 if (ntest.ge.10) then 4922 write(6,*) '---------------' 4923 write(6,*) ' optc_trnewton' 4924 write(6,*) '---------------' 4925 write(6,*) ' imode = ',imode 4926 write(6,*) ' ndim = ',ndim 4927 write(6,*) ' trrad = ',trrad 4928 end if 4929 4930 iret = 0 ! start optimistic 4931 4932 if (imode.eq.1) then 4933 if (ntest.ge.10) then 4934 write(6,*) '--------------------' 4935 write(6,*) ' trying Newton step' 4936 write(6,*) '--------------------' 4937 end if 4938 4939 ! get a copy of hred 4940 scr(1:ndim*ndim) = hred(1:ndim*ndim) 4941 4942 ! and solve for eigenvalues 4943 irg = 0 4944 call rg(ndim,ndim,scr,eigr,eigi,irg,xdum,iscr,xscr,ierr) 4945 if (ierr.ne.0) then 4946 write(6,*) 'Internal ERROR: rg gives ierr = ',ierr 4947 stop 'optc_trnewton: rg' 4948 end if 4949 4950 if (ntest.ge.100) then 4951 write(6,*) 'eigenvalues of Hred:' 4952 do ii = 1, ndim 4953 write(6,*) ii,eigr(ii),eigi(ii) 4954 end do 4955 end if 4956 4957 eig_min = 0d0 4958 do ii = 1, ndim 4959 eig_min = min(eig_min,eigr(ii)) 4960 end do 4961 4962 if (eig_min.lt.0d0) then 4963 write(6,*) 'detected negative eigenvalue in subspace' 4964 write(6,*) 'you are clearly not in the local region!' 4965 write(6,*) 'I will switch to NEO algorithm' 4966 iret = 1 4967 end if 4968 4969 scr(1:ndim*ndim) = hred(1:ndim*ndim) 4970 xscr(1:ndim) = -gred(1:ndim) 4971 ! solve reduced LEQ 4972 call dgefa(scr,ndim,ndim,iscr,ierr) 4973 if (ierr.ne.0) then 4974 write(6,*) 'Internal ERROR: dgefa gives ierr = ',ierr 4975 stop 'optc_trnewton: dgefa' 4976 end if 4977 call dgesl(scr,ndim,ndim,iscr,xscr,0) 4978 4979 xnrm = sqrt(inprod(xscr,xscr,ndim)) 4980 4981 if (ntest.ge.100) then 4982 write(6,*) 'the Newton step: ' 4983 call wrtmat(xscr,ndim,1,ndim,1) 4984 write(6,*) ' norm = ',xnrm 4985 write(6,*) ' trrad = ',trrad 4986 end if 4987 4988 if (xnrm.gt.trrad) then 4989 write(6,*) 'detected too long step!' 4990 write(6,*) 'you are clearly not in the local region!' 4991 write(6,*) 'I will switch to NEO algorithm' 4992 iret = 1 4993 end if 4994 4995 if (iret.eq.0) then 4996 cred(1:ndim) = xscr(1:ndim) 4997 ! energy prediction 4998 call matvcb(hred,cred,xscr, 4999 & ndim,ndim,0) 5000 de_pred = inprod(gred,cred,ndim) + 5001 & 0.5d0*inprod(xscr,cred,ndim) 5002 end if 5003 xlamb = 0d0 5004 5005 end if 5006 5007 if (imode.eq.2) then 5008 5009 if (ntest.ge.10) then 5010 write(6,*) '--------------------------------------' 5011 write(6,*) ' solving Newton eigenvector equations' 5012 write(6,*) '--------------------------------------' 5013 end if 5014 5015 ! gamma is set outside 5016c TEST 5017c gamma = 20d0 5018c TEST 5019 ! signal that no previous gamma is available 5020 gamma_prev = -1d0 5021 gamma_low = gamma_min 5022 gamma_high = 1d20 5023 iter = 0 5024 do 5025 iter = iter + 1 5026 5027 do ! trial loop: 1 w/o symmetrization, 2 w/ symmetrization 5028 ! set up augmented Hessian (g scaled with gamma) 5029 do jj = 1, ndim 5030 do ii = 1, ndim 5031 scr((jj-1)*(ndim+1)+ii) = hred((jj-1)*ndim+ii) 5032 end do 5033 scr((jj-1)*(ndim+1)+ndim+1) = gamma*gred(jj) 5034 end do 5035 do ii = 1, ndim 5036 scr(ndim*(ndim+1)+ii) = gamma*gred(ii) 5037 end do 5038 scr((ndim+1)*(ndim+1))=0d0 5039 5040 if (isymmet.eq.1) then 5041 do ii = 1, ndim 5042 do jj = 1, ii-1 5043 xel = 0.5d0*(scr((ii-1)*(ndim+1)+jj)+ 5044 & scr((jj-1)*(ndim+1)+ii)) 5045 scr((ii-1)*(ndim+1)+jj)=xel 5046 scr((jj-1)*(ndim+1)+ii)=xel 5047 end do 5048 end do 5049 end if 5050 5051 if (ntest.ge.100) then 5052 write(6,*) 'gradient scaled augmented Hessian (gamma = ', 5053 & gamma,')' 5054 call wrtmat2(scr,ndim+1,ndim+1,ndim+1,ndim+1) 5055 end if 5056 5057 irg = 1 5058c call test_symmat(scr,ndim+1,ndim+1) 5059 call rg(ndim+1,ndim+1,scr,eigr,eigi,irg,scr2,iscr,xscr,ierr) 5060 if (ierr.ne.0) then 5061 write(6,*) 'Internal ERROR: rg gives ierr = ',ierr 5062 stop 'optc_trnewton: rg' 5063 end if 5064 5065 if (ntest.ge.100) then 5066 write(6,*) 'eigenvalues of Hred:' 5067 do ii = 1, ndim+1 5068 write(6,*) ii,eigr(ii),eigi(ii) 5069 end do 5070 end if 5071 5072 eig_min = 100d0 5073 idx_min = 1 5074 do ii = 1, ndim+1 5075 if (eig_min.gt.eigr(ii)) then 5076 eig_min = eigr(ii) 5077 idx_min = ii 5078 end if 5079 end do 5080c if (abs(eigi(idx_min)).gt.epsilon(1d0)) then 5081c write(6,*) 'Help, imaginary lowest eigenvalue!' 5082c stop 'trnewton' 5083c end if 5084 if (abs(eigi(idx_min)).lt.epsilon(1d0)) exit 5085 5086 if (abs(eigi(idx_min)).gt.epsilon(1d0) 5087 & .and.isymmet.eq.1) then 5088 ! then something is REALLY wrong 5089 write(6,*) 5090 & 'Strange, imaginary lowest eigenvalue persists!' 5091 stop 'trnewton' 5092 end if 5093 5094 write(6,*) '>> symmetrizing trick (auweia)' 5095 isymmet = 1 ! OK, so we try to symmetrize (desperate trick) 5096 5097 end do 5098 isymmet = 0 5099 5100 ! rescale eigenvector 5101 xscr(1:ndim) = (1d0/gamma)*scr2((idx_min-1)*(ndim+1)+1: 5102 & (idx_min-1)*(ndim+1)+ndim) 5103 xscr(ndim+1) = scr2((idx_min-1)*(ndim+1)+ndim+1) 5104 if (ntest.ge.100) then 5105 write(6,*) 'the raw eigenvector to ',eig_min 5106 call wrtmat(xscr,1,ndim+1,1,ndim+1) 5107 end if 5108 5109 ! normalize such that the last element is 1 5110 fac = xscr(ndim+1) 5111 call scalve(xscr,1d0/fac,ndim+1) 5112 5113 if (ntest.ge.100) then 5114 write(6,*) 'the renormalized eigenvector to ',eig_min 5115 call wrtmat(xscr,1,ndim+1,1,ndim+1) 5116 end if 5117 5118 ! get step length (only the first ndim elements) 5119 xnrm = sqrt(inprod(xscr,xscr,ndim)) 5120 5121 if (ntest.ge.100) then 5122 write(6,'(x,a,i4,2(x,e12.6))') 5123 & ' iter, gamma, step length: ', iter, gamma, xnrm 5124 write(6,'(x,a,18x,e12.6)') 5125 & ' trrad = ', trrad 5126 end if 5127 5128 if (abs(xnrm-trrad).gt.thrsh) then 5129 if (xnrm.gt.trrad) then 5130 ! save highest gamma for that so far xnrm.gt.trrad 5131 ! i.e. the low limit on gamma 5132 gamma_low = max(gamma,gamma_low) 5133 if (gamma_prev.eq.-1d0) then 5134 gamma_new = gamma + 0.1d0 5135 else 5136 dlt = -(gamma-gamma_prev)/(xnrm-xnrm_prev)*(xnrm-trrad) 5137 if (dlt.lt.0d0) dlt = xmxstp 5138 gamma_new = gamma + min(xmxstp,dlt) 5139 end if 5140 else 5141 ! save lowest gamma for that so far xnrm.lt.trrad 5142 ! i.e. the high limit on gamma 5143 gamma_high = min(gamma,gamma_high) 5144 if (gamma.eq.gamma_min) then 5145 write(6,*) 'arrived at lowest possible gamma;' 5146 write(6,*) 'suggesting to take newton step instead!' 5147 iret = 1 5148 exit 5149 end if 5150c if (abs(xnrm-trrad).lt.0.1d0.and.gamma_prev.ne.-1d0) then 5151 if (gamma_prev.ne.-1d0) then 5152 dlt = -(gamma-gamma_prev)/(xnrm-xnrm_prev)*(xnrm-trrad) 5153 if (dlt.gt.0d0) dlt = -xmxstp 5154 dlt = max(-xmxstp,dlt) 5155 else 5156 dlt = -0.1d0 5157 end if 5158 gamma_new = max(gamma_min,gamma + dlt) 5159 end if 5160 ! stabilize by bracketing 5161 if (gamma_new.gt.gamma_high.or. 5162 & gamma_new.lt.gamma_low) then 5163 gamma_new = 0.5*(gamma_high+gamma_low) 5164 end if 5165 if (ntest.ge.100) then 5166 write(6,*) 'optimization of gamma:' 5167 write(6,*) ' gamma_high, gamma_low: ',gamma_high,gamma_low 5168 write(6,*) ' previous gamma, step: ',gamma_prev,xnrm_prev 5169 write(6,*) ' current gamma, step : ',gamma,xnrm 5170 write(6,*) ' new gamma: ',gamma_new 5171 end if 5172 gamma_prev = gamma 5173 gamma = gamma_new 5174 xnrm_prev = xnrm 5175 else 5176 if (ntest.ge.100) write(6,*) 'exited optim. loop' 5177 exit 5178 end if 5179 5180 if (iter.gt.maxit) then 5181 write(6,*) ' problem in optimizing gamma' 5182 stop 'trnewton' 5183 end if 5184 5185 end do 5186 5187 cred(1:ndim) = xscr(1:ndim) 5188c de_pred = 0.5d0*eig_min/(gamma*gamma) 5189c print *,' energy prediction 1: ',0.5d0*eig_min/(gamma*gamma) 5190 call matvcb(hred,cred,xscr, 5191 & ndim,ndim,0) 5192 de_pred = inprod(gred,cred,ndim) + 5193 & 0.5d0*inprod(xscr,cred,ndim) 5194c print *,' energy prediction 2: ',de_pred 5195 5196 xlamb = eig_min!*gamma*gamma 5197 5198 end if 5199 5200 if (ntest.ge.100) then 5201 write(6,*) 'final solution for lambda = ',xlamb 5202 call wrtmat(cred,ndim,1,ndim,1) 5203 write(6,*) ' norm = ',xnrm 5204 write(6,*) ' trrad = ',trrad 5205 end if 5206 5207 return 5208 end 5209*----------------------------------------------------------------------* 5210 subroutine optc_trn_resid(cred,ndim,xlamb, 5211 & lures,xresnrm, 5212 & luvec_sbsp,lumv_sbsp,lugrvf, 5213 & vec1,vec2) 5214*----------------------------------------------------------------------* 5215* set up residual vector in full space 5216*----------------------------------------------------------------------* 5217 implicit none 5218 5219 integer, parameter :: 5220 & ntest = 00 5221 5222 integer, intent(in) :: 5223 & ndim, lures, luvec_sbsp, lumv_sbsp, lugrvf 5224 real(8), intent(in) :: 5225 & xlamb, cred(ndim) 5226 real(8), intent(out) :: 5227 & xresnrm 5228 real(8), intent(inout) :: 5229 & vec1(*), vec2(*) 5230 5231 integer :: 5232 & lblk, luscr1, luscr2 5233 5234 integer, external :: 5235 & iopen_nus 5236 real(8), external :: 5237 & inprdd 5238 5239 lblk = -1 5240 5241 luscr1 = iopen_nus('RES_SCR1') 5242 5243 call mvcsmd(lumv_sbsp,cred,luscr1,lures,vec1,vec2,ndim,1,lblk) 5244 5245 if (abs(xlamb).gt.epsilon(1d0)) then 5246 luscr2 = iopen_nus('RES_SCR2') 5247 call mvcsmd(luvec_sbsp,cred,lures,luscr2,vec1,vec2,ndim,1,lblk) 5248 call vecsmd(vec1,vec2,-xlamb,1d0,lures,luscr1,luscr2,1,lblk) 5249 call vecsmd(vec1,vec2,1d0,1d0,luscr2,lugrvf,lures,1,lblk) 5250 call relunit(luscr2,'delete') 5251 else 5252 call vecsmd(vec1,vec2,1d0,1d0,luscr1,lugrvf,lures,1,lblk) 5253 end if 5254 5255 xresnrm = sqrt(inprdd(vec1,vec1,lures,lures,1,lblk)) 5256 5257 call relunit(luscr1,'delete') 5258 5259 return 5260 end 5261 5262 subroutine optc_prjout(nrdvec,lurdvec,luvec, 5263 & vec1,vec2,nwfpar,lincore) 5264 5265 implicit none 5266 5267 integer, parameter :: 5268 & ntest = 00 5269 5270 logical, intent(in) :: 5271 & lincore 5272 integer, intent(in) :: 5273 & nrdvec, lurdvec, luvec, nwfpar 5274 real(8), intent(inout) :: 5275 & vec1(*), vec2(*) 5276 5277 integer :: 5278 & irdvec, luscr1, luscr2 5279 real(8) :: 5280 & ovl(nrdvec), xnrm 5281 real(8), external :: 5282 & inprod, inprdd 5283 integer, external :: 5284 & iopen_nus 5285 5286 if (lincore) then 5287 call vec_from_disc(vec1,nwfpar,1,-1,luvec) 5288 end if 5289 if (ntest.ge.100) then 5290 if (lincore) then 5291 xnrm = sqrt(inprod(vec1,vec1,nwfpar)) 5292 else 5293 xnrm = sqrt(inprdd(vec1,vec1,luvec,luvec,1,-1)) 5294 end if 5295 write(6,*) ' norm of unprojected vector: ',xnrm 5296 end if 5297 5298 call rewino(lurdvec) 5299 do irdvec = 1, nrdvec 5300 if (lincore) then 5301 call vec_from_disc(vec2,nwfpar,0,-1,lurdvec) 5302 ovl(irdvec) = inprod(vec1,vec2,nwfpar) 5303 else 5304 call rewino(luvec) 5305 ovl(irdvec) = inprdd(vec1,vec2,luvec,lurdvec,0,-1) 5306 end if 5307 if (ntest.ge.100) 5308 & write(6,*) ' overlap with vec ',irdvec,' :',ovl(irdvec) 5309 if (lincore) 5310 & vec1(1:nwfpar) = 5311 & vec1(1:nwfpar)-ovl(irdvec)*vec2(1:nwfpar) 5312 end do 5313 5314 if (lincore) then 5315 call vec_to_disc(vec1,nwfpar,1,-1,luvec) 5316 else 5317 luscr1 = iopen_nus('PRJOUT_SCR1') 5318 luscr2 = iopen_nus('PRJOUT_SCR2') 5319 5320 call mvcsmd(lurdvec,ovl,luscr1,luscr2,vec1,vec2,nrdvec,1,-1) 5321 call vecsmd(vec1,vec2,1d0,-1d0,luvec,luscr1,luscr2,1,-1) 5322 call copvcd(luscr2,luvec,vec1,1,-1) 5323 5324 call relunit(luscr1,'delete') 5325 call relunit(luscr2,'delete') 5326 end if 5327 5328 if (ntest.ge.100) then 5329 if (lincore) then 5330 xnrm = sqrt(inprod(vec1,vec1,nwfpar)) 5331 else 5332 xnrm = sqrt(inprdd(vec1,vec1,luvec,luvec,1,-1)) 5333 end if 5334 write(6,*) ' norm of projected vector: ',xnrm 5335 5336 end if 5337 5338 return 5339 end 5340c $Id$ 5341