C E X A M P L E 12
C
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_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 FMSPSH (MYNAME)
CALL FMSIST ('MDATAU', 1)
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
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
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), LUB(4))
ELSE
CALL FMSWRT (LUB(1), LDISK, CMD(LB), LUB(4))
END IF
LDISK = LDISK + LUB(4)
END DO
C
C (4) Perform matrix algebra:
C
C 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 FMSPOP (MYNAME)
CALL FMSEND
END
C=======================================================================
REAL*8 FUNCTION AIJNUM(I,J)
C=======================================================================
INTEGER I, J, IJ
REAL*8 AI, AJ, AN
COMMON /MYDATA/MOD, N, NRHS
C Generate a matrix term for row I, column J that is
C reproducible and will cause pivoting.
AI = DFLOAT(I)
AJ = DFLOAT(J)
AN = DFLOAT(N)
AIJNUM = SQRT(2.0*AI*AJ)
AIJNUM = SQRT(10.0*(AIJNUM - INT(AIJNUM)))
RETURN
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(*)
LOGICAL IDIAG
C
D PRINT 2000, IROW1, IROW2, JCOL1, JCOL2
IF(IROW1 .EQ. JCOL1) THEN
IDIAG = .TRUE.
ELSE
IDIAG = .FALSE.
END IF
DO I = IROW1,IROW2
J1 = MAX0(LOWEQ(I),JCOL1)
J2 = MIN0(I-1,JCOL2)
DO J = J1,J2
A(LOCEQ(I)+IJSTEP*J) = AIJNUM(I,J)
END DO
IF(IDIAG) THEN
D(I) = AIJNUM(I,I)
END IF
END DO
RETURN
D2000 FORMAT('RSUBLK filling A(',I6,':',I6,',',I6,':',I6,')')
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(*)
C
D PRINT 2000, IROW1, IROW2, JCOL1, JCOL2
DO J = JCOL1,JCOL2
I1 = MAX0(LOWEQ(J),IROW1)
I2 = MIN0(J-1,IROW2)
DO I = I1,I2
A(LOCEQ(J)+IJSTEP*I) = AIJNUM(I,J)
END DO
END DO
RETURN
D2000 FORMAT('RNUBLK filling A(',I6,':',I6,',',I6,':',I6,')')
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_RE, AIJ_IM
COMPLEX*16 A(0:*)
LOGICAL IDIAG
C
D PRINT 2000, IROW1, IROW2, JCOL1, JCOL2
IF(IROW1 .EQ. JCOL1) THEN
IDIAG = .TRUE.
ELSE
IDIAG = .FALSE.
END IF
DO I = IROW1,IROW2
J1 = MAX0(LOWEQ(I),JCOL1)
J2 = MIN0(I-1,JCOL2)
DO J = J1,J2
AIJ_RE = AIJNUM(I,J)
AIJ_IM = - AIJ_RE
A(LOCEQ(I)+IJSTEP*J) = DCMPLX(AIJ_RE,AIJ_IM)
END DO
IF(IDIAG) THEN
D(I) = AIJNUM(I,I)
END IF
END DO
RETURN
D2000 FORMAT('CHUBLK filling A(',I6,':',I6,',',I6,':',I6,')')
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_RE, AIJ_IM
COMPLEX*16 A(0:*), D(*)
LOGICAL IDIAG
C
D PRINT 2000, IROW1, IROW2, JCOL1, JCOL2
IF(IROW1 .EQ. JCOL1) THEN
IDIAG = .TRUE.
ELSE
IDIAG = .FALSE.
END IF
DO I = IROW1,IROW2
J1 = MAX0(LOWEQ(I),JCOL1)
J2 = MIN0(I-1,JCOL2)
DO J = J1,J2
AIJ_RE = AIJNUM(I,J)
AIJ_IM = - AIJ
A(LOCEQ(I)+IJSTEP*J) = DCMPLX(AIJ_RE,AIJ_IM)
END DO
IF(IDIAG) THEN
AIJ_RE = AIJNUM(I,I)
AIJ_IM = - AIJ_RE
D(I) = DCMPLX(AIJ_RE,AIJ_IM)
END IF
END DO
RETURN
D2000 FORMAT('CSUBLK filling A(',I6,':',I6,',',I6,':',I6,')')
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_RE, AIJ_IM
COMPLEX*16 A(0:*), D(*)
C
D PRINT 2000, IROW1, IROW2, JCOL1, JCOL2
DO J = JCOL1,JCOL2
I1 = MAX0(LOWEQ(J),IROW1)
I2 = MIN0(J-1,IROW2)
DO I = I1,I2
AIJ_RE = AIJNUM(I,J)
AIJ_IM = - AIJ_RE
A(LOCEQ(J)+IJSTEP*I) = DCMPLX(AIJ_RE,AIJ_IM)
END DO
END DO
RETURN
D2000 FORMAT('CNUBLK filling A(',I6,':',I6,',',I6,':',I6,')')
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:*)
C
D PRINT 2000, NUMEQ, JEQN1, JEQN2
DO J=JEQN1,JEQN2
DO I=1,NUMEQ
LIJ = LOCI(I) + LOCJ(J,LUFLAG(I))
A(LIJ) = AIJNUM(I,J)
END DO
END DO
RETURN
D2000 FORMAT('RNUSLB filling A(1:',I6,',',I6,':',I6,')')
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_RE, AIJ_IM
C
D PRINT 2000, NUMEQ, JEQN1, JEQN2
DO J=JEQN1,JEQN2
DO I=1,NUMEQ
LIJ = LOCI(I) + LOCJ(J,LUFLAG(I))
AIJ_RE = AIJNUM(I,J)
AIJ_IM = - AIJ_RE
A(LIJ) = DCMPLX(AIJ_RE,AIJ_IM)
END DO
END DO
RETURN
D2000 FORMAT('CNUSLB filling A(1:',I6,',',I6,':',I6,')')
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