1 SUBROUTINE check3d_symmetry(a,b,c,ap,bp,cp,mblk,mblkc,nblocks, &
3 USE descriptor_mod, mm_loc=>mm, csrc_loc=>csrc
18 INTEGER,
INTENT(in) :: nblocks, mblk, mblkc
19 REAL(dp),
TARGET,
DIMENSION(mblk,mblkc,nblocks),
INTENT(inout) :: &
21 REAL(dp),
TARGET,
DIMENSION(Locq,Locp,nblocks),
INTENT(inout) :: &
23 REAL(dp),
TARGET,
DIMENSION(Locq,Locp,nblocks-1),
INTENT(inout) ::&
26 REAL(dp),
INTENT(out) :: asym_index
35 REAL(dp),
parameter :: tol = 1.0d-6
37 REAL(dp),
POINTER :: amat(:,:), bmat(:,:), cmat(:,:)
38 REAL(dp),
POINTER :: amatp(:,:), bmatp(:,:), cmatp(:,:)
39 REAL(dp),
DIMENSION(:,:),
POINTER :: Asrc, Bsrc, Csrc
40 REAL(dp),
DIMENSION(:,:),
ALLOCATABLE :: Atmp
41 INTEGER,
DIMENSION(DLEN_) :: descAtmp
46 LOGICAL,
PARAMETER :: check_matp = .false.
49 REAL(dp) :: dnorm, anorm,bnorm,cnorm, alpha,beta
50 REAL(dp),
EXTERNAL :: pdlange
51 REAL(dp),
DIMENSION(1) :: work
53 IF (.NOT.lscalapack)
RETURN
61 ALLOCATE( atmp(locq,locp) )
64 ALLOCATE( atmp(mblk,mblkc) )
90 CALL pdgeadd(
'Notranspose', mm,nn,
91 alpha, asrc,1,1,descatmp,
92 beta, atmp,1,1,descatmp )
96 anorm = pdlange(
'F',mm,nn,atmp,1,1,descatmp,work)
105 CALL pdgeadd(
'Transpose', mm,nn,
106 alpha, asrc,1,1,descatmp,
107 beta, atmp,1,1,descatmp )
112 dnorm = pdlange(
'F', mm,nn, atmp,1,1,descatmp, work)
113 w1 = w1 + dnorm*dnorm
114 w2 = w2 + anorm*anorm
121 IF (k < nblocks)
THEN
138 bnorm = pdlange(
'F',mm,nn,bsrc,1,1,descatmp,work)
139 cnorm = pdlange(
'F',mm,nn,csrc,1,1,descatmp,work)
150 CALL pdgeadd(
'Notranspose', mm,nn,
151 alpha,bsrc,1,1,descatmp,
152 beta,atmp,1,1,descatmp )
162 CALL pdgeadd(
'Transpose', mm,nn,
163 alpha,csrc,1,1,descatmp,
164 beta,atmp,1,1,descatmp )
169 dnorm = pdlange(
'F', mm,nn, atmp,1,1,descatmp, work)
171 w1 = w1 + dnorm*dnorm
172 w2 = w2 + bnorm*bnorm + cnorm*cnorm
181 asym_index = sqrt(w1/w2)
189 END SUBROUTINE check3d_symmetry
192 SUBROUTINE check3d_symmetry_serial(dblk,lblk,ublk,mpol,ntor, &
193 mblk,nblocks,asym_index)
200 INTEGER,
INTENT(in) :: mblk, nblocks, mpol, ntor
201 REAL(dp),
INTENT(out) :: asym_index
202 REAL(dp),
INTENT(in),
DIMENSION(mblk,mblk,nblocks) :: &
207 INTEGER :: ntype, m, n, icol
220 IF (m.EQ.0 .AND. n.LT.0) cycle
222 sum((ublk(icol,:,1:nblocks-1)
223 - lblk(:,icol,2:nblocks))**2)
224 w2 = w2 + sum(ublk(icol,:,1:nblocks-1)**2
225 + lblk(:,icol,2:nblocks)**2)
226 w1 = w1 + sum((dblk(icol,:,:) -
228 w2 = w2 + sum(dblk(icol,:,:)**2)
234 asym_index = sqrt(w1/w2)
239 END SUBROUTINE check3d_symmetry_serial
242 SUBROUTINE checkforces(xc, gc)
243 USE stel_kinds,
ONLY: dp
244 USE stel_constants,
ONLY: zero, twopi
245 USE descriptor_mod,
ONLY: iam
251 USE quantities,
ONLY: fsubsmncf, fsubumnsf, fsubvmnsf, &
252 mpol, ntor, ns, wp0, hs_i
253 USE hessian,
ONLY: l_compute_hessian
254 USE nscalingtools,
ONLY: startglobrow, endglobrow
255 USE v3_utilities,
ONLY:
assert
260 REAL(dp) :: xc(ndims*(1+mpol)*(2*ntor+1),ns), &
261 gc(ndims*(1+mpol)*(2*ntor+1),ns)
265 REAL(dp),
ALLOCATABLE :: gc1(:,:)
266 INTEGER :: ntype, n, m, icol, js, mblk, nt1, icount
267 REAL(dp) :: wplus, wmins, force, eps, fnorm
268 LOGICAL :: lsave, lsave2, lwrite, l_edge, l_s, l_u, l_v
269 CHARACTER*(10) :: force_array(3) = (/
' FORCE-S: ',
' FORCE-U: ',
' FORCE-V: '
277 lsave2= l_compute_hessian
283 l_compute_hessian = .false.
293 ALLOCATE (gc1(mblk,ns))
295 fnorm = (twopi*twopi)*hs_i
298 IF (iam .EQ. 0) print *,
'WRITING CHECK-FORCES FILES FORT.(3000+iam)'
300 WRITE (3000+iam, *)
'fnorm: ', fnorm,
' iam: ', iam
309 radial:
DO icount = 1, 4
311 IF (icount .EQ. 1) js = 1
312 IF (icount .EQ. 2) js = 2
313 IF (icount .EQ. 3) js = ns/2
314 IF (icount .EQ. 4) js = ns
316 lwrite = (startglobrow.LE.js .AND. endglobrow.GE.js)
318 WRITE (3000+iam, *)
'POINT: ', js,
' START: ',
319 startglobrow,
' END: ', endglobrow
335 force = -(wplus-wmins)/(2*eps)
343 IF (ntype .GT. 3) nt1 = ntype-3
344 WRITE (3000+iam,100) m, n, ntype, force, force_array(nt1),
352 IF (lwrite)
WRITE (3000+iam,*)
358 100
FORMAT(1x,
'M: ',i4,
' N: ',i4,
' NTYPE: ',i4,
' GRAD-W: ',1pe14.4, a,
360 CALL assert(icol.EQ.mblk,
'icol != mblk in CheckForces')
362 l_compute_hessian = lsave2
368 END SUBROUTINE checkforces
371 SUBROUTINE tridiag_test (xc, gc, eps)
375 USE stel_kinds,
ONLY: dp
376 USE stel_constants,
ONLY: zero
377 USE quantities,
ONLY: fsubsmncf, fsubumnsf, fsubvmnsf, mpol, ntor, ns
378 USE hessian,
ONLY: l_compute_hessian
382 USE descriptor_mod,
ONLY: iam
383 USE v3_utilities,
ONLY:
assert
388 REAL(dp),
DIMENSION(ndims*(1+mpol)*(2*ntor+1),ns) :: &
390 REAL(dp),
INTENT(IN) :: eps
394 INTEGER :: js, icol, ntype, n, m
395 REAL(dp),
ALLOCATABLE :: gc0(:,:)
396 LOGICAL :: l_save, l_AppSave
407 ALLOCATE (gc0(
SIZE(gc,1),
SIZE(gc,2)))
418 IF (m .eq. 0 .and. n.lt.0) cycle
423 CALL assert(all(gc(:,1:js-2) .EQ. zero),
424 'TRI_DIAGONAL FAILED FOR LOWER BANDS')
425 CALL assert(all(gc(:,js+2:ns) .EQ. zero),
426 'TRI_DIAGONAL FAILED FOR UPPER BANDS')
431 IF (iam .EQ. 0) print *,
'TRI_DIAGONAL TEST PASSED'
436 END SUBROUTINE tridiag_test