home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / prog_oth / prep_065.lzh / PREP / VECDEM.P < prev   
Text File  |  1991-08-16  |  2KB  |  72 lines
  1. c Demo to demonstrate some PREP facilities.  This program is a demo
  2. c only and will not compile without a lot of variable definitions.
  3.  
  4. #include "vecdem.h"
  5.  
  6.         subroutine w_accel_l(psi, lin_fac, source, omega)
  7.         include "ellipdim"
  8.  
  9.         if (w_bypass) return
  10.         w_error = FALSE
  11.  
  12. c Set up the basis consisting of past iterates
  13. [    basis(#,#,1) = psi(#,#)
  14.     basis(#,#,2) = psi(#,#) - psi_alt(#,#,1)
  15.     basis(#,#,3) = psi(#,#) - 2*psi_alt(#,#,1) + psi_alt(#,#,2)
  16.     basis(#,#,4) = 1      ]
  17.     PERIODIC( basis1 )
  18.     PERIODIC( basis2 )
  19.     PERIODIC( basis3 )
  20.     PERIODIC( basis4 )
  21.  
  22. c Calculate the matrix and the source vector
  23.         do i = 1, w_dim
  24.     ii = i
  25.     do j = i, w_dim
  26.     jj = j
  27.            call make_mat_l(psi, lin_fac, source, omega, i, j)
  28.         end_do
  29.     end_do
  30.  
  31.     do i = 1, w_dim
  32.            w_source(i) = 0
  33.            w_source(i) = source(#,#)*basis(#,#,i) + w_source(i)
  34.         end_do
  35.  
  36. c invert the symmetric matrix
  37.         call linsys(w_matrix, w_dim, w_dim, w_source, w_coeff, ising, lfirst,
  38.      *              lprint, work)
  39.         if (ising == 1) then
  40.            write(*,*) ' WARNING:  W_matrix is singular '
  41.            w_error = TRUE
  42.            return
  43.         endif
  44.  
  45. c calculate the improved solution
  46.         psi(#,#) = 0
  47.         do i = 1, w_dim
  48.            psi(#,#) = psi(#,#) + w_coeff(i)*basis(#,#,i)
  49.         end_do
  50.  
  51. c output section for error checking
  52.         do i = 1, w_dim
  53.            write(*,100) i, .5*w_matrix(i,i) - w_source(i),
  54.      *                  i, w_coeff(i)
  55.         end_do
  56.  
  57.     do_limits = { w_dim }
  58.         action = 0
  59.         do i = 1, w_dim
  60.            action = action + w_matrix(i,#)*w_coeff(i)*w_coeff(#)
  61.         end_do
  62.         action = action/2
  63.         action = action - w_source(#)*w_coeff(#)
  64.         write(*,*) ' new action = ',action
  65.  
  66.         return
  67.  
  68.  
  69. 100     format(' action(',i1')= ',g16.9,'    w_coeff(',i1,')= ', g16.9)
  70.  
  71.         end
  72.