The flow of this EXAMPLE is as follows:
1. Populate matrix [A] with test data using FMSROW.
2. Populate vector {X} with test data using FMSWRT.
3. Factor [A] into [F].
4. Solve [A]{X}={B}.
5. Read [A] and compute average values in coarse windows.
6. Read [F] and compute average values in coarse windows.
7. Check the answer.
The following input may be specified:
- FMS Module
1 = Real Symmetric
2 = Real Nonsymmetric
3 = Complex Hermitian
4 = Complex Symmetric
5 = Complex Nonsymmetric
- Number of equations
- Number of solution vectors
- The cell size to accumulate average values.
For example, if the matrix has rank 1000 and a cell size of 100 is specified, then an array of 10x10 cells will be computed and printed. Each cell will contain the average value of the 100x100 matrix elements it represents. - Any FMS Parameter
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), LUIX(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, 0LUA, 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