Merged some serial modifications from therm sample version.

psblas3-type-indexed
Salvatore Filippone 19 years ago
parent a4262fe6b4
commit f6ea0ade44

@ -12,7 +12,7 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsdp.o psb_dcsmm.o psb_dcsmv.o \
psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \
psb_zcoins.o psb_zcsprt.o psb_zneigh.o psb_ztransp.o psb_ztransc.o\
psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspinfo.o psb_zspscal.o\
psb_getifield.o psb_setifield.o
psb_getifield.o psb_setifield.o psb_update_mod.o
INCDIRS = -I ../../lib -I .
@ -25,6 +25,7 @@ lib: auxd cood csrd jadd f77d dpd lib1
lib1: $(FOBJS)
psb_dcoins.o psb_zcoins.o: psb_update_mod.o
auxd:
(cd aux; make lib)

@ -1,76 +0,0 @@
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
SUBROUTINE DCOOZERO(M,N,DESCRA,A,IA1,IA2,
+ INFOA,IA,JA,MZ,NZ,IERROR)
C
C This subroutione performs the operation:
C
C A(IA : IA + MZ - 1, JA : JA + NZ - 1) = 0
C
C This isn't accomplished by removing elements
C from sparse matrix representation, but assigning them
C the zero value.
C Columns are supposed to be ordered
C into the same row. This subroutine will
C not work properly otherwise.
C
IMPLICIT NONE
C .. Scalar Arguments ..
INTEGER M,N,IA,JA,MZ,NZ,IERROR
C .. Array Arguments ..
INTEGER IA1(*),IA2(*),INFOA(*)
CHARACTER DESCRA*11
DOUBLE PRECISION A(*)
C .. Local scalars ..
INTEGER I, J, JBEGIN, JEND, AUX, NNZ
DOUBLE PRECISION
IERROR=0
IF (((JA + NZ - 1) .GT. N) .OR.
+ ((IA + MZ - 1) .GT. M) .OR.
+ (IA .LT. 1) .OR. (JA .LT. 1)) THEN
IERROR = 1
GOTO 9999
ENDIF
NNZ = INFOA(1)
I = 1
DO WHILE ((IA1(I).LT.IA).AND.(I.LE.NNZ))
I = I + 1
ENDDO
DO WHILE ((IA1(I).LE.(IA+MZ-1)).AND.(I.LE.NNZ))
IF ((JA.LE.IA2(I)).AND.(IA2(I).LE.(JA+NZ-1))) THEN
A(I) = 0.0D0
ENDIF
I = I + 1
ENDDO
RETURN
END

@ -1,134 +0,0 @@
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
SUBROUTINE DCRUPDATE(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
C
C .. Matrix A to be updated is required to be stored with
C .. column indices belonging to the same row ordered.
C .. Block H to be inserted don't need to be stored in such way.
C
C Flag = 0: put elements to 0.0D0;
C Flag = 1: replace elements with new value;
C Flag = 2: sum block value to elements;
C
IMPLICIT NONE
C .. Scalar Arguments ..
INTEGER IA, JA, IH, JH, M, N,
+ IERROR, FLAG, LIWORK
C .. Array Arguments ..
INTEGER IA1(*),IA2(*),IH1(*),IH2(*),
+ INFOA(*),INFOH(*),IWORK(*),
+ GLOB_TO_LOC(*)
CHARACTER DESCRA*11,DESCRH*11
DOUBLE PRECISION A(*),H(*)
C .. Local scalars ..
INTEGER I, J, K, XBLCK, XMATR,
+ AUX, AUX1, AUX2, AUX3,
+ LOC_COLUMN, LOC_POINTER, SHIFT1,
+ SHIFT2
C .. Local arrays ..
IERROR = 0
DO I = 1, M
XBLCK = IH + I - 1
XMATR = IA + I - 1
SHIFT1 = IA2(XMATR + 1) - IA2(XMATR)
SHIFT2 = 2 * SHIFT1
C If columns are already sorted, return point is 100
C CALL MRGSRT(IA2(XMATR + 1) - IA2(XMATR),
C + IA1(IA2(XMATR)),
C + IWORK(SHIFT2 + 1),
C + *100)
GOTO 100
K = IWORK(SHIFT2 + 1)
C If columns have been sorted by mrgsrt
DO J = 1, IA2(XMATR + 1) - IA2(XMATR)
IWORK(J) = IA1(IA2(XMATR) - 1 + K)
IWORK(SHIFT1 + J) = IA2(XMATR) - 1 + K
K = IWORK(SHIFT2 + 1 + K)
ENDDO
GOTO 101
C Else
100 CONTINUE
DO J = IA2(XMATR), IA2(XMATR + 1) - 1
AUX = J - IA2(XMATR) + 1
IWORK(AUX) = IA1(J)
IWORK(SHIFT1 + AUX) = J
ENDDO
C End If
C Now IWORK(1: .. ) contains ordered column indices
C and IWORK(SHIFT1 + 1: .. ) contains position of those
C indices in the stored matrix data structures.
101 CONTINUE
DO J = IH2(XBLCK), IH2(XBLCK + 1) - 1
IF ((JH .LE. IH1(J)) .AND.
+ (IH1(J) .LE. (JH + N - 1))) THEN
LOC_COLUMN = GLOB_TO_LOC(JA - JH + IH1(J))
C Binary search
AUX1 = 1
AUX2 = IA2(XMATR + 1) - IA2(XMATR)
DO K = 1, IA2(XMATR + 1) - IA2(XMATR)
AUX = (AUX1 + AUX2) / 2
IF (LOC_COLUMN .GT. IWORK(AUX)) THEN
AUX1 = AUX + 1
ELSE
AUX2 = AUX - 1
ENDIF
IF ((LOC_COLUMN .EQ. IWORK(AUX)) .OR.
+ (AUX1 .GT. AUX2))
+ EXIT
ENDDO
IF (LOC_COLUMN .EQ. IWORK(AUX)) THEN
LOC_POINTER = IWORK(SHIFT1 + AUX)
ELSE
IERROR = 1
GOTO 9999
ENDIF
IF (FLAG .EQ. 0) THEN
A(LOC_POINTER) = 0.0D0
ELSE IF (FLAG .EQ. 1) THEN
A(LOC_POINTER) = H(J)
ELSE IF (FLAG .EQ. 2) THEN
A(LOC_POINTER) = A(LOC_POINTER) + H(J)
ELSE
IERROR = 1
ENDIF
ENDIF
ENDDO
ENDDO
9999 RETURN
END

@ -5,8 +5,7 @@ include ../../../Make.inc
#
FOBJS = dcsrck.o dcsrmm.o dcsrsm.o dcsrmv.o dcsrsv.o dcrnrmi.o \
dcrcrupd.o dcocrupd.o dcsrprt.o dcsrmv4.o dcsrmv2.o dcsrmv3.o \
dcsrrws.o\
dcsrprt.o dcsrmv4.o dcsrmv2.o dcsrmv3.o dcsrrws.o\
zcrnrmi.o zcsrmm.o zcsrrws.o zcsrsm.o zsrmv.o zsrsv.o zcsrck.o
OBJS=$(FOBJS)

@ -1,80 +0,0 @@
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
SUBROUTINE DCOCRUPD(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
C
C .. Matrix A to be updated is required to be stored with
C .. column indices belonging to the same row ordered.
C .. Block H to be inserted don't need to be stored in such a way.
C
C Flag = 0: put elements to 0.0D0;
C Flag = 1: replace elements with new value;
C Flag = 2: sum block value to elements;
C
IMPLICIT NONE
include 'psb_const.fh'
C .. Scalar Arguments ..
INTEGER IA, JA, IH, JH, M, N,
+ IERROR, FLAG, LIWORK
C .. Array Arguments ..
INTEGER IA1(*),IA2(*),IH1(*),IH2(*),
+ INFOA(*),INFOH(*),IWORK(*),
+ GLOB_TO_LOC(*)
CHARACTER DESCRA*11,DESCRH*11
DOUBLE PRECISION A(*),H(*)
C .. Local scalars ..
INTEGER J, NNZ, IP1, NNZI
C .. Local arrays ..
IERROR = 0
IF (IBITS(INFOA(PSB_UPD_),2,1).EQ.1) THEN
C
C Smart update capability
C
IP1 = INFOA(PSB_UPD_PNT_)
NNZ = IA2(IP1+PSB_NNZ_)
NNZI = INFOH(1)
DO J = 1, NNZI
NNZ = NNZ + 1
A(NNZ) = H(J)
ENDDO
IA2(IP1+PSB_NNZ_) = NNZ
ELSE
IERROR = 2
ENDIF
9999 CONTINUE
RETURN
END

@ -1,181 +0,0 @@
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
SUBROUTINE DCRCRUPD(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
C
C .. Matrix A to be updated is required to be stored with
C .. column indices belonging to the same row ordered.
C .. Block H to be inserted don't need to be stored in such a way.
C
C Flag = 0: put elements to 0.0D0;
C Flag = 1: replace elements with new value;
C Flag = 2: sum block value to elements;
C
IMPLICIT NONE
include 'psb_const.fh'
C .. Scalar Arguments ..
INTEGER IA, JA, IH, JH, M, N,
+ IERROR, FLAG, LIWORK
C .. Array Arguments ..
INTEGER IA1(*),IA2(*),IH1(*),IH2(*),
+ INFOA(*),INFOH(*),IWORK(*),
+ GLOB_TO_LOC(*)
CHARACTER DESCRA*11,DESCRH*11
DOUBLE PRECISION A(*),H(*)
C .. Local scalars ..
INTEGER I, J, XBLCK, XMATR,
+ NRC, IPH, JPH, JPA, LPA, IRET, LNK, NNZ, IP1
C .. Local arrays ..
IERROR = 0
IF (IBITS(INFOA(PSB_UPD_),2,1).EQ.1) THEN
C
C Smart update capability
C
IP1 = INFOA(PSB_UPD_PNT_)
NNZ = IA2(IP1+PSB_NNZ_)
DO I = 1, M
XBLCK = IH + I - 1
DO J = IH2(XBLCK),IH2(XBLCK+1) - 1
NNZ = NNZ + 1
A(NNZ) = H(J)
ENDDO
ENDDO
IA2(IP1+PSB_NNZ_) = NNZ
ELSE
IF (FLAG.EQ.0) THEN
DO I = 1, M
XBLCK = IH + I - 1
XMATR = IA + I - 1
NRC = IH2(XBLCK+1) - IH2(XBLCK)
IPH = IH2(XBLCK)
IF (LIWORK.LT.2*NRC+2) THEN
IERROR = 10
RETURN
ENDIF
DO J = 1, NRC
IWORK(J) = GLOB_TO_LOC(JA - JH + IH1(IPH+J-1))
ENDDO
CALL MRGSRT(NRC,IWORK(1),IWORK(NRC+1),IRET)
JPA = IA2(XMATR)
LPA = IA2(XMATR+1)
LNK = IWORK(NRC+1)
DO J = 1, NRC
JPH = IWORK(LNK)
DO WHILE ((IA1(JPA).NE.JPH).AND.(JPA.LT.LPA))
JPA = JPA + 1
ENDDO
IF (IA1(JPA).EQ.JPH) THEN
A(JPA) = 0.0D0
LNK = IWORK(NRC+1+LNK)
ELSE
IERROR = 1
GOTO 9999
ENDIF
enddo
ENDDO
ELSE IF (FLAG.EQ.1) THEN
DO I = 1, M
XBLCK = IH + I - 1
XMATR = IA + I - 1
NRC = IH2(XBLCK+1) - IH2(XBLCK)
IPH = IH2(XBLCK)
IF (LIWORK.LT.2*NRC+2) THEN
IERROR = 10
RETURN
ENDIF
DO J = 1, NRC
IWORK(J) = GLOB_TO_LOC(JA - JH + IH1(IPH+J-1))
ENDDO
CALL MRGSRT(NRC,IWORK(1),IWORK(NRC+1),IRET)
JPA = IA2(XMATR)
LPA = IA2(XMATR+1)
LNK = IWORK(NRC+1)
DO J = 1, NRC
JPH = IWORK(LNK)
DO WHILE((IA1(JPA).NE.JPH).AND.(JPA.LT.LPA))
JPA = JPA + 1
ENDDO
IF (IA1(JPA).EQ.JPH) THEN
A(JPA) = H(IPH+LNK-1)
LNK = IWORK(NRC+1+LNK)
ELSE
IERROR = 1
GOTO 9999
ENDIF
ENDDO
enddo
ELSE IF (FLAG.EQ.2) THEN
DO I = 1, M
XBLCK = IH + I - 1
XMATR = IA + I - 1
NRC = IH2(XBLCK+1) - IH2(XBLCK)
IPH = IH2(XBLCK)
IF (LIWORK.LT.2*NRC+2) THEN
IERROR = 10
RETURN
ENDIF
DO J = 1, NRC
IWORK(J) = GLOB_TO_LOC(JA - JH + IH1(IPH+J-1))
ENDDO
CALL MRGSRT(NRC,IWORK(1),IWORK(NRC+1),IRET)
JPA = IA2(XMATR)
LPA = IA2(XMATR+1)
LNK = IWORK(NRC+1)
DO J = 1, NRC
JPH = IWORK(LNK)
DO WHILE((IA1(JPA).NE.JPH).AND.(JPA.LT.LPA))
JPA = JPA + 1
ENDDO
IF (IA1(JPA).EQ.JPH) THEN
A(JPA) = A(JPA) + H(IPH+LNK-1)
LNK = IWORK(NRC+1+LNK)
ELSE
IERROR = 1
GOTO 9999
ENDIF
ENDDO
enddo
ELSE
IERROR = 2
ENDIF
ENDIF
9999 CONTINUE
RETURN
END

@ -1,85 +0,0 @@
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
SUBROUTINE DCRZERO(M,N,DESCRA,A,IA1,IA2,
+ INFOA,IA,JA,MZ,NZ,IERROR)
C
C This subroutione performs the operation:
C
C A(IA : IA + MZ - 1, JA : JA + NZ - 1) = 0
C
C This isn't accomplished by removing elements
C from sparse matrix representation, but assigning them
C the zero value.
C Columns are supposed to be ordered
C into the same row. This subroutine will
C not work properly if this hypotesis is not
C verified.
C
IMPLICIT NONE
C .. Scalar Arguments ..
INTEGER M,N,IA,JA,MZ,NZ,IERROR
C .. Array Arguments ..
INTEGER IA1(*),IA2(*),INFOA(*)
CHARACTER DESCRA*11
DOUBLE PRECISION A(*)
C .. Local scalars ..
INTEGER I, J, JBEGIN, JEND, AUX
DOUBLE PRECISION
IERROR=0
IF (((JA + NZ - 1) .GT. N) .OR.
+ ((IA + MZ - 1) .GT. M) .OR.
+ (IA .LT. 1) .OR. (JA .LT. 1)) THEN
IERROR = 1
GOTO 9999
ENDIF
DO I = IA, IA + M - 1
C Scan current line until found first element
DO JBEGIN = IA2(I), IA2(I + 1) - 1
IF (IA1(JBEGIN) .GE. JA) EXIT
ENDDO
C If reached last column end not yet
C encountered proper column, skip this row
IF ((JBEGIN .EQ. IA2(I + 1) - 1) .AND.
+ (IA1(JBEGIN) .LT. JA)) CYCLE
C Now I'm sure there's at least one element
C to process: scan until found last element
AUX = JA + N - 1
DO JEND = JBEGIN, IA2(I + 1) - 1
IF (IA1(JEND) .GE. AUX) EXIT
ENDDO
IF (IA1(JEND) .GT. AUX) JEND = JEND - 1
DO J = JBEGIN, JEND
A(J) = 0.0D0
ENDDO
ENDDO
9999 RETURN
END

@ -4,7 +4,7 @@ include ../../../Make.inc
# The object files
#
FOBJS = daxpby.o dcsmm.o dcsnmi.o dcsrp.o dcssm.o \
dcsupd.o dgelp.o dlpupd.o dswmm.o dswprt.o \
dgelp.o dlpupd.o dswmm.o dswprt.o \
dswsm.o smmp.o dcsrws.o \
zcsnmi.o zaxpby.o zcsmm.o zcssm.o zswmm.o zswsm.o\
zcsrws.o zgelp.o zlpupd.o

@ -1,149 +0,0 @@
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
SUBROUTINE DCSUPD(M, N, FIDA, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, FIDH, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
IMPLICIT NONE
C .. Scalar Arguments ..
INTEGER IA, JA, IH, JH, M, N,
+ IERROR, FLAG, LIWORK
C .. Array Arguments ..
INTEGER IA1(*),IA2(*),IH1(*),IH2(*),
+ INFOA(*),INFOH(*),IWORK(*),
+ GLOB_TO_LOC(*)
CHARACTER DESCRA*11,DESCRH*11, FIDA*5, FIDH*5
DOUBLE PRECISION A(*),H(*)
C .. Local Scalars..
INTEGER ERR_ACT
C .. Local Array..
integer int_val(5)
double precision real_val(5)
character*20 name, strings(2)
C .. Executable Statements ..
C
C
C Check parameters
C
IERROR = 0
NAME = 'DCSUPD\0'
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
IF (M.LT.0) THEN
IERROR = 10
INT_VAL(1) = 1
INT_VAL(2) = M
ELSE IF (N.LT.0) THEN
IERROR = 10
INT_VAL(1) = 2
INT_VAL(2) = N
ENDIF
C
C Error handling
C
IF(IERROR.NE.0) THEN
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ENDIF
C
C Check for M, N
C
IF(M.LE.0 .OR. N.LE.0) THEN
GOTO 9999
ENDIF
C
C Switching on FIDA
C
IF (FIDA(1:3).EQ.'CSR') THEN
IF (FIDH(1:3).EQ.'CSR') THEN
C
C Submatrix H in CSR format into A matrix in CSR format
C
CALL DCRCRUPD(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
ELSE IF (FIDH(1:3).EQ.'COO') THEN
C
C Submatrix H in COO format into A matrix in CSR format
C
CALL DCOCRUPD(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
ELSE
write(*,*) 'FIDA', FIDA(1:3), 'FIDH', FIDH(1:3)
IERROR = 3010
STRINGS(1) = FIDH(1:3)
NAME = 'DCSUPD\0'
ENDIF
ELSE IF (FIDA(1:3).EQ.'COO') THEN
IF (FIDH(1:3).EQ.'COO')THEN
CALL DCOCRUPD(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
ENDIF
ELSE IF (FIDA(1:3).EQ.'JAD') THEN
IF (FIDH(1:3).EQ.'COO')THEN
CALL DCOJDUPD(M, N, DESCRA, A, IA1,
+ IA2, INFOA, IA, JA, DESCRH, H, IH1, IH2,
+ INFOH, IH, JH, FLAG, GLOB_TO_LOC,
+ IWORK, LIWORK, IERROR)
ENDIF
C IERROR = 3010
C STRINGS(1) = FIDH(1:3)
C NAME = 'DCSUPD\0'
ENDIF
IF(IERROR.NE.0) THEN
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
END IF
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
RETURN
9999 CONTINUE
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
IF ( ERR_ACT .NE. 0 ) THEN
CALL FCPSB_SERROR()
RETURN
ENDIF
RETURN
END

@ -38,7 +38,8 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
use psb_realloc_mod
use psb_string_mod
use psb_error_mod
use psb_serial_mod, only : psb_spinfo
use psb_serial_mod, only : psb_spinfo, psb_csdp
use psb_update_mod
implicit none
integer, intent(in) :: nz, imin,imax,jmin,jmax
@ -49,9 +50,11 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
character(len=5) :: ufida
integer :: i,j,ir,ic,nr,nc, ng, nza, isza,spstate, nnz,&
& ip1, nzl, err_act, int_err(5), iupd
& ip1, nzl, err_act, int_err(5), iupd, irst
logical, parameter :: debug=.false.
logical :: rebuild_
character(len=20) :: name, ch_err
type(psb_dspmat_type) :: tmp
name='psb_dcoins'
info = 0
@ -84,8 +87,12 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
goto 9999
end if
!!$ if (present(rebuild)) then
!!$ rebuild_ = rebuild
!!$ else
rebuild_ = .false.
!!$ end if
!!$ ufida = toupper(a%fida)
call touppers(a%fida,ufida)
ng = size(gtl)
spstate = psb_sp_getifld(psb_state_,a,info)
@ -101,7 +108,7 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
end if
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info.ne.izero) then
if(info /= izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
@ -110,7 +117,7 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info.ne.izero) then
if(info /= izero) then
info=4010
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
@ -119,7 +126,7 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,&
& imin,imax,jmin,jmax,info)
if(info.ne.izero) then
if(info /= izero) then
info=4010
ch_err='psb_inner_ins'
call psb_errpush(info,name,a_err=ch_err)
@ -150,7 +157,7 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
if(info.ne.izero) then
if (info /= izero) then
info=4010
ch_err='psb_inner_upd'
call psb_errpush(info,name,a_err=ch_err)
@ -167,22 +174,51 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
case (psb_upd_dflt_, psb_upd_srch_)
select case(toupper(a%fida))
case ('CSR')
!!$ write(0,*) 'Calling csr_inner_upd'
call csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
!!$ write(0,*) 'From csr_inner_upd:',info
case ('COO')
call coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
call psb_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
if (info > 0) then
if (rebuild_) then
write(0,*) 'COINS: Going through rebuild_ fingers crossed!'
irst = info
call psb_nullify_sp(tmp)
tmp%fida='COO'
call psb_csdp(a,tmp,info)
call psb_setifield(psb_spmat_bld_,psb_state_,tmp,info)
call psb_sp_transfer(tmp,a,info)
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info /= izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
case default
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info /= izero) then
info=4010
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
endif
if (irst <= nz) then
call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),&
& nza,a%ia1,a%ia2,a%aspk,gtl,ng,imin,imax,jmin,jmax,info)
end if
else
info = 2231
call psb_errpush(info,name)
goto 9999
end if
else if (info < 0) then
info = 2230
call psb_errpush(info,name)
goto 9999
end select
end if
case default
@ -278,214 +314,5 @@ contains
end subroutine psb_inner_ins
subroutine csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
real(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc=i2-i1
if (.true.) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = -1
return
end if
else
!!$
ip = -1
lb = i1
ub = i2-1
do
if (lb > ub) exit
m = (lb+ub)/2
!!$ write(0,*) 'Debug: ',lb,m,ub
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = -1
return
end if
end if
end if
end do
case(psb_dupl_add_)
! Add
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = -2
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine csr_inner_upd
subroutine coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
real(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl,isrt
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then
info = -4
return
end if
ilr = -1
ilc = -1
nnz = psb_sp_getifld(psb_nnz_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
info = -2
return
end if
end if
end do
case(psb_dupl_add_)
! Add
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = -2
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine coo_inner_upd
end subroutine psb_dcoins

@ -0,0 +1,538 @@
!!$
!!$ Parallel Sparse BLAS v2.0
!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
module psb_update_mod
use psb_spmat_type
use psb_const_mod
use psb_realloc_mod
use psb_string_mod
interface psb_srch_upd
module procedure psb_d_srch_upd, psb_z_srch_upd
end interface
interface coo_srch_upd
module procedure d_coo_srch_upd, z_coo_srch_upd
end interface
interface csr_srch_upd
module procedure d_csr_srch_upd, z_csr_srch_upd
end interface
contains
subroutine psb_d_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
real(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
info = 0
select case(toupper(a%fida))
case ('CSR')
!!$ write(0,*) 'Calling csr_srch_upd'
call csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
!!$ write(0,*) 'From csr_srch_upd:',info
case ('COO')
call coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
case default
info = -9
end select
end subroutine psb_d_srch_upd
subroutine psb_z_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
complex(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
info = 0
select case(toupper(a%fida))
case ('CSR')
!!$ write(0,*) 'Calling csr_srch_upd'
call csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
!!$ write(0,*) 'From csr_srch_upd:',info
case ('COO')
call coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
case default
info = -9
end select
end subroutine psb_z_srch_upd
subroutine d_csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
real(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc=i2-i1
if (.true.) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = i
return
end if
else
!!$
ip = -1
lb = i1
ub = i2-1
do
if (lb > ub) exit
m = (lb+ub)/2
!!$ write(0,*) 'Debug: ',lb,m,ub
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = i
return
end if
end if
end if
end do
case(psb_dupl_add_)
! Add
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = i
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine d_csr_srch_upd
subroutine d_coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
real(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl,isrt
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then
info = -4
return
end if
ilr = -1
ilc = -1
nnz = psb_sp_getifld(psb_nnz_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
info = i
return
end if
end if
end do
case(psb_dupl_add_)
! Add
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = i
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine d_coo_srch_upd
subroutine z_csr_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
complex(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc=i2-i1
if (.true.) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = i
return
end if
else
!!$
ip = -1
lb = i1
ub = i2-1
do
if (lb > ub) exit
m = (lb+ub)/2
!!$ write(0,*) 'Debug: ',lb,m,ub
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = i
return
end if
end if
end if
end do
case(psb_dupl_add_)
! Add
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = i
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine z_csr_srch_upd
subroutine z_coo_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
complex(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl,isrt
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then
info = -4
return
end if
ilr = -1
ilc = -1
nnz = psb_sp_getifld(psb_nnz_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
info = i
return
end if
end if
end do
case(psb_dupl_add_)
! Add
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = i
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine z_coo_srch_upd
end module psb_update_mod

@ -38,7 +38,8 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
use psb_realloc_mod
use psb_string_mod
use psb_error_mod
use psb_serial_mod, only : psb_spinfo
use psb_serial_mod, only : psb_spinfo, psb_csdp
use psb_update_mod
implicit none
integer, intent(in) :: nz, imin,imax,jmin,jmax
@ -49,9 +50,11 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
character(len=5) :: ufida
integer :: i,j,ir,ic,nr,nc, ng, nza, isza,spstate, nnz,&
& ip1, nzl, err_act, int_err(5), iupd
& ip1, nzl, err_act, int_err(5), iupd, irst
logical, parameter :: debug=.false.
logical :: rebuild_
character(len=20) :: name, ch_err
type(psb_zspmat_type) :: tmp
name='psb_zcoins'
info = 0
@ -84,8 +87,12 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
goto 9999
end if
!!$ if (present(rebuild)) then
!!$ rebuild_ = rebuild
!!$ else
rebuild_ = .false.
!!$ end if
!!$ ufida = toupper(a%fida)
call touppers(a%fida,ufida)
ng = size(gtl)
spstate = psb_sp_getifld(psb_state_,a,info)
@ -101,7 +108,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
end if
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info.ne.izero) then
if(info /= izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
@ -110,7 +117,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info.ne.izero) then
if(info /= izero) then
info=4010
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
@ -119,7 +126,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,&
& imin,imax,jmin,jmax,info)
if(info.ne.izero) then
if(info /= izero) then
info=4010
ch_err='psb_inner_ins'
call psb_errpush(info,name,a_err=ch_err)
@ -150,7 +157,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
if(info.ne.izero) then
if (info /= izero) then
info=4010
ch_err='psb_inner_upd'
call psb_errpush(info,name,a_err=ch_err)
@ -167,22 +174,51 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
case (psb_upd_dflt_, psb_upd_srch_)
select case(toupper(a%fida))
case ('CSR')
!!$ write(0,*) 'Calling csr_inner_upd'
call csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
!!$ write(0,*) 'From csr_inner_upd:',info
case ('COO')
call coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
call psb_srch_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
if (info > 0) then
if (rebuild_) then
write(0,*) 'COINS: Going through rebuild_ fingers crossed!'
irst = info
call psb_nullify_sp(tmp)
tmp%fida='COO'
call psb_csdp(a,tmp,info)
call psb_setifield(psb_spmat_bld_,psb_state_,tmp,info)
call psb_sp_transfer(tmp,a,info)
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info /= izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
case default
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info /= izero) then
info=4010
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
endif
if (irst <= nz) then
call psb_inner_ins((nz-irst+1),ia(irst:nz),ja(irst:nz),val(irst:nz),&
& nza,a%ia1,a%ia2,a%aspk,gtl,ng,imin,imax,jmin,jmax,info)
end if
else
info = 2231
call psb_errpush(info,name)
goto 9999
end if
else if (info < 0) then
info = 2230
call psb_errpush(info,name)
goto 9999
end select
end if
case default
@ -278,214 +314,5 @@ contains
end subroutine psb_inner_ins
subroutine csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
complex(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc=i2-i1
if (.true.) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = -1
return
end if
else
!!$
ip = -1
lb = i1
ub = i2-1
do
if (lb > ub) exit
m = (lb+ub)/2
!!$ write(0,*) 'Debug: ',lb,m,ub
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1)
info = -1
return
end if
end if
end if
end do
case(psb_dupl_add_)
! Add
ilr = -1
ilc = -1
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
i1 = a%ia2(ir)
i2 = a%ia2(ir+1)
nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = -2
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine csr_inner_upd
subroutine coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,&
& imin,imax,jmin,jmax,nzl,info)
implicit none
type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng
integer, intent(in) :: ia(*),ja(*),gtl(*)
integer, intent(inout) :: nza
complex(kind(1.d0)), intent(in) :: val(*)
integer, intent(out) :: info
integer :: i,ir,ic,check_flag, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,nnz,dupl,isrt
info = 0
dupl = psb_sp_getifld(psb_dupl_,a,info)
if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then
info = -4
return
end if
ilr = -1
ilc = -1
nnz = psb_sp_getifld(psb_nnz_,a,info)
select case(dupl)
case(psb_dupl_ovwrt_,psb_dupl_err_)
! Overwrite.
! Cannot test for error, should have been caught earlier.
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
info = -2
return
end if
end if
end do
case(psb_dupl_add_)
! Add
do i=1, nz
ir = ia(i)
ic = ja(i)
if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
ir = gtl(ir)
ic = gtl(ic)
if (ir /= ilr) then
call ibsrch(i1,ir,nnz,a%ia1)
i2 = i1
do
if (i2+1 > nnz) exit
if (a%ia1(i2+1) /= a%ia1(i2)) exit
i2 = i2 + 1
end do
do
if (i1-1 < 1) exit
if (a%ia1(i1-1) /= a%ia1(i1)) exit
i1 = i1 - 1
end do
ilr = ir
end if
nc = i2-i1+1
call issrch(ip,ic,nc,a%ia2(i1:i2))
if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else
info = -2
return
end if
end if
end do
case default
info = -3
write(0,*) 'Duplicate handling: ',dupl
end select
end subroutine coo_inner_upd
end subroutine psb_zcoins

Loading…
Cancel
Save