      subroutine ax2y(n,la,a,x,y)
      implicit double precision (a-h,o-z)
      dimension a(la,n),x(n),y(n)

      character*120 rcsID
      data rcsID
     &/'$Id: matrix2.f,v 1.5 2005/10/08 01:58:55 xpw Exp $'/

      do 100 i=1,n
      y(i)=0.d0
      do 100 j=1,n
  100 y(i)=y(i)+a(i,j)*x(j)
      return
      end
c------------------------------------------
      subroutine ax2y2(n,m,la,a,x,y)
      implicit double precision (a-h,o-z)
      dimension a(la,m),x(m),y(n)
      do 100 i=1,n
      y(i)=0.d0
      do 100 j=1,m
  100 y(i)=y(i)+a(i,j)*x(j)
      return
      end
c------------------------------------------
c
      subroutine atx2y(m,n,la,a,x,y)
      implicit double precision (a-h,o-z)
      dimension a(la,m),x(n),y(m)
      do 100 i=1,m
      y(i) = 0.d0
      do 100 j=1,n
 100  y(i) = y(i) + a(j,i)*x(j)
      return
      end
c------------------------------------------
c
	subroutine ab2cn(l,m,n,la,lb,lc,a,b,c)
	implicit double precision (a-h,o-z)
	dimension a(la,m),b(lb,n),c(lc,n)
	do 10 i=1,l
	do 10 j=1,n
	c(i,j)=0.d0
	do 10 k=1,m
 10	c(i,j)=c(i,j)+a(i,k)*b(k,j)
	return
	end
c-------------------------------------------
c
	subroutine atb2cn(l,m,n,la,lb,lc,a,b,c)
	implicit double precision (a-h,o-z)
	dimension a(la,l),b(lb,n),c(lc,n)
	do 10 i=1,l
	do 10 j=1,n
	c(i,j)=0.d0
	do 10 k=1,m
 10	c(i,j)=c(i,j)+a(k,i)*b(k,j)
	return
	end
c-------------------------------------------
c
	subroutine abt2cn(l,m,n,la,lb,lc,a,b,c)
	implicit double precision (a-h,o-z)
	dimension a(la,m),b(lb,m),c(lc,n)
	do 10 i=1,l
	do 10 j=1,n
	c(i,j)=0.d0
	do 10 k=1,m
 10	c(i,j)=c(i,j)+a(i,k)*b(j,k)
	return
	end
c-------------------------------------------
c
	subroutine rtr2u(n,r,u)
c...purpose:
c     to compute R^t R
c     the partial(!) upper triangular part of the results are vector 
c     stored in u
c     n is the matrix dimension
c     R is vector stored upper triangular matrix
	implicit double precision (a-h,o-z)
	dimension r(n*(n+1)/2), u(n*(n+1)/2)
	l=0
	do j=1,n
	do i=1,j
	l=l+1
	u(l)=0.d0
	do k=1,i
	ki=i*(i-1)/2+k
	kj=j*(j-1)/2+k
	u(l)=u(l)+r(ki)*r(kj)
	end do  ! k
	end do   ! i
	end do    ! j
	end
c-------------------------------------------
c
	subroutine rtr2c(n,r,lc,c)
c...purpose:
c     to compute R^t R
c     the whole results are stored in matrix c
c     n is the matrix dimension
c     R is vector stored upper triangular matrix
c     lc is the leading dimension of c in the calling routine
	implicit double precision (a-h,o-z)
	dimension r(n*(n+1)/2), c(lc,n)
	do j=1,n
	do i=1,j
	c(i,j)=0.d0
	do k=1,i
	ki=i*(i-1)/2+k
	kj=j*(j-1)/2+k
	c(i,j)=c(i,j)+r(ki)*r(kj)
	end do  ! k
	end do   ! i
	end do    ! j

	do i=1,n
	do j=1,i
	c(i,j)=c(j,i)
	end do
	end do
	end
c-------------------------------------------
	subroutine RAC(R,N,A,IMAXA,IRA,JCA,W,IMAXW,IRW)
C
C         TO COMPUTE  IRW ROWS OF THE PRODUCT W=R*A
C
C     WHERE R IS UPPER TRIANGULAR VECTOR SUBSCRIPTED AND OF
C     DIMENSION N, A HAS JCA COLUMNS AND IRA ROWS. 
C     WHEN IRA .LT. N THE BOTTOM N-IRA ROWS OF A ARE ASSUMED TO BE ZERO
C
C *** INPUTS:
C
C     R(N*(N+1)/2)  UPPER TRIANGULAR VECTOR STORED INPUT MATRIX
C
C     N             DIMENSION OF R.
C
C     A(IRA,JCA)    INPUT RECTANGULAR MATRIX
C
C     IMAXA         ROW DIMENSION OF A
C
C     IRA           NUMBER OF ROWS IN THE A MATRIX
C                   (IRA.LE.IMAXA AND IRA .LE. N)
C
C     JCA           NUMBER OF COLUMNS IN THE A MATRIX
C
C     IMAXW         ROW DIMENSION OF W
C
C     IRW           NUMBER OF ROWS OF THE PRODUCT R*A THAT ARE
C                   COMPUTED (IRW.LE.IMAXW AND IRW .LE. N)
C
C *** OUTPUT:
C
C     W(IRW,N)      RESULTING RECTANGULAR MATRIX, = R*A
C                   IF THE COMPILER ALLOWS IT W CAN OVERWRITE A

	implicit double precision (a-h,o-z)
      DIMENSION  R(*),A(IMAXA,JCA),W(IMAXW,JCA)
      DOUBLE PRECISION SUM
C
      DO 40 J=1,JCA
        II=0
        DO 30 I=1,IRW
          II=II+I
          SUM=0.D0
          IF (I.GT.IRA) GO TO 15
              IK=II
              DO 10 K=I,IRA
                 SUM=SUM+R(IK)*A(K,J)
                 IK=IK+K
  10          CONTINUE
  15      CONTINUE
          W(I,J)=SUM
  30    CONTINUE
  40  CONTINUE
      RETURN
      END

