This example solves a system of equations [A]{X}={B} and then computes the error {E}=[A]{X}-{B}, where the following may be specified as input: This example is designed to support pre-sales benchmarking of machines.
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, LB,  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, LB,  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