C E X A M P L E 18
C
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_18')
C
C Number of vectors to reduce during factoring:
PARAMETER (NUMRED = 0)
C
C Skip operations during solving (no):
PARAMETER (ISKIP = 0)
C
C FMS matrix and vector file attributes:
INTEGER LUF(25)
INTEGER LUB(25), LUC(25), LUX(25), LUY(25), LUE(25)
INTEGER LUS0(25)
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)
C
C Input parameters:
LOGICAL ASK
INTEGER ASK_I
C
C Local variables:
INTEGER IDTYPE, I, LB, LC, LX, LY, LE, LDISK, NV
INTEGER LENB, LENC, LENX, LENY, LENE
REAL*8 EI, ERROR, ERROR2
REAL*8 R_VCX, R_VYB, RZERO, RONE
COMPLEX*16 C_VCX, C_VYB, CZERO, CONE
C
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C Common block to communicate with fill routines:
COMMON /MYDATA/MOD, N, NRHS
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH (MYNAME)
CALL FMSIST ('MDATAU', 1)
RZERO = 0.0D0
RONE = 1.0D0
CZERO = DCMPLX(0.0D0, 0.0D0)
CONE = DCMPLX(1.0D0, 1.0D0)
1 CONTINUE
WRITE (6,*) 'Enter the FMS module number (2 or 5)'
WRITE (6,*) ' 2 = Real Nonsymmetric'
WRITE (6,*) ' 5 = Complex Nonsymmetric'
READ (5,*) MOD
IF(MOD .EQ. 2) THEN
IDTYPE = 1
ELSE IF (MOD .EQ. 5) THEN
IDTYPE = 2
ELSE
GO TO 1
END IF
N = ASK_I('Enter the number of equations')
NRHS = ASK_I('Enter the number of vectors')
WRITE (6,*) 'You may now alter any FMS parameter.'
WRITE (6,*) 'When you are finished, type the letters RETURN'
CALL FMSSET
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C (2) Open FMS files:
IF(MOD.EQ.2) CALL RNDI (LOWEQ, N, 'LUF', LUF)
IF(MOD.EQ.5) CALL CNDI (LOWEQ, N, 'LUF', LUF)
CALL FMSOV (N, IDTYPE, NRHS, 'LUB', LUB)
CALL FMSOV (N, IDTYPE, NRHS, 'LUC', LUC)
CALL FMSOV (N, IDTYPE, NRHS, 'LUX', LUX)
CALL FMSOV (N, IDTYPE, NRHS, 'LUY', LUY)
CALL FMSOV (N, IDTYPE, NRHS, 'LUE', LUE)
C
C Allocate storage for one {B}, {C}, {X}, {Y} and {E}:
C Initialize {B} and {C}:
IF(IDTYPE .EQ. 1) THEN
LENB = LUB(4)
CALL FMSRMG (RMD, LB, LENB)
LENC = LUC(4)
CALL FMSRMG (RMD, LC, LENC)
LENX = LUX(4)
CALL FMSRMG (RMD, LX, LENX)
LENY = LUY(4)
CALL FMSRMG (RMD, LY, LENY)
LENE = LUE(4)
CALL FMSRMG (RMD, LE, LENE)
DO I = 0,(N-1)
RMD(LB+I) = RONE
RMD(LC+I) = RONE
END DO
ELSE
LENB = LUB(4)/2
CALL FMSCMG (CMD, LB, LENB)
LENC = LUC(4)/2
CALL FMSCMG (CMD, LC, LENC)
LENX = LUX(4)/2
CALL FMSCMG (CMD, LX, LENX)
LENY = LUY(4)/2
CALL FMSCMG (CMD, LY, LENY)
LENE = LUE(4)/2
CALL FMSCMG (CMD, LE, LENE)
DO I = 0,(N-1)
CMD(LB+I) = CONE
CMD(LC+I) = CONE
END DO
END IF
C
C (3) Write data to FMS files:
C
C Matrix file:
IF(MOD.EQ.2) CALL RNDA (LUS0, 0, LUF)
IF(MOD.EQ.5) CALL CNDA (LUS0, 0, LUF)
C
C Vector files:
LDISK = 1
DO NV = 1,NRHS
IF(IDTYPE .EQ. 1) THEN
CALL FMSWRT (LUB(1), LDISK, RMD(LB), LUB(4))
CALL FMSWRT (LUC(1), LDISK, RMD(LC), LUC(4))
ELSE
CALL FMSWRT (LUB(1), LDISK, CMD(LB), LUB(4))
CALL FMSWRT (LUC(1), LDISK, CMD(LC), LUC(4))
END IF
LDISK = LDISK + LUB(4)
END DO
C
C (4) Perform matrix algebra:
C
C Populate and factor matrix [A]:
IF(MOD.EQ.2) CALL RNDF (LUF, LUF, LUB, LUX, NUMRED)
IF(MOD.EQ.5) CALL CNDF (LUF, LUF, LUB, LUX, NUMRED)
C
C Solve [A]{X} = {B}:
IF(MOD.EQ.2) CALL RNDS (LUF, LUB, LUX, NRHS, ISKIP)
IF(MOD.EQ.5) CALL CNDS (LUF, LUB, LUX, NRHS, ISKIP)
C
C Solve [A]T{Y} = {C}:
C Set the FMS Parameter LUTRAN to solve with [A]T:
CALL FMSIST ('LUTRAN', 1)
IF(MOD.EQ.2) CALL RNDS (LUF, LUC, LUY, NRHS, ISKIP)
IF(MOD.EQ.5) CALL CNDS (LUF, LUC, LUY, NRHS, ISKIP)
C Reset the FMS Parameter LUTRAN for a normal solve:
CALL FMSIST ('LUTRAN', 0)
C
C Repopulate matrix [A]:
CALL FMSIST ('MDATAU', 1)
IF(MOD.EQ.2) CALL RNDA (LUS0, 0, LUF)
IF(MOD.EQ.5) CALL CNDA (LUS0, 0, LUF)
C
C Multiply {E} = {B} - [A]{X}:
CALL FMSMM (LUE, LUB, -1, 'N', LUF, LUX)
C
C (5) Read data from FMS files:
C Compute VCX and VYB:
C Check the answer:
ERROR = RZERO
ERROR2 = RZERO
LDISK = 1
NWREAD = IDTYPE*N
DO NV = 1,NRHS
IF(IDTYPE .EQ. 1) THEN
C Read in {B}nv, {C}nv, {X}nv, {Y}nv and {E}nv
CALL FMSRED (LUB(1), LDISK, RMD(LB), NWREAD)
CALL FMSRED (LUC(1), LDISK, RMD(LC), NWREAD)
CALL FMSRED (LUX(1), LDISK, RMD(LX), NWREAD)
CALL FMSRED (LUY(1), LDISK, RMD(LY), NWREAD)
CALL FMSRED (LUE(1), LDISK, RMD(LE), NWREAD)
LDISK = LDISK + LUB(4)
R_VCX = RZERO
R_VYB = RZERO
DO I = 0,(N-1)
R_VCX = R_VCX + RMD(LC+I) * RMD(LX+I)
R_VYB = R_VYB + RMD(LY+I) * RMD(LB+I)
EI = ABS(RMD(LE+I))
IF(EI .GT. ERROR) ERROR = EI
END DO
EI = ABS(R_VCX - R_VYB)
IF(EI .GT. ERROR2) ERROR2 = EI
ELSE
C Read in {B}nv, {C}nv, {X}nv, {Y}nv and {E}nv
CALL FMSRED (LUB(1), LDISK, CMD(LB), NWREAD)
CALL FMSRED (LUC(1), LDISK, CMD(LC), NWREAD)
CALL FMSRED (LUX(1), LDISK, CMD(LX), NWREAD)
CALL FMSRED (LUY(1), LDISK, CMD(LY), NWREAD)
CALL FMSRED (LUE(1), LDISK, CMD(LE), NWREAD)
LDISK = LDISK + LUB(4)
C_VCX = CZERO
C_VYB = CZERO
DO I = 0,(N-1)
C_VCX = C_VCX + CMD(LC+I) * CMD(LX+I)
C_VYB = C_VYB + CMD(LY+I) * CMD(LB+I)
EI = ABS(CMD(LE+I))
IF(EI .GT. ERROR) ERROR = EI
END DO
EI = ABS(C_VCX - C_VYB)
IF(EI .GT. ERROR2) ERROR2 = EI
END IF
END DO
WRITE(6,*) 'MAXIMUM ERROR IN {B}-[A]{X} =', ERROR
WRITE(6,*) 'MAXIMUM ERROR IN {C}T{X}-{Y}T{B} =', ERROR2
C
C (6) Close FMS files:
CALL FMSCV (LUE)
CALL FMSCV (LUY)
CALL FMSCV (LUX)
CALL FMSCV (LUC)
CALL FMSCV (LUB)
CALL FMSCM (LUF)
C
C Return the memory:
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMR (RMD, LE, LENE)
CALL FMSRMR (RMD, LY, LENY)
CALL FMSRMR (RMD, LX, LENX)
CALL FMSRMR (RMD, LC, LENC)
CALL FMSRMR (RMD, LB, LENB)
ELSE
CALL FMSCMR (CMD, LE, LENE)
CALL FMSCMR (CMD, LY, LENY)
CALL FMSCMR (CMD, LX, LENX)
CALL FMSCMR (CMD, LC, LENC)
CALL FMSCMR (CMD, LB, LENB)
END IF
IF(ASK('Do you want another solution?')) GO TO 1
CALL FMSPOP (MYNAME)
CALL FMSEND
END
C=======================================================================
LOGICAL FUNCTION ASK(QUESTION)
C=======================================================================
CHARACTER* (*) QUESTION
CHARACTER*1 IYN
WRITE(6,2000) QUESTION
READ (5,1000) IYN
IF( (IYN .EQ. 'Y') .OR. (IYN .EQ. 'y') ) THEN
ASK = .TRUE.
ELSE
ASK = .FALSE.
END IF
RETURN
1000 FORMAT (A)
2000 FORMAT (1X,A,' (y,n)>')
END
C=======================================================================
INTEGER FUNCTION ASK_I(STRING)
C=======================================================================
CHARACTER* (*) STRING
WRITE(6,2000) STRING
READ (5,*) ASK_I
RETURN
2000 FORMAT (1X,A,'>')
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(*), AIJ
LOGICAL IDIAG
C
IF(IROW1 .EQ. JCOL1) THEN
IDIAG = .TRUE.
ELSE
IDIAG = .FALSE.
END IF
AIJ = SQRT(2.0)
DO I = IROW1,IROW2
J1 = MAX0(LOWEQ(I),JCOL1)
J2 = MIN0(I-1,JCOL2)
DO J = J1,J2
A(LOCEQ(I)+IJSTEP*J) = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
IF(IDIAG) THEN
D(I) = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END IF
END DO
RETURN
END
C=======================================================================
SUBROUTINE RNUBLK (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(*), AIJ
C
AIJ = SQRT(2.0)
DO J = JCOL1,JCOL2
I1 = MAX0(LOWEQ(J),IROW1)
I2 = MIN0(J-1,IROW2)
DO I = I1,I2
A(LOCEQ(J)+IJSTEP*I) = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE CSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C=======================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
REAL*8 AIJ, AIJ_RE, AIJ_IM
COMPLEX*16 A(0:*), D(*)
LOGICAL IDIAG
C
IF(IROW1 .EQ. JCOL1) THEN
IDIAG = .TRUE.
ELSE
IDIAG = .FALSE.
END IF
AIJ = SQRT(2.0)
DO I = IROW1,IROW2
J1 = MAX0(LOWEQ(I),JCOL1)
J2 = MIN0(I-1,JCOL2)
DO J = J1,J2
AIJ_RE = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
AIJ_IM = AIJ
A(LOCEQ(I)+IJSTEP*J) = DCMPLX(AIJ_RE,AIJ_IM)
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
IF(IDIAG) THEN
AIJ_RE = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
AIJ_IM = AIJ
D(I) = DCMPLX(AIJ_RE,AIJ_IM)
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END IF
END DO
RETURN
END
C=======================================================================
SUBROUTINE CNUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C=======================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
REAL*8 AIJ, AIJ_RE, AIJ_IM
COMPLEX*16 A(0:*), D(*)
C
AIJ = SQRT(2.0)
DO J = JCOL1,JCOL2
I1 = MAX0(LOWEQ(J),IROW1)
I2 = MIN0(J-1,IROW2)
DO I = I1,I2
AIJ_RE = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
AIJ_IM = AIJ
A(LOCEQ(J)+IJSTEP*I) = DCMPLX(AIJ_RE,AIJ_IM)
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE RNUSLB (A, LOCI, LOCJ, LUFLAG, JEQN1, JEQN2, NUMEQ)
C=======================================================================
INTEGER JEQN1, JEQN2, NUMEQ
INTEGER LOCI(NUMEQ), LOCJ(JEQN1:JEQN2,2), LUFLAG(NUMEQ)
REAL*8 A(0:*), AIJ
C
AIJ = SQRT(2.0)
DO I=1,NUMEQ
DO J=JEQN1,JEQN2
LIJ = LOCI(I) + LOCJ(J,LUFLAG(I))
A(LIJ) = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE CNUSLB (A, LOCI, LOCJ, LUFLAG, JEQN1, JEQN2, NUMEQ)
C=======================================================================
INTEGER JEQN1, JEQN2, NUMEQ
INTEGER LOCI(NUMEQ), LOCJ(JEQN1:JEQN2,2), LUFLAG(NUMEQ)
COMPLEX*16 A(0:*)
REAL*8 AIJ, AIJ_RE, AIJ_IM
C
AIJ = SQRT(2.0)
DO I=1,NUMEQ
DO J=JEQN1,JEQN2
LIJ = LOCI(I) + LOCJ(J,LUFLAG(I))
AIJ_RE = AIJ
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
AIJ_IM = AIJ
A(LIJ) = DCMPLX(AIJ_RE,AIJ_IM)
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
END DO
RETURN
END