Skip to content

GPEC, DCON, RDCON, STRIDE - BUGFIX - Fix LAPACK and ZVODE type and size errors #299

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
matt-pharr merged 2 commits into from
Oct 9, 2026
Merged

GPEC, DCON, RDCON, STRIDE - BUGFIX - Fix LAPACK and ZVODE type and size errors #299

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b426793
GPEC, DCON, RDCON, STRIDE - BUGFIX - Fix LAPACK and ZVODE type and si…
matt-pharr Oct 9, 2026
5efba83
DCON, RDCON - MINOR - Address review comments (make output happen bef…
matt-pharr Oct 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions dcon/free.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out)
ev(ipert)=wvt(ipert,ipert)
ENDDO
plasma1=REAL(ep(1))
vacuum1=REAL(ev(1))
vacuum1=ev(1)
total1=REAL(et(1))
c-----------------------------------------------------------------------
c write data for ahb and deallocate.
Expand All @@ -252,15 +252,15 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out)
c-----------------------------------------------------------------------
c write to screen and copy to output.
c-----------------------------------------------------------------------
IF(verbose) WRITE(*,10) REAL(ep(1)),REAL(ev(1)),
IF(verbose) WRITE(*,10) REAL(ep(1)),ev(1),
$ REAL(et(1)),AIMAG(et(1))
c-----------------------------------------------------------------------
c write eigenvalues to file.
c-----------------------------------------------------------------------
message=""
WRITE(out_unit,'(/1x,a)')"Total Energy Eigenvalues:"
WRITE(out_unit,20)
WRITE(out_unit,30)(isol,REAL(ep(isol)),REAL(ev(isol)),
WRITE(out_unit,30)(isol,REAL(ep(isol)),ev(isol),
$ REAL(et(isol)),AIMAG(et(isol)),
$ TRIM(message(isol)),isol=1,mpert)
WRITE(out_unit,20)
Expand All @@ -271,7 +271,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out)
m=mlow+(/(isol,isol=0,mpert-1)/)
DO isol=1,mpert
WRITE(out_unit,40)
WRITE(out_unit,50)isol,imax(1),REAL(ep(isol)),REAL(ev(isol)),
WRITE(out_unit,50)isol,imax(1),REAL(ep(isol)),ev(isol),
$ REAL(et(isol)),AIMAG(et(isol)),TRIM(message(isol))
WRITE(out_unit,60)
WRITE(out_unit,70)(ipert,m(ipert),wt(ipert,isol),
Expand Down Expand Up @@ -302,6 +302,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out)
wp(:,ipert)=vr(:,eindex(mpert+1-ipert))
ep(ipert)=tt(eindex(mpert+1-ipert))
ENDDO
lwork=2*mpert-1
CALL zheev('V','U',mpert,wv,mpert,ev,work,lwork,rwork,info)
c-----------------------------------------------------------------------
c optionally write netcdf file.
Expand Down
20 changes: 20 additions & 0 deletions dcon/ode_output.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -273,13 +273,20 @@ SUBROUTINE ode_output_get_evals
10 FORMAT(/a,i4,1p,2(a,es10.3))
20 FORMAT(/5x,"i",5x,"eval",7x,"evali",7x,"err"/)
30 FORMAT(i6,1p,3e11.3)
40 FORMAT(/"istep = ",i4,", psifac = ",1p,es10.3,
$ ": skipped, ",a," is singular (zgetrf info = ",i0,")")
c-----------------------------------------------------------------------
c compute plasma response matrix.
c-----------------------------------------------------------------------
temp=CONJG(TRANSPOSE(u(:,1:mpert,1)))
wp=u(:,1:mpert,2)
wp=CONJG(TRANSPOSE(wp))
CALL zgetrf(mpert,mpert,temp,mpert,ipiv,info)
IF(info /= 0)THEN
IF(out_evals .AND. istep > 0)
$ WRITE(evals_out_unit,40)istep,psifac,"u(:,:,1)",info
RETURN
ENDIF
CALL zgetrs('N',mpert,mpert,temp,mpert,ipiv,wp,mpert,info)
wp=(wp+CONJG(TRANSPOSE(wp)))/2
c-----------------------------------------------------------------------
Expand All @@ -291,6 +298,19 @@ SUBROUTINE ode_output_get_evals
key=-ABS(1/evals)
CALL bubble(key,index,1,mpert)
c-----------------------------------------------------------------------
c compute inverse plasma response matrix (zheev destroys wp).

@logan-nc logan-nc Oct 9, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am surprised to see this needed here since I see no new zheev in the diff and even a Xtmp=X guard added. Is this fixing a bug introduced by another PR? If so, it needs to be clearly marked as urgent!

It seems like that style of guard is the better way to go rather than introducing repeat wp calculations that are hard to maintain.

Even the rdcon diff is different from the dcon one already. Inevitably someone will change the original wp calculation location and not realize they need to go fix these other locations as well.

@matt-pharr matt-pharr Oct 9, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could do instead an Xtmp=X type thing, yes, and I see your point why that one would be better. Though I am not sure I understand the rest of your comment about another PR. The ZHEEV statement that had the bug looks like it was added 10 years ago when you added RDCON.

c-----------------------------------------------------------------------
temp=CONJG(TRANSPOSE(u(:,1:mpert,2)))
wp=CONJG(TRANSPOSE(u(:,1:mpert,1)))
CALL zgetrf(mpert,mpert,temp,mpert,ipiv,info)
IF(info /= 0)THEN
IF(out_evals .AND. istep > 0)
$ WRITE(evals_out_unit,40)istep,psifac,"u(:,:,2)",info
RETURN
ENDIF
CALL zgetrs('N',mpert,mpert,temp,mpert,ipiv,wp,mpert,info)
wp=(wp+CONJG(TRANSPOSE(wp)))/2
c-----------------------------------------------------------------------
c compute and sort inverse eigenvalues.
c-----------------------------------------------------------------------
lwork=2*mpert-1
Expand Down
3 changes: 2 additions & 1 deletion gpec/gpdiag.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1821,7 +1821,8 @@ SUBROUTINE gpdiag_reluctpowout(rout,bpout,bout,rcout)
INTEGER :: i,j,lwork
REAL(r8), DIMENSION(3*mpert-2) :: rwork
COMPLEX(r8), DIMENSION(2*mpert-1) :: work
COMPLEX(r8), DIMENSION(mpert) :: temp,ev,evo
COMPLEX(r8), DIMENSION(mpert) :: temp
REAL(r8), DIMENSION(mpert) :: ev,evo
COMPLEX(r8), DIMENSION(mpert,mpert) :: sqrta,mat,mato,
$ sqrtao

Expand Down
5 changes: 3 additions & 2 deletions gpec/gpglobal.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,8 @@ MODULE gpglobal_mod
$ surf_indev,vsurf_indev,fsurf_indev,surf_indinvev
$ ,eft,efp
REAL(r8), DIMENSION(:,:), ALLOCATABLE :: xzpts,
$ chperr,chpsqr,grri,grre,griw,grrw,gdr,gdz,gdpsi,gdthe,gdphi
$ chperr,chpsqr,grri,grre,griw,grrw,gdr,gdz,gdpsi,gdthe,gdphi,
$ permeabsv
REAL(r8), DIMENSION(3,3) :: w,v


Expand All @@ -74,7 +75,7 @@ MODULE gpglobal_mod
$ bwp_mn_rmatch
COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: wt,wt0,chp_mn,kap_mn,
$ permeabev,chimats,chemats,flxmats,kaxmats,singbno_mn,
$ plas_indev,plas_indinvev,reluctev,indrelev,permeabsv,
$ plas_indev,plas_indinvev,reluctev,indrelev,
$ surf_indmats,surf_indevmats,vsurf_indmats,fsurf_indmats,
$ surf_indinvmats,surf_indinvevmats,surfet,surfep,
$ amat,bmat,cmat,fmats,gmats,kmats,coil_indmat,coilmn,
Expand Down
4 changes: 2 additions & 2 deletions gpec/gpout.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ SUBROUTINE gpout_response(rout,bpout,bout,rcout,tout,jout)
WRITE(out_unit,'(1x,I4,8(1x,es17.8))') i,
$ surf_indev(i),REAL(plas_indev(resp_index,i)),
$ REAL(permeabev(resp_index,i)),AIMAG(permeabev(resp_index,i)),
$ REAL(permeabsv(resp_index,i)),REAL(reluctev(resp_index,i)),
$ permeabsv(resp_index,i),REAL(reluctev(resp_index,i)),
$ REAL(et(i)),REAL(plas_indinvev(resp_index,i))
ENDDO
WRITE(out_unit,*)
Expand Down Expand Up @@ -5985,7 +5985,7 @@ SUBROUTINE gpout_control_filter(mode,finmn,foutmn,ftypes,fmodes,
mat = permeabmats(resp_index,:,:)
wmat=MATMUL(MATMUL(CONJG(TRANSPOSE(mat)),wmatt),mat)
! convert to bsqrtA/|sqrtA|
wmat = MATMUL(MATMUL(ptof,wvecs),ptof)
wmat = MATMUL(MATMUL(ptof,wmat),ptof)
wmatt = MATMUL(MATMUL(ptof,wmatt),ptof)
wvecs = wmat
wtvecs = wmatt
Expand Down
2 changes: 1 addition & 1 deletion rdcon/dcon.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -431,7 +431,7 @@ SUBROUTINE dcon_run
mx0=0
vx0=0
IF (netcdf_out)
$ CALL rdcon_netcdf_out(mx0,mx0,mx0,mx0,vx0,vx0,vx0)
$ CALL rdcon_netcdf_out(mx0,mx0,mx0,mx0,vx0,REAL(vx0,r8),vx0)
ENDIF
IF(mat_flag .OR. ode_flag)DEALLOCATE(amat,bmat,cmat,ipiva,jmat)
IF(bin_euler)CALL bin_close(euler_bin_unit)
Expand Down
12 changes: 7 additions & 5 deletions rdcon/free.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,8 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero)
REAL(r8) :: v1
REAL(r8), DIMENSION(mpert) :: singfac
! eigenvalues are complex for gpec
COMPLEX(r8), DIMENSION(mpert) :: ep,ev,et,tt
COMPLEX(r8), DIMENSION(mpert) :: ep,et,tt
REAL(r8), DIMENSION(mpert) :: ev

REAL(r8), DIMENSION(3*mpert-1) :: rwork
REAL(r8), DIMENSION(2*mpert) :: rwork2
Expand Down Expand Up @@ -219,7 +220,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero)
ev(ipert)=wvt(ipert,ipert)
ENDDO
plasma1=REAL(ep(1))
vacuum1=REAL(ev(1))
vacuum1=ev(1)
total1=REAL(et(1))
c-----------------------------------------------------------------------
c write data for ahb and deallocate.
Expand All @@ -242,15 +243,15 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero)
c-----------------------------------------------------------------------
c write to screen and copy to output.
c-----------------------------------------------------------------------
IF(verbose) WRITE(*,10) REAL(ep(1)),REAL(ev(1)),
IF(verbose) WRITE(*,10) REAL(ep(1)),ev(1),
$ REAL(et(1)),AIMAG(et(1))
c-----------------------------------------------------------------------
c write eigenvalues to file.
c-----------------------------------------------------------------------
message=""
WRITE(out_unit,'(/1x,a)')"Total Energy Eigenvalues:"
WRITE(out_unit,20)
WRITE(out_unit,30)(isol,REAL(ep(isol)),REAL(ev(isol)),
WRITE(out_unit,30)(isol,REAL(ep(isol)),ev(isol),
$ REAL(et(isol)),AIMAG(et(isol)),
$ TRIM(message(isol)),isol=1,mpert)
WRITE(out_unit,20)
Expand All @@ -262,7 +263,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero)
DO isol=1,mpert
WRITE(out_unit,40)
WRITE(out_unit,50)isol,imax(1,isol),REAL(ep(isol)),
$ REAL(ev(isol)),REAL(et(isol)),
$ ev(isol),REAL(et(isol)),
$ AIMAG(et(isol)),TRIM(message(isol))
WRITE(out_unit,60)
WRITE(out_unit,70)(ipert,m(ipert),wt(ipert,isol),
Expand Down Expand Up @@ -310,6 +311,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero)
wp(:,ipert)=vr(:,eindex(mpert+1-ipert))
ep(ipert)=tt(eindex(mpert+1-ipert))
ENDDO
lwork=2*mpert-1
CALL zheev('V','U',mpert,wv,mpert,ev,work,lwork,rwork,info)

c-----------------------------------------------------------------------
Expand Down
28 changes: 28 additions & 0 deletions rdcon/ode_output.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -259,13 +259,20 @@ SUBROUTINE ode_output_get_evals
10 FORMAT(/a,i4,1p,2(a,es10.3))
20 FORMAT(/5x,"i",5x,"eval",7x,"evali",7x,"err"/)
30 FORMAT(i6,1p,3e11.3)
40 FORMAT(/"istep = ",i4,", psifac = ",1p,es10.3,
$ ": skipped, ",a," is singular (zgetrf info = ",i0,")")
c-----------------------------------------------------------------------
c compute plasma response matrix.
c-----------------------------------------------------------------------
temp=CONJG(TRANSPOSE(u(:,1:mpert,1)))
wp=u(:,1:mpert,2)
wp=CONJG(TRANSPOSE(wp))
CALL zgetrf(mpert,mpert,temp,mpert,ipiv,info)
IF(info /= 0)THEN
IF(out_evals .AND. istep > 0)
$ WRITE(evals_out_unit,40)istep,psifac,"u(:,:,1)",info
RETURN
ENDIF
CALL zgetrs('N',mpert,mpert,temp,mpert,ipiv,wp,mpert,info)
wp=(wp+CONJG(TRANSPOSE(wp)))/2
c-----------------------------------------------------------------------
Expand All @@ -277,6 +284,27 @@ SUBROUTINE ode_output_get_evals
key=-ABS(1/evals)
CALL bubble(key,index,1,mpert)
c-----------------------------------------------------------------------
c compute inverse plasma response matrix (zheev destroys wp).
c-----------------------------------------------------------------------
temp=CONJG(TRANSPOSE(u(:,1:mpert,2)))
wp=CONJG(TRANSPOSE(u(:,1:mpert,1)))
CALL zgetrf(mpert,mpert,temp,mpert,ipiv,info)
IF(info /= 0)THEN
IF(out_evals .AND. istep > 0)
$ WRITE(evals_out_unit,40)istep,psifac,"u(:,:,2)",info
RETURN
ENDIF
CALL zgetrs('N',mpert,mpert,temp,mpert,ipiv,wp,mpert,info)
wp=(wp+CONJG(TRANSPOSE(wp)))/2
c-----------------------------------------------------------------------
c compute and sort inverse eigenvalues.
c-----------------------------------------------------------------------
lwork=2*mpert-1
CALL zheev('N','U',mpert,wp,mpert,evalsi,work,lwork,rwork,info)
indexi=(/(ipert,ipert=1,mpert)/)
key=-ABS(evalsi)
CALL bubble(key,indexi,1,mpert)
c-----------------------------------------------------------------------
c write ascii eigenvalues.
c-----------------------------------------------------------------------
IF(out_evals .AND. istep > 0)THEN
Expand Down
8 changes: 4 additions & 4 deletions rdcon/rdcon_netcdf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,8 @@ END SUBROUTINE check
c-----------------------------------------------------------------------
SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)

COMPLEX(r8), DIMENSION(mpert), INTENT(IN) :: ep,ev,et
COMPLEX(r8), DIMENSION(mpert), INTENT(IN) :: ep,et
REAL(r8), DIMENSION(mpert), INTENT(IN) :: ev
COMPLEX(r8), DIMENSION(mpert,mpert), INTENT(IN) :: wp,wv,wt,wt0

INTEGER :: i, ncid,
Expand Down Expand Up @@ -308,7 +309,7 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
CALL check( nf90_def_var(ncid, "W_v_eigenvector", nf90_double,
$ (/m_dim, mo_dim, i_dim/), wv_id) )
CALL check( nf90_def_var(ncid, "W_v_eigenvalue", nf90_double,
$ (/mo_dim, i_dim/), wvv_id) )
$ (/mo_dim/), wvv_id) )
CALL check( nf90_put_att(ncid,wv_id,"long_name",
$ "Vacuum Energy Eigenmodes") )
CALL check( nf90_put_att(ncid,wvv_id,"long_name",
Expand Down Expand Up @@ -459,8 +460,7 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
IF(vac_flag)THEN
CALL check( nf90_put_var(ncid,wv_id,RESHAPE((/REAL(wv),
$ AIMAG(wv)/),(/mpert,mpert,2/))) )
CALL check( nf90_put_var(ncid,wvv_id,RESHAPE((/REAL(ev),
$ AIMAG(ev)/),(/mpert,2/))) )
CALL check( nf90_put_var(ncid,wvv_id,ev) )
CALL check( nf90_put_var(ncid,wt_id,RESHAPE((/REAL(wt),
$ AIMAG(wt)/),(/mpert,mpert,2/))) )
CALL check( nf90_put_var(ncid,wtv_id,RESHAPE((/REAL(et),
Expand Down
4 changes: 3 additions & 1 deletion stride/riccati.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -458,6 +458,7 @@ SUBROUTINE riccati_eig(nstep,t,X,out_unit,bin_unit,dX)

CHARACTER(80) :: format2
REAL(r8), DIMENSION(mpert) :: eigval
COMPLEX(r8), DIMENSION(mpert,mpert) :: Xtemp
REAL(r8), DIMENSION(3*mpert-1) :: rwork_diagnose
COMPLEX(r8), DIMENSION(2*mpert-1) :: work_diagnose
INTEGER, DIMENSION(1) :: jmin
Expand All @@ -473,7 +474,8 @@ SUBROUTINE riccati_eig(nstep,t,X,out_unit,bin_unit,dX)
c preliminaries.
c-----------------------------------------------------------------------
lwork=2*mpert-1
CALL zheev('N','U',mpert,X,mpert,eigval,work_diagnose,
Xtemp=X ! zheev destroys its input, and X is INTENT(IN)
CALL zheev('N','U',mpert,Xtemp,mpert,eigval,work_diagnose,
$ lwork,rwork_diagnose,info)
jmin=MINLOC(ABS(eigval))
crit1=eigval(jmin(1)) * sq%f(3)**2
Expand Down
2 changes: 1 addition & 1 deletion stride/sing.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1118,7 +1118,7 @@ SUBROUTINE sing_derL(startPsi,L)
INTEGER :: ipert,jpert,isol,iqty,info
REAL(r8) :: q
REAL(r8), DIMENSION(mpert) :: singfac
COMPLEX(r8), PARAMETER :: one=1
COMPLEX(r8), PARAMETER :: one=1,zero=0
COMPLEX(r8), DIMENSION(mpert,mpert) :: fmatb,gmatb
COMPLEX(r8), DIMENSION(2*mband+1,mpert) :: kmatb
COMPLEX(r8), DIMENSION(mpert,mpert) :: dTemp
Expand Down
Loading