KTH framework for Nek5000 toolboxes; testing version  0.0.1
gmres.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine uzawa_gmres(res,h1,h2,h2inv,intype,iter)
3 
4 c Solve the pressure equation by right-preconditioned
5 c GMRES iteration.
6 c intype = 0 (steady)
7 c intype = 1 (explicit)
8 c intype = -1 (implicit)
9 
10  include 'SIZE'
11  include 'TOTAL'
12  include 'GMRES'
13  common /ctolpr/ divex
14  common /cprint/ ifprint
15  logical ifprint
16  real res (lx2*ly2*lz2*lelv)
17  real h1 (lx1,ly1,lz1,lelv)
18  real h2 (lx1,ly1,lz1,lelv)
19  real h2inv(lx1,ly1,lz1,lelv)
20 
21  common /scrmg/ wp(lx2,ly2,lz2,lelv)
22 
23  common /ctmp0/ wk1(lgmres),wk2(lgmres)
24  common /cgmres1/ y(lgmres)
25 
26  real alpha, l, temp
27  integer j,m
28 c
29  logical iflag
30  save iflag
31  data iflag /.false./
32  real norm_fac
33  save norm_fac
34 c
35  real*8 etime1,dnekclock
36 c
37  if(.not.iflag) then
38  iflag=.true.
39  call uzawa_gmres_split0(ml_gmres,mu_gmres,bm2,bm2inv,
40  $ lx2*ly2*lz2*nelv)
41  norm_fac = 1./sqrt(volvm2)
42  endif
43 c
44  etime1 = dnekclock()
45  etime_p = 0.
46  divex = 0.
47  iter = 0
48  m = lgmres
49 c
50  call chktcg2(tolps,res,iconv)
51  if (param(21).gt.0.and.tolps.gt.abs(param(21)))
52  $ tolps = abs(param(21))
53 c if (param(21).lt.0) tolps = abs(param(21))
54  if (istep.eq.0) tolps = 1.e-4
55  tolpss = tolps
56 c
57  ntot2 = lx2*ly2*lz2*nelv
58 c
59  iconv = 0
60  call rzero(x_gmres,ntot2)
61 
62  do while(iconv.eq.0.and.iter.lt.100)
63 
64  if(iter.eq.0) then
65  ! -1
66  call col3(r_gmres,ml_gmres,res,ntot2) ! r = L res
67 c call copy(r_gmres,res,ntot2)
68  else
69  !update residual
70  call copy(r_gmres,res,ntot2) ! r = res
71  call cdabdtp(w_gmres,x_gmres,h1,h2,h2inv,intype) ! w = A x
72  call add2s2(r_gmres,w_gmres,-1.,ntot2) ! r = r - w
73  ! -1
74  call col2(r_gmres,ml_gmres,ntot2) ! r = L r
75  endif
76  ! ______
77  gamma_gmres(1) = sqrt(glsc2(r_gmres,r_gmres,ntot2))! gamma = \/ (r,r)
78  ! 1
79  if(iter.eq.0) then
80  div0 = gamma_gmres(1)*norm_fac
81  if (param(21).lt.0) tolpss=abs(param(21))*div0
82  endif
83 
84  !check for lucky convergence
85  rnorm = 0.
86  if(gamma_gmres(1) .eq. 0.) goto 9000
87  temp = 1./gamma_gmres(1)
88  call cmult2(v_gmres(1,1),r_gmres,temp,ntot2)! v = r / gamma
89  ! 1 1
90  do j=1,m
91  iter = iter+1
92  ! -1
93  call col3(w_gmres,mu_gmres,v_gmres(1,j),ntot2) ! w = U v
94  ! j
95 
96  etime2 = dnekclock()
97  if(param(43).eq.1) then
98  call uzprec(z_gmres(1,j),w_gmres,h1,h2,intype,wp)
99  else ! -1
100  call hsmg_solve(z_gmres(1,j),w_gmres) ! z = M w
101 c call copy(z_gmres(1,j),w_gmres,ntot2) ! z = M w
102  endif
103  etime_p = etime_p + dnekclock()-etime2
104 
105  call cdabdtp(w_gmres,z_gmres(1,j), ! w = A z
106  $ h1,h2,h2inv,intype) ! j
107 
108  ! -1
109  call col2(w_gmres,ml_gmres,ntot2) ! w = L w
110 
111 c !modified Gram-Schmidt
112 c do i=1,j
113 c h_gmres(i,j)=glsc2(w_gmres,v_gmres(1,i),ntot2) ! h = (w,v )
114 c ! i,j i
115 c call add2s2(w_gmres,v_gmres(1,i),-h_gmres(i,j),ntot2) ! w = w - h v
116 c enddo ! i,j i
117 
118 
119 c 2-PASS GS, 1st pass:
120 
121  do i=1,j
122  h_gmres(i,j)=vlsc2(w_gmres,v_gmres(1,i),ntot2) ! h = (w,v )
123  enddo ! i,j i
124 
125  call gop(h_gmres(1,j),wk1,'+ ',j) ! sum over P procs
126 
127  do i=1,j
128  call add2s2(w_gmres,v_gmres(1,i),-h_gmres(i,j),ntot2) ! w = w - h v
129  enddo ! i,j i
130 
131 
132 c 2-PASS GS, 2nd pass:
133 c
134 c do i=1,j
135 c wk1(i)=vlsc2(w,v_gmres(1,i),ntot2) ! h = (w,v )
136 c enddo ! i,j i
137 c !
138 c call gop(wk1,wk2,'+ ',j) ! sum over P procs
139 c
140 c do i=1,j
141 c call add2s2(w,v_gmres(1,i),-wk1(i),ntot2) ! w = w - h v
142 c h(i,j) = h(i,j) + wk1(i) ! i,j i
143 c enddo
144 
145 
146  !apply Givens rotations to new column
147  do i=1,j-1
148  temp = h_gmres(i,j)
149  h_gmres(i ,j)= c_gmres(i)*temp
150  $ + s_gmres(i)*h_gmres(i+1,j)
151  h_gmres(i+1,j)= -s_gmres(i)*temp
152  $ + c_gmres(i)*h_gmres(i+1,j)
153  enddo
154  ! ______
155  alpha = sqrt(glsc2(w_gmres,w_gmres,ntot2)) ! alpha = \/ (w,w)
156  rnorm = 0.
157  if(alpha.eq.0.) goto 900 !converged
158  l = sqrt(h_gmres(j,j)*h_gmres(j,j)+alpha*alpha)
159  temp = 1./l
160  c_gmres(j) = h_gmres(j,j) * temp
161  s_gmres(j) = alpha * temp
162  h_gmres(j,j) = l
163  gamma_gmres(j+1) = -s_gmres(j) * gamma_gmres(j)
164  gamma_gmres(j) = c_gmres(j) * gamma_gmres(j)
165 
166 c call outmat(h,m,j,' h ',j)
167 
168  rnorm = abs(gamma_gmres(j+1))*norm_fac
169  ratio = rnorm/div0
170  if (ifprint.and.nio.eq.0)
171  $ write (6,66) iter,tolpss,rnorm,div0,ratio,istep
172  66 format(i5,1p4e12.5,i8,' Divergence')
173 
174 #ifndef FIXITER
175  if (rnorm .lt. tolpss) goto 900 !converged
176 #else
177  if (iter.gt.param(151)-1) goto 900
178 #endif
179  if (j.eq.m) goto 1000 !not converged, restart
180 
181  temp = 1./alpha
182  call cmult2(v_gmres(1,j+1),w_gmres,temp,ntot2) ! v = w / alpha
183  ! j+1
184  enddo
185  900 iconv = 1
186  1000 continue
187  !back substitution
188  ! -1
189  !c = H gamma
190  do k=j,1,-1
191  temp = gamma_gmres(k)
192  do i=j,k+1,-1
193  temp = temp - h_gmres(k,i)*c_gmres(i)
194  enddo
195  c_gmres(k) = temp/h_gmres(k,k)
196  enddo
197  !sum up Arnoldi vectors
198  do i=1,j
199  call add2s2(x_gmres,z_gmres(1,i),c_gmres(i),ntot2)
200  ! x = x + c z
201  ! i i
202  enddo
203 c if(iconv.eq.1) call dbg_write(x,lx2,ly2,lz2,nelv,'esol',3)
204  enddo
205  9000 continue
206 c
207  divex = rnorm
208 c iter = iter - 1
209 c
210 c DIAGNOSTICS
211 c call copy (w,x,ntot2)
212  call ortho (w_gmres) ! Orthogonalize wrt null space, if present
213 c call copy(r,res,ntot2) !r = res
214 c call cdabdtp(r,w,h1,h2,h2inv,intype) ! r = A w
215 c do i=1,ntot2
216 c r(i) = res(i) - r(i) ! r = res - r
217 c enddo
218 c call uzawa_gmres_temp(r,bm2inv,ntot2)
219 c ! ______
220 c gamma(1) = sqrt(glsc2(r,r,ntot2)/volvm2) ! gamma = \/ (r,r)
221 c ! 1
222 c print *, 'GMRES end resid:',gamma(1)
223 c END DIAGNOSTICS
224  call copy(res,x_gmres,ntot2)
225 
226  call ortho (res) ! Orthogonalize wrt null space, if present
227 
228  etime1 = dnekclock()-etime1
229  if (nio.eq.0) write(6,9999) istep,' U-PRES gmres ',
230  & iter,divex,div0,tolpss,etime_p,etime1
231 c call flush_hack
232  9999 format(i11,a,i6,1p5e13.4)
233 
234  return
235  end
236 
237 c-----------------------------------------------------------------------
238 
239  subroutine uzawa_gmres_split0(l,u,b,binv,n)
240  integer n
241  real l(n),u(n),b(n),binv(n)
242  integer i
243  do i=1,n
244  l(i)=sqrt(binv(i))
245  u(i)=sqrt(b(i))
246  if(abs(u(i)*l(i)-1.0).gt.1e-13) print *, i, u(i)*l(i)
247  enddo
248  return
249  end
250 
251 c-----------------------------------------------------------------------
252  subroutine uzawa_gmres_split(l,u,b,binv,n)
253  integer n
254  real l(n),u(n),b(n),binv(n)
255  integer i
256  do i=1,n
257 c l(i)=sqrt(binv(i))
258 c u(i)=sqrt(b(i))
259 
260 c u(i)=sqrt(b(i))
261 c l(i)=1./u(i)
262 
263 c l(i)=sqrt(binv(i))
264  l(i)=1.
265  u(i)=1./l(i)
266 
267 
268 c if(abs(u(i)*l(i)-1.0).gt.1e-13)write(6,*) i,u(i)*l(i),' gmr_sp'
269  enddo
270  return
271  end
272 
273 c-----------------------------------------------------------------------
274  subroutine uzawa_gmres_temp(a,b,n)
275  integer n
276  real a(n),b(n)
277  integer i
278  do i=1,n
279  a(i)=sqrt(b(i))*a(i)
280  enddo
281  return
282  end
283 c-----------------------------------------------------------------------
284  subroutine ax(w,x,h1,h2,n)
285  include 'SIZE'
286  include 'TOTAL'
287 
288 c
289 c w = A*x for pressure iteration
290 c
291 
292  integer n
293  real w(n),x(n),h1(n),h2(n)
294 
295  imsh = 1
296  isd = 1
297  call axhelm (w,x,h1,h2,imsh,isd)
298  call dssum (w,lx1,ly1,lz1)
299  call col2 (w,pmask,n)
300 
301  return
302  end
303 c-----------------------------------------------------------------------
304  subroutine hmh_gmres(res,h1,h2,wt,iter)
305 
306 c Solve the Helmholtz equation by right-preconditioned
307 c GMRES iteration.
308 
309 
310  include 'SIZE'
311  include 'TOTAL'
312  include 'FDMH1'
313  include 'GMRES'
314  common /ctolpr/ divex
315  common /cprint/ ifprint
316  logical ifprint
317  real res (lx1*ly1*lz1*lelv)
318  real h1 (lx1,ly1,lz1,lelv)
319  real h2 (lx1,ly1,lz1,lelv)
320  real wt (lx1,ly1,lz1,lelv)
321 
322  common /scrcg/ d(lx1*ly1*lz1*lelv),wk(lx1*ly1*lz1*lelv)
323 
324  common /cgmres1/ y(lgmres)
325  common /ctmp0/ wk1(lgmres),wk2(lgmres)
326  real alpha, l, temp
327  integer outer
328 
329  logical iflag,if_hyb
330  save iflag,if_hyb
331 c data iflag,if_hyb /.false. , .true. /
332  data iflag,if_hyb /.false. , .false. /
333  real norm_fac
334  save norm_fac
335 
336  real*8 etime1,dnekclock
337 
338 
339  n = lx1*ly1*lz1*nelv
340 
341  etime1 = dnekclock()
342  etime_p = 0.
343  divex = 0.
344  maxit = iter
345  iter = 0
346  m = lgmres
347 
348  if(.not.iflag) then
349  iflag=.true.
350  call uzawa_gmres_split(ml_gmres,mu_gmres,bm1,binvm1,
351  $ lx1*ly1*lz1*nelv)
352  norm_fac = 1./sqrt(volvm1)
353  endif
354 
355  if (param(100).ne.2) call set_fdm_prec_h1b(d,h1,h2,nelv)
356 
357  call chktcg1(tolps,res,h1,h2,pmask,vmult,1,1)
358  if (param(21).gt.0.and.tolps.gt.abs(param(21)))
359  $ tolps = abs(param(21))
360  if (istep.eq.0) tolps = 1.e-4
361  tolpss = tolps
362 c
363  iconv = 0
364  call rzero(x_gmres,n)
365 
366  outer = 0
367  do while (iconv.eq.0)
368 
369  outer = outer+1
370 
371  if(iter.eq.0) then ! -1
372  call col3(r_gmres,ml_gmres,res,n) ! r = L res
373 c call copy(r,res,n)
374  else
375  !update residual
376  call copy (r_gmres,res,n) ! r = res
377  call ax (w_gmres,x_gmres,h1,h2,n) ! w = A x
378  call add2s2(r_gmres,w_gmres,-1.,n) ! r = r - w
379  ! -1
380  call col2(r_gmres,ml_gmres,n) ! r = L r
381  endif
382  ! ______
383  gamma_gmres(1) = sqrt(glsc3(r_gmres,r_gmres,wt,n)) ! gamma = \/ (r,r)
384  ! 1
385  if(iter.eq.0) then
386  div0 = gamma_gmres(1)*norm_fac
387  if (param(21).lt.0) tolpss=abs(param(21))*div0
388  endif
389 
390  !check for lucky convergence
391  rnorm = 0.
392  if(gamma_gmres(1) .eq. 0.) goto 9000
393  temp = 1./gamma_gmres(1)
394  call cmult2(v_gmres(1,1),r_gmres,temp,n) ! v = r / gamma
395  ! 1 1
396  do j=1,m
397  iter = iter+1
398  ! -1
399  call col3(w_gmres,mu_gmres,v_gmres(1,j),n) ! w = U v
400  ! j
401 
402 c . . . . . Overlapping Schwarz + coarse-grid . . . . . . .
403 
404  etime2 = dnekclock()
405 
406 c if (outer.gt.2) if_hyb = .true. ! Slow outer convergence
407  if (ifmgrid) then
408  if (param(40).ge.0 .and. param(40).le.2) then
409  call h1mg_solve(z_gmres(1,j),w_gmres,if_hyb) ! z = M w
410  else if (param(40).eq.3) then
411  call fem_amg_solve(z_gmres(1,j),w_gmres)
412  endif
413  else ! j
414  kfldfdm = ldim+1
415  if (param(100).eq.2) then
416  call h1_overlap_2 (z_gmres(1,j),w_gmres,pmask)
417  else
418  call fdm_h1
419  $ (z_gmres(1,j),w_gmres,d,pmask,vmult,nelv,
420  $ ktype(1,1,kfldfdm),wk)
421  endif
422  call crs_solve_h1 (wk,w_gmres) ! z = M w
423  call add2 (z_gmres(1,j),wk,n) ! j
424  endif
425 
426 
427  call ortho (z_gmres(1,j)) ! Orthogonalize wrt null space, if present
428  etime_p = etime_p + dnekclock()-etime2
429 c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
430 
431 
432  call ax (w_gmres,z_gmres(1,j),h1,h2,n) ! w = A z
433  ! j
434 
435  ! -1
436  call col2(w_gmres,ml_gmres,n) ! w = L w
437 
438 c !modified Gram-Schmidt
439 
440 c do i=1,j
441 c h_gmres(i,j)=glsc3(w_gmres,v_gmres(1,i),wt,n) ! h = (w,v )
442 c ! i,j i
443 
444 c call add2s2(w_gmres,v_gmres(1,i),-h_gmres(i,j),n) ! w = w - h v
445 c enddo ! i,j i
446 
447 c 2-PASS GS, 1st pass:
448 
449  do i=1,j
450  h_gmres(i,j)=vlsc3(w_gmres,v_gmres(1,i),wt,n) ! h = (w,v )
451  enddo ! i,j i
452 
453  call gop(h_gmres(1,j),wk1,'+ ',j) ! sum over P procs
454 
455  do i=1,j
456  call add2s2(w_gmres,v_gmres(1,i),-h_gmres(i,j),n) ! w = w - h v
457  enddo ! i,j i
458 
459 
460 c 2-PASS GS, 2nd pass:
461 c
462 c do i=1,j
463 c wk1(i)=vlsc3(w_gmres,v_gmres(1,i),wt,n) ! h = (w,v )
464 c enddo ! i,j i
465 c !
466 c call gop(wk1,wk2,'+ ',j) ! sum over P procs
467 c
468 c do i=1,j
469 c call add2s2(w_gmres,v_gmres(1,i),-wk1(i),n) ! w = w - h v
470 c h_gmres(i,j) = h_gmres(i,j) + wk1(i) ! i,j i
471 c enddo
472 
473  !apply Givens rotations to new column
474  do i=1,j-1
475  temp = h_gmres(i,j)
476  h_gmres(i ,j)= c_gmres(i)*temp
477  $ + s_gmres(i)*h_gmres(i+1,j)
478  h_gmres(i+1,j)= -s_gmres(i)*temp
479  $ + c_gmres(i)*h_gmres(i+1,j)
480  enddo
481  ! ______
482  alpha = sqrt(glsc3(w_gmres,w_gmres,wt,n)) ! alpha = \/ (w,w)
483  rnorm = 0.
484  if(alpha.eq.0.) goto 900 !converged
485  l = sqrt(h_gmres(j,j)*h_gmres(j,j)+alpha*alpha)
486  temp = 1./l
487  c_gmres(j) = h_gmres(j,j) * temp
488  s_gmres(j) = alpha * temp
489  h_gmres(j,j) = l
490  gamma_gmres(j+1) = -s_gmres(j) * gamma_gmres(j)
491  gamma_gmres(j) = c_gmres(j) * gamma_gmres(j)
492 
493  rnorm = abs(gamma_gmres(j+1))*norm_fac
494  ratio = rnorm/div0
495  if (ifprint.and.nio.eq.0)
496  $ write (6,66) iter,tolpss,rnorm,div0,ratio,istep
497  66 format(i5,1p4e12.5,i8,' Divergence')
498 
499 #ifndef FIXITER
500  if (iter+1.gt.maxit) goto 900
501  if (rnorm .lt. tolpss) goto 900 !converged
502 #else
503  if (iter.gt.param(151)-1) goto 900
504 #endif
505  if (j.eq.m) goto 1000 !not converged, restart
506 
507  temp = 1./alpha
508  call cmult2(v_gmres(1,j+1),w_gmres,temp,n) ! v = w / alpha
509  ! j+1
510  enddo
511  900 iconv = 1
512  1000 continue
513  !back substitution
514  ! -1
515  !c = H gamma
516  do k=j,1,-1
517  temp = gamma_gmres(k)
518  do i=j,k+1,-1
519  temp = temp - h_gmres(k,i)*c_gmres(i)
520  enddo
521  c_gmres(k) = temp/h_gmres(k,k)
522  enddo
523  !sum up Arnoldi vectors
524  do i=1,j
525  call add2s2(x_gmres,z_gmres(1,i),c_gmres(i),n) ! x = x + c z
526  enddo ! i i
527 c if(iconv.eq.1) call dbg_write(x,lx1,ly1,lz1,nelv,'esol',3)
528  enddo
529  9000 continue
530 
531  divex = rnorm
532  call copy(res,x_gmres,n)
533 
534  call ortho (res) ! Orthogonalize wrt null space, if present
535 
536  etime1 = dnekclock()-etime1
537  if (nio.eq.0) write(6,9999) istep,iter,divex,div0,tolpss,etime_p,
538  & etime1,if_hyb
539 c call flush_hack
540  9999 format(4x,i7,' PRES gmres ',4x,i5,1p5e13.4,1x,l4)
541 
542  if (outer.le.2) if_hyb = .false.
543 
544  return
545  end
546 c-----------------------------------------------------------------------
547  subroutine set_overlap2
548 c
549 c Sets up the gather scatter and the SEM operators
550 c
551  include 'SIZE'
552  include 'TOTAL'
553 
554  common /c_is1/ glo_num(lxs*lys*lzs*lelv)
555  integer*8 glo_num
556  common /ivrtx/ vertex((2**ldim)*lelt)
557  common /handle/ gsh_dd
558  integer vertex,gsh_dd
559 
560  mz = ldim-2
561  nx = lx1+2
562  ny = ly1+2
563  nz = lz1+2*mz
564  call setupds_no_crn(gsh_dd,nx,ny,nz,nelv,nelgv,vertex,glo_num)
565  call swap_lengths ! Set up Matrices for FDM
566  call gen_fast_g
567 
568  return
569  end
570 c-----------------------------------------------------------------------
571  subroutine h1_overlap_2(u,v,mask)
572 c
573 c Local overlapping Schwarz solves with overlap of 2
574 c
575  include 'SIZE'
576  include 'TOTAL'
577 
578 
579  common /cwork1/ v1(lxs,lys,lzs,lelt)
580  common /handle/ gsh_dd
581  integer gsh_dd
582 
583  real u(lx1,lx1,lz1,1),v(1),mask(1)
584  integer e
585 
586  n = lx1*ly1*lz1*nelfld(ifield)
587 
588  call dd_swap_vals(v1,v,gsh_dd)
589 
590  iz1 = 0
591  if (if3d) iz1=1
592  do ie=1,nelfld(ifield)
593  do iz=1,lz1
594  do iy=1,ly1
595  do ix=1,lx1
596  u(ix,iy,iz,ie) = v1(ix+1,iy+1,iz+iz1,ie)
597  enddo
598  enddo
599  enddo
600  enddo
601 
602  call dssum (u,lx1,ly1,lz1)
603  call col2 (u,mask,n)
604 
605 
606  return
607  end
608 c-----------------------------------------------------------------------
609  subroutine dd_swap_vals(v1,v0,gsh_dd)
610 
611  include 'SIZE'
612  include 'TOTAL'
613 
614  common /work1/ w1(lxs,lys,lzs) ! work arrarys for locals
615  $ ,w2(lxs,lys,lzs)
616  integer gsh_dd
617 
618  real v1(lxs,lys,lzs,lelt)
619  real v0(lx1,ly1,lz1,lelt)
620 
621  integer e
622 
623  mz = ldim-2
624 
625  nx = lx1+2
626  ny = ly1+2
627  nz = lz1+2*mz
628 
629  n = lx1*ly1*lz1*nelv
630  m = (lx1+2)*(ly1+2)*(lz1+2*mz)*nelv
631 
632  do e=1,nelv
633  call rzero_g (v1,e,nx,ny,nz)
634  call fill_interior_g(v1,v0,e,lx1,lz1,mz,nelv) ! v0 --> v1(int)
635  call dface_ext_g (v1,2,e,nx,nz) ! v1(int) --> v1(face)
636  enddo
637 c
638 c ~ ~T
639 c This is the Q Q part
640 c
641  call fgslib_gs_op(gsh_dd,v1,1,1,0) ! 1 ==> + ! swap v1 & add vals
642 
643  do e =1,nelv
644  call dface_add1si_g (v1,-1.,2,e,nx,nz)
645  call fastdm1_g (v1(1,1,1,e),e,w1,w2)
646  call s_face_to_int2_g(v1,-1.,2,e,nx,nz)
647  enddo
648 c
649 c Exchange/add elemental solutions
650 c
651  call fgslib_gs_op (gsh_dd,v1,1,1,0)
652  do e =1,nelv
653  call s_face_to_int2_g(v1,1.,2,e,nx,nz)
654  enddo
655 
656  return
657  end
658 c-----------------------------------------------------------------------
659  subroutine gen_fast_g
660 c
661 c Generate fast diagonalization matrices for each element
662 c
663  include 'SIZE'
664  include 'INPUT'
665  include 'PARALLEL'
666  include 'SOLN'
667  include 'WZ'
668 
669  parameter(lxss=lxs*lxs)
670  common /fastg/ sr(lxss,2,lelv),ss(lxss,2,lelv),st(lxss,2,lelv)
671  $ , df(lxs*lys*lzs,lelv)
672 
673  common /ctmpf/ lr(2*lx1+4),ls(2*lx1+4),lt(2*lx1+4)
674  $ , llr(lelt),lls(lelt),llt(lelt)
675  $ , lmr(lelt),lms(lelt),lmt(lelt)
676  $ , lrr(lelt),lrs(lelt),lrt(lelt)
677  real lr ,ls ,lt
678  real llr,lls,llt
679  real lmr,lms,lmt
680  real lrr,lrs,lrt
681 
682  integer lbr,rbr,lbs,rbs,lbt,rbt
683 
684 
685  call load_semhat_weighted ! Fills the SEMHAT arrays
686 
687 
688  ierr = 0
689  do ie=1,nelv
690 
691  call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,3,ierr)
692 c
693 c Set up matrices for each element.
694 c
695  call set_up_fast_1d_sem_g( sr(1,1,ie),lr,nr ,lbr,rbr
696  $ ,llr(ie),lmr(ie),lrr(ie),ie)
697  call set_up_fast_1d_sem_g( ss(1,1,ie),ls,ns ,lbs,rbs
698  $ ,lls(ie),lms(ie),lrs(ie),ie)
699  if (if3d) then
700  call set_up_fast_1d_sem_g( st(1,1,ie),lt,nt ,lbt,rbt
701  $ ,llt(ie),lmt(ie),lrt(ie),ie)
702  endif
703 c
704 c Set up diagonal inverse
705 c
706  if (if3d) then
707  eps = 1.e-5 * (vlmax(lr(2),nr-2)
708  $ + vlmax(ls(2),ns-2) + vlmax(lt(2),nt-2))
709  l = 1
710  do k=1,nt
711  do j=1,ns
712  do i=1,nr
713  diag = lr(i) + ls(j) + lt(k)
714  if (diag.gt.eps) then
715  df(l,ie) = 1.0/diag
716  else
717 c write(6,3) ie,'Reset Eig in gen fast:',i,j,k,l
718 c $ ,eps,diag,lr(i),ls(j),lt(k)
719 c 3 format(i6,1x,a21,4i5,1p5e12.4)
720  df(l,ie) = 0.0
721  endif
722  l = l+1
723  enddo
724  enddo
725  enddo
726  else
727  eps = 1.e-5*(vlmax(lr(2),nr-2) + vlmax(ls(2),ns-2))
728  l = 1
729  do j=1,ns
730  do i=1,nr
731  diag = lr(i) + ls(j)
732 
733  if (diag.gt.eps) then
734  df(l,ie) = 1.0/diag
735 
736  else
737 c write(6,2) ie,'Reset Eig in gen fast:',i,j,l
738 c $ ,eps,diag,lr(i),ls(j)
739 c 2 format(i6,1x,a21,3i5,1p4e12.4)
740  df(l,ie) = 0.0
741  endif
742  l = l+1
743  enddo
744  enddo
745  endif
746 c
747 c Next element ....
748 c
749  enddo
750 
751  return
752  end
753 c-----------------------------------------------------------------------
754  subroutine set_up_fast_1d_sem_g(s,lam,n,lbc,rbc,ll,lm,lr,ie)
755  include 'SIZE'
756  include 'SEMHAT'
757 
758  parameter(lr3=2*lxs*lxs)
759  common /fast1dsem/ g(lr3),w(lr3)
760 
761  real g,w
762  real s(1),lam(1),ll,lm,lr
763  integer lbc,rbc
764 
765  integer bb0,bb1,eb0,eb1,n,n1
766  logical l,r
767 
768  n=lx1-1
769  !bcs on E are from normal vel component
770  if(lbc.eq.2 .or. lbc.eq.3) then !wall,sym - Neumann
771  eb0=0
772  bb0=0
773  else !outflow,element - Dirichlet
774  eb0=1
775  bb0=1
776  endif
777  if(rbc.eq.2 .or. rbc.eq.3) then !wall,sym - Neumann
778  eb1=n
779  bb1=n
780  else !outflow,element - Dirichlet
781  eb1=n-1
782  bb1=n-1
783  endif
784  l = (lbc.eq.0)
785  r = (rbc.eq.0)
786 
787 c calculate A tilde operator
788  call set_up_fast_1d_sem_op_a(s,eb0,eb1,l,r,ll,lm,lr,ah)
789 c calculate B tilde operator
790  call set_up_fast_1d_sem_op_b(g,bb0,bb1,l,r,ll,lm,lr,bh)
791 
792  n=n+3
793 c call outmat (s,n,n,' A ',ie)
794 c call outmat (g,n,n,' B ',ie)
795  call generalev(s,g,lam,n,w)
796  if(.not.l) call row_zero(s,n,n,1)
797  if(.not.r) call row_zero(s,n,n,n)
798  call transpose(s(n*n+1),n,s,n) ! compute the transpose of s
799 c call outmat (s,n,n,' S ',ie)
800 c call outmat (s(n*n+1),n,n,' St ',1)
801 c call exitt
802 
803 
804  return
805  end
806 c-----------------------------------------------------------------------
807  subroutine set_up_fast_1d_sem_op_a(g,b0,b1,l
808  $ ,r,ll,lm,lr,ah)
809 c -1 T
810 c G = J B J
811 c
812 c gives the inexact restriction of this matrix to
813 c an element plus two node on either side
814 c
815 c g - the output matrix
816 c b0, b1 - the range for Bhat indices for the element
817 c (enforces boundary conditions)
818 c l, r - whether there is a left or right neighbor
819 c ll,lm,lr - lengths of left, middle, and right elements
820 c ah - hat matrix for A or B
821 c
822 c result is inexact because:
823 c neighbor's boundary condition at far end unknown
824 c length of neighbor's neighbor unknown
825 c (these contribs should be small for large N and
826 c elements of nearly equal size)
827 c
828  include 'SIZE'
829  real g(0:lx1+1,0:lx1+1)
830  real ah(0:lx1-1,0:lx1-1)
831  real ll,lm,lr
832  integer b0,b1
833  logical l,r
834 
835  real bl,bm,br
836  integer n
837 
838  n =lx1-1
839 
840 c
841 c compute the weight of A hat
842 c
843  if (l) bl = 2. /ll
844  bm = 2. /lm
845  if (r) br = 2. /lr
846 
847  call rzero(g,(n+3)*(n+3))
848  do j=1,n+1
849  do i=1,n+1
850  g(i,j) = ah(i-1,j-1)*bm
851  enddo
852  enddo
853 
854  if (l) then
855  g(0,0) = g(0,0) + bl*ah(n-1,n-1)
856  g(0,1) = g(0,1) + bl*ah(n-1,n )
857  g(1,0) = g(1,0) + bl*ah(n ,n-1)
858  g(1,1) = g(1,1) + bl*ah(n ,n )
859  elseif (b0.eq.0) then !Neumann BC
860  g(0,0) = 1.
861  else !Dirichlet BC
862  g(0,0) = 1.
863  g(1,1) = 1.
864  do i=2,n+1
865  g(i,1) = 0.
866  g(1,i) = 0.
867  enddo
868  endif
869 
870  if (r) then
871  g(n+1,n+1) = g(n+1,n+1) + br*ah(0,0)
872  g(n+1,n+2) = g(n+1,n+2) + br*ah(0,1)
873  g(n+2,n+1) = g(n+2,n+1) + br*ah(0,1)
874  g(n+2,n+2) = g(n+2,n+2) + br*ah(1,1)
875  elseif (b1.eq.n) then !Neumann BC
876  g(n+2,n+2)=1.
877  else !Dirichlet BC
878  g(n+2,n+2) = 1.
879  g(n+1,n+1) = 1.
880  do i=1,n
881  g(i,n+1) = 0.
882  g(n+1,i) = 0.
883  enddo
884  endif
885 
886  return
887  end
888 c-----------------------------------------------------------------------
889  subroutine set_up_fast_1d_sem_op_b(g,b0,b1,l
890  $ ,r,ll,lm,lr,bh)
891 c -1 T
892 c G = J B J
893 c
894 c gives the inexact restriction of this matrix to
895 c an element plus two node on either side
896 c
897 c g - the output matrix
898 c b0, b1 - the range for Bhat indices for the element
899 c (enforces boundary conditions)
900 c l, r - whether there is a left or right neighbor
901 c ll,lm,lr - lengths of left, middle, and right elements
902 c bh - hat matrix for B
903 c
904 c result is inexact because:
905 c neighbor's boundary condition at far end unknown
906 c length of neighbor's neighbor unknown
907 c (these contribs should be small for large N and
908 c elements of nearly equal size)
909 c
910  include 'SIZE'
911  real g(0:lx1+1,0:lx1+1)
912  real bh(0:lx1-1)
913  real ll,lm,lr
914  integer b0,b1
915  logical l,r
916 
917  real bl,bm,br
918  integer n
919 
920  n =lx1-1
921 c
922 c compute the weight of B hat
923 c
924  if (l) bl = ll / 2.
925  bm = lm / 2.
926  if (r) br = lr / 2.
927 
928  call rzero(g,(n+3)*(n+3))
929  do i=1,n+1
930  g(i,i) = bh(i-1)*bm
931  enddo
932 
933  if (l) then
934  g(0,0) = g(0,0) + bl*bh(n-1)
935  g(1,1) = g(1,1) + bl*bh(n )
936  elseif (b0.eq.0) then !Neumann BC
937  g(0,0) = 1.
938  else !Dirichlet BC
939  g(0,0) = 1.
940  g(1,1) = 1.
941  do i=2,n+1
942  g(i,1) = 0.
943  g(1,i) = 0.
944  enddo
945  endif
946 
947  if (r) then
948  g(n+1,n+1) = g(n+1,n+1) + br*bh(0)
949  g(n+2,n+2) = g(n+2,n+2) + br*bh(1)
950  elseif (b1.eq.n) then !Neumann BC
951  g(n+2,n+2)=1.
952  else !Dirichlet BC
953  g(n+2,n+2) = 1.
954  g(n+1,n+1) = 1.
955  do i=1,n
956  g(i,n+1) = 0.
957  g(n+1,i) = 0.
958  enddo
959  endif
960 
961  return
962  end
963 c-----------------------------------------------------------------------
964  subroutine fill_interior_g(v1,v,e,nx,nz,iz1,nel)
965 
966  real v1(nx+2,nx+2,nz+2*iz1,nel) ! iz1=ldim-2
967  real v (nx ,nx ,nz ,nel)
968  integer e
969 
970  ny = nx
971 
972  do iz=1,nz
973  do iy=1,ny
974  do ix=1,nx
975  v1(ix+1,iy+1,iz+iz1,e) = v(ix,iy,iz,e)
976  enddo
977  enddo
978  enddo
979 
980  return
981  end
982 c-----------------------------------------------------------------------
983  subroutine dface_ext_g(x,t,e,nx,nz)
984 c Extend interior to face of element
985 
986  include 'SIZE'
987  include 'INPUT'
988  real x(nx,nx,nz,1)
989  integer e,t
990 
991  ny = nx
992 
993  if (if3d) then
994  do iz=2,nz-1
995  do ix=2,nx-1
996  x(ix,1 ,iz,e) = x(ix, 1+t,iz,e)
997  x(ix,ny,iz,e) = x(ix,ny-t,iz,e)
998  enddo
999  enddo
1000 
1001  do iz=2,nz-1
1002  do iy=2,ny-1
1003  x(1 ,iy,iz,e) = x( 1+t,iy,iz,e)
1004  x(nx,iy,iz,e) = x(nx-t,iy,iz,e)
1005  enddo
1006  enddo
1007 
1008  do iy=2,ny-1
1009  do ix=2,nx-1
1010  x(ix,iy,1 ,e) = x(ix,iy, 1+t,e)
1011  x(ix,iy,nz,e) = x(ix,iy,nz-t,e)
1012  enddo
1013  enddo
1014  else
1015  do ix=2,nx-1
1016  x(ix,1 ,1,e) = x(ix, 1+t,1,e)
1017  x(ix,ny,1,e) = x(ix,ny-t,1,e)
1018  enddo
1019  do iy=2,ny-1
1020  x(1 ,iy,1,e) = x( 1+t,iy,1,e)
1021  x(nx,iy,1,e) = x(nx-t,iy,1,e)
1022  enddo
1023  endif
1024 
1025  return
1026  end
1027 c-----------------------------------------------------------------------
1028  subroutine dface_add1si_g(x,c,t,e,nx,nz)
1029 c Scale interior and add to face of element
1030 
1031  include 'SIZE'
1032  include 'INPUT'
1033  real x(nx,nx,nz,1)
1034 
1035  integer e,t
1036 
1037  ny = nx
1038 
1039  if (if3d) then
1040 
1041  do iz=2,nz-1
1042  do ix=2,nx-1
1043  x(ix,1 ,iz,e) = x(ix,1 ,iz,e) + c*x(ix, 1+t,iz,e)
1044  x(ix,ny,iz,e) = x(ix,ny,iz,e) + c*x(ix,ny-t,iz,e)
1045  enddo
1046  enddo
1047 
1048  do iz=2,nz-1
1049  do iy=2,ny-1
1050  x(1 ,iy,iz,e) = x(1 ,iy,iz,e) + c*x( 1+t,iy,iz,e)
1051  x(nx,iy,iz,e) = x(nx,iy,iz,e) + c*x(nx-t,iy,iz,e)
1052  enddo
1053  enddo
1054 
1055  do iy=2,ny-1
1056  do ix=2,nx-1
1057  x(ix,iy,1 ,e) = x(ix,iy,1 ,e) + c*x(ix,iy, 1+t,e)
1058  x(ix,iy,nz,e) = x(ix,iy,nz,e) + c*x(ix,iy,nz-t,e)
1059  enddo
1060  enddo
1061 
1062  else ! 2D
1063 
1064  do ix=2,nx-1
1065  x(ix,1 ,1,e) = x(ix,1 ,1,e) + c*x(ix, 1+t,1,e)
1066  x(ix,ny,1,e) = x(ix,ny,1,e) + c*x(ix,ny-t,1,e)
1067  enddo
1068  do iy=2,ny-1
1069  x(1 ,iy,1,e) = x(1 ,iy,1,e) + c*x( 1+t,iy,1,e)
1070  x(nx ,iy,1,e) = x(nx,iy,1,e) + c*x(nx-t,iy,1,e)
1071  enddo
1072 
1073  endif
1074 
1075  return
1076  end
1077 c-----------------------------------------------------------------------
1078  subroutine fastdm1_g(R,ie,w1,w2)
1079 c
1080 c Fast diagonalization solver for FEM on mesh 1
1081 c
1082  include 'SIZE'
1083 
1084  parameter(lxss=lxs*lxs)
1085  common /fastg/ sr(lxss,2,lelv),ss(lxss,2,lelv),st(lxss,2,lelv)
1086  $ , df(lxs*lys*lzs,lelv)
1087 
1088  parameter(lxyz = lxs*lys*lzs)
1089 
1090  real r(1),w1(1),w2(1)
1091 
1092  nx = lx1+2
1093 c
1094 c T
1095 c S r
1096  call tensr3 (w1,nx,r ,nx,sr(1,2,ie),ss(1,1,ie),st(1,1,ie),w2)
1097 c
1098 c
1099 c -1 T
1100 c D S r
1101 c
1102  call col2 (w1,df(1,ie),lxyz)
1103 c
1104 c
1105 c -1 T
1106 c S D S r
1107 c
1108  call tensr3 (r ,nx,w1,nx,sr(1,1,ie),ss(1,2,ie),st(1,2,ie),w2)
1109 
1110  return
1111  end
1112 c-----------------------------------------------------------------------
1113  subroutine s_face_to_int2_g(x,c,t,e,nx,nz)
1114 c
1115 c Scale face and add to interior of element
1116 c
1117  include 'SIZE'
1118  include 'INPUT'
1119  real x(nx,nx,nz,1)
1120  integer t,e
1121 
1122  ny=nx
1123 
1124 
1125  if (if3d) then
1126 
1127  do iz=2,nz-1
1128  do ix=2,nx-1
1129  x(ix, 1+t,iz,e) = c*x(ix,1 ,iz,e) + x(ix, 1+t,iz,e)
1130  x(ix,ny-t,iz,e) = c*x(ix,ny,iz,e) + x(ix,ny-t,iz,e)
1131  enddo
1132  enddo
1133 
1134  do iz=2,nz-1
1135  do iy=2,ny-1
1136  x( 1+t,iy,iz,e) = c*x(1 ,iy,iz,e) + x( 1+t,iy,iz,e)
1137  x(nx-t,iy,iz,e) = c*x(nx,iy,iz,e) + x(nx-t,iy,iz,e)
1138  enddo
1139  enddo
1140 
1141  do iy=2,ny-1
1142  do ix=2,nx-1
1143  x(ix,iy, 1+t,e) = c*x(ix,iy,1 ,e) + x(ix,iy, 1+t,e)
1144  x(ix,iy,nz-t,e) = c*x(ix,iy,nz,e) + x(ix,iy,nz-t,e)
1145  enddo
1146  enddo
1147 
1148  else
1149 c 2D
1150  do ix=2,nx-1
1151  x(ix, 1+t,1,e) = c*x(ix,1 ,1,e) + x(ix, 1+t,1,e)
1152  x(ix,ny-t,1,e) = c*x(ix,ny,1,e) + x(ix,ny-t,1,e)
1153  enddo
1154  do iy=2,ny-1
1155  x( 1+t,iy,1,e) = c*x(1 ,iy,1,e) + x( 1+t,iy,1,e)
1156  x(nx-t,iy,1,e) = c*x(nx,iy,1,e) + x(nx-t,iy,1,e)
1157  enddo
1158  endif
1159 
1160  return
1161  end
1162 c-----------------------------------------------------------------------
1163  subroutine outfldr_g(x,txt10,nx,nz,ichk)
1164  include 'SIZE'
1165  include 'TSTEP'
1166  real x(nx,nx,nz,lelt)
1167  character*10 txt10
1168 
1169  integer idum,e
1170  save idum
1171  data idum /3/
1172 
1173  if (idum.lt.0) return
1174 
1175  ny = nx
1176 
1177 
1178  mtot = nx*ny*nz*nelv
1179  if (nx.gt.8.or.nelv.gt.16) return
1180  xmin = glmin(x,mtot)
1181  xmax = glmax(x,mtot)
1182 
1183  nell = nelt
1184  rnel = nell
1185  snel = sqrt(rnel)+.1
1186  ne = snel
1187  ne1 = nell-ne+1
1188  k = 1
1189  do ie=1,1
1190  ne = 0
1191  write(6,116) txt10,k,ie,xmin,xmax,istep,time,nx
1192  do l=12,0,-4
1193  write(6,117)
1194  do j=ny,1,-1
1195  if (nx.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1196  if (nx.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1197  if (nx.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1198  if (nx.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1199  if (nx.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1200  if (nx.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1201  if (nx.eq.8) write(6,118) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1202  enddo
1203  enddo
1204  enddo
1205 
1206  102 FORMAT(4(2f9.5,2x))
1207  103 FORMAT(4(3f9.5,2x))
1208  104 FORMAT(4(4f7.3,2x))
1209  105 FORMAT(5f9.5,10x,5f9.5)
1210  106 FORMAT(6f7.1,5x,6f7.1,5x,6f7.1,5x,6f7.1)
1211  107 FORMAT(7f8.4,5x,7f8.4)
1212  108 FORMAT(8f8.4,4x,8f8.4)
1213  118 FORMAT(8f12.9)
1214 
1215  116 FORMAT( /,5x,' ^ ',/,
1216  $ 5x,' Y | ',/,
1217  $ 5x,' | ',a10,/,
1218  $ 5x,' +----> ','Plane = ',i2,'/',i2,2x,2e12.4,/,
1219  $ 5x,' X ','Step =',i9,f15.5,i2)
1220  117 FORMAT(' ')
1221 
1222  if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
1223  return
1224  end
1225 c-----------------------------------------------------------------------
1226  subroutine outfldi_g(x,txt10,nx,nz,ichk)
1227  include 'SIZE'
1228  include 'TSTEP'
1229  integer x(nx,nx,nz,lelt)
1230  character*10 txt10
1231 
1232  integer idum,e,xmin,xmax
1233  save idum
1234  data idum /3/
1235 
1236  if (idum.lt.0) return
1237 
1238  ny = nx
1239 
1240 
1241  mtot = nx*ny*nz*nelv
1242  if (nx.gt.8.or.nelv.gt.16) return
1243  xmin = iglmin(x,mtot)
1244  xmax = iglmax(x,mtot)
1245 
1246  nell = nelt
1247  rnel = nell
1248  snel = sqrt(rnel)+.1
1249  ne = snel
1250  ne1 = nell-ne+1
1251  k = 1
1252  do ie=1,1
1253  ne = 0
1254  write(6,116) txt10,k,ie,xmin,xmax,istep,time,nx
1255  do l=12,0,-4
1256  write(6,117)
1257  do j=ny,1,-1
1258  if (nx.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1259  if (nx.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1260  if (nx.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1261  if (nx.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1262  if (nx.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1263  if (nx.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1264  if (nx.eq.8) write(6,118) ((x(i,j,k,e+l),i=1,nx),e=1,4)
1265  enddo
1266  enddo
1267  enddo
1268 
1269  102 FORMAT(4(2i9,2x))
1270  103 FORMAT(4(3i9,2x))
1271  104 FORMAT(4(4i7,2x))
1272  105 FORMAT(5i9,10x,5i9)
1273  106 FORMAT(6i7,5x,6i7,5x,6i7,5x,6i7)
1274  107 FORMAT(7i8,5x,7i8)
1275  108 FORMAT(8i8,4x,8i8)
1276  118 FORMAT(8i12)
1277 
1278  116 FORMAT( /,5x,' ^ ',/,
1279  $ 5x,' Y | ',/,
1280  $ 5x,' | ',a10,/,
1281  $ 5x,' +----> ','Plane = ',i2,'/',i2,2x,2i12,/,
1282  $ 5x,' X ','Step =',i9,f15.5,i2)
1283  117 FORMAT(' ')
1284 
1285  if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
1286  return
1287  end
1288 c-----------------------------------------------------------------------
1289  subroutine setupds_no_crn(gs_h,nx,ny,nz,nel,melg,vertex,glo_num)
1290  include 'SIZE'
1291  include 'INPUT'
1292  include 'PARALLEL'
1293  include 'NONCON'
1294  integer gs_h,vertex(1),e
1295  integer*8 ngv,glo_num(nx,ny,nz,nel)
1296 
1297 
1298  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1299 
1300 c set up the global numbering
1301  call set_vert(glo_num,ngv,nx,nel,vertex,.false.)
1302 
1303 c zero out corners
1304  mz1 = max(1,nz-1)
1305  do e=1,nel
1306  do k=1,nz,mz1
1307  do j=1,ny,ny-1
1308  do i=1,nx,nx-1
1309  glo_num(i,j,k,e) = 0
1310  enddo
1311  enddo
1312  enddo
1313  enddo
1314 
1315 
1316  ntot = nx*ny*nz*nel
1317 
1318  t0 = dnekclock()
1319  call fgslib_gs_setup(gs_h,glo_num,ntot,nekcomm,mp) ! initialize gs code
1320  t1 = dnekclock()
1321 
1322  et = t1-t0
1323 
1324 c call gs_chkr(glo_num)
1325 
1326  if (nio.eq.0) then
1327  write(6,1) et,nx,nel,ntot,ngv,gs_h
1328  1 format(' gs_init time',1pe11.4,' seconds ',i3,4i10)
1329  endif
1330 c
1331  return
1332  end
1333 c-----------------------------------------------------------------------
1334  subroutine rzero_g(a,e,nx,ny,nz)
1335 
1336  real a(nx,ny,nz,1)
1337  integer e
1338 
1339  do i=1,nx
1340  do j=1,ny
1341  do k=1,nz
1342  a(i,j,k,e) = 0.0
1343  enddo
1344  enddo
1345  enddo
1346 
1347  return
1348  END
1349 c-----------------------------------------------------------------------
1350 
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
checkit
subroutine checkit(idum)
Definition: connect1.f:1263
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
row_zero
subroutine row_zero(a, m, n, e)
Definition: fast3d.f:1617
swap_lengths
subroutine swap_lengths
Definition: fast3d.f:1541
load_semhat_weighted
subroutine load_semhat_weighted
Definition: fast3d.f:1182
get_fast_bc
subroutine get_fast_bc(lbr, rbr, lbs, rbs, lbt, rbt, e, bsym, ierr)
Definition: fast3d.f:803
tensr3
subroutine tensr3(v, nv, u, nu, A, Bt, Ct, w)
Definition: fasts.f:126
fem_amg_solve
#define fem_amg_solve
Definition: fem_amg_preco.c:9
fill_interior_g
subroutine fill_interior_g(v1, v, e, nx, nz, iz1, nel)
Definition: gmres.f:965
set_up_fast_1d_sem_op_b
subroutine set_up_fast_1d_sem_op_b(g, b0, b1, l, r, ll, lm, lr, bh)
Definition: gmres.f:891
dface_ext_g
subroutine dface_ext_g(x, t, e, nx, nz)
Definition: gmres.f:984
h1_overlap_2
subroutine h1_overlap_2(u, v, mask)
Definition: gmres.f:572
gen_fast_g
subroutine gen_fast_g
Definition: gmres.f:660
outfldr_g
subroutine outfldr_g(x, txt10, nx, nz, ichk)
Definition: gmres.f:1164
hmh_gmres
subroutine hmh_gmres(res, h1, h2, wt, iter)
Definition: gmres.f:305
dface_add1si_g
subroutine dface_add1si_g(x, c, t, e, nx, nz)
Definition: gmres.f:1029
s_face_to_int2_g
subroutine s_face_to_int2_g(x, c, t, e, nx, nz)
Definition: gmres.f:1114
rzero_g
subroutine rzero_g(a, e, nx, ny, nz)
Definition: gmres.f:1335
dd_swap_vals
subroutine dd_swap_vals(v1, v0, gsh_dd)
Definition: gmres.f:610
uzawa_gmres
subroutine uzawa_gmres(res, h1, h2, h2inv, intype, iter)
Definition: gmres.f:3
uzawa_gmres_split0
subroutine uzawa_gmres_split0(l, u, b, binv, n)
Definition: gmres.f:240
setupds_no_crn
subroutine setupds_no_crn(gs_h, nx, ny, nz, nel, melg, vertex, glo_num)
Definition: gmres.f:1290
fastdm1_g
subroutine fastdm1_g(R, ie, w1, w2)
Definition: gmres.f:1079
set_up_fast_1d_sem_op_a
subroutine set_up_fast_1d_sem_op_a(g, b0, b1, l, r, ll, lm, lr, ah)
Definition: gmres.f:809
set_up_fast_1d_sem_g
subroutine set_up_fast_1d_sem_g(s, lam, n, lbc, rbc, ll, lm, lr, ie)
Definition: gmres.f:755
ax
subroutine ax(w, x, h1, h2, n)
Definition: gmres.f:285
set_overlap2
subroutine set_overlap2
Definition: gmres.f:548
uzawa_gmres_temp
subroutine uzawa_gmres_temp(a, b, n)
Definition: gmres.f:275
uzawa_gmres_split
subroutine uzawa_gmres_split(l, u, b, binv, n)
Definition: gmres.f:253
outfldi_g
subroutine outfldi_g(x, txt10, nx, nz, ichk)
Definition: gmres.f:1227
fdm_h1
subroutine fdm_h1(z, r, d, mask, mult, nel, kt, rr)
Definition: hmholtz.f:937
chktcg1
subroutine chktcg1(tol, res, h1, h2, mask, mult, imesh, isd)
Definition: hmholtz.f:528
generalev
subroutine generalev(a, b, lam, n, w)
Definition: hmholtz.f:1369
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
set_fdm_prec_h1b
subroutine set_fdm_prec_h1b(d, h1, h2, nel)
Definition: hmholtz.f:1274
h1mg_solve
subroutine h1mg_solve(z, rhs, if_hybrid)
Definition: hsmg.f:1856
hsmg_solve
subroutine hsmg_solve(e, r)
Definition: hsmg.f:1377
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
glmin
function glmin(a, n)
Definition: math.f:973
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
glsc2
function glsc2(x, y, n)
Definition: math.f:794
glsc3
function glsc3(a, b, mult, n)
Definition: math.f:776
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
iglmax
function iglmax(a, n)
Definition: math.f:913
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
iglmin
function iglmin(a, n)
Definition: math.f:900
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
uzprec
subroutine uzprec(rpcg, rcg, h1m1, h2m1, intype, wp)
Definition: navier1.f:899
ortho
subroutine ortho(respr)
Definition: navier1.f:224
cdabdtp
subroutine cdabdtp(ap, wp, h1, h2, h2inv, intype)
Definition: navier1.f:259
chktcg2
subroutine chktcg2(tol, res, iconv)
Definition: navier1.f:1090
vlsc3
function vlsc3(X, Y, B, N)
Definition: navier4.f:324
crs_solve_h1
subroutine crs_solve_h1(uf, vf)
Definition: navier8.f:1341
set_vert
subroutine set_vert(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:4
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341