c#> matvec_product.dc2
c Subroutine:	MATVEC_PRODUCT
c
c Purpose:      3-dimensional matrix - vector product                   
c
c Category:	MATH
c
c File:		irco.shl
c
c Author:	Do Kester
c
c Use:	call MATVEC_PRODUCT( 
c		RMAT, 	I	double precision array(3,3)
c		X, 	I	double precision array(3,N)
c		Y, 	O	double precision array(3,N)
c		N )	I	integer
c       N     Number of X and Y vectors                     
c       X     an array with N 3-dimensional (column) vectors
c       RMAT  a 3 x 3 matrix                                
c       Y     an array with N (column) vectors, Y = RMAT times X
c                                                                       
c Description:
c       Mainly meant as a shorthand for                                 
c         Y(i,j) = RMAT(i,1) * X(1,j)                                   
c                + RMAT(i,2) * X(2,j)                                   
c                + RMAT(i,3) * X(3,j)                                   
c                                                                       
c       Use it with N=1 for simple 3-dimensional vectors,               
c              with N=3 for 3 x 3 matrix multiplication,                
c              with arbitrary N for e.g.                                
c                   large volume coordinate transformations.            
c                                                                       
c       It is allowed that X and Y are the same array.                  
c                                                                       
c       Peculiar things will happen if you use it for                   
c       matrix multiplication, and RMAT and Y are the same array.       
c                                                                      
c Update:       880112 Timo Prusti & Peter Arendz : standarized         
c               870518 Albrecht de Jonge documentation improved         
c               850625 Albrecht de Jonge                                
c#<
                                                                        
        subroutine matvec_product( rmat, x, y, n )

        integer           n
        double precision  x(3,n), y(3,n), rmat( 3, 3 )

        double precision  t1, t2, t3
        integer           j

        for j = 1, n
          t1 = x(1,j)
          t2 = x(2,j)
          t3 = x(3,j)

          y(1,j) = rmat(1,1) * t1 + rmat(1,2) * t2 + rmat(1,3) * t3
          y(2,j) = rmat(2,1) * t1 + rmat(2,2) * t2 + rmat(2,3) * t3
          y(3,j) = rmat(3,1) * t1 + rmat(3,2) * t2 + rmat(3,3) * t3

        cfor

        return
        end