      SUBROUTINE STRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
     $                   B, LDB )
****************************************************************************
*                                                                          *
*   DATA PARALLEL BLAS based on MPL                                        *
*                                                                          *
*   Version 1.0   1/9-92 ,                                                 *
*   For MasPar MP-1 computers                                              *
*                                                                          *
*   para//ab, University of Bergen, NORWAY                                 *
*                                                                          *
*   These programs must be called using F90 style array syntax.            *
*   Note that the F77 style calling sequence has been retained             *
*   in this version for compatibility reasons, be aware that               *
*   parameters related to the array dimensions and shape therefore may     *
*   be redundant and without any influence.                                *
*   The calling sequence may be changed in a future version.               *
*   Please report any BUGs, ideas for improvement or other                 *
*   comments to                                                            *
*                    adm@parallab.uib.no                                   *
*                                                                          *
*   Future versions may then reflect your suggestions.                     *
*   The most current version of this software is available                 *
*   from netlib@nac.no , send the message `send index from maspar'         *
*                                                                          *
*   REVISIONS:                                                             *
*                                                                          *
****************************************************************************
      implicit none
*     .. Scalar Arguments ..
      CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
      INTEGER            M, N, LDA, LDB
      REAL               ALPHA
*     .. Array Arguments ..
       REAL, array(:,:) :: a,b
      intent(in)    :: a
      intent(inout) :: b
*     ..
*
*  Purpose
*  =======
*
*  STRMM  performs one of the matrix-matrix operations
*
*     B := alpha*op( A )*B,   or   B := alpha*B*op( A ),
*
*  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
*  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
*
*     op( A ) = A   or   op( A ) = A'.
*
*  Parameters
*  ==========
*
*  SIDE   - CHARACTER*1.
*           On entry,  SIDE specifies whether  op( A ) multiplies B from
*           the left or right as follows:
*
*              SIDE = 'L' or 'l'   B := alpha*op( A )*B.
*
*              SIDE = 'R' or 'r'   B := alpha*B*op( A ).
*
*           Unchanged on exit.
*
*  UPLO   - CHARACTER*1.
*           On entry, UPLO specifies whether the matrix A is an upper or
*           lower triangular matrix as follows:
*
*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
*
*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
*
*           Unchanged on exit.
*
*  TRANSA - CHARACTER*1.
*           On entry, TRANSA specifies the form of op( A ) to be used in
*           the matrix multiplication as follows:
*
*              TRANSA = 'N' or 'n'   op( A ) = A.
*
*              TRANSA = 'T' or 't'   op( A ) = A'.
*
*              TRANSA = 'C' or 'c'   op( A ) = A'.
*
*           Unchanged on exit.
*
*  DIAG   - CHARACTER*1.
*           On entry, DIAG specifies whether or not A is unit triangular
*           as follows:
*
*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
*
*              DIAG = 'N' or 'n'   A is not assumed to be unit
*                                  triangular.
*
*           Unchanged on exit.
*
*  M      - INTEGER.
*           On entry, M specifies the number of rows of B. M must be at
*           least zero.
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the number of columns of B.  N must be
*           at least zero.
*           Unchanged on exit.
*
*  ALPHA  - REAL            .
*           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
*           zero then  A is not referenced and  B need not be set before
*           entry.
*           Unchanged on exit.
*
*  A      - REAL array of DIMENSION ( LDA, k ), where k is m
*           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
*           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
*           upper triangular part of the array  A must contain the upper
*           triangular matrix  and the strictly lower triangular part of
*           A is not referenced.
*           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
*           lower triangular part of the array  A must contain the lower
*           triangular matrix  and the strictly upper triangular part of
*           A is not referenced.
*           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
*           A  are not referenced either,  but are assumed to be  unity.
*           Unchanged on exit.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
*           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
*           then LDA must be at least max( 1, n ).
*           Unchanged on exit.
*
*  B      - REAL array of DIMENSION ( LDB, n ).
*           Before entry,  the leading  m by n part of the array  B must
*           contain the matrix  B,  and  on exit  is overwritten  by the
*           transformed matrix.
*
*  LDB    - INTEGER.
*           On entry, LDB specifies the first dimension of B as declared
*           in  the  calling  (sub)  program.   LDB  must  be  at  least
*           max( 1, m ).
*           Unchanged on exit.
*
*
*  Level 3 Blas routine.
*
*  -- Written on 8-February-1989.
*     Jack Dongarra, Argonne National Laboratory.
*     Iain Duff, AERE Harwell.
*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
*     Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. External Subroutines ..
      EXTERNAL           XERBLA
      EXTERNAL           SGEMM
*     .. Intrinsic Functions
      INTRINSIC          MAX
      INTRINSIC          matmul
      INTRINSIC          transpose
*     .. Local Scalars ..
      LOGICAL            LSIDE, NOUNIT, UPPER
      INTEGER            I, INFO, J, K, NROWA
      REAL TEMP
*     .. Local Arrays ..
      integer, array(max(m,n),max(m,n)) :: maska
       REAL,    array(max(m,n),max(m,n)) :: aloc
       REAL,    array(m,n)               :: btmp
*     .. Parameters ..
      REAL               ONE         , ZERO
      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Executable Statements ..
*
cts  cmpf mpl mpl_itrmsk
cmpf ondpu a,b
*     Test the input parameters.
*
      LSIDE  = LSAME( SIDE  , 'L' )
      IF( LSIDE )THEN
         NROWA = M
      ELSE
         NROWA = N
      END IF
      NOUNIT = LSAME( DIAG  , 'N' )
      UPPER  = LSAME( UPLO  , 'U' )
*
      INFO   = 0
      IF(      ( .NOT.LSIDE                ).AND.
     $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
         INFO = 1
      ELSE IF( ( .NOT.UPPER                ).AND.
     $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
         INFO = 2
      ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
         INFO = 3
      ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
     $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
         INFO = 4
      ELSE IF( M  .LT.0               )THEN
         INFO = 5
      ELSE IF( N  .LT.0               )THEN
         INFO = 6
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
         INFO = 9
      ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
         INFO = 11
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'STRMM ', INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( N.EQ.0 ) RETURN
      IF( m.EQ.0 ) RETURN
*
*     And when  alpha.eq.zero.
*
      IF( ALPHA.EQ.ZERO )THEN
         B (1:m,1:n) = ZERO
         RETURN
      END IF
*
*     Start the operations.
*
      IF( LSIDE )THEN
*
*           Form  B := alpha*op(A)*B.
*
      forall (i = 1:m, j = 1:m) maska(i,j) = j-i
cts   CALL MPL_ITRMSK(maska(1:m,1:m),m)
c
        aloc(1:m,1:m) = zero
c
        if (upper) then
           where ( maska(1:m,1:m) .ge. 0 ) aloc(1:m,1:m) = a(1:m,1:m)
	else
           where ( maska(1:m,1:m) .le. 0 ) aloc(1:m,1:m) = a(1:m,1:m)
   	endif
c
	if (.not. nounit) then
	    where ( maska(1:m,1:m) .eq. 0 ) aloc(1:m,1:m) = one
   	endif
c
        IF( LSAME( TRANSA, 'N' ) )THEN
*          b(1:m,1:n)=alpha*matmul(aloc(1:m,1:m),b(1:m,1:n))
cts        CALL SGEMM('N','N',m,n,m,alpha,aloc(1:m,1:m),m
cts  $        ,b(1:m,1:n),m,ZERO,btmp(1:m,1:n),m)
           CALL SGEMM('N','N',m,n,m,alpha,aloc,m
     $        ,b,m,ZERO,btmp,m)
           b = btmp
        ELSE
*         b(1:m,1:n)=alpha*matmul(transpose(aloc(1:m,1:m)),b(1:m,1:n))
cts        CALL SGEMM('T','N',m,n,m,alpha,aloc(1:m,1:m),m
cts  $        ,b(1:m,1:n),m,ZERO,btmp(1:m,1:n),m)
           CALL SGEMM('T','N',m,n,m,alpha,aloc,m
     $        ,b,m,ZERO,btmp,m)
           b = btmp
        ENDIF
      ELSE
*
*           Form  B := alpha*B*op(A).
*
        forall (i = 1:n, j = 1:n) maska(i,j) = j-i
cts   CALL MPL_ITRMSK(maska(1:n,1:n),n)
c
        aloc(1:n,1:n) = zero
c
        if (upper) then
           where ( maska(1:n,1:n) .ge. 0 ) aloc(1:n,1:n) = a(1:n,1:n)
	else
           where ( maska(1:n,1:n) .le. 0 ) aloc(1:n,1:n) = a(1:n,1:n)
   	endif
c
	if (.not. nounit) then
	    where ( maska(1:n,1:n) .eq. 0 ) aloc(1:n,1:n) = one
   	endif
c
        IF( LSAME( TRANSA, 'N' ) )THEN
*          b(1:m,1:n)=alpha*matmul(b(1:m,1:n),aloc(1:n,1:n))
cts        CALL SGEMM('N','N',m,n,n,alpha,b(1:m,1:n),m
cts  $        ,aloc(1:n,1:n),n,ZERO,btmp(1:m,1:n),m)
           CALL SGEMM('N','N',m,n,n,alpha,b,m
     $        ,aloc,n,ZERO,btmp,m)
           b = btmp
        ELSE
*          b(1:m,1:n)=alpha*matmul(b(1:m,1:n),transpose(aloc(1:n,1:n)))
cts        CALL SGEMM('N','T',m,n,n,alpha,b(1:m,1:n),m
cts  $        ,aloc(1:n,1:n),n,ZERO,btmp(1:m,1:n),m)
           CALL SGEMM('N','T',m,n,n,alpha,b,m
     $        ,aloc,n,ZERO,btmp,m)
           b = btmp
        ENDIF
      ENDIF
*
      RETURN
*
*     End of STRMM .
*
      END