fd2d.f

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

	 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