APPENDIX

Below are reprinted subroutines from the current version of PVI. Since the paper version of PVI was published, I modified ident() and diag() so that their outputs can be overlaid on their inputs.

subroutine null(     xx, n)
integer i, n;	real xx( n)
do i= 1, n
	xx(i) = 0.
return; end

subroutine conjnull( conj, add, x, nx,  y, ny )
integer ix, iy,      conj, add,    nx,     ny
real  				x( nx), y( ny )
if( add == 0 )
	if( conj == 0 )
		do iy= 1, ny 
			y(iy) = 0.
	else
		do ix= 1, nx 
			x(ix) = 0. 
return; end

subroutine ident( conj, add, epsilon, n, pp,    qq   )
integer i,        conj, add,          n
real 			     epsilon,    pp(n), qq(n)  # equivalence (pp,qq) OK
if( conj == 0 ) {
	if( add == 0 ) { do i=1,n {  qq(i) =         epsilon * pp(i) } }
	else           { do i=1,n {  qq(i) = qq(i) + epsilon * pp(i) } }
	}
else {  if( add == 0 ) { do i=1,n {  pp(i) =         epsilon * qq(i) } }
	else           { do i=1,n {  pp(i) = pp(i) + epsilon * qq(i) } }
	}
return; end

subroutine diag( conj, add, lambda,n,  pp,    qq)
integer i,       conj, add,        n 		      # equivalence (pp,qq) OK
real 			    lambda(n), pp(n), qq(n)
if( conj == 0 ) {
	if( add == 0 ) { do i=1,n {  qq(i) =         lambda(i) * pp(i) } }
	else           { do i=1,n {  qq(i) = qq(i) + lambda(i) * pp(i) } }
	}
else {  if( add == 0 ) { do i=1,n {  pp(i) =         lambda(i) * qq(i) } }
	else           { do i=1,n {  pp(i) = pp(i) + lambda(i) * qq(i) } }
	}
return; end

# CINLOI -- Convolution INternal with Lags.  Output is conjugate with INPUT.
#
subroutine cinloi( conj, add, lag1,lag2,nb1,nb2,bb,  n1,n2, xx, yy)
integer            conj, add, lag1,lag2,nb1,nb2,     n1,n2        # lag=1 causal
real	                                bb(nb1,nb2), xx(n1,n2), yy(n1,n2)
integer y1,y2, x1,x2, b1,b2
call conjnull(     conj, add,                        xx,n1*n2,  yy,n1*n2 )
if( conj == 0 )
	do b2=1,nb2 { do y2= 1+nb2-lag2, n2-lag2+1 {  x2= y2 - b2 + lag2      
	do b1=1,nb1 { do y1= 1+nb1-lag1, n1-lag1+1 {  x1= y1 - b1 + lag1      
				yy(y1,y2) = yy(y1,y2) + bb(b1,b2) * xx(x1,x2)
				}} }}
else
	do b2=1,nb2 { do y2= 1+nb2-lag2, n2-lag2+1 {  x2= y2 - b2 + lag2      
	do b1=1,nb1 { do y1= 1+nb1-lag1, n1-lag1+1 {  x1= y1 - b1 + lag1      
				xx(x1,x2) = xx(x1,x2) + bb(b1,b2) * yy(y1,y2)
				}} }}
return; end