DataMuseum.dk

top - download

    Length: 3336 (0xd08)
    Types: TextFile
    Names: »vecdemo.p«

└─⟦a0efdde77⟧ Bits:30001252 EUUGD11 Tape, 1987 Spring Conference Helsinki
    └─ ⟦this⟧ »EUUGD11/euug-87hel/sec1/prep/vecdemo.p« 

TextFile

c***********************************************************************
c                                                                      *
c                    subroutine W_ACCEL_LIN                            *
c                                                                      *
c Do the Wachspress accelleration.                                     *
c   The solution is expressed as a linear combination of the previous  *
c iterate and the lowest order fourier modes and the coefficients      *
c are found so as to minimize the error.                               *
c                                                                      *
c P.R.OVE  7/6/85                                                      *
c***********************************************************************

        subroutine w_accel_l(psi, lin_fac, source, omega)
        use ellipdim
        do_limits = { mx, my }

        if (w_bypass) return
        w_error = FALSE

c**********************************************************************
c Set up the basis consisting of past iterates                        *
c**********************************************************************
[	basis(#,#,1) = psi(#,#)
	basis(#,#,2) = psi(#,#) - psi_alt(#,#,1)
	basis(#,#,3) = psi(#,#) - 2*psi_alt(#,#,1) + psi_alt(#,#,2)
	basis(#,#,4) = 1      ]
	call periodic( mx, my, basis1 )
	call periodic( mx, my, basis2 )
	call periodic( mx, my, basis3 )
	call periodic( mx, my, basis4 )

c**********************************************************************
c Calculate the Wachspress matrix and the source vector               *
c**********************************************************************
        do i = 1, w_dim
	ii = i
	do j = i, w_dim
	   jj = j
           call make_mat_l(psi, lin_fac, source, omega, i, j)
        end_do
	end_do

	do i = 1, w_dim
           w_source(i) = 0
           w_source(i) = source(#,#)*basis(#,#,i) + w_source(i)
        end_do

c**********************************************************************
c invert the symmetric matrix and improve the solution psi.           *
c**********************************************************************
        call linsys(w_matrix, w_dim, w_dim, w_source, w_coeff,
     *              ising, lfirst, lprint, work)
        if (ising == 1) then
c          write(*,*) ' WARNING:  W_matrix is singular '
           w_error = TRUE
           goto 99
        endif

c calculate the improved solution
        psi(#,#) = 0
        do i = 1, w_dim
           psi(#,#) = psi(#,#) + w_coeff(i)*basis(#,#,i)
        end_do

c**********************************************************************
c output section for error checking  (optional)                       *
c**********************************************************************
	go to 99
        do i = 1, w_dim
           write(*,100) i, .5*w_matrix(i,i) - w_source(i),
     *                  i, w_coeff(i)
100        format(' action(',i1')= ',g16.9,'    w_coeff(',i1,')= ',
     *               g16.9)
        end_do

	do_limits = { w_dim }
        action = 0
        do i = 1, w_dim
           action = action + w_matrix(i,#)*w_coeff(i)*w_coeff(#)
        end_do
        action = action/2
        action = action - w_source(#)*w_coeff(#)
        write(*,*) ' new action = ',action

99      return
        end