diff --git a/src/serial/Makefile b/src/serial/Makefile index 3da00891..da7377d2 100644 --- a/src/serial/Makefile +++ b/src/serial/Makefile @@ -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) diff --git a/src/serial/coo/Makefile b/src/serial/coo/Makefile index 7c6eaae8..fe25ca68 100644 --- a/src/serial/coo/Makefile +++ b/src/serial/coo/Makefile @@ -3,7 +3,7 @@ include ../../../Make.inc # # The object files # -FOBJS = dcooprt.o dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o\ +FOBJS = dcooprt.o dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o\ zcoomm.o zcoomv.o zcoonrmi.o zcoorws.o zcoosm.o zcoosv.o diff --git a/src/serial/coo/dcoozero.f b/src/serial/coo/dcoozero.f deleted file mode 100644 index 85bfd74a..00000000 --- a/src/serial/coo/dcoozero.f +++ /dev/null @@ -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 diff --git a/src/serial/coo/dcrupdate.f b/src/serial/coo/dcrupdate.f deleted file mode 100644 index bf7e1b96..00000000 --- a/src/serial/coo/dcrupdate.f +++ /dev/null @@ -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 - - - - diff --git a/src/serial/csr/Makefile b/src/serial/csr/Makefile index bd0edd1b..f13c0ac7 100644 --- a/src/serial/csr/Makefile +++ b/src/serial/csr/Makefile @@ -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) diff --git a/src/serial/csr/dcocrupd.f b/src/serial/csr/dcocrupd.f deleted file mode 100644 index 8389ea38..00000000 --- a/src/serial/csr/dcocrupd.f +++ /dev/null @@ -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 - - - - diff --git a/src/serial/csr/dcrcrupd.f b/src/serial/csr/dcrcrupd.f deleted file mode 100644 index 6e5c7af3..00000000 --- a/src/serial/csr/dcrcrupd.f +++ /dev/null @@ -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 - - - - diff --git a/src/serial/csr/dcrzero.f b/src/serial/csr/dcrzero.f deleted file mode 100644 index 00fcb264..00000000 --- a/src/serial/csr/dcrzero.f +++ /dev/null @@ -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 diff --git a/src/serial/f77/Makefile b/src/serial/f77/Makefile index ac825df5..ba897019 100644 --- a/src/serial/f77/Makefile +++ b/src/serial/f77/Makefile @@ -3,8 +3,8 @@ 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 \ +FOBJS = daxpby.o dcsmm.o dcsnmi.o dcsrp.o dcssm.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 diff --git a/src/serial/f77/dcsupd.f b/src/serial/f77/dcsupd.f deleted file mode 100644 index c5e5ad76..00000000 --- a/src/serial/f77/dcsupd.f +++ /dev/null @@ -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 - - - diff --git a/src/serial/psb_dcoins.f90 b/src/serial/psb_dcoins.f90 index 2e1d855d..164537df 100644 --- a/src/serial/psb_dcoins.f90 +++ b/src/serial/psb_dcoins.f90 @@ -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,& + call psb_srch_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,& - & imin,imax,jmin,jmax,nzl,info) - - case default - + 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 + + 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 diff --git a/src/serial/psb_update_mod.f90 b/src/serial/psb_update_mod.f90 new file mode 100644 index 00000000..2201f609 --- /dev/null +++ b/src/serial/psb_update_mod.f90 @@ -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 + diff --git a/src/serial/psb_zcoins.f90 b/src/serial/psb_zcoins.f90 index d852673d..0f38cead 100644 --- a/src/serial/psb_zcoins.f90 +++ b/src/serial/psb_zcoins.f90 @@ -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,& + call psb_srch_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,& - & imin,imax,jmin,jmax,nzl,info) - - case default - + 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 + + 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