C E X A M P L E 12
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 LUX(25), LUB(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 data:
LOGICAL ASK
INTEGER ASK_I
C
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
C
C Local variables:
INTEGER IDTYPE, ISTYPE, I, LB, LAX, LDISK, NV
INTEGER LENX, LENF
REAL*8 EI, ERROR
C
C Common block to communicate with fill routines:
COMMON /MYDATA/MOD, N, NRHS
DATA LOWEQ/-1/
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSIST ('MDATAU', 1)
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
1 CONTINUE
WRITE (6,*) 'Enter the FMS module number (1 to 5)'
WRITE (6,*) ' 1 = Real Symmetric'
WRITE (6,*) ' 2 = Real Nonsymmetric'
WRITE (6,*) ' 3 = Complex Hermitian'
WRITE (6,*) ' 4 = Complex Symmetric'
WRITE (6,*) ' 5 = Complex Nonsymmetric'
READ (5,*) MOD
IF( (MOD.LT.1) .OR. (MOD.GT.5) ) GO TO 1
N = ASK_I('Enter the number of equations')
NRHS = ASK_I('Enter the number of solution vectors')
WRITE (6,*) 'You may now alter any FMS parameter.'
WRITE (6,*) 'When you are finished, type the letters RETURN'
CALL FMSSET
IDTYPE = 1
ISTYPE = 1
IF(MOD .GE. 3) IDTYPE = 2
IF(MOD .EQ. 2) ISTYPE = 2
IF(MOD .EQ. 3) ISTYPE = 3
IF(MOD .EQ. 5) ISTYPE = 2
C
C (2) Open FMS files:
IF(MOD.EQ.1) CALL RSDI (LOWEQ, N, 'LUF', LUF)
IF(MOD.EQ.2) CALL RNDI (LOWEQ, N, 'LUF', LUF)
IF(MOD.EQ.3) CALL CHDI (LOWEQ, N, 'LUF', LUF)
IF(MOD.EQ.4) CALL CSDI (LOWEQ, N, 'LUF', LUF)
IF(MOD.EQ.5) CALL CNDI (LOWEQ, N, 'LUF', LUF)
CALL FMSOV (N, IDTYPE, NRHS, 'LUX', LUX)
CALL FMSOV (N, IDTYPE, NRHS, 'LUB', LUB)
C
C Populate test vector:
IF(IDTYPE .EQ. 1) THEN
LENX = LUX(4)
CALL FMSRMG (RMD, LB, LENX)
CALL FMSRMG (RMD, LAX, LENX)
DO I = 0,(N-1)
RMD(LB+I) = 0.0D0
END DO
RMD(LB) = 1.0D0
ELSE
LENX = LUX(4)/2
CALL FMSCMG (CMD, LB, LENX)
CALL FMSCMG (CMD, LAX, LENX)
DO I = 0,(N-1)
CMD(LB+I) = (0.0D0,0.0D0)
END DO
IF(MOD.EQ.3) THEN
CMD(LB) = (1.0D0,0.0D0)
ELSE
CMD(LB) = (1.0D0,1.0D0)
END IF
END IF
C
C (3) Write data to FMS files:
C
C Matrix file:
IF(MOD.EQ.1) CALL RSDA (LUS0, 0, LUF)
IF(MOD.EQ.2) CALL RNDA (LUS0, 0, LUF)
IF(MOD.EQ.3) CALL CHDA (LUS0, 0, LUF)
IF(MOD.EQ.4) CALL CSDA (LUS0, 0, LUF)
IF(MOD.EQ.5) CALL CNDA (LUS0, 0, LUF)
C
C Vector file:
LDISK = 1
DO NV = 1,NRHS
IF(IDTYPE .EQ. 1) THEN
CALL FMSWRT (LUB(1), LDISK, RMD(LB), N)
ELSE
CALL FMSWRT (LUB(1), LDISK, CMD(LB), 2*N)
END IF
LDISK = LDISK + LUB(4)
END DO
C
C (4) Perform matrix algebra:
C
C Populate and factor matrix [A]:
IF(MOD.EQ.1) CALL RSDF (LUF, LUF, LUB, LUX, NUMRED)
IF(MOD.EQ.2) CALL RNDF (LUF, LUF, LUB, LUX, NUMRED)
IF(MOD.EQ.3) CALL CHDF (LUF, LUF, LUB, LUX, NUMRED)
IF(MOD.EQ.4) CALL CSDF (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.1) CALL RSDS (LUF, LUB, LUX, NRHS, ISKIP)
IF(MOD.EQ.2) CALL RNDS (LUF, LUB, LUX, NRHS, ISKIP)
IF(MOD.EQ.3) CALL CHDS (LUF, LUB, LUX, NRHS, ISKIP)
IF(MOD.EQ.4) CALL CSDS (LUF, LUB, LUX, NRHS, ISKIP)
IF(MOD.EQ.5) CALL CNDS (LUF, LUB, LUX, NRHS, ISKIP)
C
C Compute the enengy [F] = {X}T{B}:
IF(ISTYPE .EQ. 1) THEN
LENF = NRHS*(NRHS+1)/2
ELSE
LENF = NRHS*NRHS
END IF
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMG (RMD, LF, LENF)
ELSE
CALL FMSCMG (CMD, LF, LENF)
END IF
IF(MOD.EQ.1) CALL RSDVVM (LUX, LUB, RMD(LF), NRHS)
IF(MOD.EQ.2) CALL RNDVVM (LUX, LUB, RMD(LF), NRHS, NRHS)
IF(MOD.EQ.3) CALL CHDVVM (LUX, LUB, CMD(LF), NRHS)
IF(MOD.EQ.4) CALL CSDVVM (LUX, LUB, CMD(LF), NRHS)
IF(MOD.EQ.5) CALL CNDVVM (LUX, LUB, CMD(LF), NRHS, NRHS)
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMR (RMD, LF, LENF)
ELSE
CALL FMSCMR (CMD, LF, LENF)
END IF
C
C Repopulate matrix [A]:
CALL FMSIST ('MDATAU', 1)
IF(MOD.EQ.1) CALL RSDA (LUS0, 0, LUF)
IF(MOD.EQ.2) CALL RNDA (LUS0, 0, LUF)
IF(MOD.EQ.3) CALL CHDA (LUS0, 0, LUF)
IF(MOD.EQ.4) CALL CSDA (LUS0, 0, LUF)
IF(MOD.EQ.5) CALL CNDA (LUS0, 0, LUF)
C
C Multiply {B} = [A]{X}:
IF(MOD.EQ.1) CALL RSDMVM (LUF, LUX, LUB, NRHS)
IF(MOD.EQ.2) CALL RNDMVM (LUF, LUX, LUB, NRHS)
IF(MOD.EQ.3) CALL CHDMVM (LUF, LUX, LUB, NRHS)
IF(MOD.EQ.4) CALL CSDMVM (LUF, LUX, LUB, NRHS)
IF(MOD.EQ.5) CALL CNDMVM (LUF, LUX, LUB, NRHS)
C
C (5) Read data from FMS files:
C Check the answer:
ERROR = 0.0D0
LDISK = 1
DO NV = 1,NRHS
IF(IDTYPE .EQ. 1) THEN
CALL FMSRED (LUB(1), LDISK, RMD(LAX), N)
DO I = 0,(N-1)
EI = ABS(RMD(LAX+I) - RMD(LB+I))
IF(EI .GT. ERROR) ERROR = EI
END DO
ELSE
CALL FMSRED (LUB(1), LDISK, CMD(LAX), 2*N)
DO I = 0,(N-1)
EI = ABS(CMD(LAX+I) - CMD(LB+I))
IF(EI .GT. ERROR) ERROR = EI
END DO
END IF
LDISK = LDISK + LUB(4)
END DO
WRITE(6,*) 'MAXIMUM ERROR =', ERROR
C
C (6) Close FMS files:
CALL FMSCM (LUF)
CALL FMSCV (LUX)
CALL FMSCV (LUB)
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMR (RMD, LB, LENX)
CALL FMSRMR (RMD, LAX, LENX)
ELSE
CALL FMSCMR (CMD, LB, LENX)
CALL FMSCMR (CMD, LAX, LENX)
END IF
IF(ASK('Do you want another solution?')) GO TO 1
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 CHUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C ================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
REAL*8 D(*), AIJ, AIJ_RE, AIJ_IM
COMPLEX*16 A(0:*)
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) = CMPLX(AIJ_RE,AIJ_IM)
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 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) = CMPLX(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) = CMPLX(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) = CMPLX(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) = CMPLX(AIJ_RE,AIJ_IM)
AIJ = SQRT(10.0*(AIJ - INT(AIJ)))
END DO
END DO
RETURN
END