fd2d.f90

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

  real ( kind = 8 ) a4(n)
  integer i
  integer indx
  integer isgn
  integer j
  integer temp
  real ( kind = 8 ) rtemp

  if ( n <= 1 ) then
    return
  end if
!
!  Initialize.
!
  i = 0
  indx = 0
  isgn = 0
  j = 0
!
!  Call the external heap sorter.
!
  do

    call sort_heap_external ( n, indx, i, j, isgn )
!
!  Interchange the I and J objects and "carry" K.
!
    if ( 0 < indx ) then

      temp = a1(i)
      a1(i) = a1(j)
      a1(j) = temp

      temp = a2(i)
      a2(i) = a2(j)
      a2(j) = temp

      rtemp = a3(i)
      a3(i) = a3(j)
      a3(j) = rtemp

      rtemp = a4(i)
      a4(i) = a4(j)
      a4(j) = rtemp
!
!  Compare the I and J objects.
!
    else if ( indx < 0 ) then

      if ( a1(i) < a1(j) ) then
        isgn = -1
      else if ( a1(i) == a1(j) .and. a2(i) <= a2(j) ) then
        isgn = -1
      else
        isgn = +1
      end if

    else if ( indx == 0 ) then

      exit

    end if

  end do

  return
end
subroutine matvec_triad ( n, x, y, nelt, ia, ja, a, isym )

!***********************************************************************
!
!! MATVEC_TRIAD computes A*X for a matrix A stored in SLAP Triad form.
!
!  Modified:
!
!    21 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of elements in the vectors.
!
!    Input, real ( kind = 8 ) X(N), the vector to be multiplied by A.
!
!    Output, real ( kind = 8 ) Y(N), the product A * X.
!
!    Input, integer NELT, the number of nonzero entries in A.
!
!    Input, integer IA(NELT), JA(NELT), real ( kind = 8 ) A(NELT), the data
!    structure storing the sparse matrix.
!
!    Input, integer ISYM, is 0 if all nonzero entries of the matrix
!    are stored, and 1 if only the diagonal and upper or lower triangle
!    are stored.
!
  implicit none

  integer n
  integer nelt

  real ( kind = 8 ) a(nelt)
  integer ia(nelt)
  integer isym
  integer ja(nelt)
  integer k
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) y(n)

  y(1:n) = 0.0D+00

  do k = 1, nelt

    if ( ia(k) < 1 .or. n < ia(k) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'MATVEC_TRIAD - Fatal error!'
      write ( *, '(a,i6,a,i6)' ) '  IA(',k,') = ', ia(k)
      stop
    end if

    if ( ja(k) < 1 .or. n < ja(k) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'MATVEC_TRIAD - Fatal error!'
      write ( *, '(a,i6,a,i6)' ) '  JA(',k,') = ', ja(k)
      stop
    end if

    y(ia(k)) = y(ia(k)) + a(k) * x(ja(k))

  end do

  return
end
subroutine msolve_identity ( n, r, z, nelt, ia, ja, a, isym, rwork, iwork )

!***********************************************************************
!
!! MSOLVE_IDENTITY applies the identity matrix preconditioner.
!
!  Discussion:
!
!    Most SLAP solver routines require a preconditioner routine
!    that can solve M * Z = R.  If no preconditioning is required,
!    then you can simply behave as though the preconditioning matrix
!    M was the identity matrix.
!
!  Modified:
!
!    21 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of elements in the vectors.
!
!    Input, real ( kind = 8 ) R(N), the right hand side.
!
!    Output, real ( kind = 8 ) Z(N), the solution of M * Z = R.
!
!    Input, integer NELT, the number of nonzero entries in A.
!
!    Input, integer IA(NELT), JA(NELT), real ( kind = 8 ) A(NELT), the data
!    structure storing the sparse matrix.
!
!    Input, integer ISYM, is 0 if all nonzero entries of the matrix
!    are stored, and 1 if only the diagonal and upper or lower triangle
!    are stored.
!
!    Input, real ( kind = 8 ) RWORK(*), integer IWORK(*), arrays that can be
!    used to pass information to the preconditioner.
!
  implicit none

  integer n
  integer nelt

  real ( kind = 8 ) a(nelt)
  integer ia(nelt)
  integer isym
  integer iwork(*)
  integer ja(nelt)
  real ( kind = 8 ) r(n)
  real ( kind = 8 ) rwork(*)
  real ( kind = 8 ) z(n)

  z(1:n) = r(1:n)

  return
end
subroutine r8_swap ( x, y )

!*****************************************************************************80
!
!! R8_SWAP swaps two R8's.
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none

  real ( kind = 8 ) x
  real ( kind = 8 ) y
  real ( kind = 8 ) z

  z = x
  x = y
  y = z

  return
end
subroutine sparse ( i, j, v, m, n, mode, nz_num, nz, b_i, b_j, b_v )

!*****************************************************************************80
!
!! SPARSE manages the storage of sparse matrix information.
!
!  Discussion:
!
!    To use this routine, first call it with MODE = 0 to initialize.
!
!    Then, for every nonzero entry of the matrix, call with MODE = 1,
!    just to count the entries.  The value of NZ_NUM at the end of
!    this step can be used to allocate storage for B_I, B_J, and B_V.
!
!    Then, for every nonzero entry, call again with MODE = 2, and the
!    I, J, and V values set.  These will be transferred into the
!    B_I, B_J and B_V arrays.
!
!    At the moment, the data still needs to be sorted and merged, but
!    that must be done externally to this routine.
!
!  Modified:
!
!    30 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, the row and column index of the entry.
!
!    Input, real ( kind = 8 ) V, the value.
!
!    Input, integer M, N, the number of rows and columns in the matrix.
!
!    Input, integer MODE
!    0, initialize NZ and NZ_NUM.
!    1, "report" each entry, return a count in NZ_NUM.
!    2, store each entry, and keep a running count in NZ.
!    3, sort all the entries (not set up yet).
!    4, compress a sorted set of entries (not set up yet).
!
!    Input/output, integer NZ_NUM, is used in MODE = 0 and 1 only.  
!    Mode 0 initializes NZ_NUM to 0.  Then, every call with MODE = 1
!    increments NZ_NUM by 1.  This results in a crude count of the maximum
!    number of storage entries needed.
!
!    Input/output, integer NZ, is used only in MODE = 0 and 2.
!    A call with MODE = 0 initializes NZ to 0.  Then, each call with
!    MODE = 2 presumably is providing another set of (I,J,V) data
!    to be stored in the sparse matrix array.  NZ is incremented
!    by 1, and the data is stored at that index.
!
!    Input/output, integer B_I(*), B_J(*), real ( kind = 8 ) B_V(*).  These
!    vectors are used to store the sparse matrix data.  Presumably, a dimension
!    of NZ_NUM (determined with MODE = 1) would be enough to store 
!    the entries (assigned with MODE = 2).  Initially, the data is
!    stored as the user reports it, which means that the data is unsorted,
!    and that the same row and column indices may appear multiple times.
!
  implicit none

  integer b_i(*)
  integer b_j(*)
  real ( kind = 8 ) b_v(*)
  integer i
  integer j
  integer m
  integer mode
  integer n
  integer nz
  integer nz_num
  real ( kind = 8 ) v
!
!  Initialize.
!
  if ( mode == 0 ) then

    nz_num = 0
    nz = 0
!
!  Just count entries one by one.
!
  else if ( mode == 1 ) then

    nz_num = nz_num + 1
!
!  Store.
!
  else if ( mode == 2 ) then

    nz = nz + 1
    b_i(nz) = i
    b_j(nz) = j
    b_v(nz) = v

  else if ( mode == 3 ) then

  else if ( mode == 4 ) then

  end if

  return
end
subroutine sort_heap_external ( n, indx, i, j, isgn )

!*****************************************************************************80
!
!! SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.
!
!  Discussion:
!
!    The actual list of data is not passed to the routine.  Hence this
!    routine may be used to sort integers, reals, numbers, names,
!    dates, shoe sizes, and so on.  After each call, the routine asks
!    the user to compare or interchange two items, until a special
!    return value signals that the sorting is completed.
!
!  Modified:
!
!    05 February 2004
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Parameters:
!
!    Input, integer N, the number of items to be sorted.
!
!    Input/output, integer INDX, the main communication signal.
!
!    The user must set INDX to 0 before the first call.
!    Thereafter, the user should not change the value of INDX until
!    the sorting is done.
!
!    On return, if INDX is
!
!      greater than 0,
!      * interchange items I and J;
!      * call again.
!
!      less than 0,
!      * compare items I and J;
!      * set ISGN = -1 if I < J, ISGN = +1 if J < I;
!      * call again.
!
!      equal to 0, the sorting is done.
!
!    Output, integer I, J, the indices of two items.
!    On return with INDX positive, elements I and J should be interchanged.
!    On return with INDX negative, elements I and J should be compared, and
!    the result reported in ISGN on the next call.
!
!    Input, integer ISGN, results of comparison of elements I and J.
!    (Used only when the previous call returned INDX less than 0).
!    ISGN <= 0 means I is less than or equal to J;
!    0 <= ISGN means I is greater than or equal to J.
!
  implicit none

  integer i
  integer, save :: i_save = 0
  integer indx
  integer isgn
  integer j
  integer, save :: j_save = 0
  integer, save :: k = 0
  integer, save :: k1 = 0
  integer n
  integer, save :: n1 = 0
!
!  INDX = 0: This is the first call.
!
  if ( indx == 0 ) then

    i_save = 0
    j_save = 0
    k = n / 2
    k1 = k
    n1 = n
!
!  INDX < 0: The user is returning the results of a comparison.
!
  else if ( indx < 0 ) then

    if ( indx == -2 ) then

      if ( isgn < 0 ) then
        i_save = i_save + 1
      end if

      j_save = k1
      k1 = i_save
      indx = -1
      i = i_save
      j = j_save
      return

    end if

    if ( 0 < isgn ) then
      indx = 2
      i = i_save
      j = j_save
      return
    end if

    if ( k <= 1 ) then

      if ( n1 == 1 ) then
        i_save = 0
        j_save = 0
        indx = 0
      else
        i_save = n1
        n1 = n1 - 1
        j_save = 1
        indx = 1
      end if

      i = i_save
      j = j_save
      return

    end if

    k = k - 1
    k1 = k
!
!  0 < INDX, the user was asked to make an interchange.
!
  else if ( indx == 1 ) then

    k1 = k

  end if

  do

    i_save = 2 * k1

    if ( i_save == n1 ) then
      j_save = k1
      k1 = i_save
      indx = -1
      i = i_save
      j = j_save
      return
    else if ( i_save <= n1 ) then
      j_save = i_save + 1
      indx = -2
      i = i_save
      j = j_save
      return
    end if

    if ( k <= 1 ) then
      exit
    end if

    k = k - 1
    k1 = k

  end do

  if ( n1 == 1 ) then
    i_save = 0
    j_save = 0
    indx = 0
    i = i_save
    j = j_save
  else
    i_save = n1
    n1 = n1 - 1
    j_save = 1
    indx = 1
    i = i_save
    j = j_save
  end if

  return
end