*----------------------------------------------------------------------------*
* The PROMAT routine will first compute Zbar from the Z matrix and then *
* multiply Zbar by a vector V, the resultant vector is put in RV *
*----------------------------------------------------------------------------*
c
subroutine promat(n,ro,r,v,rv)
c
implicit real(8) (a-h,o-z)
real(8) v(n),rv(n),r(n*(n+1)/2),ro,a1
c
n1=n-1
n2=n-2
i1=1
j1=1
a1=r(1)*v(2)+ro*v(1)
do 2 j=1,n2
j1=j1+j
2 a1=a1+r(j1)*v(j+2)
rv(1)=a1
do 3 i=2,n2
i1=i1+i
j1=i1
a1=r(i1)*v(i+1)+ro*v(i)
do 4 j=i,n2
j1=j1+j
4 a1=a1+r(j1)*v(j+2)
j1=i1-i+1
i2=i-1
do 5 k=1,i2
j=i-k
j1=j1-1
5 a1=a1+r(j1)*v(j)
3 rv(i)=a1
i1=i1+n1
a1=r(i1)*v(n)+ro*v(n1)
j1=i1-n2
do 6 k=1,n2
j=n1-k
j1=j1-1
6 a1=a1+r(j1)*v(j)
rv(n1)=a1
i1=i1-n1
a1=ro*v(n)
do 7 j=1,n1
7 a1=a1+r(i1+j)*v(j)
rv(n)=a1
return
end
c
*----------------------------------------------------------------------------*
* PROJEC determine the projection of vector S onto the space defined by the *
* 3 orthogonal vectors V1,V2,V3 *
*----------------------------------------------------------------------------*
c
subroutine projec(n,v1,v2,v3,s)
c
implicit real(8) (a-h,o-z)
real(8) s(n),v1(n),v2(n),v3(n),a1,a2,a3
c
a1=provec(n,v1,s)/dble(n)
a2=provec(n,v2,s)
a3=provec(n,v3,s)
do 1 i=1,n
1 s(i)=s(i)-(a1+a2*v2(i)+a3*v3(i))
return
end
c
*----------------------------------------------------------------------------*
* PROVEC is a double precision routine to compute the dot product of 2 *
* double precision vectors X and Y *
*----------------------------------------------------------------------------*
c
real(8) function provec (n,x,y)
c
implicit real(8) (a-h,o-z)
real(8) x(n),y(n)
c
provec=0.d0
do 1 i=1,n
1 provec=provec+x(i)*y(i)
return
end
c
*----------------------------------------------------------------------------*
* DPROVEC is a double precision routine to compute the dot product of a *
* single precision X and double precision Y vectors *
*----------------------------------------------------------------------------*
real(8) function dprovec (n,x,y)
c
real(4) x(n)
real(8) y(n)
c
dprovec=0.d0
do 1 i=1,n
1 dprovec=dprovec+dble(x(i))*y(i)
return
end