Skip to content

Commit

Permalink
240411.105737.HKT revoke: try using r2update instead of outprod i…
Browse files Browse the repository at this point in the history
…n `shiftbase`
  • Loading branch information
zaikunzhang committed Apr 11, 2024
1 parent 42c84f1 commit 0e547e8
Showing 1 changed file with 6 additions and 7 deletions.
13 changes: 6 additions & 7 deletions fortran/common/shiftbase.f90
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ module shiftbase_mod
!
! Started: July 2020
!
! Last Modified: Thursday, April 11, 2024 AM10:50:53
! Last Modified: Thursday, April 11, 2024 AM10:57:22
!--------------------------------------------------------------------------------------------------!

implicit none
Expand Down Expand Up @@ -42,10 +42,10 @@ subroutine shiftbase_lfqint(kopt, xbase, xpt, zmat, bmat, pq, hq, idz)
!--------------------------------------------------------------------------------------------------!

! Common modules
use, non_intrinsic :: consts_mod, only : RP, IK, ONE, ZERO, HALF, QUART, DEBUGGING
use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, HALF, QUART, DEBUGGING
use, non_intrinsic :: debug_mod, only : assert
use, non_intrinsic :: infnan_mod, only : is_finite
use, non_intrinsic :: linalg_mod, only : inprod, matprod, outprod, issymmetric, r2update
use, non_intrinsic :: linalg_mod, only : inprod, matprod, outprod, issymmetric

implicit none

Expand All @@ -72,7 +72,7 @@ subroutine shiftbase_lfqint(kopt, xbase, xpt, zmat, bmat, pq, hq, idz)
real(RP) :: qxoptq
real(RP) :: sxpt(size(xpt, 2))
real(RP) :: v(size(xbase))
!real(RP) :: vxopt(size(xbase), size(xbase))
real(RP) :: vxopt(size(xbase), size(xbase))
real(RP) :: xopt(size(xbase))
real(RP) :: xoptsq
real(RP) :: xptxav(size(xpt, 1), size(xpt, 2))
Expand Down Expand Up @@ -142,9 +142,8 @@ subroutine shiftbase_lfqint(kopt, xbase, xpt, zmat, bmat, pq, hq, idz)
! Update the quadratic model. Note that PQ remains unchanged. For HQ, see (7.14) of the NEWUOA paper.
!v = matprod(xptxav, pq) ! Vector V in (7.14) of the NEWUOA paper
v = matprod(xpt, pq) - HALF * sum(pq) * xopt ! This one seems to work better numerically.
!vxopt = outprod(v, xopt) !!MATLAB: vxopt = v * xopt'; % v and xopt should be both columns
!hq = (vxopt + transpose(vxopt)) + hq
call r2update(hq, ONE, xopt, v)
vxopt = outprod(v, xopt) !!MATLAB: vxopt = v * xopt'; % v and xopt should be both columns
hq = (vxopt + transpose(vxopt)) + hq !call r2update(hq, ONE, xopt, v)
!call symmetrize(hq) ! Do this if the update above does not ensure symmetry.

! The following instructions complete the shift of XBASE.
Expand Down

0 comments on commit 0e547e8

Please sign in to comment.