You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
psblas3/base/serial/dp/dvtfg.f

203 lines
14 KiB
Fortran

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
C
C Compute level numbers for the triangular matrix.
C
SUBROUTINE DVTFG (UPLO,M,JA,IA,NG,IPA,IPAT,KLEN,IWORK1,IWORK2,
* IWORK3)
use psb_string_mod
use psb_error_mod
implicit none
C .. Scalar Arguments ..
INTEGER M, NG
CHARACTER UPLO
C .. Array Arguments ..
INTEGER IA(*), IPA(*), IPAT(*), IWORK3(*), JA(*),
* KLEN(*), IWORK2(*), IWORK1(*)
C .. Local Scalars ..
INTEGER I, J, L, L0, L1, LEV, NP, iret
C .. Intrinsic Functions ..
INTRINSIC MAX
integer :: debug_level, debug_unit
character(len=20) :: name='DVTFG'
C .. Executable Statements ..
C
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
NG = 0
C
C CHECK ON THE NUMBERS OF THE ELEMENTS OF THE MATRIX
C
IF ((IA(M+1)-1).EQ.0) THEN
C
C THE MATRIX HASN'T ELEMENTS
C THE OUTPUT PERMUTATIONS ARE POSED TO THE IDENTITY MATRIX
C
DO 20 I = 1, M
IPA(I) = I
IPAT(I) = I
KLEN(I) = 0
20 CONTINUE
ELSE
C
C COMPUTE LEVEL NUMBER FOR EACH ROWS
C
C IWORK1: AUXILIARY VECTOR WHICH CONTAINS
C LEVEL NUMBERS OF EACH ROW
C
DO 40 I = 1, M
IWORK1(I) = 0
IWORK3(I) = 0
40 CONTINUE
IF (toupper(UPLO).EQ.'L') THEN
C
C LOWER TRIANGULAR SPARSE MATRIX
C
DO 80 I = 1, M
IWORK1(I) = 1
DO 60 J = IA(I), IA(I+1) - 1
IWORK1(I) = MAX(IWORK1(I),IWORK1(JA(J))+1)
60 CONTINUE
80 CONTINUE
ELSE IF (toupper(UPLO).EQ.'U') THEN
C
C UPPER TRIANGULAR SPARSE MATRIX
C
DO 120 I = M, 1, -1
IWORK1(I) = 1
DO 100 J = IA(I), IA(I+1) - 1
IWORK1(I) = MAX(IWORK1(I),IWORK1(JA(J))+1)
100 CONTINUE
120 CONTINUE
END IF
C
C COUNT NUMBER OF ROWS IN EACH EQUIVALENCE GROUPS
C
C NOTE: GROUP = SET OF ROWS HAVING THE SAME LEVEL NUMBER
C
C IWORK3: AUXILIARY VECTOR WHICH CONTAINS
C THE NUMBER OF ROWS FOR EACH GROUPS
C
DO 140 I = 1, M
IWORK3(IWORK1(I)) = IWORK3(IWORK1(I)) + 1
140 CONTINUE
C
C SET UP IWORK2:
C IWORK2(I) POINTS TO THE BEGINNING OF I-TH GROUP
C
IWORK2(1) = 1
DO 160 I = 2, M + 1
IWORK2(I) = IWORK3(I-1) + IWORK2(I-1)
IF (IWORK2(I).EQ.M+1) THEN
NG = I - 1
GO TO 180
END IF
160 CONTINUE
180 CONTINUE
C
C NG : TOTAL NUMBER OF LEVELS
C IWORK3: VECTOR CONTAINING THE NUMBER OF THE ROWS SORTED BY
C EQUIVALENCE GROUPS.
C
DO 200 I = 1, M
IWORK3(IWORK2(IWORK1(I))) = I
IWORK2(IWORK1(I)) = IWORK2(IWORK1(I)) + 1
200 CONTINUE
C
C REGENERATE IWORK2: POINTER INTO IWORK3 FOR EQUIVALENCE GROUPS
C
DO 220 I = NG + 1, 2, -1
IWORK2(I) = IWORK2(I-1)
220 CONTINUE
IWORK2(1) = 1
C
C IWORK1: ROWS LENGTH IN NEW ORDERING
C
DO 240 L = 1, M
IWORK1(L) = IA(IWORK3(L)) - IA(IWORK3(L)+1)
240 CONTINUE
C
C SORT ROWS BY DECREASING NUMBER OF NONZERO ELEMENTS.
C
C IPA: VECTOR OF NUMBER OF ROWS SORTED BY EQUIVALENCE
C GROUPS AND BY ROW LENGTH
C IPAT: AUXILIARY VECTOR NEED TO THE SORTER ROUTINES
C
L1 = IWORK2(2) - IWORK2(1)
DO 260 L = 1, L1
IPA(L) = IWORK3(L)
260 CONTINUE
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Group ',1,':',(ipa(l),l=1,l1)
DO 360 LEV = 2, NG
C
C LOOP ON GROUPS
C L1: LENGTH OF CURRENT GROUP
C L0: POINTER TO IPA TO THE FIRST LOCATIONS RESERVED
C FOR CURRENT GROUP
C
L1 = IWORK2(LEV+1) - IWORK2(LEV)
L0 = IWORK2(LEV) - 1
CALL MSORT_UP(L1,IWORK1(IWORK2(LEV)),IPAT,IRET)
IF (IRET.EQ.0) THEN
NP = IPAT(1)
DO 280 L = 1, L1
IPA(L0+L) = IWORK3(L0+NP)
NP = IPAT(1+NP)
280 CONTINUE
ELSE
C
C VECTOR ALREADY SORTED. NO CHANGE IS NEED
C
DO 320 L = 1, L1
IPA(L0+L) = IWORK3(L0+L)
320 CONTINUE
ENDIF
if (debug_level >= psb_debug_serial_)
+ write(debug_unit,*) trim(name),
+ ': Group ',lev,':',(ipa(l0+l),l=1,l1)
360 CONTINUE
C
C IPAT = IPA-1
C
DO 380 I = 1, M
IPAT(IPA(I)) = I
380 CONTINUE
C DO 400 I = 1, NG
C KLEN(I) = IWORK2(I+1) - IWORK2(I)
DO 400 I = 1, NG+1
KLEN(I) = IWORK2(I)
400 CONTINUE
END IF
RETURN
END