fd2d.f90

是一个接收永磁体磁场计算的程序

!
    if ( time_steps * ( ( 10 * time_step ) / time_steps ) &
      == 10 * time_step ) then
      write ( *, * ) '  TIME_STEP = ', time_step
!      call timestamp ( )
    end if

    hhat(1:n) = u(1:n) / ( alpha + abs ( u(1:n) ) )
    f(1:n) = u(1:n) - u(1:n) * abs ( u(1:n) ) - v(1:n) * hhat(1:n)
    g(1:n) = beta * v(1:n) * hhat(1:n) - gamma * v(1:n)
    y1(1:n) = u(1:n) + delt * f(1:n)

    if ( .false. ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Right hand side Y1'
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, * ) i, y1(i)
      end do
    end if

    y2(1:n) = v(1:n) + delt * g(1:n)
!
!  Solve BU * u_new = y1
!        BV * v_new = y2
!
    if ( solve == 0 ) then

      call dgmres ( n, y1, u, nz_num, b_i, b_j, b_u, isym, matvec_triad, &
        msolve_identity, itol, tol, it_max, it, err, ierr, iunit, dummy, &
        dummy, rgwk, lrgw, igwk, ligw, dummy, idummy )

      if ( time_steps * ( ( 10 * time_step ) / time_steps ) &
        == 10 * time_step ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  Number of iterations:  ', it
        write ( *, '(a,g14.6)' ) '  Convergence measure is ', rgwk(1)
        write ( *, '(a,g14.6)' ) '  Error estimate ', err
        write ( *, '(a,i6)' ) '  Error code is ', ierr
      end if

    else

      call ds3_jac_sl ( n, nz_num, isym, b_i, b_j, b_u, y1, u, tol, it_max, &
        job, it, diff )

      if ( time_steps * ( ( 10 * time_step ) / time_steps ) &
        == 10 * time_step ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  Number of iterations:  ', it
        write ( *, '(a,g14.6)' ) '  Last change in solution ', diff
      end if

    end if

    if ( .false. ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution U'
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, * ) i, u(i)
      end do
    end if

    if ( solve == 0 ) then

      call dgmres ( n, y2, v, nz_num, b_i, b_j, b_v, isym, matvec_triad, &
        msolve_identity, itol, tol, it_max, it, err, ierr, iunit, dummy, &
        dummy, rgwk, lrgw, igwk, ligw, dummy, idummy )

      if ( .false. ) then
        write ( *, '(a)'       ) ' '
        write ( *, '(a,i6)'    ) '  Number of iterations:  ', it
        write ( *, '(a,g14.6)' ) '  Convergence measure is ', rgwk(1)
        write ( *, '(a,g14.6)' ) '  Error estimate ', err
        write ( *, '(a,i6)'    ) '  Error code is ', ierr
      end if

    else

      call ds3_jac_sl ( n, nz_num, isym, b_i, b_j, b_v, y2, v, tol, it_max, &
        job, it, diff )

      if ( time_steps * ( ( 10 * time_step ) / time_steps ) &
        == 10 * time_step ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  Number of iterations:  ', it
        write ( *, '(a,g14.6)' ) '  Last change in solution ', diff
      end if

    end if

    if ( .false. ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution V'
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, * ) i, v(i)
      end do
    end if

  end do
!
!  Make a PPM plot of the solution components U and V at the final time.
!
  allocate ( b_ppm(1:n) )
  allocate ( g_ppm(1:n) )
  allocate ( r_ppm(1:n) )

  u_max = maxval ( u(1:n) )
  r_ppm(1:n) = nint ( 255 *           u(1:n)   / u_max )
  g_ppm(1:n) = 0
  b_ppm(1:n) = nint ( 255 * ( u_max - u(1:n) ) / u_max )

  call ppma_write ( 'u2d.ppma', dimj, dimj, r_ppm, g_ppm, b_ppm, ierror )

  v_max = maxval ( v(1:n) )
  r_ppm(1:n) = nint ( 255 *           v(1:n)   / v_max )
  g_ppm(1:n) = 0
  b_ppm(1:n) = nint ( 255 * ( v_max - v(1:n) ) / v_max )

  call ppma_write ( 'v2d.ppma', dimj, dimj, r_ppm, g_ppm, b_ppm, ierror )
!
!  Free up the allocated memory.
!
  deallocate ( b_i )
  deallocate ( b_j )
  deallocate ( b_ppm )
  deallocate ( b_u )
  deallocate ( b_v )
  deallocate ( f )
  deallocate ( g )
  deallocate ( g_ppm )
  deallocate ( hhat )
  deallocate ( r_ppm )
  if ( solve == 0 ) then
    deallocate ( rgwk )
  end if
  deallocate ( u )
  deallocate ( u_grid )
  deallocate ( u0 )
  deallocate ( v )
  deallocate ( v_grid )
  deallocate ( v0 )
  deallocate ( x )
  deallocate ( x_grid )
  deallocate ( y )
  deallocate ( y_grid )
  deallocate ( y1 )
  deallocate ( y2 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FD2D:'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
!  call timestamp ( )

  stop
end
subroutine ds3_diagonal2 ( n, nz_num, isym, row, col, a, b )

!*****************************************************************************80
!
!! DS3_DIAGONAL2 reorders two square DS3 matrices so diagonal entries are first.
!
!  Discussion:
!
!    The DS3 storage format corresponds to the SLAP Triad format.
!
!    The DS3 storage format stores the row, column and value of each nonzero
!    entry of a sparse matrix.  The entries may be given in any order.  No
!    check is made for the erroneous case in which a given matrix entry is
!    specified more than once.
!
!    There is a symmetry option for square matrices.  If the symmetric storage
!    option is used, the format specifies that only nonzeroes on the diagonal
!    and lower triangle are stored.  However, this routine makes no attempt
!    to enforce this.  The only thing it does is to "reflect" any nonzero
!    offdiagonal value.  Moreover, no check is made for the erroneous case
!    in which both A(I,J) and A(J,I) are specified, but with different values.
!
!    This routine reorders the entries of A so that the first N entries
!    are exactly the diagonal entries of the matrix, in order.
!
!  Modified:
!
!    24 November 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer NZ_NUM, the number of nonzero elements in the matrix.
!
!    Input, integer ISYM, is 0 if the matrix is not symmetric, and 1
!    if the matrix is symmetric.  If the matrix is symmetric, then
!    only the nonzeroes on the diagonal and in the lower triangle are stored.
!
!    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
!    of the nonzero elements.
!
!    Input/output, real ( kind = 8 ) A(NZ_NUM), the nonzero elements
!    of the matrix.
!
!    Input/output, real ( kind = 8 ) B(NZ_NUM), the nonzero elements
!    of the second matrix.
!
  implicit none

  integer n
  integer nz_num

  real ( kind = 8 ) a(nz_num)
  real ( kind = 8 ) b(nz_num)
  integer col(nz_num)
  integer found
  integer i
  integer isym
  integer k
  integer row(nz_num)

  found = 0

  do k = 1, nz_num

    do while ( row(k) == col(k) )

      if ( row(k) == k ) then
        found = found + 1
        exit
      end if

      i = row(k)

      call i4_swap ( row(i), row(k) )
      call i4_swap ( col(i), col(k) )
      call r8_swap ( a(i), a(k) )
      call r8_swap ( b(i), b(k) )

      found = found + 1

      if ( n <= found ) then
        exit
      end if

    end do

    if ( n <= found ) then
      exit
    end if

  end do

  if ( found < n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DS3_DIAGONAL - Warning!'
    write ( *, '(a,i6)' ) '  Number of diagonal entries expected was ', n
    write ( *, '(a,i6)' ) '  Number found was ', found
    stop
  end if

  return
end
subroutine ds3_jac_sl ( n, nz_num, isym, row, col, a, b, x, tol, it_max, &
  job, it, diff )

!*****************************************************************************80
!
!! DS3_JAC_SL solves a DS3 system using Jacobi iteration.
!
!  Discussion:
!
!    The DS3 storage format corresponds to the SLAP Triad format.
!
!    The DS3 storage format stores the row, column and value of each nonzero
!    entry of a sparse matrix.  The entries may be given in any order.  No
!    check is made for the erroneous case in which a given matrix entry is
!    specified more than once.
!
!    There is a symmetry option for square matrices.  If the symmetric storage
!    option is used, the format specifies that only nonzeroes on the diagonal
!    and lower triangle are stored.  However, this routine makes no attempt
!    to enforce this.  The only thing it does is to "reflect" any nonzero
!    offdiagonal value.  Moreover, no check is made for the erroneous case
!    in which both A(I,J) and A(J,I) are specified, but with different values.
!
!    This routine REQUIRES that the matrix be square, that the matrix
!    have nonzero diagonal entries, and that the first N entries of
!    the array A be exactly the diagonal entries of the matrix, in order.
!
!  Modified:
!
!    27 November 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!    Input, integer NZ_NUM, the number of nonzero elements in the matrix.
!
!    Input, integer ISYM, is 0 if the matrix is not symmetric, and 1
!    if the matrix is symmetric.  If the matrix is symmetric, then
!    only the nonzeroes on the diagonal and in the lower triangle are stored.
!
!    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
!    of the nonzero elements.
!
!    Input, real ( kind = 8 ) A(NZ_NUM), the nonzero elements of the matrix.
!
!    Input, real ( kind = 8 ) B(N), the right hand side of the linear system.
!
!    Input/output, real ( kind = 8 ) X(N), an approximate solution
!    to the system.
!
!    Input, real ( kind = 8 ) TOL, a tolerance.  If the maximum change in
!    the solution is less than TOL, the iteration is terminated early.
!
!    Input, integer IT_MAX, the maximum number of iterations.
!
!    Input, integer JOB, specifies the system to solve.
!    0, solve A * x = b.
!    nonzero, solve A' * x = b.
!
!    Output, integer IT, the number of iterations taken.
!
!    Output, real ( kind = 8 ) DIFF, the maximum change in the solution
!    on the last iteration.
!
  implicit none

  integer n
  integer nz_num

  real ( kind = 8 ) a(nz_num)
  real ( kind = 8 ) b(n)
  integer col(nz_num)
  real ( kind = 8 ) diff
  integer i
  integer isym
  integer it
  integer it_max
  integer it_num
  integer j
  integer job
  integer k
  integer row(nz_num)
  real ( kind = 8 ) tol
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xnew(n)
  real ( kind = 8 ) x_norm

  if ( job == 0 ) then

    do it_num = 1, it_max

      it = it_num
!
!  Initialize to right hand side.
!
      xnew(1:n) = b(1:n)
!
!  Subtract off-diagonal terms.
!
      do k = n+1, nz_num
        i = row(k)
        j = col(k)
        xnew(i) = xnew(i) - a(k) * x(j)
        if ( isym == 1 ) then
          xnew(j) = xnew(j) - a(k) * x(i)
        end if
      end do
!
!  Divide by diagonal terms.
!
      xnew(1:n) = xnew(1:n) / a(1:n)
!
!  Measure change:
!
      x_norm = maxval ( abs ( x(1:n) ) )
      diff = maxval ( abs ( xnew(1:n) - x(1:n) ) )
!
!  Update.
!
      x(1:n) = xnew(1:n)
!
!  Test for early termination.
!
      if ( diff <= tol  * ( x_norm + 1.0D+00 ) ) then
        exit
      end if

    end do

  else

    do it_num = 1, it_max

      it = it_num
!
!  Initialize to right hand side.
!
      xnew(1:n) = b(1:n)
!
!  Subtract off-diagonal terms.
!
      do k = n+1, nz_num
        i = col(k)
        j = row(k)
        xnew(i) = xnew(i) - a(k) * x(j)
        if ( isym == 1 ) then
          xnew(j) = xnew(j) - a(k) * x(i)
        end if
      end do
!
!  Divide by diagonal terms.
!
      xnew(1:n) = xnew(1:n) / a(1:n)
!
!  Measure change:
!
      x_norm = maxval ( abs ( x(1:n) ) )
      diff = maxval ( abs ( xnew(1:n) - x(1:n) ) )
!
!  Update.
!
      x(1:n) = xnew(1:n)
!
!  Test for early termination.
!
      if ( diff <= tol * ( x_norm + 1.0D+00 ) ) then
        exit
      end if

    end do

  end if

  return
end
subroutine i4_swap ( i, j )

!*****************************************************************************80
!
!! I4_SWAP swaps two I4's.
!
!  Modified:
!
!    30 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On output, the values of I and
!    J have been interchanged.
!
  implicit none

  integer i
  integer j
  integer k

  k = i
  i = j
  j = k

  return
end
subroutine i4vec2_sort_a_plus2 ( n, a1, a2, a3, a4 )

!*****************************************************************************80
!
!! I4VEC2_SORT_A_PLUS2 ascending sorts integer pairs, and adjusts real vectors.
!
!  Discussion:
!
!    The data items have the form (I,J,K,L), where I and J are integers
!    and K and L are real values.  The data is to be ascending sorted on the
!    I and J fields.
!
!  Modified:
!
!    05 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of items of data.
!
!    Input/output, integer A1(N), integer A2(N), the data to be sorted.
!
!    Input/output, real ( kind = 8 ) A3(N), A4(N), the data to be rearranged
!    in accordance with the sorting of A1 and A2.
!
  implicit none

  integer n

  integer a1(n)
  integer a2(n)
  real ( kind = 8 ) a3(n)