diff --git a/dcon/free.f b/dcon/free.f index 6dd86eff3..ceefbb051 100644 --- a/dcon/free.f +++ b/dcon/free.f @@ -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. @@ -252,7 +252,7 @@ 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. @@ -260,7 +260,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out) 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) @@ -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), @@ -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. diff --git a/dcon/ode_output.f b/dcon/ode_output.f index 35b7adf51..60e01414e 100644 --- a/dcon/ode_output.f +++ b/dcon/ode_output.f @@ -273,6 +273,8 @@ 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----------------------------------------------------------------------- @@ -280,6 +282,11 @@ SUBROUTINE ode_output_get_evals 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----------------------------------------------------------------------- @@ -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). +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 diff --git a/gpec/gpdiag.f b/gpec/gpdiag.f index 3d63f52bf..d5b6825ef 100644 --- a/gpec/gpdiag.f +++ b/gpec/gpdiag.f @@ -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 diff --git a/gpec/gpglobal.F b/gpec/gpglobal.F index ce75712dd..a61eeffc5 100644 --- a/gpec/gpglobal.F +++ b/gpec/gpglobal.F @@ -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 @@ -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, diff --git a/gpec/gpout.f b/gpec/gpout.f index fa5ccf325..69036acaa 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -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,*) @@ -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 diff --git a/rdcon/dcon.f b/rdcon/dcon.f index a8e056246..e1874bc50 100644 --- a/rdcon/dcon.f +++ b/rdcon/dcon.f @@ -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) diff --git a/rdcon/free.f b/rdcon/free.f index 2a20399d0..a2e169913 100644 --- a/rdcon/free.f +++ b/rdcon/free.f @@ -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 @@ -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. @@ -242,7 +243,7 @@ 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. @@ -250,7 +251,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero) 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) @@ -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), @@ -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----------------------------------------------------------------------- diff --git a/rdcon/ode_output.f b/rdcon/ode_output.f index 00193f135..75711a88d 100644 --- a/rdcon/ode_output.f +++ b/rdcon/ode_output.f @@ -259,6 +259,8 @@ 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----------------------------------------------------------------------- @@ -266,6 +268,11 @@ SUBROUTINE ode_output_get_evals 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----------------------------------------------------------------------- @@ -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 diff --git a/rdcon/rdcon_netcdf.f b/rdcon/rdcon_netcdf.f index 1b2c52d34..e605886ed 100644 --- a/rdcon/rdcon_netcdf.f +++ b/rdcon/rdcon_netcdf.f @@ -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, @@ -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", @@ -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), diff --git a/stride/riccati.f b/stride/riccati.f index 678a086c0..a77286e4e 100755 --- a/stride/riccati.f +++ b/stride/riccati.f @@ -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 @@ -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 diff --git a/stride/sing.F b/stride/sing.F index 9aaaa4710..b7c4b9c06 100644 --- a/stride/sing.F +++ b/stride/sing.F @@ -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