C E X A M P L E 1 9
C
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_19')
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:
C Unfactored matrix:
INTEGER LUA(25)
C Factored matrix:
INTEGER LUF(25)
C R.H.S. and solution Vectors:
INTEGER LUX(25)
INTEGER LUA0(25)
DATA LUA0(1)/0/
INTEGER LUX0(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 Local variables:
C FMS Module number:
INTEGER MOD
C Number of equations:
INTEGER N
REAL*8 R8_N
C Number of R.H.S. vectors:
INTEGER NRHS
C Data type:
INTEGER IDTYPE
C Loop counter:
INTEGER I
C Vector starting address:
INTEGER LX
C Vector storage length:
INTEGER LENX
C Matrix starting address:
INTEGER LA
C Disk location:
INTEGER LDISK
C Vector number:
INTEGER NV
C Current error:
REAL*8 EI
C Maximum overall error:
REAL*8 ERROR
C Input function:
LOGICAL ASK
INTEGER ASK_I
C Diagonal value:
COMPLEX*16 C16_DIA
C Constants:
REAL*8 R8_ZERO
DATA R8_ZERO/0.0D0/
COMPLEX*16 C16_ZERO
DATA C16_ZERO/(0.0D0,0.0D0)/
REAL*8 R8_ONE
DATA R8_ONE/1.0D0/
COMPLEX*16 C16_ONE(1)
DATA C16_ONE/(1.0D0,-0.0D0)/
REAL*8 R8_ALPHA(1)
DATA R8_ALPHA/1.0D0/
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C Common block to communicate with user supplied routines:
COMMON /MYDATA/NCELL, NCELLS, LA
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH (MYNAME)
CALL FMSIST ('MFMAT' , 2)
CALL FMSIST ('LOGTIM', 3)
CALL FMSIST ('IPRF' , 1026)
CALL FMSIST ('IPRS' , 1026)
100 CONTINUE
1 CONTINUE
WRITE (6,*) 'The FMS modules are numbered as follows:'
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'
MOD = ASK_I('Enter the FMS module number (1 to 5)')
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')
NCELL = ASK_I('Enter the cell size to average and print')
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 Check MFMAT:
CALL FMSIGT ('MFMAT', MFMAT)
IF(MFMAT .EQ. 3) THEN
WRITE(6,*) 'Cannot save [A] while pivoting'
WRITE(6,*) 'Building [A] as a Block Matrix'
MFMAT = 2
END IF
IDTYPE = 1
R8_N = DFLOAT(N)
IF(MOD .EQ. 3) IDTYPE = 2
IF(MOD .EQ. 4) IDTYPE = 2
IF(MOD .EQ. 5) IDTYPE = 2
NCELLS = (N + NCELL - 1)/NCELL
C
C (2) Open FMS files:
IF(MOD.EQ.1) CALL RSDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.2) CALL RNDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.3) CALL CHDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.4) CALL CSDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.5) CALL CNDI (LOWEQ, N, 'LUA', LUA)
CALL FMSOM (LUA, 'LUF', LUF)
CALL FMSOV (N, IDTYPE, NRHS, 'LUX', LUX)
C
C (3) Write data to FMS files:
C
C Matrix File:
IF(IDTYPE .EQ. 1) THEN
C Real data:
C Allocate a scratch vector:
LENX = LUX(4)
CALL FMSRMG (RMD, LX, LENX)
C Initialize it to zero:
DO I = 1,LENX
RMD(LX-1+I) = R8_ZERO
END DO
C
C Write the R.H.S. vector file:
RMD(LX) = R8_ONE
LDISK = 1
DO NV = 1,NRHS
CALL FMSWRT (LUX(1), LDISK, RMD(LX), LUX(4))
LDISK = LDISK + LUX(4)
END DO
RMD(LX) = R8_ZERO
C
C Write the matrix file:
CALL FMSROW (0, RMD(LX), LUA)
RMD(LX) = R8_N
DO I=2,N
RMD(LX-1+I) = -R8_ONE
END DO
CALL FMSROW (1, RMD(LX), LUA)
RMD(LX) = -R8_ONE
DO I=2,N
RMD(LX-1+I) = R8_ZERO
END DO
DO I=2,N
RMD(LX-1+I) = R8_ONE
CALL FMSROW (I, RMD(LX), LUA)
RMD(LX-1+I) = R8_ZERO
END DO
CALL FMSROW (N+1, RMD(LX), LUA)
ELSE
C Complex data:
C Allocate a scratch vector:
LENX = LUX(4)/2
CALL FMSCMG (CMD, LX, LENX)
DO I = 1,LENX
CMD(LX-1+I) = C16_ZERO
END DO
C
C Populate the vector:
IF(MOD.EQ.3) THEN
C16_DIA = (1.0D0,0.0D0)
ELSE
C16_DIA = (1.0D0,1.0D0)
END IF
CMD(LX) = C16_DIA
LDISK = 1
DO NV = 1,NRHS
CALL FMSWRT (LUX(1), LDISK, CMD(LX), LUX(4))
LDISK = LDISK + LUX(4)
END DO
C
C Populate the matrix:
CALL FMSROW (0, CMD(LX), LUA)
IF(MOD .EQ. 3) THEN
CMD(LX) = DCMPLX(R8_N,0.0D0)
ELSE
CMD(LX) = DCMPLX(R8_N,R8_N)
END IF
DO I=2,N
CMD(LX-1+I) = -C16_DIA
END DO
CALL FMSROW (1, CMD(LX), LUA)
CMD(LX) = -C16_DIA
DO I=2,N
CMD(LX-1+I) = C16_ZERO
END DO
DO I=2,N
CMD(LX-1+I) = C16_DIA
CALL FMSROW (I, CMD(LX), LUA)
CMD(LX-1+I) = C16_ZERO
END DO
CALL FMSROW (N+1, CMD(LX), LUA)
END IF
C
C (4) Perform matrix algebra:
C Assemble and factor the matrix, saving [A]:
CALL FMSIST ('MDATAU', 0)
IF(MOD.EQ.1) CALL RSDAF (LUA, R8_ALPHA,
1 1, LUS0, 0, LUA, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.2) CALL RNDAF (LUA, R8_ALPHA,
1 1, LUS0, 0, LUA, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.3) CALL CHDAF (LUA, R8_ALPHA,
1 1, LUS0, 0, LUA, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.4) CALL CSDAF (LUA, C16_ONE,
1 1, LUS0, 0, LUA, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.5) CALL CNDAF (LUA, C16_ONE,
1 1, LUS0, 0, LUA, LUF, LUX, LUX, NUMRED)
C
C Perform forward reduction and back-substitution:
IF(MOD.EQ.1) CALL RSDS (LUF, LUX, LUX, NRHS, ISKIP)
IF(MOD.EQ.2) CALL RNDS (LUF, LUX, LUX, NRHS, ISKIP)
IF(MOD.EQ.3) CALL CHDS (LUF, LUX, LUX, NRHS, ISKIP)
IF(MOD.EQ.4) CALL CSDS (LUF, LUX, LUX, NRHS, ISKIP)
IF(MOD.EQ.5) CALL CNDS (LUF, LUX, LUX, NRHS, ISKIP)
C
C Accumulate and print average windows of [A]:
C Allocate storage for an incore array of cells:
NN = NCELLS*NCELLS
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMG (RMD, LA, NN)
DO I=1,NN
RMD(LA-1+I) = R8_ZERO
END DO
ELSE
CALL FMSCMG (CMD, LA, NN)
DO I=1,NN
CMD(LA-1+I) = C16_ZERO
END DO
END IF
CALL FMSIST ('MDATAU', 2)
IF(MOD.EQ.1) CALL RSDAF (LUA, R8_ALPHA, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.2) CALL RNDAF (LUA, R8_ALPHA, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.3) CALL CHDAF (LUA, R8_ALPHA, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.4) CALL CSDAF (LUA, C16_ONE, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.5) CALL CNDAF (LUA, C16_ONE, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
C
C Fill upper triangle of symmetric matrices:
IF(MOD.EQ.1) CALL RFILL (RMD(LA), NCELLS)
IF(MOD.EQ.3) CALL CFILL (CMD(LA), NCELLS)
IF(MOD.EQ.4) CALL CFILL (CMD(LA), NCELLS)
C
C You can change the format of the output.
C These are the defaults:
CALL FMSCST ('RFMAT', '(E12.4)')
CALL FMSCST ('CFMAT', '(E12.4,1H,,E11.4)')
C
IF(IDTYPE .EQ. 1) THEN
WRITE(6,2000)
CALL FMSPRF (RMD(LA), NCELLS, NCELLS, 2, IDTYPE)
CALL RAVG (RMD(LA), NCELLS, NCELL, N)
WRITE(6,2001)
CALL FMSPRF (RMD(LA), NCELLS, NCELLS, 2, IDTYPE)
ELSE
WRITE(6,2000)
CALL FMSPRF (CMD(LA), NCELLS, NCELLS, 2, IDTYPE)
CALL CAVG (CMD(LA), NCELLS, NCELL, N)
WRITE(6,2001)
CALL FMSPRF (CMD(LA), NCELLS, NCELLS, 2, IDTYPE)
END IF
C
C Accumulate and print average windows of [F]:
IF(IDTYPE .EQ. 1) THEN
DO I=1,NN
RMD(LA-1+I) = R8_ZERO
END DO
ELSE
DO I=1,NN
CMD(LA-1+I) = C16_ZERO
END DO
END IF
IF(MOD.EQ.1) CALL RSDAF (LUF, R8_ALPHA, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.2) CALL RNDAF (LUF, R8_ALPHA, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.3) CALL CHDAF (LUF, R8_ALPHA, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.4) CALL CSDAF (LUF, C16_ONE, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
IF(MOD.EQ.5) CALL CNDAF (LUF, C16_ONE, 1, LUS0, 0, LUA0, LUA0,
1 LUX0, LUX0, 0)
C
C Fill upper triangle of symmetric matrices:
IF(MOD.EQ.1) CALL RFILL (RMD(LA), NCELLS)
IF(MOD.EQ.3) CALL CFILL (CMD(LA), NCELLS)
IF(MOD.EQ.4) CALL CFILL (CMD(LA), NCELLS)
C
IF(IDTYPE .EQ. 1) THEN
WRITE(6,2002)
CALL FMSPRF (RMD(LA), NCELLS, NCELLS, 2, IDTYPE)
CALL RAVG (RMD(LA), NCELLS, NCELL, N)
WRITE(6,2003)
CALL FMSPRF (RMD(LA), NCELLS, NCELLS, 2, IDTYPE)
ELSE
WRITE(6,2002)
CALL FMSPRF (CMD(LA), NCELLS, NCELLS, 2, IDTYPE)
CALL CAVG (CMD(LA), NCELLS, NCELL, N)
WRITE(6,2003)
CALL FMSPRF (CMD(LA), NCELLS, NCELLS, 2, IDTYPE)
END IF
C
C (5) Read data from FMS files:
C Check the answer:
ERROR = 0.0D0
LDISK = 1
DO 60 NV = 1,NRHS
IF(IDTYPE .EQ. 1) THEN
CALL FMSRED (LUX(1), LDISK, RMD(LX), N)
DO 50 I = 1,N
EI = ABS(RMD(LX-1+I) - 1.0D0)
IF(EI .GT. ERROR) ERROR = EI
50 CONTINUE
ELSE
CALL FMSRED (LUX(1), LDISK, CMD(LX), 2*N)
DO 51 I = 1,N
EI = ABS(CMD(LX-1+I) - 1.0D0)
IF(EI .GT. ERROR) ERROR = EI
51 CONTINUE
END IF
LDISK = LDISK + LUX(4)
60 CONTINUE
WRITE(6,*) 'MAXIMUM ERROR =', ERROR
C
C (6) Close FMS files:
CALL FMSCM (LUF)
CALL FMSCM (LUA)
CALL FMSCV (LUX)
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMR (RMD, LA, NN)
CALL FMSRMR (RMD, LX, LENX)
ELSE
CALL FMSCMR (CMD, LA, NN)
CALL FMSCMR (CMD, LX, LENX)
END IF
IF(ASK('Do you want another solution?')) GO TO 100
CALL FMSPOP (MYNAME)
CALL FMSEND
2000 FORMAT (/'Accumulated values of [A] in cells.')
2001 FORMAT (/'Average values of [A] in cells.')
2002 FORMAT (/'Accumulated values of [F] in cells.')
2003 FORMAT (/'Average values of [F] in cells.')
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(*)
COMMON /MYDATA/NCELL, NCELLS, LA
C
C FMS memory management requires the following arrays:
POINTER (RMD_PTR, RMD)
REAL*8 RMD(0:1)
WRITE(6,2000) IROW1, IROW2, JCOL1, JCOL2, IJSTEP
CALL FMSIGT ('MEMPTR', RMD_PTR)
C
C Add this data to the appropriate cells:
DO I = IROW1, IROW2
IC = (I+NCELL-1)/NCELL
JMAX = MIN0(I, JCOL2)
DO J = JCOL1, JMAX
JC = (J+NCELL-1)/NCELL
LCIJ = LA - 1 + IC + NCELLS*(JC-1)
IF(J .LT. I) THEN
C This is an off-diagonal term:
LIJ = LOCEQ(I) + IJSTEP*J
RMD(LCIJ) = RMD(LCIJ) + A(LIJ)
ELSE
C This is a diagonal term:
RMD(LCIJ) = RMD(LCIJ) + D(I)
END IF
END DO
END DO
WRITE(6,2001)
RETURN
2000 FORMAT (
1 'RSUBLK filling A(',I6,':',I6,',',I6,':',I6,'), IJSTEP=',I6)
2001 FORMAT ('RSUBLK ending.')
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(*)
COMMON /MYDATA/NCELL, NCELLS, LA
C
C FMS memory management requires the following arrays:
POINTER (RMD_PTR, RMD)
REAL*8 RMD(0:1)
WRITE(6,2000) IROW1, IROW2, JCOL1, JCOL2, IJSTEP
CALL FMSIGT ('MEMPTR', RMD_PTR)
C
C Add this data to the appropriate cells:
DO J = JCOL1, JCOL2
JC = (J+NCELL-1)/NCELL
IMAX = MIN0(J-1, IROW2)
DO I = IROW1, IMAX
IC = (I+NCELL-1)/NCELL
LCIJ = LA - 1 + IC + NCELLS*(JC-1)
LIJ = LOCEQ(J) + IJSTEP*I
RMD(LCIJ) = RMD(LCIJ) + A(LIJ)
END DO
END DO
WRITE(6,2001)
RETURN
2000 FORMAT (
1 'RNUBLK filling A(',I6,':',I6,',',I6,':',I6,'), IJSTEP=',I6)
2001 FORMAT ('RNUBLK ending.')
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(*)
COMPLEX*16 A(0:*)
COMMON /MYDATA/NCELL, NCELLS, LA
C
C FMS memory management requires the following arrays:
POINTER (CMD_PTR, CMD)
COMPLEX*16 CMD(0:1)
WRITE(6,2000) IROW1, IROW2, JCOL1, JCOL2, IJSTEP
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C Add this data to the appropriate cells:
DO I = IROW1, IROW2
IC = (I+NCELL-1)/NCELL
JMAX = MIN0(I, JCOL2)
DO J = JCOL1, JMAX
JC = (J+NCELL-1)/NCELL
LCIJ = LA - 1 + IC + NCELLS*(JC-1)
IF(J .LT. I) THEN
C This is an off-diagonal term:
LIJ = LOCEQ(I) + IJSTEP*J
CMD(LCIJ) = CMD(LCIJ) + A(LIJ)
ELSE
C This is a diagonal term:
CMD(LCIJ) = CMD(LCIJ) + DCMPLX(D(I),0.0D0)
END IF
END DO
END DO
WRITE(6,2001)
RETURN
2000 FORMAT (
1 'CHUBLK filling A(',I6,':',I6,',',I6,':',I6,'), IJSTEP=',I6)
2001 FORMAT ('CHUBLK ending.')
END
C ================================================================
SUBROUTINE CSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C ================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
COMPLEX*16 A(0:*), D(*)
COMMON /MYDATA/NCELL, NCELLS, LA
C
C FMS memory management requires the following arrays:
POINTER (CMD_PTR, CMD)
COMPLEX*16 CMD(0:1)
WRITE(6,2000) IROW1, IROW2, JCOL1, JCOL2, IJSTEP
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C Add this data to the appropriate cells:
DO I = IROW1, IROW2
IC = (I+NCELL-1)/NCELL
JMAX = MIN0(I, JCOL2)
DO J = JCOL1, JMAX
JC = (J+NCELL-1)/NCELL
LCIJ = LA - 1 + IC + NCELLS*(JC-1)
IF(J .LT. I) THEN
C This is an off-diagonal term:
LIJ = LOCEQ(I) + IJSTEP*J
CMD(LCIJ) = CMD(LCIJ) + A(LIJ)
ELSE
C This is a diagonal term:
CMD(LCIJ) = CMD(LCIJ) + D(I)
END IF
END DO
END DO
WRITE(6,2001)
RETURN
2000 FORMAT (
1 'CSUBLK filling A(',I6,':',I6,',',I6,':',I6,'), IJSTEP=',I6)
2001 FORMAT ('CSUBLK ending.')
END
C ================================================================
SUBROUTINE CNUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C ================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
COMPLEX*16 A(0:*), D(*)
COMMON /MYDATA/NCELL, NCELLS, LA
C
C FMS memory management requires the following arrays:
POINTER (CMD_PTR, CMD)
COMPLEX*16 CMD(0:1)
WRITE(6,2000) IROW1, IROW2, JCOL1, JCOL2, IJSTEP
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C Add this data to the appropriate cells:
DO J = JCOL1, JCOL2
JC = (J+NCELL-1)/NCELL
IMAX = MIN0(J-1, IROW2)
DO I = IROW1, IMAX
IC = (I+NCELL-1)/NCELL
LCIJ = LA - 1 + IC + NCELLS*(JC-1)
LIJ = LOCEQ(J) + IJSTEP*I
CMD(LCIJ) = CMD(LCIJ) + A(LIJ)
END DO
END DO
WRITE(6,2001)
RETURN
2000 FORMAT (
1 'CNUBLK filling A(',I6,':',I6,',',I6,':',I6,'), IJSTEP=',I6)
2001 FORMAT ('CNUBLK ending.')
END
C=======================================================================
LOGICAL FUNCTION ASK(QUESTION)
C=======================================================================
CHARACTER* (*) QUESTION
CHARACTER*1 IYN
WRITE(6,2000) QUESTION
READ (5,1000) IYN
WRITE(6,2001) 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)>')
2001 FORMAT (4X,'You entered ',A)
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 RFILL (A, N)
C=======================================================================
C Fill upper triangle of [A]
REAL*8 A(N,N)
INTEGER N
INTEGER I, J
DO J=2,N
DO I=1,(J-1)
A(I,J) = A(J,I)
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE CFILL (A, N)
C=======================================================================
C Fill upper triangle of [A]
COMPLEX*16 A(N,N)
INTEGER N
INTEGER I, J
DO J=2,N
DO I=1,(J-1)
A(I,J) = A(J,I)
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE RAVG (A, NCELLS, NCELL, N)
C=======================================================================
REAL*8 A(NCELLS,NCELLS)
INTEGER NCELLS, NCELL, N
INTEGER IC, JC, NI, NJ
NJ = NCELL
DO JC=1,NCELLS
IF(JC .EQ. NCELLS) NJ = N - NCELL*(JC-1)
NI = NCELL
DO IC=1,NCELLS
IF(IC .EQ. NCELLS) NI = N - NCELL*(IC-1)
A(IC,JC) = A(IC,JC)/DFLOAT(NI*NJ)
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE CAVG (A, NCELLS, NCELL, N)
C=======================================================================
COMPLEX*16 A(NCELLS,NCELLS)
INTEGER NCELLS, NCELL, N
INTEGER IC, JC, NI, NJ
NJ = NCELL
DO JC=1,NCELLS
IF(JC .EQ. NCELLS) NJ = N - NCELL*(JC-1)
NI = NCELL
DO IC=1,NCELLS
IF(IC .EQ. NCELLS) NI = N - NCELL*(IC-1)
A(IC,JC) = A(IC,JC)/DFLOAT(NI*NJ)
END DO
END DO
RETURN
END