psblas3:
. Make.inc.in base/comm/Makefile base/internals/Makefile base/internals/psi_exist_ovr_elem.f base/internals/psi_exist_ovr_elem.f90 base/internals/psi_list_search.f base/internals/psi_list_search.f90 base/internals/psi_srtlist.f90 base/internals/srtlist.f base/modules/Makefile base/serial/Makefile base/serial/f77 base/serial/f77/Makefile base/serial/f77/caxpby.f base/serial/f77/cnumbmm.f base/serial/f77/daxpby.f base/serial/f77/dnumbmm.f base/serial/f77/iaxpby.f base/serial/f77/saxpby.f base/serial/f77/snumbmm.f base/serial/f77/symbmm.f base/serial/f77/zaxpby.f base/serial/f77/znumbmm.f base/serial/psb_aspxpby.f90 base/serial/psi_c_serial_impl.f90 base/serial/psi_d_serial_impl.f90 base/serial/psi_i_serial_impl.f90 base/serial/psi_s_serial_impl.f90 base/serial/psi_z_serial_impl.f90 base/serial/smmp.f90 base/tools/Makefile configure.ac configure prec/psb_d_bjacprec.f90 prec/psb_d_diagprec.f90 prec/psb_d_nullprec.f90 test/fileread/Makefile test/hello/Makefile test/kernel/Makefile test/pargen/Makefile test/serial/Makefile test/torture/Makefile test/util/Makefile Merged changes from fixprec branch.trunk
parent
5473b6e783
commit
077998c1bd
@ -1,74 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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
|
||||
INTEGER FUNCTION PSI_EXIST_OVR_ELEM(OVR_ELEM,
|
||||
+ DIM_LIST,ELEM_SEARCHED)
|
||||
use psb_const_mod
|
||||
C PURPOSE:
|
||||
C == = ====
|
||||
C
|
||||
C If ELEM_SEARCHED exist in the list OVR_ELEM returns its position in
|
||||
C the list, else returns -1
|
||||
C
|
||||
C
|
||||
C INPUT
|
||||
C == = ===
|
||||
C OVRLAP_ELEMENT_D.: Contains for all overlap points belonging to
|
||||
C the current process:
|
||||
C 1. overlap point index
|
||||
C 2. Number of domains sharing that overlap point
|
||||
C the end is marked by a -1...............................
|
||||
C
|
||||
C DIM_LIST..........: Dimension of list OVRLAP_ELEMENT_D
|
||||
C
|
||||
C ELEM_SEARCHED.....:point's Local index identifier to be searched.
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
C ....Scalars parameters....
|
||||
INTEGER(psb_ipk_) :: DIM_LIST,ELEM_SEARCHED
|
||||
C ...Array Parameters....
|
||||
INTEGER(psb_ipk_) :: OVR_ELEM(DIM_LIST,*)
|
||||
|
||||
C ...Local Scalars....
|
||||
INTEGER(psb_ipk_) :: I
|
||||
|
||||
I=1
|
||||
DO WHILE ((I.LE.DIM_LIST).AND.(OVR_ELEM(I,1).NE.ELEM_SEARCHED))
|
||||
I=I+1
|
||||
ENDDO
|
||||
IF ((I.LE.DIM_LIST).AND.(OVR_ELEM(I,1).EQ.ELEM_SEARCHED)) THEN
|
||||
PSI_EXIST_OVR_ELEM=I
|
||||
ELSE
|
||||
PSI_EXIST_OVR_ELEM=-1
|
||||
ENDIF
|
||||
END
|
||||
|
@ -0,0 +1,73 @@
|
||||
!
|
||||
! Parallel Sparse BLAS version 3.5
|
||||
! (C) Copyright 2006, 2010, 2015, 2017
|
||||
! Salvatore Filippone Cranfield University
|
||||
! Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
integer function psi_exist_ovr_elem(ovr_elem, dim_list,elem_searched)
|
||||
use psb_const_mod
|
||||
! PURPOSE:
|
||||
! == = ====
|
||||
!
|
||||
! If ELEM_SEARCHED exist in the list OVR_ELEM returns its position in
|
||||
! the list, else returns -1
|
||||
!
|
||||
!
|
||||
! INPUT
|
||||
! == = ===
|
||||
! OVRLAP_ELEMENT_D.: Contains for all overlap points belonging to
|
||||
! the current process:
|
||||
! 1. overlap point index
|
||||
! 2. Number of domains sharing that overlap point
|
||||
! the end is marked by a -1...............................
|
||||
!
|
||||
! DIM_LIST..........: Dimension of list OVRLAP_ELEMENT_D
|
||||
!
|
||||
! ELEM_SEARCHED.....:point's Local index identifier to be searched.
|
||||
|
||||
implicit none
|
||||
|
||||
! ....Scalars parameters....
|
||||
integer(psb_ipk_) :: dim_list,elem_searched
|
||||
! ...array parameters....
|
||||
integer(psb_ipk_) :: ovr_elem(dim_list,*)
|
||||
|
||||
! ...local scalars....
|
||||
integer(psb_ipk_) :: i
|
||||
|
||||
i=1
|
||||
do while ((i.le.dim_list).and.(ovr_elem(i,1).ne.elem_searched))
|
||||
i=i+1
|
||||
enddo
|
||||
if ((i.le.dim_list).and.(ovr_elem(i,1).eq.elem_searched)) then
|
||||
psi_exist_ovr_elem=i
|
||||
else
|
||||
psi_exist_ovr_elem=-1
|
||||
endif
|
||||
end function psi_exist_ovr_elem
|
||||
|
@ -1,57 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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
|
||||
INTEGER FUNCTION PSI_LIST_SEARCH(LIST,LENGHT_LIST,ELEM)
|
||||
use psb_const_mod
|
||||
C !RETURNS POSITION OF ELEM IN A ARRAY LIST
|
||||
C !OF LENGHT LENGHT_LIST, IF THIS ELEMENT NOT EXISTS
|
||||
C !RETURNS -1
|
||||
INTEGER(psb_ipk_) :: LIST(*)
|
||||
INTEGER(psb_ipk_) :: LENGHT_LIST
|
||||
INTEGER(psb_ipk_) :: ELEM
|
||||
|
||||
INTEGER(psb_ipk_) :: I
|
||||
|
||||
I=1
|
||||
DO WHILE ((I.LE.LENGHT_LIST).AND.(LIST(I).NE.ELEM))
|
||||
I=I+1
|
||||
ENDDO
|
||||
IF (I.LE.LENGHT_LIST) THEN
|
||||
IF (LIST(I).EQ.ELEM) THEN
|
||||
PSI_LIST_SEARCH=I
|
||||
ELSE
|
||||
PSI_LIST_SEARCH=-1
|
||||
ENDIF
|
||||
ELSE
|
||||
PSI_LIST_SEARCH=-1
|
||||
ENDIF
|
||||
END
|
||||
|
@ -0,0 +1,58 @@
|
||||
!
|
||||
! Parallel Sparse BLAS version 3.5
|
||||
! (C) Copyright 2006, 2010, 2015, 2017
|
||||
! Salvatore Filippone Cranfield University
|
||||
! Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
integer function psi_list_search(list,lenght_list,elem)
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
!returns position of elem in a array list
|
||||
!of lenght lenght_list, if this element does not exist
|
||||
!returns -1
|
||||
integer(psb_ipk_) :: list(*)
|
||||
integer(psb_ipk_) :: lenght_list
|
||||
integer(psb_ipk_) :: elem
|
||||
|
||||
integer(psb_ipk_) :: i
|
||||
|
||||
i=1
|
||||
do while ((i.le.lenght_list).and.(list(i).ne.elem))
|
||||
i=i+1
|
||||
enddo
|
||||
if (i.le.lenght_list) then
|
||||
if (list(i).eq.elem) then
|
||||
psi_list_search=i
|
||||
else
|
||||
psi_list_search=-1
|
||||
endif
|
||||
else
|
||||
psi_list_search=-1
|
||||
endif
|
||||
end function psi_list_search
|
||||
|
@ -0,0 +1,203 @@
|
||||
!
|
||||
! Parallel Sparse BLAS version 3.5
|
||||
! (C) Copyright 2006, 2010, 2015, 2017
|
||||
! Salvatore Filippone Cranfield University
|
||||
! Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
!**********************************************************************
|
||||
! *
|
||||
! The communication step among processors at each *
|
||||
! matrix-vector product is a variable all-to-all *
|
||||
! collective communication that we reimplement *
|
||||
! in terms of point-to-point communications. *
|
||||
! The data in input is a list of dependencies: *
|
||||
! for each node a list of all the nodes it has to *
|
||||
! communicate with. The lists are guaranteed to be *
|
||||
! symmetric, i.e. for each pair (I,J) there is a *
|
||||
! pair (J,I). The idea is to organize the ordering *
|
||||
! so that at each communication step as many *
|
||||
! processors as possible are communicating at the *
|
||||
! same time, i.e. a step is defined by the fact *
|
||||
! that all edges (I,J) in it have no common node. *
|
||||
! *
|
||||
! Formulation of the problem is: *
|
||||
! Given an undirected graph (forest): *
|
||||
! Find the shortest series of steps to cancel all *
|
||||
! graph edges, where at each step all edges belonging *
|
||||
! to a matching in the graph are canceled. *
|
||||
! *
|
||||
! An obvious lower bound to the optimum number of steps *
|
||||
! is the largest degree of any node in the graph. *
|
||||
! *
|
||||
! The algorithm proceeds as follows: *
|
||||
! 1. Build a list of all edges, e.g. copy the *
|
||||
! dependencies lists keeping only (I,J) with I<J *
|
||||
! 2. Compute an auxiliary vector with the degree of *
|
||||
! each node of the graph. *
|
||||
! 3. While there are edges in the graph do: *
|
||||
! 4. Weight the edges with the sum of the degrees *
|
||||
! of their nodes and sort them into descending order *
|
||||
! 5. Scan the list of edges; if neither node of the *
|
||||
! edge has been marked yet, cancel the edge and mark *
|
||||
! the two nodes *
|
||||
! 6. If no edge was chosen but the graph is nonempty *
|
||||
! raise an error condition *
|
||||
! 7. Queue the edges in the matchin to the output *
|
||||
! sequence; *
|
||||
! 8. Decrease by 1 the degree of all marked nodes, *
|
||||
! then clear all marks *
|
||||
! 9. Cycle to 3. *
|
||||
! 10. For each node: scan the edge sequence; if an *
|
||||
! edge has the node as an endpoint, queue the other *
|
||||
! node in the dependency list for the current one *
|
||||
! *
|
||||
!**********************************************************************
|
||||
subroutine srtlist(dep_list,dl_lda,ldl,np,dg,dgp,upd, edges,idx,ich,info)
|
||||
use psb_serial_mod
|
||||
implicit none
|
||||
integer(psb_ipk_) :: np, dl_lda, info
|
||||
integer(psb_ipk_) :: dep_list(dl_lda,*), ldl(*),dg(*), dgp(*),&
|
||||
& idx(*), upd(*),edges(2,*),ich(*)
|
||||
integer(psb_ipk_) :: i,j, nedges,ip1,ip2,nch,ip,iedge,&
|
||||
& i1,ix,ist,iswap(2)
|
||||
integer(psb_ipk_) :: no_comm
|
||||
parameter (no_comm=-1)
|
||||
|
||||
|
||||
if (np .lt. 0) then
|
||||
info = 1
|
||||
return
|
||||
endif
|
||||
|
||||
!
|
||||
! dg contains number of communications
|
||||
!
|
||||
do i=1, np
|
||||
dg(i)=ldl(i)
|
||||
enddo
|
||||
|
||||
|
||||
nedges = 0
|
||||
do i=1, np
|
||||
do j=1, dg(i)
|
||||
ip = dep_list(j,i) + 1
|
||||
if (ip.gt.i) nedges = nedges + 1
|
||||
enddo
|
||||
enddo
|
||||
|
||||
iedge = 0
|
||||
do i=1, np
|
||||
do j=1, dg(i)
|
||||
ip = dep_list(j,i) + 1
|
||||
if (ip.gt.i) then
|
||||
iedge = iedge + 1
|
||||
edges(1,iedge) = i
|
||||
edges(2,iedge) = ip
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
ist = 1
|
||||
do while (ist.le.nedges)
|
||||
|
||||
do i=1, np
|
||||
upd(i) = 0
|
||||
enddo
|
||||
do i=ist, nedges
|
||||
dgp(i) = -(dg(edges(1,i))+dg(edges(2,i)))
|
||||
enddo
|
||||
|
||||
call psb_msort(dgp(ist:nedges),ix=idx(ist:nedges))
|
||||
i1 = ist
|
||||
nch = 0
|
||||
do i = ist, nedges
|
||||
ix = idx(i)+ist-1
|
||||
ip1 = edges(1,ix)
|
||||
ip2 = edges(2,ix)
|
||||
if ((upd(ip1).eq.0).and.(upd(ip2).eq.0)) then
|
||||
upd(ip1) = -1
|
||||
upd(ip2) = -1
|
||||
nch = nch + 1
|
||||
ich(nch) = ix
|
||||
endif
|
||||
enddo
|
||||
if (nch.eq.0) then
|
||||
write(psb_err_unit,*)&
|
||||
& 'srtlist ?????? impossible error !!!!!?????',&
|
||||
& nedges,ist
|
||||
do i=ist, nedges
|
||||
ix = idx(i)+ist-1
|
||||
write(psb_err_unit,*)&
|
||||
& 'SRTLIST: Edge:',ix,edges(1,ix),&
|
||||
& edges(2,ix),dgp(ix)
|
||||
enddo
|
||||
info = psb_err_input_value_invalid_i_
|
||||
return
|
||||
endif
|
||||
call psb_msort(ich(1:nch))
|
||||
do i=1, nch
|
||||
iswap(1) = edges(1,ist)
|
||||
iswap(2) = edges(2,ist)
|
||||
edges(1,ist) = edges(1,ich(i))
|
||||
edges(2,ist) = edges(2,ich(i))
|
||||
edges(1,ich(i)) = iswap(1)
|
||||
edges(2,ich(i)) = iswap(2)
|
||||
ist = ist + 1
|
||||
enddo
|
||||
do i=1, np
|
||||
dg(i) = dg(i) + upd(i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i=1, np
|
||||
if (dg(i).ne.0) then
|
||||
write(psb_err_unit,*)&
|
||||
& 'SRTLIST Error on exit:',i,dg(i)
|
||||
endif
|
||||
dg(i) = 0
|
||||
enddo
|
||||
do j=1,nedges
|
||||
i = edges(1,j)
|
||||
dg(i) = dg(i)+1
|
||||
dep_list(dg(i),i) = edges(2,j)-1
|
||||
i = edges(2,j)
|
||||
dg(i) = dg(i)+1
|
||||
dep_list(dg(i),i) = edges(1,j)-1
|
||||
enddo
|
||||
do i=1, np
|
||||
if (dg(i).ne.ldl(i)) then
|
||||
write(psb_err_unit,*) &
|
||||
& 'SRTLIST Mismatch on output',i,dg(i),ldl(i)
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
return
|
||||
end subroutine srtlist
|
||||
|
||||
|
@ -1,215 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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
|
||||
***********************************************************************
|
||||
* *
|
||||
* The communication step among processors at each *
|
||||
* matrix-vector product is a variable all-to-all *
|
||||
* collective communication that we reimplement *
|
||||
* in terms of point-to-point communications. *
|
||||
* The data in input is a list of dependencies: *
|
||||
* for each node a list of all the nodes it has to *
|
||||
* communicate with. The lists are guaranteed to be *
|
||||
* symmetric, i.e. for each pair (I,J) there is a *
|
||||
* pair (J,I). The idea is to organize the ordering *
|
||||
* so that at each communication step as many *
|
||||
* processors as possible are communicating at the *
|
||||
* same time, i.e. a step is defined by the fact *
|
||||
* that all edges (I,J) in it have no common node. *
|
||||
* *
|
||||
* Formulation of the problem is: *
|
||||
* Given an undirected graph (forest): *
|
||||
* Find the shortest series of steps to cancel all *
|
||||
* graph edges, where at each step all edges belonging *
|
||||
* to a matching in the graph are canceled. *
|
||||
* *
|
||||
* An obvious lower bound to the optimum number of steps *
|
||||
* is the largest degree of any node in the graph. *
|
||||
* *
|
||||
* The algorithm proceeds as follows: *
|
||||
* 1. Build a list of all edges, e.g. copy the *
|
||||
* dependencies lists keeping only (I,J) with I<J *
|
||||
* 2. Compute an auxiliary vector with the degree of *
|
||||
* each node of the graph. *
|
||||
* 3. While there are edges in the graph do: *
|
||||
* 4. Weight the edges with the sum of the degrees *
|
||||
* of their nodes and sort them into descending order *
|
||||
* 5. Scan the list of edges; if neither node of the *
|
||||
* edge has been marked yet, cancel the edge and mark *
|
||||
* the two nodes *
|
||||
* 6. If no edge was chosen but the graph is nonempty *
|
||||
* raise an error condition *
|
||||
* 7. Queue the edges in the matchin to the output *
|
||||
* sequence; *
|
||||
* 8. Decrease by 1 the degree of all marked nodes, *
|
||||
* then clear all marks *
|
||||
* 9. Cycle to 3. *
|
||||
* 10. For each node: scan the edge sequence; if an *
|
||||
* edge has the node as an endpoint, queue the other *
|
||||
* node in the dependency list for the current one *
|
||||
* *
|
||||
***********************************************************************
|
||||
SUBROUTINE SRTLIST(DEP_LIST,DL_LDA,LDL,NP,dg,dgp,upd,
|
||||
+ edges,idx,ich,INFO)
|
||||
use psb_serial_mod
|
||||
IMPLICIT NONE
|
||||
INTEGER(psb_ipk_) :: NP, DL_LDA, INFO
|
||||
INTEGER(psb_ipk_) :: DEP_LIST(DL_LDA,*), LDL(*),DG(*), DGP(*),
|
||||
+ IDX(*), UPD(*),EDGES(2,*),ICH(*)
|
||||
INTEGER(psb_ipk_) :: I,J, NEDGES,IP1,IP2,NCH,IP,IEDGE,
|
||||
+ I1,IX,IST,ISWAP(2)
|
||||
INTEGER(psb_ipk_) :: NO_COMM
|
||||
PARAMETER (NO_COMM=-1)
|
||||
|
||||
|
||||
IF (NP .LT. 0) THEN
|
||||
INFO = 1
|
||||
RETURN
|
||||
ENDIF
|
||||
|
||||
C
|
||||
C dg contains number of communications
|
||||
C
|
||||
DO I=1, NP
|
||||
DG(I)=LDL(I)
|
||||
ENDDO
|
||||
|
||||
|
||||
NEDGES = 0
|
||||
DO I=1, NP
|
||||
DO J=1, DG(I)
|
||||
IP = DEP_LIST(J,I) + 1
|
||||
c$$$ write(psb_err_unit,*)
|
||||
c$$$ 'SRTLIST Input :',i,ip
|
||||
IF (IP.GT.I)
|
||||
+ NEDGES = NEDGES + 1
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
IEDGE = 0
|
||||
DO I=1, NP
|
||||
DO J=1, DG(I)
|
||||
IP = DEP_LIST(J,I) + 1
|
||||
IF (IP.GT.I) THEN
|
||||
IEDGE = IEDGE + 1
|
||||
EDGES(1,IEDGE) = I
|
||||
EDGES(2,IEDGE) = IP
|
||||
ENDIF
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
IST = 1
|
||||
DO WHILE (IST.LE.NEDGES)
|
||||
|
||||
DO I=1, NP
|
||||
UPD(I) = 0
|
||||
ENDDO
|
||||
DO I=IST, NEDGES
|
||||
DGP(I) = -(DG(EDGES(1,I))+DG(EDGES(2,I)))
|
||||
ENDDO
|
||||
|
||||
call psb_msort(dgp(ist:nedges),ix=idx(ist:nedges))
|
||||
I1 = IST
|
||||
NCH = 0
|
||||
DO I = IST, NEDGES
|
||||
IX = IDX(I)+IST-1
|
||||
IP1 = EDGES(1,IX)
|
||||
IP2 = EDGES(2,IX)
|
||||
IF ((UPD(IP1).eq.0).AND.(UPD(IP2).eq.0)) THEN
|
||||
UPD(IP1) = -1
|
||||
UPD(IP2) = -1
|
||||
NCH = NCH + 1
|
||||
ICH(NCH) = IX
|
||||
ENDIF
|
||||
ENDDO
|
||||
IF (NCH.eq.0) THEN
|
||||
write(psb_err_unit,*)
|
||||
+ 'SRTLIST ?????? Impossible error !!!!!?????',
|
||||
+ nedges,ist
|
||||
do i=ist, nedges
|
||||
IX = IDX(I)+IST-1
|
||||
write(psb_err_unit,*)
|
||||
+ 'SRTLIST: Edge:',ix,edges(1,ix),
|
||||
+ edges(2,ix),dgp(ix)
|
||||
enddo
|
||||
info = psb_err_input_value_invalid_i_
|
||||
return
|
||||
ENDIF
|
||||
call psb_msort(ich(1:nch))
|
||||
DO I=1, NCH
|
||||
ISWAP(1) = EDGES(1,IST)
|
||||
ISWAP(2) = EDGES(2,IST)
|
||||
EDGES(1,IST) = EDGES(1,ICH(I))
|
||||
EDGES(2,IST) = EDGES(2,ICH(I))
|
||||
EDGES(1,ICH(I)) = ISWAP(1)
|
||||
EDGES(2,ICH(I)) = ISWAP(2)
|
||||
IST = IST + 1
|
||||
ENDDO
|
||||
DO I=1, NP
|
||||
DG(I) = DG(I) + UPD(I)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
DO I=1, NP
|
||||
IF (DG(I).NE.0) THEN
|
||||
write(psb_err_unit,*)
|
||||
+ 'SRTLIST Error on exit:',i,dg(i)
|
||||
ENDIF
|
||||
DG(I) = 0
|
||||
ENDDO
|
||||
DO J=1,NEDGES
|
||||
I = EDGES(1,J)
|
||||
DG(I) = DG(I)+1
|
||||
DEP_LIST(DG(I),I) = EDGES(2,J)-1
|
||||
I = EDGES(2,J)
|
||||
DG(I) = DG(I)+1
|
||||
DEP_LIST(DG(I),I) = EDGES(1,J)-1
|
||||
ENDDO
|
||||
DO I=1, NP
|
||||
IF (DG(I).NE.LDL(I)) THEN
|
||||
write(psb_err_unit,*)
|
||||
+ 'SRTLIST Mismatch on output',i,dg(i),ldl(i)
|
||||
ENDIF
|
||||
ENDDO
|
||||
|
||||
c$$$ write(psb_err_unit,*)
|
||||
c$$$ 'Output communication:',t2-t1
|
||||
c$$$ do i=1,np
|
||||
c$$$ do j=1,ldl(i)
|
||||
c$$$ write(psb_err_unit,*)
|
||||
c$$$ 'SRTLIST', i,dep_list(j,i)+1
|
||||
c$$$ enddo
|
||||
c$$$ enddo
|
||||
|
||||
RETURN
|
||||
END
|
||||
|
||||
|
@ -1,39 +0,0 @@
|
||||
include ../../../Make.inc
|
||||
|
||||
#
|
||||
# The object files
|
||||
#
|
||||
FOBJS = iaxpby.o daxpby.o saxpby.o \
|
||||
caxpby.o zaxpby.o symbmm.o \
|
||||
cnumbmm.o dnumbmm.o snumbmm.o znumbmm.o
|
||||
|
||||
|
||||
OBJS=$(FOBJS)
|
||||
|
||||
#
|
||||
# Where the library should go, and how it is called.
|
||||
# Note that we are regenerating most of libsparker.a on the fly.
|
||||
LIBDIR=../..
|
||||
INCDIR=../..
|
||||
MODDIR=../../modules
|
||||
FINCLUDES=$(FMFLAG). $(FMFLAG)$(MODDIR) $(FMFLAG)$(INCDIR)
|
||||
LIBFILE=$(LIBDIR)/$(LIBNAME)
|
||||
|
||||
#
|
||||
# No change should be needed below
|
||||
#
|
||||
|
||||
|
||||
default: lib
|
||||
|
||||
lib: $(OBJS)
|
||||
$(AR) $(LIBDIR)/$(LIBNAME) $(OBJS)
|
||||
$(RANLIB) $(LIBDIR)/$(LIBNAME)
|
||||
|
||||
clean: cleanobjs
|
||||
|
||||
veryclean: cleanobjs
|
||||
|
||||
cleanobjs:
|
||||
/bin/rm -f $(OBJS)
|
||||
|
@ -1,200 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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 caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
|
||||
use psb_const_mod
|
||||
complex(psb_spk_), parameter :: one=(1.0,0.0)
|
||||
complex(psb_spk_), parameter :: zero=(0.0,0.0)
|
||||
integer(psb_ipk_) :: n, m, lldx, lldy, info
|
||||
complex(psb_spk_) X(lldx,*), Y(lldy,*)
|
||||
complex(psb_spk_) alpha, beta
|
||||
integer(psb_ipk_) :: i, j
|
||||
integer(psb_ipk_) :: int_err(5)
|
||||
character name*20
|
||||
name='caxpby'
|
||||
|
||||
|
||||
C
|
||||
C Error handling
|
||||
C
|
||||
info = psb_success_
|
||||
if (m.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (n.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=n
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldx.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=5
|
||||
int_err(2)=1
|
||||
int_err(3)=lldx
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldy.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=8
|
||||
int_err(2)=1
|
||||
int_err(3)=lldy
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (alpha.eq.zero) then
|
||||
if (beta.eq.zero) then
|
||||
do j=1, n
|
||||
do i=1,m
|
||||
y(i,j) = zero
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
c$$$
|
||||
c$$$ Do nothing!
|
||||
c$$$
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.one) then
|
||||
|
||||
if (beta.eq.zero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.-one) then
|
||||
|
||||
if (beta.eq.zero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if (beta.eq.zero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call fcpsb_serror()
|
||||
return
|
||||
|
||||
end
|
@ -1,85 +0,0 @@
|
||||
c == =====================================================================
|
||||
c Sparse Matrix Multiplication Package
|
||||
c
|
||||
c Randolph E. Bank and Craig C. Douglas
|
||||
c
|
||||
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
c
|
||||
c Compile this with the following command (or a similar one):
|
||||
c
|
||||
c f77 -c -O smmp.f
|
||||
c
|
||||
c == =====================================================================
|
||||
subroutine cnumbmm(n, m, l,
|
||||
* ia, ja, diaga, a,
|
||||
* ib, jb, diagb, b,
|
||||
* ic, jc, diagc, c,
|
||||
* temp)
|
||||
c
|
||||
use psb_const_mod
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga,
|
||||
* ib(*), jb(*), diagb,
|
||||
* ic(*), jc(*), diagc
|
||||
c
|
||||
complex(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
c
|
||||
c numeric matrix multiply c=a*b
|
||||
c
|
||||
maxlmn = max(l,m,n)
|
||||
do 10 i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
10 continue
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
c
|
||||
c c = a*b
|
||||
c
|
||||
do 50 i = 1,n
|
||||
do 30 jj = ia(i),ia(i+1)
|
||||
c a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) goto 30
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
c b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm)
|
||||
* temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do 20 k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
20 continue
|
||||
30 continue
|
||||
c c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
c$$$ if (mod(i,100) == 1)
|
||||
c$$$ + write(psb_err_unit,*)
|
||||
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do 40 j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
40 continue
|
||||
50 continue
|
||||
return
|
||||
end
|
@ -1,198 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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 daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
|
||||
use psb_const_mod
|
||||
integer(psb_ipk_) :: n, m, lldx, lldy, info
|
||||
real(psb_dpk_) X(lldx,*), Y(lldy,*)
|
||||
real(psb_dpk_) alpha, beta
|
||||
integer(psb_ipk_) :: i, j
|
||||
integer(psb_ipk_) :: int_err(5)
|
||||
character name*20
|
||||
name='daxpby'
|
||||
|
||||
|
||||
C
|
||||
C Error handling
|
||||
C
|
||||
info = psb_success_
|
||||
if (m.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (n.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=n
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldx.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=5
|
||||
int_err(2)=1
|
||||
int_err(3)=lldx
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldy.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=8
|
||||
int_err(2)=1
|
||||
int_err(3)=lldy
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (alpha.eq.dzero) then
|
||||
if (beta.eq.dzero) then
|
||||
do j=1, n
|
||||
do i=1,m
|
||||
y(i,j) = dzero
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.done) then
|
||||
c$$$
|
||||
c$$$ Do nothing!
|
||||
c$$$
|
||||
|
||||
else if (beta.eq.-done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.done) then
|
||||
|
||||
if (beta.eq.dzero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.-done) then
|
||||
|
||||
if (beta.eq.dzero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if (beta.eq.dzero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-done) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call fcpsb_serror()
|
||||
return
|
||||
|
||||
end
|
@ -1,85 +0,0 @@
|
||||
c == =====================================================================
|
||||
c Sparse Matrix Multiplication Package
|
||||
c
|
||||
c Randolph E. Bank and Craig C. Douglas
|
||||
c
|
||||
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
c
|
||||
c Compile this with the following command (or a similar one):
|
||||
c
|
||||
c f77 -c -O smmp.f
|
||||
c
|
||||
c == =====================================================================
|
||||
subroutine dnumbmm(n, m, l,
|
||||
* ia, ja, diaga, a,
|
||||
* ib, jb, diagb, b,
|
||||
* ic, jc, diagc, c,
|
||||
* temp)
|
||||
use psb_const_mod
|
||||
c
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga,
|
||||
* ib(*), jb(*), diagb,
|
||||
* ic(*), jc(*), diagc
|
||||
c
|
||||
real(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
c
|
||||
c numeric matrix multiply c=a*b
|
||||
c
|
||||
maxlmn = max(l,m,n)
|
||||
do 10 i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
10 continue
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
c
|
||||
c c = a*b
|
||||
c
|
||||
do 50 i = 1,n
|
||||
do 30 jj = ia(i),ia(i+1)
|
||||
c a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) goto 30
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
c b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm)
|
||||
* temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do 20 k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
20 continue
|
||||
30 continue
|
||||
c c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
c$$$ if (mod(i,100) == 1)
|
||||
c$$$ + write(psb_err_unit,*)
|
||||
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do 40 j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
40 continue
|
||||
50 continue
|
||||
return
|
||||
end
|
@ -1,198 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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 iaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
|
||||
use psb_const_mod
|
||||
integer n, m, lldx, lldy, info
|
||||
integer(psb_ipk_) X(lldx,*), Y(lldy,*)
|
||||
integer(psb_ipk_) alpha, beta
|
||||
integer(psb_ipk_) :: i, j
|
||||
integer(psb_ipk_) :: int_err(5)
|
||||
character name*20
|
||||
name='saxpby'
|
||||
|
||||
|
||||
C
|
||||
C Error handling
|
||||
C
|
||||
info = psb_success_
|
||||
if (m.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (n.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=n
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldx.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=5
|
||||
int_err(2)=1
|
||||
int_err(3)=lldx
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldy.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=8
|
||||
int_err(2)=1
|
||||
int_err(3)=lldy
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (alpha.eq.izero) then
|
||||
if (beta.eq.izero) then
|
||||
do j=1, n
|
||||
do i=1,m
|
||||
y(i,j) = izero
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.ione) then
|
||||
c$$$
|
||||
c$$$ Do nothing!
|
||||
c$$$
|
||||
|
||||
else if (beta.eq.-ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.ione) then
|
||||
|
||||
if (beta.eq.izero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.-ione) then
|
||||
|
||||
if (beta.eq.izero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if (beta.eq.izero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-ione) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call fcpsb_serror()
|
||||
return
|
||||
|
||||
end
|
@ -1,198 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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 saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
|
||||
use psb_const_mod
|
||||
integer n, m, lldx, lldy, info
|
||||
real(psb_spk_) X(lldx,*), Y(lldy,*)
|
||||
real(psb_spk_) alpha, beta
|
||||
integer(psb_ipk_) :: i, j
|
||||
integer(psb_ipk_) :: int_err(5)
|
||||
character name*20
|
||||
name='saxpby'
|
||||
|
||||
|
||||
C
|
||||
C Error handling
|
||||
C
|
||||
info = psb_success_
|
||||
if (m.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (n.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=n
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldx.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=5
|
||||
int_err(2)=1
|
||||
int_err(3)=lldx
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldy.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=8
|
||||
int_err(2)=1
|
||||
int_err(3)=lldy
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (alpha.eq.szero) then
|
||||
if (beta.eq.szero) then
|
||||
do j=1, n
|
||||
do i=1,m
|
||||
y(i,j) = szero
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.sone) then
|
||||
c$$$
|
||||
c$$$ Do nothing!
|
||||
c$$$
|
||||
|
||||
else if (beta.eq.-sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.sone) then
|
||||
|
||||
if (beta.eq.szero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.-sone) then
|
||||
|
||||
if (beta.eq.szero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if (beta.eq.szero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-sone) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call fcpsb_serror()
|
||||
return
|
||||
|
||||
end
|
@ -1,85 +0,0 @@
|
||||
c == =====================================================================
|
||||
c Sparse Matrix Multiplication Package
|
||||
c
|
||||
c Randolph E. Bank and Craig C. Douglas
|
||||
c
|
||||
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
c
|
||||
c Compile this with the following command (or a similar one):
|
||||
c
|
||||
c f77 -c -O smmp.f
|
||||
c
|
||||
c == =====================================================================
|
||||
subroutine snumbmm(n, m, l,
|
||||
* ia, ja, diaga, a,
|
||||
* ib, jb, diagb, b,
|
||||
* ic, jc, diagc, c,
|
||||
* temp)
|
||||
use psb_const_mod
|
||||
c
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga,
|
||||
* ib(*), jb(*), diagb,
|
||||
* ic(*), jc(*), diagc
|
||||
c
|
||||
real(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
c
|
||||
c numeric matrix multiply c=a*b
|
||||
c
|
||||
maxlmn = max(l,m,n)
|
||||
do 10 i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
10 continue
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
c
|
||||
c c = a*b
|
||||
c
|
||||
do 50 i = 1,n
|
||||
do 30 jj = ia(i),ia(i+1)
|
||||
c a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) goto 30
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
c b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm)
|
||||
* temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do 20 k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
20 continue
|
||||
30 continue
|
||||
c c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
c$$$ if (mod(i,100) == 1)
|
||||
c$$$ + write(psb_err_unit,*)
|
||||
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do 40 j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
40 continue
|
||||
50 continue
|
||||
return
|
||||
end
|
@ -1,119 +0,0 @@
|
||||
c == =====================================================================
|
||||
c Sparse Matrix Multiplication Package
|
||||
c
|
||||
c Randolph E. Bank and Craig C. Douglas
|
||||
c
|
||||
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
c
|
||||
c Compile this with the following command (or a similar one):
|
||||
c
|
||||
c f77 -c -O smmp.f
|
||||
c
|
||||
c == =====================================================================
|
||||
subroutine symbmm
|
||||
* (n, m, l,
|
||||
* ia, ja, diaga,
|
||||
* ib, jb, diagb,
|
||||
* ic, jc, diagc,
|
||||
* index)
|
||||
use psb_const_mod
|
||||
use psb_realloc_mod
|
||||
use psb_sort_mod, only: psb_msort
|
||||
c
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga,
|
||||
* ib(*), jb(*), diagb, diagc, index(*)
|
||||
integer(psb_ipk_), allocatable :: ic(:),jc(:)
|
||||
integer(psb_ipk_) :: nze, info
|
||||
|
||||
c
|
||||
c symbolic matrix multiply c=a*b
|
||||
c
|
||||
if (size(ic) < n+1) then
|
||||
write(psb_err_unit,*)
|
||||
+ 'Called realloc in SYMBMM '
|
||||
call psb_realloc(n+1,ic,info)
|
||||
if (info /= psb_success_) then
|
||||
write(psb_err_unit,*)
|
||||
+ 'realloc failed in SYMBMM ',info
|
||||
end if
|
||||
endif
|
||||
maxlmn = max(l,m,n)
|
||||
do 10 i=1,maxlmn
|
||||
index(i)=0
|
||||
10 continue
|
||||
if (diagc.eq.0) then
|
||||
ic(1)=1
|
||||
else
|
||||
ic(1)=n+2
|
||||
endif
|
||||
minlm = min(l,m)
|
||||
minmn = min(m,n)
|
||||
c
|
||||
c main loop
|
||||
c
|
||||
do 50 i=1,n
|
||||
istart=-1
|
||||
length=0
|
||||
c
|
||||
c merge row lists
|
||||
c
|
||||
do 30 jj=ia(i),ia(i+1)
|
||||
c a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) goto 30
|
||||
j = i
|
||||
else
|
||||
j=ja(jj)
|
||||
endif
|
||||
c b = d + ...
|
||||
if (index(j).eq.0 .and. diagb.eq.1 .and. j.le.minlm)then
|
||||
index(j)=istart
|
||||
istart=j
|
||||
length=length+1
|
||||
endif
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' SymbMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
do 20 k=ib(j),ib(j+1)-1
|
||||
if ((jb(k)<1).or.(jb(k)>maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ 'Problem in SYMBMM 1:',j,k,jb(k),maxlmn
|
||||
else
|
||||
if(index(jb(k)).eq.0) then
|
||||
index(jb(k))=istart
|
||||
istart=jb(k)
|
||||
length=length+1
|
||||
endif
|
||||
endif
|
||||
20 continue
|
||||
30 continue
|
||||
c
|
||||
c row i of jc
|
||||
c
|
||||
if (diagc.eq.1 .and. index(i).ne.0) length = length - 1
|
||||
ic(i+1)=ic(i)+length
|
||||
|
||||
if (ic(i+1) > size(jc)) then
|
||||
if (n > (2*i)) then
|
||||
nze = max(ic(i+1), ic(i)*((n+i-1)/i))
|
||||
else
|
||||
nze = max(ic(i+1), nint((dble(ic(i))*(dble(n)/i))) )
|
||||
endif
|
||||
call psb_realloc(nze,jc,info)
|
||||
end if
|
||||
|
||||
do 40 j= ic(i),ic(i+1)-1
|
||||
if (diagc.eq.1 .and. istart.eq.i) then
|
||||
istart = index(istart)
|
||||
index(i) = 0
|
||||
endif
|
||||
jc(j)=istart
|
||||
istart=index(istart)
|
||||
index(jc(j))=0
|
||||
40 continue
|
||||
call psb_msort(jc(ic(i):ic(i)+length -1))
|
||||
index(i) = 0
|
||||
50 continue
|
||||
return
|
||||
end
|
@ -1,200 +0,0 @@
|
||||
C
|
||||
C Parallel Sparse BLAS version 3.5
|
||||
C (C) Copyright 2006, 2010, 2015, 2017
|
||||
C Salvatore Filippone Cranfield University
|
||||
C Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
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 zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info)
|
||||
use psb_const_mod
|
||||
complex(psb_dpk_), parameter :: one=(1.0d0,0.0d0)
|
||||
complex(psb_dpk_), parameter :: zero=(0.0d0,0.0d0)
|
||||
integer(psb_ipk_) :: n, m, lldx, lldy, info
|
||||
complex(psb_dpk_) X(lldx,*), Y(lldy,*)
|
||||
complex(psb_dpk_) alpha, beta
|
||||
integer(psb_ipk_) :: i, j
|
||||
integer(psb_ipk_) :: int_err(5)
|
||||
character name*20
|
||||
name='zaxpby'
|
||||
|
||||
|
||||
C
|
||||
C Error handling
|
||||
C
|
||||
info = psb_success_
|
||||
if (m.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (n.lt.0) then
|
||||
info=psb_err_iarg_neg_
|
||||
int_err(1)=1
|
||||
int_err(2)=n
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldx.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=5
|
||||
int_err(2)=1
|
||||
int_err(3)=lldx
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
else if (lldy.lt.max(1,m)) then
|
||||
info=psb_err_iarg_not_gtia_ii_
|
||||
int_err(1)=8
|
||||
int_err(2)=1
|
||||
int_err(3)=lldy
|
||||
int_err(4)=m
|
||||
call fcpsb_errpush(info,name,int_err)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (alpha.eq.zero) then
|
||||
if (beta.eq.zero) then
|
||||
do j=1, n
|
||||
do i=1,m
|
||||
y(i,j) = zero
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
c$$$
|
||||
c$$$ Do nothing!
|
||||
c$$$
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.one) then
|
||||
|
||||
if (beta.eq.zero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else if (alpha.eq.-one) then
|
||||
|
||||
if (beta.eq.zero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = -x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
if (beta.eq.zero) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else if (beta.eq.one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if (beta.eq.-one) then
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) - y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
y(i,j) = alpha*x(i,j) + beta*y(i,j)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
|
||||
9999 continue
|
||||
call fcpsb_serror()
|
||||
return
|
||||
|
||||
end
|
@ -1,85 +0,0 @@
|
||||
c == =====================================================================
|
||||
c Sparse Matrix Multiplication Package
|
||||
c
|
||||
c Randolph E. Bank and Craig C. Douglas
|
||||
c
|
||||
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
c
|
||||
c Compile this with the following command (or a similar one):
|
||||
c
|
||||
c f77 -c -O smmp.f
|
||||
c
|
||||
c == =====================================================================
|
||||
subroutine znumbmm(n, m, l,
|
||||
* ia, ja, diaga, a,
|
||||
* ib, jb, diagb, b,
|
||||
* ic, jc, diagc, c,
|
||||
* temp)
|
||||
c
|
||||
use psb_const_mod
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga,
|
||||
* ib(*), jb(*), diagb,
|
||||
* ic(*), jc(*), diagc
|
||||
c
|
||||
complex(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
c
|
||||
c numeric matrix multiply c=a*b
|
||||
c
|
||||
maxlmn = max(l,m,n)
|
||||
do 10 i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
10 continue
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
c
|
||||
c c = a*b
|
||||
c
|
||||
do 50 i = 1,n
|
||||
do 30 jj = ia(i),ia(i+1)
|
||||
c a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) goto 30
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
c b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm)
|
||||
* temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do 20 k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
20 continue
|
||||
30 continue
|
||||
c c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
c$$$ if (mod(i,100) == 1)
|
||||
c$$$ + write(psb_err_unit,*)
|
||||
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do 40 j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)
|
||||
+ ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
40 continue
|
||||
50 continue
|
||||
return
|
||||
end
|
@ -0,0 +1,477 @@
|
||||
!
|
||||
! Parallel Sparse BLAS version 3.5
|
||||
! (C) Copyright 2006, 2010, 2015, 2017
|
||||
! Salvatore Filippone Cranfield University
|
||||
! Alfredo Buttari CNRS-IRIT, Toulouse
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
! Original code adapted from:
|
||||
! == =====================================================================
|
||||
! Sparse Matrix Multiplication Package
|
||||
!
|
||||
! Randolph E. Bank and Craig C. Douglas
|
||||
!
|
||||
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
!
|
||||
! Compile this with the following command (or a similar one):
|
||||
!
|
||||
! f77 -c -O smmp.f
|
||||
!
|
||||
! == =====================================================================
|
||||
subroutine symbmm(n, m, l, ia, ja, diaga, ib, jb, diagb,&
|
||||
& ic, jc, diagc, index)
|
||||
use psb_const_mod
|
||||
use psb_realloc_mod
|
||||
use psb_sort_mod, only: psb_msort
|
||||
!
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga, &
|
||||
& ib(*), jb(*), diagb, diagc, index(*)
|
||||
integer(psb_ipk_), allocatable :: ic(:),jc(:)
|
||||
integer(psb_ipk_) :: nze, info
|
||||
|
||||
!
|
||||
! symbolic matrix multiply c=a*b
|
||||
!
|
||||
if (size(ic) < n+1) then
|
||||
write(psb_err_unit,*)&
|
||||
& 'Called realloc in SYMBMM '
|
||||
call psb_realloc(n+1,ic,info)
|
||||
if (info /= psb_success_) then
|
||||
write(psb_err_unit,*)&
|
||||
& 'realloc failed in SYMBMM ',info
|
||||
end if
|
||||
endif
|
||||
maxlmn = max(l,m,n)
|
||||
do i=1,maxlmn
|
||||
index(i)=0
|
||||
end do
|
||||
if (diagc.eq.0) then
|
||||
ic(1)=1
|
||||
else
|
||||
ic(1)=n+2
|
||||
endif
|
||||
minlm = min(l,m)
|
||||
minmn = min(m,n)
|
||||
!
|
||||
! main loop
|
||||
!
|
||||
do i=1,n
|
||||
istart=-1
|
||||
length=0
|
||||
!
|
||||
! merge row lists
|
||||
!
|
||||
rowi: do jj=ia(i),ia(i+1)
|
||||
! a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
|
||||
j = i
|
||||
else
|
||||
j=ja(jj)
|
||||
endif
|
||||
! b = d + ...
|
||||
if (index(j).eq.0 .and. diagb.eq.1 .and. j.le.minlm)then
|
||||
index(j)=istart
|
||||
istart=j
|
||||
length=length+1
|
||||
endif
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' SymbMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
do k=ib(j),ib(j+1)-1
|
||||
if ((jb(k)<1).or.(jb(k)>maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& 'Problem in SYMBMM 1:',j,k,jb(k),maxlmn
|
||||
else
|
||||
if(index(jb(k)).eq.0) then
|
||||
index(jb(k))=istart
|
||||
istart=jb(k)
|
||||
length=length+1
|
||||
endif
|
||||
endif
|
||||
end do
|
||||
end do rowi
|
||||
|
||||
!
|
||||
! row i of jc
|
||||
!
|
||||
if (diagc.eq.1 .and. index(i).ne.0) length = length - 1
|
||||
ic(i+1)=ic(i)+length
|
||||
|
||||
if (ic(i+1) > size(jc)) then
|
||||
if (n > (2*i)) then
|
||||
nze = max(ic(i+1), ic(i)*((n+i-1)/i))
|
||||
else
|
||||
nze = max(ic(i+1), nint((dble(ic(i))*(dble(n)/i))) )
|
||||
endif
|
||||
call psb_realloc(nze,jc,info)
|
||||
end if
|
||||
|
||||
do j= ic(i),ic(i+1)-1
|
||||
if (diagc.eq.1 .and. istart.eq.i) then
|
||||
istart = index(istart)
|
||||
index(i) = 0
|
||||
endif
|
||||
jc(j)=istart
|
||||
istart=index(istart)
|
||||
index(jc(j))=0
|
||||
end do
|
||||
call psb_msort(jc(ic(i):ic(i)+length -1))
|
||||
index(i) = 0
|
||||
end do
|
||||
return
|
||||
end subroutine symbmm
|
||||
! == =====================================================================
|
||||
! Sparse Matrix Multiplication Package
|
||||
!
|
||||
! Randolph E. Bank and Craig C. Douglas
|
||||
!
|
||||
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
!
|
||||
! Compile this with the following command (or a similar one):
|
||||
!
|
||||
! f77 -c -O smmp.f
|
||||
!
|
||||
! == =====================================================================
|
||||
subroutine cnumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
|
||||
& ic, jc, diagc, c, temp)
|
||||
!
|
||||
use psb_const_mod
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga,&
|
||||
& ib(*), jb(*), diagb, ic(*), jc(*), diagc
|
||||
!
|
||||
complex(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
!
|
||||
! numeric matrix multiply c=a*b
|
||||
!
|
||||
maxlmn = max(l,m,n)
|
||||
do i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
end do
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
!
|
||||
! c = a*b
|
||||
!
|
||||
do i = 1,n
|
||||
rowi: do jj = ia(i),ia(i+1)
|
||||
! a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
! b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm) &
|
||||
& temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
end do
|
||||
end do rowi
|
||||
|
||||
! c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
!$$$ if (mod(i,100) == 1)
|
||||
!$$$ + write(psb_err_unit,*)
|
||||
!$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
end do
|
||||
end do
|
||||
|
||||
return
|
||||
end subroutine cnumbmm
|
||||
! == =====================================================================
|
||||
! Sparse Matrix Multiplication Package
|
||||
!
|
||||
! Randolph E. Bank and Craig C. Douglas
|
||||
!
|
||||
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
!
|
||||
! Compile this with the following command (or a similar one):
|
||||
!
|
||||
! f77 -c -O smmp.f
|
||||
!
|
||||
! == =====================================================================
|
||||
subroutine dnumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
|
||||
& ic, jc, diagc, c, temp)
|
||||
use psb_const_mod
|
||||
!
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
|
||||
& ic(*), jc(*), diagc
|
||||
!
|
||||
real(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
!
|
||||
! numeric matrix multiply c=a*b
|
||||
!
|
||||
maxlmn = max(l,m,n)
|
||||
do i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
end do
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
!
|
||||
! c = a*b
|
||||
!
|
||||
do i = 1,n
|
||||
rowi: do jj = ia(i),ia(i+1)
|
||||
! a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
! b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm) &
|
||||
& temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
end do
|
||||
end do rowi
|
||||
|
||||
! c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
!$$$ if (mod(i,100) == 1)
|
||||
!$$$ + write(psb_err_unit,*)
|
||||
!$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
end do
|
||||
end do
|
||||
|
||||
return
|
||||
end subroutine dnumbmm
|
||||
! == =====================================================================
|
||||
! Sparse Matrix Multiplication Package
|
||||
!
|
||||
! Randolph E. Bank and Craig C. Douglas
|
||||
!
|
||||
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
!
|
||||
! Compile this with the following command (or a similar one):
|
||||
!
|
||||
! f77 -c -O smmp.f
|
||||
!
|
||||
! == =====================================================================
|
||||
subroutine snumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
|
||||
& ic, jc, diagc, c, temp)
|
||||
use psb_const_mod
|
||||
!
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
|
||||
& ic(*), jc(*), diagc
|
||||
!
|
||||
real(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
!
|
||||
! numeric matrix multiply c=a*b
|
||||
!
|
||||
maxlmn = max(l,m,n)
|
||||
do i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
end do
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
!
|
||||
! c = a*b
|
||||
!
|
||||
do i = 1,n
|
||||
rowi: do jj = ia(i),ia(i+1)
|
||||
! a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
! b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm) &
|
||||
& temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
end do
|
||||
end do rowi
|
||||
|
||||
! c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
!$$$ if (mod(i,100) == 1)
|
||||
!$$$ + write(psb_err_unit,*)
|
||||
!$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
|
||||
do j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
end do
|
||||
end do
|
||||
|
||||
return
|
||||
end subroutine snumbmm
|
||||
! == =====================================================================
|
||||
! Sparse Matrix Multiplication Package
|
||||
!
|
||||
! Randolph E. Bank and Craig C. Douglas
|
||||
!
|
||||
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
|
||||
!
|
||||
! Compile this with the following command (or a similar one):
|
||||
!
|
||||
! f77 -c -O smmp.f
|
||||
!
|
||||
! == =====================================================================
|
||||
subroutine znumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
|
||||
& ic, jc, diagc, c, temp)
|
||||
!
|
||||
use psb_const_mod
|
||||
integer(psb_ipk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
|
||||
& ic(*), jc(*), diagc
|
||||
!
|
||||
complex(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
|
||||
!
|
||||
! numeric matrix multiply c=a*b
|
||||
!
|
||||
maxlmn = max(l,m,n)
|
||||
do i = 1,maxlmn
|
||||
temp(i) = 0.
|
||||
end do
|
||||
minlm = min(l,m)
|
||||
minln = min(l,n)
|
||||
minmn = min(m,n)
|
||||
!
|
||||
! c = a*b
|
||||
!
|
||||
do i = 1,n
|
||||
rowi: do jj = ia(i),ia(i+1)
|
||||
! a = d + ...
|
||||
if (jj.eq.ia(i+1)) then
|
||||
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
|
||||
j = i
|
||||
ajj = a(i)
|
||||
else
|
||||
j=ja(jj)
|
||||
ajj = a(jj)
|
||||
endif
|
||||
! b = d + ...
|
||||
if (diagb.eq.1 .and. j.le.minlm) &
|
||||
& temp(j) = temp(j) + ajj * b(j)
|
||||
if ((j<1).or.(j>m)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: Problem with A ',i,jj,j,m
|
||||
endif
|
||||
|
||||
do k = ib(j),ib(j+1)-1
|
||||
if((jb(k)<1).or. (jb(k) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
|
||||
else
|
||||
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
|
||||
endif
|
||||
end do
|
||||
end do rowi
|
||||
|
||||
! c = d + ...
|
||||
if (diagc.eq.1 .and. i.le.minln) then
|
||||
c(i) = temp(i)
|
||||
temp(i) = 0.
|
||||
endif
|
||||
do j = ic(i),ic(i+1)-1
|
||||
if((jc(j)<1).or. (jc(j) > maxlmn)) then
|
||||
write(psb_err_unit,*)&
|
||||
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
|
||||
else
|
||||
c(j) = temp(jc(j))
|
||||
temp(jc(j)) = 0.
|
||||
endif
|
||||
end do
|
||||
end do
|
||||
|
||||
return
|
||||
end subroutine znumbmm
|
Loading…
Reference in New Issue