C E X A M P L E 6
C
C Program name:
CHARACTER*9 MYNAME
PARAMETER (MYNAME='EXAMPLE_6')
C
C Problem size parameters:
C
C Number of elements wide and long:
PARAMETER (NEWIDE = 20, NELONG = 100)
C
C Number of nodes wide and long:
PARAMETER (NPWIDE = NEWIDE+1, NPLONG = NELONG+1)
C
C Number of elements and equations:
PARAMETER (NUMEL = NEWIDE*NELONG, MAXEQ = NPWIDE*NPLONG)
C
C Data type = real:
PARAMETER (IDTYPE = 1)
C
C Number of RHS vectors:
PARAMETER (NUMRHS = 1)
C
C Number of element files:
PARAMETER (NUMSF = 1)
C
C Length of element real record:
PARAMETER (LENR = 4*4)
C
C Length of element integer record:
PARAMETER (LENI = 2 + 4)
C
C Length of element vector record:
PARAMETER (LENV = 2)
C
C Number of element vectors:
PARAMETER (NUMELV = 0)
C
C Skip forward reduction during call to RSDS:
C (Performed during call to RSDF)
PARAMETER (ISKIP = 1)
C
C FMS matrix and vector file attributes:
INTEGER LUA(25), LUA0(25)
INTEGER LUX(25), LUAX(25)
INTEGER LUS(25)
C
C Solution vector {X}
C RHS vector {B}
C Product [A]{X}, {AX}
REAL*8 RDUM(1), X(MAXEQ), B(MAXEQ), AX(MAXEQ)
COMMON /DATA/ RDUM, X, B, AX
C
C Big diagonal value to impose constraints:
REAL*8 BIGNUM
C
C Local variables:
REAL*8 EI, ERROR
C
C Common block to communicate with RSUBLK:
COMMON /EXAMP6/ BIGNUM, IBC
C
C FMS memory management requires the following arrays:
POINTER (CMD_PTR, CMD)
POINTER (RMD_PTR, RMD)
POINTER (IMD_PTR, IMD)
COMPLEX*16 CMD(0:1)
REAL*8 RMD(0:1)
INTEGER IMD(0:1)
DATA LUA0(1)/0/
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH (MYNAME)
CALL FMSIST ('MPOSDF',1)
CALL FMSIST ('INCORE',2)
CALL FMSIST ('MDATAU',2)
CALL FMSIGT ('MEMPTR',IMD_PTR)
CALL FMSIGT ('MEMPTR',RMD_PTR)
CALL FMSIGT ('MEMPTR',CMD_PTR)
WRITE(6,*) 'NUMBER OF ELEMENTS = ', NUMEL
WRITE(6,*) 'INITIAL NUMBER OF EQUATIONS = ', MAXEQ
BIGNUM = 1.0D30
C
C (2) Open FMS files:
CALL FMSOS (LENR, LENI, LENV, NUMELV, NUMEL, 'LUS', LUS)
C Generate element data:
CALL FMSIMG (IMD, NEWNUM, MAXEQ)
CALL ELGEN (NEWIDE,NELONG,NUMEL,IMD(NEWNUM),LUS,B,MAXEQ)
C Compute profile vector; eliminate zero constraints:
CALL FMSIMG (IMD, LOWEQ, MAXEQ)
CALL FMSWF (LUS, IMD(NEWNUM), MAXEQ, IMD(LOWEQ), NUMEQ)
WRITE(6,*) 'ACTUAL NUMBER OF EQUATIONS = ', NUMEQ
CALL RSDI (IMD(LOWEQ), NUMEQ, 'LUA', LUA)
CALL FMSIMR (IMD, LOWEQ, MAXEQ)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, X, NUMEQ, LUX)
C
C (3) Write data to FMS files:
C Elements written during ELGEN.
C Assembly during factoring.
C Impose boundary conditions on RHS vector:
C Flag IBC array for constrained diagonals:
C Reorder RHS terms for new numbering system:
CALL FMSIMG (IMD, IBC, NUMEQ)
DO 31 I = 1,MAXEQ
INEW = IMD(NEWNUM-1+I)
IF(INEW .EQ. 0) GO TO 31
IF(INEW .GT. 0) THEN
IMD(IBC-1+INEW) = 0
X(INEW) = B(I)
ELSE
INEW = -INEW
IMD(IBC-1+INEW) = 1
X(INEW) = BIGNUM*B(I)
END IF
31 CONTINUE
C
C (4) Perform matrix algebra:
CALL RSDAF (LUA0, RDUM, 0, LUS, NUMSF, LUA0, LUA,
1 LUX, LUX, NUMRHS)
CALL FMSIMR (IMD, IBC, NUMEQ)
CALL RSDS (LUA, LUX, LUX, NUMRHS, ISKIP)
C
C (5) Read data from FMS files:
C Not required.
C
C Compute [A]{X} from SUM([Si]{X}) to check error:
CALL FMSRMG (RMD, LAX, NUMEQ)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, RMD(LAX), NUMEQ, LUAX)
CALL RSDSVM (LUS, NUMSF, LUX, LUAX, NUMRHS)
ERROR = 0.0D0
C Expand solution vector to include zero constrained values:
DO 50 I = 1,MAXEQ
INEW = IMD(NEWNUM-1+I)
IF(INEW .EQ. 0) THEN
C Constrained zero value:
B(I) = 0.0D0
ELSE IF(INEW .LT. 0) THEN
C Constrained nonzero value:
INEW = -INEW
B(I) = X(INEW)
ELSE
C Value which was determined:
EI = ABS(RMD(LAX-1+INEW) - B(I))
IF(EI .GT. ERROR) ERROR = EI
B(I) = X(INEW)
END IF
50 CONTINUE
WRITE(6,*) 'MAXIMUM ERROR =',ERROR
CALL FMSIMR (IMD, NEWNUM, MAXEQ)
C
C (6) Close FMS files:
CALL FMSCV (LUAX)
CALL FMSRMR (RMD, LAX, NUMEQ)
CALL FMSCM (LUA)
CALL FMSCV (LUX)
CALL FMSCS (LUS)
CALL FMSPOP (MYNAME)
CALL FMSEND
END
C=======================================================================
SUBROUTINE RSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C=======================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
REAL*8 A(0:*), D(*), BIGNUM
COMMON /EXAMP6/ BIGNUM, IBC
C
POINTER (CMD_PTR, CMD)
POINTER (RMD_PTR, RMD)
POINTER (IMD_PTR, IMD)
COMPLEX*16 CMD(0:1)
REAL*8 RMD(0:1)
INTEGER IMD(0:1)
C
CALL FMSIGT ('MEMPTR',IMD_PTR)
CALL FMSIGT ('MEMPTR',RMD_PTR)
CALL FMSIGT ('MEMPTR',CMD_PTR)
IF(IROW2 .EQ. JCOL2) THEN
C This is a diagonal block:
DO 10 IROW = IROW1,IROW2
IF(IMD(IBC-1+IROW) .NE. 0) D(IROW) = BIGNUM
10 CONTINUE
END IF
RETURN
END
C=======================================================================
SUBROUTINE ELGEN (NEWIDE,NELONG,NUMEL,NEWNUM,LUS,B,NUMEQ)
C=======================================================================
INTEGER NEWIDE, NELONG, NUMEL, NUMEQ
INTEGER LUS(25), NEWNUM(NUMEQ)
INTEGER IREC(6), NSUB, ITYPE, ID(4)
CHARACTER*5 MYNAME
PARAMETER (MYNAME='ELGEN')
REAL*8 S(4,4), B(NUMEQ)
EQUIVALENCE (IREC(1),NSUB), (IREC(2),ITYPE), (IREC(3),ID)
C
CALL FMSPSH (MYNAME)
C
C Generate a 4-node quadrilateral finite element mesh that
C is NEWIDE elements wide by NELONG elements long. Use the
C same matrix S(4,4) for all elements. Impose the following
C solution values: X=0 along the top, X=1 along the bottom.
C
C Generate element connectivity array ID(4):
C Write element integer records:
CALL FMSEEK (LUS(2), 0)
NSUB = 4
ITYPE = 1
NPWIDE = NEWIDE + 1
DO 20 I = 1,NELONG
DO 10 J = 1,NEWIDE
ID(1) = I*NPWIDE + J
ID(2) = ID(1) + 1
ID(3) = ID(2) - NPWIDE
ID(4) = ID(1) - NPWIDE
CALL FMSSWR (LUS(2), IREC, 6)
10 CONTINUE
20 CONTINUE
C
C Generate element data array S(4,4):
DO 40 I = 1,4
DO 30 J = 1,4
S(I,J) = -1.0D0
30 CONTINUE
S(I,I) = 3.0D0
40 CONTINUE
C
C Write element real records:
CALL FMSEEK (LUS(1), 0)
DO 50 NE = 1,NUMEL
CALL FMSSWR (LUS(1), S, 16)
50 CONTINUE
C
C Generate right-hand side vector B(NUMEQ)
C Generate boundary condition code vector NEWNUM(NUMEQ)
DO 60 I = 1,NUMEQ
B(I) = 0.0D0
NEWNUM(I) = 1
60 CONTINUE
C Specify X(1:NPWIDE) = 0, X(NUMEQ-NPWIDE+1:NUMEQ) = 1
I2 = NUMEQ - NPWIDE
DO 70 I1 = 1,NPWIDE
NEWNUM(I1) = 0
B(I1) = 0.0D0
I2 = I2 + 1
NEWNUM(I2) = -1
B(I2) = 1.0D0
70 CONTINUE
CALL FMSPOP (MYNAME)
RETURN
END