examples/smpi/NAS/LU/ssor.f

   1 c---------------------------------------------------------------------
   2 c---------------------------------------------------------------------
   3
   4       subroutine ssor(niter)
   5
   6 c---------------------------------------------------------------------
   7 c---------------------------------------------------------------------
   8
   9 c---------------------------------------------------------------------
  10 c   to perform pseudo-time stepping SSOR iterations
  11 c   for five nonlinear pde's.
  12 c---------------------------------------------------------------------
  13
  14       implicit none
  15       integer  niter
  16
  17       include 'mpinpb.h'
  18       include 'applu.incl'
  19
  20 c---------------------------------------------------------------------
  21 c  local variables
  22 c---------------------------------------------------------------------
  23       integer i, j, k, m
  24       integer istep
  25       double precision  tmp
  26       double precision  delunm(5), tv(5,isiz1,isiz2)
  27
  28       external timer_read
  29       double precision wtime, timer_read
  30
  31       integer IERROR
  32
  33
  34       ROOT = 0
  35
  36 c---------------------------------------------------------------------
  37 c   begin pseudo-time stepping iterations
  38 c---------------------------------------------------------------------
  39       tmp = 1.0d+00 / ( omega * ( 2.0d+00 - omega ) )
  40
  41 c---------------------------------------------------------------------
  42 c   initialize a,b,c,d to zero (guarantees that page tables have been
  43 c   formed, if applicable on given architecture, before timestepping).
  44 c---------------------------------------------------------------------
  45       do m=1,isiz2
  46          do k=1,isiz1
  47             do j=1,5
  48                do i=1,5
  49                   a(i,j,k,m) = 0.d0
  50                   b(i,j,k,m) = 0.d0
  51                   c(i,j,k,m) = 0.d0
  52                   d(i,j,k,m) = 0.d0
  53                enddo
  54             enddo
  55          enddo
  56       enddo
  57
  58 c---------------------------------------------------------------------
  59 c   compute the steady-state residuals
  60 c---------------------------------------------------------------------
  61       call rhs
  62
  63 c---------------------------------------------------------------------
  64 c   compute the L2 norms of newton iteration residuals
  65 c---------------------------------------------------------------------
  66       call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
  67      >             ist, iend, jst, jend,
  68      >             rsd, rsdnm )
  69
  70       call MPI_BARRIER( MPI_COMM_WORLD, IERROR )
  71
  72       call timer_clear(1)
  73       call timer_start(1)
  74
  75 c---------------------------------------------------------------------
  76 c   the timestep loop
  77 c---------------------------------------------------------------------
  78       do istep = 1, niter
  79
  80          if (id .eq. 0) then
  81             if (mod ( istep, 20) .eq. 0 .or.
  82      >            istep .eq. itmax .or.
  83      >            istep .eq. 1) then
  84                if (niter .gt. 1) write( *, 200) istep
  85  200           format(' Time step ', i4)
  86             endif
  87          endif
  88
  89 c---------------------------------------------------------------------
  90 c   perform SSOR iteration
  91 c---------------------------------------------------------------------
  92          do k = 2, nz - 1
  93             do j = jst, jend
  94                do i = ist, iend
  95                   do m = 1, 5
  96                      rsd(m,i,j,k) = dt * rsd(m,i,j,k)
  97                   end do
  98                end do
  99             end do
 100          end do
 101
 102          DO k = 2, nz -1
 103 c---------------------------------------------------------------------
 104 c   form the lower triangular part of the jacobian matrix
 105 c---------------------------------------------------------------------
 106             call jacld(k)
 107
 108 c---------------------------------------------------------------------
 109 c   perform the lower triangular solution
 110 c---------------------------------------------------------------------
 111             call blts( isiz1, isiz2, isiz3,
 112      >                 nx, ny, nz, k,
 113      >                 omega,
 114      >                 rsd,
 115      >                 a, b, c, d,
 116      >                 ist, iend, jst, jend,
 117      >                 nx0, ny0, ipt, jpt)
 118           END DO
 119
 120           DO k = nz - 1, 2, -1
 121 c---------------------------------------------------------------------
 122 c   form the strictly upper triangular part of the jacobian matrix
 123 c---------------------------------------------------------------------
 124             call jacu(k)
 125
 126 c---------------------------------------------------------------------
 127 c   perform the upper triangular solution
 128 c---------------------------------------------------------------------
 129             call buts( isiz1, isiz2, isiz3,
 130      >                 nx, ny, nz, k,
 131      >                 omega,
 132      >                 rsd, tv,
 133      >                 d, a, b, c,
 134      >                 ist, iend, jst, jend,
 135      >                 nx0, ny0, ipt, jpt)
 136           END DO
 137
 138 c---------------------------------------------------------------------
 139 c   update the variables
 140 c---------------------------------------------------------------------
 141
 142          do k = 2, nz-1
 143             do j = jst, jend
 144                do i = ist, iend
 145                   do m = 1, 5
 146                      u( m, i, j, k ) = u( m, i, j, k )
 147      >                    + tmp * rsd( m, i, j, k )
 148                   end do
 149                end do
 150             end do
 151          end do
 152
 153 c---------------------------------------------------------------------
 154 c   compute the max-norms of newton iteration corrections
 155 c---------------------------------------------------------------------
 156          if ( mod ( istep, inorm ) .eq. 0 ) then
 157             call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
 158      >                   ist, iend, jst, jend,
 159      >                   rsd, delunm )
 160 c            if ( ipr .eq. 1 .and. id .eq. 0 ) then
 161 c                write (*,1006) ( delunm(m), m = 1, 5 )
 162 c            else if ( ipr .eq. 2 .and. id .eq. 0 ) then
 163 c                write (*,'(i5,f15.6)') istep,delunm(5)
 164 c            end if
 165          end if
 166
 167 c---------------------------------------------------------------------
 168 c   compute the steady-state residuals
 169 c---------------------------------------------------------------------
 170          call rhs
 171
 172 c---------------------------------------------------------------------
 173 c   compute the max-norms of newton iteration residuals
 174 c---------------------------------------------------------------------
 175          if ( ( mod ( istep, inorm ) .eq. 0 ) .or.
 176      >        ( istep .eq. itmax ) ) then
 177             call l2norm( isiz1, isiz2, isiz3, nx0, ny0, nz0,
 178      >                   ist, iend, jst, jend,
 179      >                   rsd, rsdnm )
 180 c            if ( ipr .eq. 1.and.id.eq.0 ) then
 181 c                write (*,1007) ( rsdnm(m), m = 1, 5 )
 182 c            end if
 183          end if
 184
 185 c---------------------------------------------------------------------
 186 c   check the newton-iteration residuals against the tolerance levels
 187 c---------------------------------------------------------------------
 188          if ( ( rsdnm(1) .lt. tolrsd(1) ) .and.
 189      >        ( rsdnm(2) .lt. tolrsd(2) ) .and.
 190      >        ( rsdnm(3) .lt. tolrsd(3) ) .and.
 191      >        ( rsdnm(4) .lt. tolrsd(4) ) .and.
 192      >        ( rsdnm(5) .lt. tolrsd(5) ) ) then
 193 c            if (ipr .eq. 1 .and. id.eq.0) then
 194 c               write (*,1004) istep
 195 c            end if
 196             return
 197          end if
 198
 199       end do
 200
 201       call timer_stop(1)
 202       wtime = timer_read(1)
 203
 204
 205       call MPI_ALLREDUCE( wtime,
 206      >                    maxtime,
 207      >                    1,
 208      >                    MPI_DOUBLE_PRECISION,
 209      >                    MPI_MAX,
 210      >                    MPI_COMM_WORLD,
 211      >                    IERROR )
 212
 213
 214
 215       return
 216
 217  1001 format (1x/5x,'pseudo-time SSOR iteration no.=',i4/)
 218  1004 format (1x/1x,'convergence was achieved after ',i4,
 219      >   ' pseudo-time steps' )
 220  1006 format (1x/1x,'RMS-norm of SSOR-iteration correction ',
 221      > 'for first pde  = ',1pe12.5/,
 222      > 1x,'RMS-norm of SSOR-iteration correction ',
 223      > 'for second pde = ',1pe12.5/,
 224      > 1x,'RMS-norm of SSOR-iteration correction ',
 225      > 'for third pde  = ',1pe12.5/,
 226      > 1x,'RMS-norm of SSOR-iteration correction ',
 227      > 'for fourth pde = ',1pe12.5/,
 228      > 1x,'RMS-norm of SSOR-iteration correction ',
 229      > 'for fifth pde  = ',1pe12.5)
 230  1007 format (1x/1x,'RMS-norm of steady-state residual for ',
 231      > 'first pde  = ',1pe12.5/,
 232      > 1x,'RMS-norm of steady-state residual for ',
 233      > 'second pde = ',1pe12.5/,
 234      > 1x,'RMS-norm of steady-state residual for ',
 235      > 'third pde  = ',1pe12.5/,
 236      > 1x,'RMS-norm of steady-state residual for ',
 237      > 'fourth pde = ',1pe12.5/,
 238      > 1x,'RMS-norm of steady-state residual for ',
 239      > 'fifth pde  = ',1pe12.5)
 240
 241       end