fd2d.f

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

!
!  About 10 times during the iteration, print out the current time step.
!
	if ( time_steps * ( ( 10 * time_step ) / time_steps ) 
     &		== 10 * time_step ) then
		  write ( *, * ) '  TIME_STEP = ', time_step
c		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)' ) ' '
c	 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

c	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)