C
C             Parallel Sparse BLAS  version 2.2
C   (C) Copyright 2006/2007/2008
C                      Salvatore Filippone    University of Rome Tor Vergata
C                      Alfredo Buttari        University of Rome Tor Vergata
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C   1. Redistributions of source code must retain the above copyright
C      notice, this list of conditions and the following disclaimer.
C   2. Redistributions in binary form must reproduce the above copyright
C      notice, this list of conditions, and the following disclaimer in the
C      documentation and/or other materials provided with the distribution.
C   3. The name of the PSBLAS group or the names of its contributors may
C      not be used to endorse or promote products derived from this
C      software without specific written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
C POSSIBILITY OF SUCH DAMAGE.
C
***********************************************************************
*    ZSRMV modified for SPARKER                                       *
*                                                                     *
*    FUNCTION: Driver for routines performing one of the sparse       *
*              matrix vector operations                               *
*                                                                     *
*                   y = alpha*op(A)*x + beta*y                        *
*                                                                     *
*              where op(A) is one of:                                 *
*                                                                     *
*                  op(A) = A or op(A) = A'  or                        *
*                  op(A) = conjug(A') or                              *
*                  op(A) = lower or upper part of A                   *
*                                                                     *
*              alpha and beta are scalars.                            *
*              The data structure of the matrix is related            *
*              to the scalar computer.                                *
*              This is an internal routine called by:                 *
*              ZCSRMM                                                 *
*                                                                     *
*    ENTRY-POINT = ZSRMV                                              *
*                                                                     *
*                                                                     *
*    INPUT =                                                          *
*                                                                     *
*                                                                     *
*      SYMBOLIC NAME: TRANS                                           *
*      POSITION:      PARAMETER NO 1.                                 *
*      ATTRIBUTES:    CHARACTER*1                                     *
*      VALUES:        'N' 'T' 'C' 'L' 'M' 'U' 'V'                     *
*      DESCRIPTION:   Specifies the form of op(A) to be used in the   *
*                     matrix vector multiplications as follows:       *
*                                                                     *
*                     TRANS = 'N'       ,  op( A ) = A.               *
*                                                                     *
*                     TRANS = 'T'       ,  op( A ) = A'.              *
*                                                                     *
*                     TRANS = 'C'       ,  OP( A ) = conjug(A')       *
*                                                                     *
*                     TRANS = 'L' or 'U',  op( A ) = lower or         *
*                             upper part of A                         *
*                                                                     *
*                     TRANS = 'M' or 'V',  op( A ) = lower or         *
*                             upper part of conjugate of A            *
*                                                                     *
*      SYMBOLIC NAME: DIAG                                            *
*      POSITION:      PARAMETER NO 2.                                 *
*      ATTRIBUTES:    CHARACTER*1                                     *
*      VALUES:        'N' 'U'                                         *
*      DESCRIPTION:                                                   *
*                     Specifies whether or not the matrix A has       *
*                     unit diagonal as follows:                       *
*                                                                     *
*                     DIAG = 'N'  A is not assumed                    *
*                            to have unit diagonal                    *
*                                                                     *
*                     DIAG = 'U'  A is assumed                        *
*                            to have unit diagonal.                   *
*                                                                     *
*                     WARNING: it is the caller's responsibility      *
*                     to ensure that if the matrix has unit           *
*                     diagonal, there are no elements of the          *
*                     diagonal are stored in the arrays AS and JA.    *
*                                                                     *
*       SYMBOLIC NAME: M                                              *
*       POSITION:      PARAMETER NO 3.                                *
*       ATTRIBUTES:    INTEGER*4.                                     *
*       VALUES:        M >= 0                                         *
*       DESCRIPTION:   Number of rows of the matrix op(A).            *
*                                                                     *
*       SYMBOLIC NAME: N                                              *
*       POSITION:      PARAMETER NO 4.                                *
*       ATTRIBUTES:    INTEGER*4.                                     *
*       VALUES:        N >= 0                                         *
*       DESCRIPTION:   Number of columns of the matrix op(A)          *
*                                                                     *
*       SYMBOLIC NAME: ALPHA                                          *
*       POSITION:      PARAMETER NO 5.                                *
*       ATTRIBUTES:    COMPLEX*16.                                    *
*       VALUES:        any.                                           *
*       DESCRIPTION:   Specifies the scalar alpha.                    *
*                                                                     *
*                                                                     *
*       SYMBOLIC NAME: AS                                             *
*       POSITION:      PARAMETER NO 6.                                *
*       ATTRIBUTES:    COMPLEX*16: ARRAY(IA(M+1)-1)                   *
*       VALUES:        ANY                                            *
*       DESCRIPTION:   Array containing the non zero coefficients of  *
*                      the sparse matrix op(A).                       *
*                                                                     *
*       SYMBOLIC NAME: JA                                             *
*       POSITION:      PARAMETER NO 7.                                *
*       ATTRIBUTES:    INTEGER*4: ARRAY(IA(M+1)-1)                    *
*       VALUES:        0 < JA(I) <= M                                 *
*       DESCRIPTION:   Array containing the column number of the      *
*                      nonzero coefficients stored in array AS.       *
*                                                                     *
*       SYMBOLIC NAME: IA                                             *
*       POSITION:      PARAMETER NO 8.                                *
*       ATTRIBUTES:    INTEGER*4: ARRAY(*)                            *
*       VALUES:        IA(I) > 0                                      *
*       DESCRIPTION:   Contains the pointers for the beginning of     *
*                      each rows.                                     *
*                                                                     *
*                                                                     *
*       SYMBOLIC NAME: X                                              *
*       POSITION:      PARAMETER NO 9.                                *
*       ATTRIBUTES:    COMPLEX*16 ARRAY(N)                            *
*                      (or ARRAY(M) when op(A) = A')                  *
*       VALUES:        any.                                           *
*       DESCRIPTION:   Contains the values of the vector to be        *
*                      multiplied by the matrix A.                    *
*                                                                     *
*       SYMBOLIC NAME: BETA                                           *
*       POSITION:      PARAMETER NO 10.                               *
*       ATTRIBUTES:    COMPLEX*16.                                    *
*       VALUES:        any.                                           *
*       DESCRIPTION:   Specifies the scalar beta.                     *
*                                                                     *
*       SYMBOLIC NAME: Y                                              *
*       POSITION:      PARAMETER NO 11.                               *
*       ATTRIBUTES:    COMPLEX*16 ARRAY(M)                            *
*                      (or ARRAY(N) when op(A) = A')                  *
*       VALUES:        any.                                           *
*       DESCRIPTION:   Contains the values of the vector to be        *
*                      updated by the matrix-vector multiplication.   *
*                                                                     *
*       SYMBOLIC NAME: WORK                                           *
*       POSITION:      PARAMETER NO 12.                               *
*       ATTRIBUTES:    COMPLEX*16 ARRAY(M)                            *
*                      (or ARRAY(N) when op(A) = A')                  *
*       VALUES:        any.                                           *
*       DESCRIPTION:   Work area available to the program. It is used *
*                      only when TRANS = 'T'.                         *
*                                                                     *
*  OUTPUT =                                                           *
*                                                                     *
*                                                                     *
*       SYMBOLIC NAME: Y                                              *
*       POSITION:      PARAMETER NO 11.                               *
*       ATTRIBUTES:    COMPLEX*16 ARRAY(M)                            *
*                      (or ARRAY(N) when op(A) = A')                  *
*       VALUES:        any.                                           *
*       DESCRIPTION:   Contains the values of the vector              *
*                      updated by the matrix-vector multiplication.   *
*                                                                     *
*                                                                     *
***********************************************************************
C
      SUBROUTINE ZSRMV (TRANS,DIAG,M,N,ALPHA,AS,JA,IA,X,BETA,Y,WORK)
      use psb_const_mod
      use psb_string_mod
C     .. Parameters ..
      complex(psb_dpk_) ONE, ZERO
      PARAMETER (ONE=(1.0D0, 0.0D0), ZERO=(0.0D0, 0.0D0))
C     .. Scalar Arguments ..
      complex(psb_dpk_) ALPHA, BETA
      INTEGER    M, N
      CHARACTER  DIAG, TRANS
C     .. Array Arguments ..
      complex(psb_dpk_) AS(*), WORK(*), X(*), Y(*)
      INTEGER    IA(*), JA(*)
C     .. Local Scalars ..
      complex(psb_dpk_) ACC
      INTEGER    I, J, K, NCOLA, NROWA,DUM
      LOGICAL    SYM, TRA, COTRA, UNI
C     .. Executable Statements ..
C
      UNI = psb_toupper(DIAG).EQ.'U'
C
C     .. Not simmetric matrix
      TRA   = psb_toupper(TRANS).EQ.'T'
      COTRA = psb_toupper(TRANS).EQ.'C'

C     .. Symmetric matrix upper or lower 
      SYM = (psb_toupper(TRANS).EQ.'L').OR.
     +  (psb_toupper(TRANS).EQ.'U').OR.
     +  (psb_toupper(TRANS).EQ.'M').OR.
     +  (psb_toupper(TRANS).EQ.'V')
C
      IF (.NOT.(TRA.OR.COTRA)) THEN
        NROWA = M
        NCOLA = N
      ELSE
        NROWA = N
        NCOLA = M
      END IF                    !(....(CO)TRA)

      IF (ALPHA.EQ.ZERO) THEN
        IF (BETA.EQ.ZERO) THEN
          DO I = 1, M
            Y(I) = ZERO
          ENDDO
        ELSE
          DO 20 I = 1, M
            Y(I) = BETA*Y(I)
 20       CONTINUE
        ENDIF
        RETURN
      END IF

C
      IF (SYM) THEN
        IF (UNI) THEN
C
C              ......Symmetric with unitary diagonal.......
C              ....OK!!
C              To be optimized
          
          IF (BETA.NE.ZERO) THEN
            DO 40 I = 1, M
C
C                 Product for diagonal elements
C
              Y(I) = BETA*Y(I) + ALPHA*X(I)
 40         CONTINUE
          ELSE
            DO I = 1, M
              Y(I) = ALPHA*X(I)
            ENDDO
          ENDIF

C              Product for other elements
          IF ((psb_toupper(TRANS).EQ.'L').OR.
     +      (psb_toupper(TRANS).EQ.'U')) THEN
            DO 80 I = 1, M
              ACC = ZERO
              DO 60 J = IA(I), IA(I+1) - 1
                K = JA(J)
                Y(K) = Y(K) + ALPHA*AS(J)*X(I)
                ACC = ACC + AS(J)*X(K)
 60           CONTINUE
              Y(I) = Y(I) + ALPHA*ACC
 80         CONTINUE
          ELSE                  ! Perform computations on conjug(A)
            DO 82 I = 1, M
              ACC = ZERO
              DO 81 J = IA(I), IA(I+1) - 1
                K = JA(J)
                Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I)
                ACC = ACC + CONJG(AS(J)) * X(K)
 81           CONTINUE
              Y(I) = Y(I) + ALPHA*ACC
 82         CONTINUE
          ENDIF
C
        ELSE IF ( .NOT. UNI) THEN
C
C            Check if matrix is lower or upper
C
          IF ((psb_toupper(TRANS).EQ.'L').OR.
     +      (psb_toupper(TRANS).EQ.'M')) THEN
C
C               LOWER CASE: diagonal element is the last element of row
C               ....OK!

            IF (BETA.NE.ZERO) THEN
              DO 100 I = 1, M
                Y(I) = BETA*Y(I)
 100          CONTINUE
            ELSE
              DO I = 1, M
                Y(I) = ZERO
              ENDDO
            ENDIF

            IF  (psb_toupper(TRANS).EQ.'L') THEN
              DO 140 I = 1, M
                ACC = ZERO
                DO 120 J = IA(I), IA(I+1) - 1 ! it was -2
                  K = JA(J)
C   
C                    To be optimized
C   
                  IF (K.NE.I) THEN !for symmetric element 
                    Y(K) = Y(K) + ALPHA*AS(J)*X(I)
                  ENDIF
                  ACC = ACC + AS(J)*X(K)
 120            CONTINUE

                Y(I) = Y(I) + ALPHA*ACC 
 140          CONTINUE
            ELSE                ! Perform computations on conjug(A)
              DO 142 I = 1, M
                ACC = ZERO
                DO 141 J = IA(I), IA(I+1) - 1 ! it was -2
                  K = JA(J)
C   
C                    To be optimized
C   
                  IF (K.NE.I) THEN !for symmetric element 
                    Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I)
                  ENDIF
                  ACC = ACC + CONJG(AS(J)) * X(K)
 141            CONTINUE
                Y(I) = Y(I) + ALPHA * ACC 
 142          CONTINUE

            ENDIF

          ELSE                  ! ....Trans<>L, M
C
C              UPPER CASE
C              ....OK!!
C
            IF (BETA.NE.ZERO) THEN
              DO 160 I = 1, M
                Y(I) = BETA*Y(I)
 160          CONTINUE
            ELSE
              DO I = 1, M
                Y(I) = ZERO
              ENDDO
            ENDIF
            IF (psb_toupper(TRANS).EQ.'U') THEN
              DO 200 I = 1, M
                ACC = ZERO
                DO 180 J = IA(I) , IA(I+1) - 1 ! removed +1
                  K = JA(J)
C   
C                    To be optimized
C   
                  IF(K.NE.I) THEN
                    Y(K) = Y(K) + ALPHA*AS(J)*X(I)
                  ENDIF
                  ACC = ACC + AS(J)*X(K)
 180            CONTINUE
                Y(I) = Y(I) + ALPHA*ACC
 200          CONTINUE
            ELSE                ! Perform computations on conjug(A)
              DO 202 I = 1, M
                ACC = ZERO
                DO 201 J = IA(I) , IA(I+1) - 1 ! removed +1
                  K = JA(J)
C   
C                    To be optimized
C   
                  IF(K.NE.I) THEN
                    Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I)
                  ENDIF
                  ACC = ACC + CONJG(AS(J)) * X(K)
 201            CONTINUE
                Y(I) = Y(I) + ALPHA*ACC
 202          CONTINUE
            ENDIF
          END IF                ! ......TRANS=='L'
        END IF                  ! ......Not UNI
C
      ELSE IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN !......NOT SYM

        IF ( .NOT. UNI) THEN      
C
C          .......General Not Unit, No Traspose
C

          IF (BETA.NE.ZERO) THEN
            DO 240 I = 1, M
              ACC = ZERO
              DO 220 J = IA(I), IA(I+1) - 1
                ACC = ACC + AS(J)*X(JA(J))
 220          CONTINUE
              Y(I) = ALPHA*ACC + BETA*Y(I)
 240        CONTINUE
          ELSE
            DO I = 1, M
              ACC = ZERO
              DO J = IA(I), IA(I+1) - 1
                ACC = ACC + AS(J)*X(JA(J))
              ENDDO
              Y(I) = ALPHA*ACC
            ENDDO
          ENDIF
C
        ELSE IF (UNI) THEN
C
          IF (BETA.NE.ZERO) THEN
            DO 280 I = 1, M
              ACC = ZERO
              DO 260 J = IA(I), IA(I+1) - 1
                ACC = ACC + AS(J)*X(JA(J))
 260          CONTINUE
              Y(I) = ALPHA*(ACC+X(I)) + BETA*Y(I)
 280        CONTINUE
          ELSE                  !(BETA.EQ.ZERO)
            DO I = 1, M
              ACC = ZERO
              DO J = IA(I), IA(I+1) - 1
                ACC = ACC + AS(J)*X(JA(J))
              ENDDO
              Y(I) = ALPHA*(ACC+X(I))
            ENDDO
          ENDIF
        END IF                  !....End Testing on UNI
C
      ELSE IF (TRA) THEN        !....Else on SYM (swapped M and N)
C
        IF ( .NOT. UNI) THEN
C
          IF (BETA.NE.ZERO) THEN
            DO 300 I = 1, M
              Y(I) = BETA*Y(I)
 300        CONTINUE
          ELSE                  !(BETA.EQ.ZERO)
            DO I = 1, M
              Y(I) = ZERO
            ENDDO
          ENDIF
C
        ELSE IF (UNI) THEN
C

          IF (BETA.NE.ZERO) THEN
            DO 320 I = 1, M
              Y(I) = BETA*Y(I) + ALPHA*X(I)
 320        CONTINUE
          ELSE                  !(BETA.EQ.ZERO)
            DO I = 1, M
              Y(I) = ALPHA*X(I)
            ENDDO
          ENDIF

C
        END IF                  !....UNI
C
        IF (ALPHA.EQ.ONE) THEN
C
          DO 360 I = 1, N
            DO 340 J = IA(I), IA(I+1) - 1
              K = JA(J)
              Y(K) = Y(K) + AS(J)*X(I)
 340        CONTINUE
 360      CONTINUE
C
        ELSE IF (ALPHA.EQ.-ONE) THEN
C
          DO 400 I = 1, n
            DO 380 J = IA(I), IA(I+1) - 1
              K = JA(J)
              Y(K) = Y(K) - AS(J)*X(I)
 380        CONTINUE
 400      CONTINUE
C
        ELSE                    !.....Else on TRA
C
C           This work array is used for efficiency
C
          DO 420 I = 1, N
            WORK(I) = ALPHA*X(I)
 420      CONTINUE
C
          DO 460 I = 1, n
            DO 440 J = IA(I), IA(I+1) - 1
              K = JA(J)
              Y(K) = Y(K) + AS(J)*WORK(I)
 440        CONTINUE
 460      CONTINUE
C
        END IF                  !.....End testing on ALPHA

      ELSE IF (COTRA) THEN      !....Else on SYM (swapped M and N)
C
        IF ( .NOT. UNI) THEN
C
          IF (BETA.NE.ZERO) THEN
            DO 500 I = 1, M
              Y(I) = BETA*Y(I)
 500        CONTINUE
          ELSE                  !(BETA.EQ.ZERO)
            DO I = 1, M
              Y(I) = ZERO
            ENDDO
          ENDIF
C
        ELSE IF (UNI) THEN
C

          IF (BETA.NE.ZERO) THEN
            DO 520 I = 1, M
              Y(I) = BETA*Y(I) + ALPHA*X(I)
 520        CONTINUE
          ELSE                  !(BETA.EQ.ZERO)
            DO I = 1, M
              Y(I) = ALPHA*X(I)
            ENDDO
          ENDIF

C
        END IF                  !....UNI
C
        IF (ALPHA.EQ.ONE) THEN
C
          DO 560 I = 1, N
            DO 540 J = IA(I), IA(I+1) - 1
              K = JA(J)
              Y(K) = Y(K) + CONJG(AS(J))*X(I)
 540        CONTINUE
 560      CONTINUE
C
        ELSE IF (ALPHA.EQ.-ONE) THEN
C
          DO 600 I = 1, n
            DO 580 J = IA(I), IA(I+1) - 1
              K = JA(J)
              Y(K) = Y(K) - CONJG(AS(J))*X(I)
 580        CONTINUE
 600      CONTINUE
C
        ELSE                    !.....Else on TRA
C
C           This work array is used for efficiency
C
          DO 620 I = 1, N
            WORK(I) = ALPHA*X(I)
 620      CONTINUE
C
          DO 660 I = 1, n
            DO 640 J = IA(I), IA(I+1) - 1
              K = JA(J)
              Y(K) = Y(K) + CONJG(AS(J))*WORK(I)
 640        CONTINUE
 660      CONTINUE
C
        END IF                  !.....End testing on ALPHA

      END IF                    !.....End testing on SYM      
C
      RETURN
C
C     END OF ZSRMV
C
      END