C
C             Parallel Sparse BLAS  v2.0
C   (C) Copyright 2006 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
C 
***********************************************************************         
*                PROCEDURAL LOGIC SECTION                             *         
*      SUBROUTINE DJADMV (DIAG,NROW,NCOL,ALPHA,NG,A,KA,JA,IA,X,BETA,Y) *         
*      DOUBLE PRECISION  ZERO                                         *         
*      PARAMETER         (ZERO=0.0D0)                                 *         
*      DOUBLE PRECISION  ACC                                          *         
*      INTEGER           I, J, K, IPX, IPG                            *         
*      LOGICAL           UNI                                          *         
*C     .. Executable Statements ..                                    *         
*C                                                                    *         
*C                                                                    *         
*      IF (DIAG.EQ.'U') THEN                                          *         
*            DO  10 I = 1, M                                          *         
*                Y(I) = BETA*Y(I) + ALPHA*X(I)                        *         
*   10       CONTINUE                                                 *         
*      ELSE                                                           *         
*            DO 20 I = 1, M                                           *         
*               Y(I) = BETA*Y(I)                                      *         
*   20       CONTINUE                                                 *         
*      END IF                                                         *         
*                                                                     *         
*      IF (ALPHA.EQ.ZERO) THEN                                        *         
*         RETURN                                                      *         
*      END IF                                                         *         
*C                                                                    *         
*C                                                                    *         
*C        DO 200 IPG = 1, NG                                          *         
*            DO 50 K = IA(2,IPG), IA(3,IPG)-1                         *         
*               IPX = IA(1,IPG)                                       *         
*               DO 40 I = JA(K), JA(K+1) - 1                          *         
*                  Y(IPX) = Y(IPX) + ALPHA*A(I)*X(KA(I))              *         
*                  IPX = IPX + 1                                      *         
*   40          CONTINUE                                              *         
*   50       CONTINUE                                                 *         
*C                                                                    *         
*C                                                                    *         
*            IPX = IA(1,IPG)                                          *         
*            DO 70 K = IA(3,IPG), IA(2,IPG+1)-1                       *         
*               DO 60 I = JA(K), JA(K+1) - 1                          *         
*                  Y(IPX) = Y(IPX) + ALPHA*A(I)*X(KA(I))              *         
*   60          CONTINUE                                              *         
*               IPX = IPX + 1                                         *         
*   70       CONTINUE                                                 *         
*  200    CONTINUE                                                    *         
*C                                                                    *         
*      RETURN                                                         *         
*C                                                                    *         
*C                                                                    *         
*      END                                                            *         
*                                                                     *         
*                                                                     *         
***********************************************************************       
      SUBROUTINE DJADMV (DIAG,NROW,NCOL,ALPHA,NG,A,KA,JA,IA,X,
     +     BETA,Y,IERROR)
      IMPLICIT NONE
      DOUBLE PRECISION  A(*),X(*),Y(*),ALPHA,BETA,ZERO                         
      INTEGER           IA(3,*),KA(*),JA(*),NCOL,NROW,NG,IERROR
      CHARACTER         DIAG                                                   
      PARAMETER         (ZERO=0.0D0)                                           
      INTEGER           I, K, IPX, IPG, I0, IN
      INTEGER           NPG
      DOUBLE PRECISION  Y0, Y1, Y2, Y3, Y4, Y5, Y6, Y7, 
     +  Y8, Y9, Y10, Y11, Y12, Y13, Y14, Y15
c      .. Executable Statements ..                                              
c                                                                               
c                                                              
      IERROR=0                 
      IF (DIAG.EQ.'U') THEN
        IF (BETA.EQ.ZERO) THEN
          DO I = 1, NROW                                                 
            Y(I) = ALPHA*X(I)                                  
          ENDDO
        ELSE
          DO  10 I = 1, NROW                                                 
            Y(I) = BETA*Y(I) + ALPHA*X(I)                                  
 10       CONTINUE  
        ENDIF
      ELSE                                                                     
        IF (BETA.EQ.ZERO) THEN
          DO I = 1, NROW                                                  
            Y(I) = 0.D0
          ENDDO
        ELSE
          DO 20 I = 1, NROW                                                  
            Y(I) = BETA*Y(I)                                                
 20       CONTINUE                                                           
        END IF
      ENDIF

      IF (ALPHA.EQ.ZERO) THEN                                                  
        RETURN                                                                
      END IF                                                                   
c                                                                               
c                                                                               
      DO 200 IPG = 1, NG                                                    
        K   = IA(2,IPG)
        NPG = JA(K+1)-JA(K)

        IF (NPG.EQ.4) THEN
          IPX = IA(1,IPG)                                                 
          Y0  = ZERO
          Y1  = ZERO
          Y2  = ZERO
          Y3  = ZERO
          K   = IA(2,IPG)
          I0  = JA(K)
          K   = IA(3,IPG)-1
          IN  = JA(K)      
          DO I = I0, IN+3, 4
            Y0 = Y0 + A(I+0)*X(KA(I+0))  
            Y1 = Y1 + A(I+1)*X(KA(I+1))  
            Y2 = Y2 + A(I+2)*X(KA(I+2))  
            Y3 = Y3 + A(I+3)*X(KA(I+3))              
          ENDDO
          Y(IPX+0)  = Y(IPX+0) + ALPHA*Y0 
          Y(IPX+1)  = Y(IPX+1) + ALPHA*Y1 
          Y(IPX+2)  = Y(IPX+2) + ALPHA*Y2  
          Y(IPX+3)  = Y(IPX+3) + ALPHA*Y3 
          
        ELSE IF (NPG.EQ.5) THEN
          
          IPX = IA(1,IPG)                                                 
          Y0  = ZERO
          Y1  = ZERO
          Y2  = ZERO
          Y3  = ZERO
          Y4  = ZERO
          K   = IA(2,IPG)
          I0  = JA(K)
          K   = IA(3,IPG)-1
          IN  = JA(K)
          DO I = I0, IN+4, 5
            Y0 = Y0 + A(I+0)*X(KA(I+0))  
            Y1 = Y1 + A(I+1)*X(KA(I+1))  
            Y2 = Y2 + A(I+2)*X(KA(I+2))  
            Y3 = Y3 + A(I+3)*X(KA(I+3))  
            Y4 = Y4 + A(I+4)*X(KA(I+4))  
          ENDDO
          Y(IPX+0)  = Y(IPX+0) + ALPHA*Y0 
          Y(IPX+1)  = Y(IPX+1) + ALPHA*Y1 
          Y(IPX+2)  = Y(IPX+2) + ALPHA*Y2  
          Y(IPX+3)  = Y(IPX+3) + ALPHA*Y3 
          Y(IPX+4)  = Y(IPX+4) + ALPHA*Y4 

        ELSE IF (NPG.EQ.6) THEN

          IPX = IA(1,IPG)                                                 
          Y0  = ZERO
          Y1  = ZERO
          Y2  = ZERO
          Y3  = ZERO
          Y4  = ZERO
          Y5  = ZERO
          K   = IA(2,IPG)
          I0  = JA(K)
          K   = IA(3,IPG)-1
          IN  = JA(K)
          DO I = I0, IN+5, 6
            Y0 = Y0 + A(I+0)*X(KA(I+0))  
            Y1 = Y1 + A(I+1)*X(KA(I+1))  
            Y2 = Y2 + A(I+2)*X(KA(I+2))  
            Y3 = Y3 + A(I+3)*X(KA(I+3))  
            Y4 = Y4 + A(I+4)*X(KA(I+4))  
            Y5 = Y5 + A(I+5)*X(KA(I+5))  
          ENDDO
          Y(IPX+0)  = Y(IPX+0) + ALPHA*Y0 
          Y(IPX+1)  = Y(IPX+1) + ALPHA*Y1 
          Y(IPX+2)  = Y(IPX+2) + ALPHA*Y2  
          Y(IPX+3)  = Y(IPX+3) + ALPHA*Y3 
          Y(IPX+4)  = Y(IPX+4) + ALPHA*Y4 
          Y(IPX+5)  = Y(IPX+5) + ALPHA*Y5 

        ELSE IF (NPG.EQ.7) THEN 

          IPX = IA(1,IPG)                                                 
          Y0  = ZERO
          Y1  = ZERO
          Y2  = ZERO
          Y3  = ZERO
          Y4  = ZERO
          Y5  = ZERO
          Y6  = ZERO
          K   = IA(2,IPG)
          I0  = JA(K)
          K   = IA(3,IPG)-1
          IN  = JA(K)
          DO I = I0, IN+6, 7
            Y0 = Y0 + A(I+0)*X(KA(I+0))  
            Y1 = Y1 + A(I+1)*X(KA(I+1))  
            Y2 = Y2 + A(I+2)*X(KA(I+2))  
            Y3 = Y3 + A(I+3)*X(KA(I+3))  
            Y4 = Y4 + A(I+4)*X(KA(I+4))  
            Y5 = Y5 + A(I+5)*X(KA(I+5))  
            Y6 = Y6 + A(I+6)*X(KA(I+6))  
          ENDDO
          Y(IPX+0)  = Y(IPX+0) + ALPHA*Y0 
          Y(IPX+1)  = Y(IPX+1) + ALPHA*Y1 
          Y(IPX+2)  = Y(IPX+2) + ALPHA*Y2  
          Y(IPX+3)  = Y(IPX+3) + ALPHA*Y3 
          Y(IPX+4)  = Y(IPX+4) + ALPHA*Y4 
          Y(IPX+5)  = Y(IPX+5) + ALPHA*Y5 
          Y(IPX+6)  = Y(IPX+6) + ALPHA*Y6  

        ELSE IF (NPG.EQ.8) THEN

          IPX = IA(1,IPG)                                                 
          Y0  = ZERO
          Y1  = ZERO
          Y2  = ZERO
          Y3  = ZERO
          Y4  = ZERO
          Y5  = ZERO
          Y6  = ZERO
          Y7  = ZERO
          K   = IA(2,IPG)
          I0  = JA(K)
          K   = IA(3,IPG)-1
          IN  = JA(K)
          DO I = I0, IN+7, 8
            Y0 = Y0 + A(I+0)*X(KA(I+0))  
            Y1 = Y1 + A(I+1)*X(KA(I+1))  
            Y2 = Y2 + A(I+2)*X(KA(I+2))  
            Y3 = Y3 + A(I+3)*X(KA(I+3))  
            Y4 = Y4 + A(I+4)*X(KA(I+4))  
            Y5 = Y5 + A(I+5)*X(KA(I+5))  
            Y6 = Y6 + A(I+6)*X(KA(I+6))  
            Y7 = Y7 + A(I+7)*X(KA(I+7))  
          ENDDO           
          Y(IPX+0)  = Y(IPX+0) + ALPHA*Y0 
          Y(IPX+1)  = Y(IPX+1) + ALPHA*Y1 
          Y(IPX+2)  = Y(IPX+2) + ALPHA*Y2  
          Y(IPX+3)  = Y(IPX+3) + ALPHA*Y3 
          Y(IPX+4)  = Y(IPX+4) + ALPHA*Y4 
          Y(IPX+5)  = Y(IPX+5) + ALPHA*Y5 
          Y(IPX+6)  = Y(IPX+6) + ALPHA*Y6  
          Y(IPX+7)  = Y(IPX+7) + ALPHA*Y7 

        ELSE IF (NPG.EQ.16) THEN
          
          IPX = IA(1,IPG)                                                 
          Y0  = ZERO
          Y1  = ZERO
          Y2  = ZERO
          Y3  = ZERO
          Y4  = ZERO
          Y5  = ZERO
          Y6  = ZERO
          Y7  = ZERO
          Y8  = ZERO
          Y9  = ZERO
          Y10  = ZERO
          Y11  = ZERO
          Y12  = ZERO
          Y13  = ZERO
          Y14  = ZERO
          Y15  = ZERO
          K   = IA(2,IPG)
          I0  = JA(K)
          K   = IA(3,IPG)-1
          IN  = JA(K)
          DO I = I0, IN+15, 16
            Y0 = Y0 + A(I+0)*X(KA(I+0))  
            Y1 = Y1 + A(I+1)*X(KA(I+1))  
            Y2 = Y2 + A(I+2)*X(KA(I+2))  
            Y3 = Y3 + A(I+3)*X(KA(I+3))  
            Y4 = Y4 + A(I+4)*X(KA(I+4))  
            Y5 = Y5 + A(I+5)*X(KA(I+5))  
            Y6 = Y6 + A(I+6)*X(KA(I+6))  
            Y7 = Y7 + A(I+7)*X(KA(I+7))  
            Y8 = Y8 + A(I+8)*X(KA(I+8))  
            Y9 = Y9 + A(I+9)*X(KA(I+9))  
            Y10 = Y10 + A(I+10)*X(KA(I+10))  
            Y11 = Y11 + A(I+11)*X(KA(I+11))  
            Y12 = Y12 + A(I+12)*X(KA(I+12))  
            Y13 = Y13 + A(I+13)*X(KA(I+13))  
            Y14 = Y14 + A(I+14)*X(KA(I+14))  
            Y15 = Y15 + A(I+15)*X(KA(I+15))  
          ENDDO
          Y(IPX+0)  = Y(IPX+0) + ALPHA*Y0 
          Y(IPX+1)  = Y(IPX+1) + ALPHA*Y1 
          Y(IPX+2)  = Y(IPX+2) + ALPHA*Y2  
          Y(IPX+3)  = Y(IPX+3) + ALPHA*Y3 
          Y(IPX+4)  = Y(IPX+4) + ALPHA*Y4 
          Y(IPX+5)  = Y(IPX+5) + ALPHA*Y5 
          Y(IPX+6)  = Y(IPX+6) + ALPHA*Y6  
          Y(IPX+7)  = Y(IPX+7) + ALPHA*Y7 
          Y(IPX+8)  = Y(IPX+8) + ALPHA*Y8 
          Y(IPX+9)  = Y(IPX+9) + ALPHA*Y9 
          Y(IPX+10)  = Y(IPX+10) + ALPHA*Y10  
          Y(IPX+11)  = Y(IPX+11) + ALPHA*Y11 
          Y(IPX+12)  = Y(IPX+12) + ALPHA*Y12 
          Y(IPX+13)  = Y(IPX+13) + ALPHA*Y13 
          Y(IPX+14)  = Y(IPX+14) + ALPHA*Y14  
          Y(IPX+15)  = Y(IPX+15) + ALPHA*Y15

        ELSE

          DO  K = IA(2,IPG), IA(3,IPG)-1                                   
            IPX = IA(1,IPG)                                                 
            DO  I = JA(K), JA(K+1) - 1                                    
              Y(IPX) = Y(IPX) + ALPHA*A(I)*X(KA(I))                        
              IPX = IPX + 1                                                
            ENDDO
          ENDDO
        END IF

c        CSR Product

        IPX = IA(1,IPG)                            
        DO 70 K = IA(3,IPG), IA(2,IPG+1)-1         
          DO 60 I = JA(K), JA(K+1) - 1
            Y(IPX) = Y(IPX) + ALPHA*A(I)*X(KA(I))
 60       CONTINUE                                
          IPX = IPX + 1                           
 70     CONTINUE                                   
 200  CONTINUE                                                              
c                                                                               
      RETURN                                                                   
      END