diff options
Diffstat (limited to 'Trivac/src/VECPER.f')
| -rwxr-xr-x | Trivac/src/VECPER.f | 204 |
1 files changed, 204 insertions, 0 deletions
diff --git a/Trivac/src/VECPER.f b/Trivac/src/VECPER.f new file mode 100755 index 0000000..1916e03 --- /dev/null +++ b/Trivac/src/VECPER.f @@ -0,0 +1,204 @@ +*DECK VECPER + SUBROUTINE VECPER(HNAME,IMPV,ISEG,L4,MUIN,LON,LTSW,NBL,LBL,MUV, + 1 IPV) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Ordering of matrix elements for supervectorial operations on a matrix +* in compressed diagonal storage mode. +* +*Copyright: +* Copyright (C) 2002 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): A. Hebert +* +*Parameters: input +* HNAME name of the matrix (for edition purpose only). +* IMPV print parameter for statistics (equal to zero for no print). +* ISEG number of elements in a vector register. +* L4 matrix order. +* MUIN position of each diagonal element in non-ordered matrix. +* +*Parameters: output +* LON number of groups of linear systems. +* LTSW maximum bandwidth (=2 for tridiagonal systems). +* NBL number of linear systems in each group. +* LBL number of unknowns in each group. +* MUV position of each diagonal element in ordered matrix. +* IPV permutation vector for the ordered unknowns. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + CHARACTER HNAME*4 + INTEGER IMPV,ISEG,L4,MUIN(L4),LON,LTSW,NBL(LON),LBL(LON),MUV(L4), + 1 IPV(L4) +*---- +* LOCAL VARIABLES +*---- + INTEGER, DIMENSION(:), ALLOCATABLE :: ISET,ISORT,IOFSET,IORD +*---- +* DETERMINE THE TOTAL NUMBER OF LINEAR SYSTEMS AND COMPUTE THE ORDER +* OF EACH OF THEM +*---- + ALLOCATE(ISET(L4)) + ISET(1)=0 + K1=MUIN(1)+1 + DO 10 I=2,L4 + ISET(I)=0 + K2=MUIN(I) + DO 5 J=I-K2+K1,I-1 + ISET(J)=1 + 5 CONTINUE + K1=K2+1 + 10 CONTINUE + NSYS=0 + DO 15 I=1,L4 + IPV(I)=0 + MUV(I)=0 + IF(ISET(I).EQ.0) NSYS=NSYS+1 + 15 CONTINUE + LON=1+(NSYS-1)/ISEG + IF(IMPV.GE.2) WRITE (6,'(/35H VECPER: NUMBER OF INDEPENDANT LINE, + 1 22HAR SYSTEMS IN MATRIX '',A4,3H'' =,I7/9X,17HNUMBER OF GROUPS , + 1 19HOF LINEAR SYSTEMS =,I6)') HNAME,NSYS,LON + ALLOCATE(IORD(NSYS),IOFSET(NSYS)) + ISYS=0 + IORD0=0 + IOFSET(1)=1 + DO 20 I=1,L4 + IF(ISET(I).EQ.0) THEN + ISYS=ISYS+1 + IORD(ISYS)=I-IORD0 + IF(I.NE.L4) IOFSET(ISYS+1)=I+1 + IORD0=I + ENDIF + 20 CONTINUE + DEALLOCATE(ISET) +*---- +* SORT THE LINEAR SYSTEMS BY DECREASING ORDER +*---- + ALLOCATE(ISORT(NSYS)) + JNEW=NSYS + DO 25 ISYS=NSYS,1,-1 + IF(IORD(ISYS).EQ.1) THEN + ISORT(JNEW)=ISYS + JNEW=JNEW-1 + ENDIF + 25 CONTINUE + INEW=0 + 30 IBIG=0 + DO 50 ISYS=1,NSYS + IF(IORD(ISYS).EQ.1) GO TO 50 + DO 40 KSYS=1,INEW + IF(ISORT(KSYS).EQ.ISYS) GO TO 50 + 40 CONTINUE + IBIG=MAX(IBIG,IORD(ISYS)) + 50 CONTINUE + IF(IBIG.LE.1) GO TO 70 + DO 60 ISYS=1,NSYS + IF(IORD(ISYS).EQ.IBIG) THEN + INEW=INEW+1 + ISORT(INEW)=ISYS + ENDIF + 60 CONTINUE + GO TO 30 + 70 IF(INEW.NE.JNEW) CALL XABORT('VECPER: ALGORITHM FAILURE 1') + DO 80 I=1,LON + ISYS=ISORT((I-1)*ISEG+1) + NBL(I)=ISEG + LBL(I)=IORD(ISYS) + 80 CONTINUE + NBL(LON)=NSYS-(LON-1)*ISEG + IF(IMPV.GE.2) WRITE (6,'(9X,33HMAXIMUM ORDER OF AN INDEPENDANT L, + 1 14HINEAR SYSTEM =,I9)') LBL(1) + IF(IMPV.GE.3) THEN + I1=1 + DO 90 I=1,(LON-1)/8+1 + I2=I1+7 + IF(I2.GT.LON) I2=LON + WRITE (6,200) (J,J=I1,I2) + WRITE (6,210) (NBL(J),J=I1,I2) + WRITE (6,220) (LBL(J),J=I1,I2) + I1=I1+8 + 90 CONTINUE + ENDIF +*---- +* COMPUTE THE PERMUTATION MATRIX +*---- + LBL0=0 + KSYS=0 + DO 105 J=1,LON + DO 101 K=1,NBL(J) + KSYS=KSYS+1 + ISYS=ISORT(KSYS) + IOF0=IOFSET(ISYS) + IOF1=IOF0+IORD(ISYS)-1 + IF(IOF1.GT.L4) CALL XABORT('VECPER: ALGORITHM FAILURE 2') + DO 100 I=IOF0,IOF1 + IPV(I)=(LBL0+I-IOF0)*ISEG+K + 100 CONTINUE + 101 CONTINUE + LBL0=LBL0+LBL(J) + 105 CONTINUE + DO 110 I=1,L4 + IF(IPV(I).LE.0) CALL XABORT('VECPER: ALGORITHM FAILURE 3') + IF(IPV(I).GT.LBL0*ISEG) CALL XABORT('VECPER: ALGORITHM FAILURE 4') + 110 CONTINUE + L4NEW=0 + DO 115 J=1,LON + L4NEW=L4NEW+LBL(J)*NBL(J) + 115 CONTINUE + IF(IMPV.GE.2) WRITE (6,'(/35H VECPER: INCREASING NUMBER OF UNKNO, + 1 8HWNS FROM,I7,3H TO,I7,11H. FILL-IN =,F7.2,3H %.)') L4,L4NEW, + 2 100.0*(REAL(L4NEW)/REAL(L4)-1.0) +*---- +* COMPUTE THE VECTORIAL BANDWIDTH +*---- + LBL0=0 + KSYS=0 + IIMAX=0 + LTSW=0 + MAXNEW=0 + MUVOLD=0 + DO 150 J=1,LON + DO 120 I=1,LBL(J) + MUV(LBL0+I)=1 + 120 CONTINUE + DO 131 K=1,NBL(J) + KSYS=KSYS+1 + ISYS=ISORT(KSYS) + IOF0=IOFSET(ISYS)-1 + DO 130 I=2,IORD(ISYS) + IBIG=MUIN(IOF0+I)-MUIN(IOF0+I-1) + IF(IBIG.GT.MUV(LBL0+I)) MUV(LBL0+I)=IBIG + 130 CONTINUE + 131 CONTINUE + DO 140 I=1,LBL(J) + LTSW=MAX(LTSW,MUV(LBL0+I)) + IIMAX=IIMAX+MUV(LBL0+I) + MUV(LBL0+I)=IIMAX + 140 CONTINUE + LBL0=LBL0+LBL(J) + MAXNEW=MAXNEW+(MUV(LBL0)-MUVOLD)*NBL(J) + MUVOLD=MUV(LBL0) + 150 CONTINUE + IF(IMPV.GE.2) WRITE (6,'(/35H VECPER: INCREASING NUMBER OF TERMS, + 1 17H IN MATRICES FROM,I9,3H TO,I9,11H. FILL-IN =,F7.2,3H %./9X, + 2 19HMAXIMUM BANDWIDTH =,I4)') MUIN(L4),MAXNEW, + 3 100.0*(REAL(MAXNEW)/REAL(MUIN(L4))-1.0),LTSW +* + DEALLOCATE(ISORT,IOFSET,IORD) + RETURN +* + 200 FORMAT (//13H GROUP ,8(I8,5X,1HI)) + 210 FORMAT ( 13H NB. SYSTEMS ,8(I8,5X,1HI)) + 220 FORMAT ( 13H NB. UNKNOWNS,8(I8,5X,1HI)) + END |