PROGRAM RCSFMS
C-----------------------------------------------------------------------
C This is a sample program that illustrates the following:
C
C * How to use the FMS out-of-core solver in a RCS program,
C
C * How to generate matrix elements efficiently,
C
C * How to use the FMS parallel processing utilities to generate
C matrix elements in parallel.
C
C It is assumed that the model consists of a mesh composed of
C triangles, and the unknowns (equations) are the values on the
C triangle edges. It is also assumed that each triangle impacts
C every other triangle, producing a full matrix.
C
C It is also assumed that the data is double precision complex,
C (COMPLEX*16).
C
C The test matrix produced by this example has a solution vector
C of all (0.5,-0.5)'s.
C
C-----------------------------------------------------------------------
C DECLARATIONS
C-----------------------------------------------------------------------
C
C Program name:
CHARACTER*9 MYNAME
PARAMETER (MYNAME='EXAMPLE_9')
C
C Input data functions:
INTEGER ASK_I
C
C Subroutines called in parallel (passed to FMSPAR) must be
C declared external.
EXTERNAL FILL
C
C FMS file object handles (attribute lists)
C Matrix file:
INTEGER LUA(25)
C Vector file:
INTEGER LUX(25)
C
C Data type:
INTEGER IDTYPE
C Complex*16:
PARAMETER (IDTYPE=2)
C
C FMS profile vector:
C LOWEQ(1)=-1 flags a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C Constants used for this test matrix:
COMPLEX*16 CZERO, CONE, ANSWER
PARAMETER (CZERO = (0.0D0, 0.0D0))
PARAMETER (CONE = (1.0D0, 0.0D0))
PARAMETER (ANSWER = (0.5D0,-0.5D0))
C
C Variables used to check answer:
REAL*8 ERROR, ETEST
C
C Local variables:
INTEGER L_X, LENX
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 (1) Initialize FMS:
C-----------------------------------------------------------------------
CALL FMSINI
C This must be the first statement in your program.
C It performs the following functions:
C * Initializes FMS,
C * Sets FMS Parameters to default values,
C * If children processes start at the beginning of the program,
C they are put to sleep.
C
CALL FMSPSH (MYNAME)
C
C The following calls to subroutine FMSIGT obtain the pointers
C to the FMS memory management arrays IMD, RMD and CMD:
C
C FMS contains a list of environment variables called FMS
C Parameters. These parameters control special features.
C You may modify the default value of the FMS parameters in one
C of 3 ways:
C * Call FMSSET to interactively modify parameters at run time.
C * Include FMSSET, Parameter statements and RETURN in the license
C file to modify the default values,
C * Call FMSIST, FMSRST or FMSCST in the program.
C
C Here are some examples:
C
C If you want to process this matrix with partial pivoting,
C include the following:
C CALL FMSIST ('MFMAT',3)
C
C or
C
C If you want to average the upper and lower triangle parts to
C produce a symmetric matrix, include the following:
C CALL FMSIST ('LUOK', 1)
C NOTE: For symmetric processing, you must also change the
C complex nonsymmetric FMS calls (beginning with CN) to
C complex symmetric (beginning with CS).
C
C Loop back to here to do next problem:
100 CONTINUE
C
C Allow for interactive modification of the FMS parameters:
C If you do not want to do this, remove this call.
PRINT *
PRINT *, 'You may now set any FMS parameter.'
PRINT *, 'When you have finished, type the letters RETURN'
CALL FMSSET
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C You may also obtain the value of any FMS parameter by calling
C FMSIGT, FMSRGT or FMSCGT. Here we get the number of processors.
CALL FMSIGT ('MAXCPU', MAXCPU)
C-----------------------------------------------------------------------
C Application Specific
C-----------------------------------------------------------------------
C
C Generate a model
C ================
C For this sample problem, we will generate a mesh of a cylinder.
C
C Find out how many equations and RHS vectors to use:
PRINT *
PRINT *, 'Enter the model parameters below.'
PRINT *, 'If you want to stop, enter a value of 0'
PRINT *
MAXEQU = ASK_I ('Enter the number of equations')
MAXEQU = MAX0(MAXEQU,6)
IF(MAXEQU .LE. 0) GO TO 200
NUMVEC = ASK_I ('Enter the number of vectors')
IF(NUMVEC .LE. 0) GO TO 200
C
PRINT *
PRINT *,'Start of generating model'
C
C FMS contains memory management routines. Subroutines FMSIMG,
C FMSRMG and FMSCMG allocate INTEGER, REAL*8 and COMPLEX*16
C memory respectively. The arguments include a reference
C array, the returned offset, and the number of words requested.
C Allocate memory for edge list:
CALL FMSIMG (IMD, L_IDEDGE, 3*MAXEQU)
C CALL FMSCST ('SHOW','MEMORY')
C
C Generate the mesh:
NEQEST = MAXEQU
CALL MGEN (NEQEST, NFACES, NEDGES, IMD(L_IDEDGE))
C
PRINT *,'Number of faces....................=', NFACES
PRINT *,'Number of edges....................=', NEDGES
PRINT *,'End of generating model'
C
C For this problem, the number of equations is the number of
C edges.
NUMEQ = NEDGES
C
C Determine order for processing faces
C ====================================
PRINT *,'Start of computing face order'
C
C The amount of storage required to accumulate matrix data before
C transferring to FMS can be reduced by changing the order of
C processing faces. The following code determines the optimum
C order, IORDER(NFACES) and the number of vectors required
C for accumulation, MAXVEC.
C
C Allocate storage for sorting faces:
CALL FMSIMG (IMD, L_IORDER, NFACES)
CALL FMSIMG (IMD, L_MAXUSE, NEDGES)
CALL FMSIMG (IMD, L_NLEFT, NEDGES)
CALL FMSIMG (IMD, L_IUSE , NFACES)
CALL FMSIMG (IMD, L_IAGE , NEDGES)
C
CALL ORDERF (
1 NFACES , NEDGES , IMD(L_IDEDGE),
2 MAXVEC , IMD(L_IORDER), IMD(L_MAXUSE),
3 IMD(L_NLEFT), IMD(L_IUSE ), IMD(L_IAGE ))
C
PRINT *,'Maximum number of active vectors...=', MAXVEC
C
C Return scratch storage:
CALL FMSIMR (IMD, L_IAGE , NEDGES)
CALL FMSIMR (IMD, L_IUSE , NFACES)
PRINT *,'End of computing face order'
C-----------------------------------------------------------------------
C (2) Open FMS files:
C-----------------------------------------------------------------------
CALL CNDI (LOWEQ, NUMEQ, 'Matrix', LUA)
CALL FMSOV (NUMEQ, IDTYPE, NUMVEC, 'Vectors', LUX)
C-----------------------------------------------------------------------
C (3) Write data to FMS files:
C-----------------------------------------------------------------------
PRINT *
PRINT *,'Starting fill'
C
C Allocate storage for computing matrix elements:
C ===============================================
C
C MAXVEC shared vectors for accumulators:
C Increase MAXVEC to allow processes to operate independently:
MAXVEC = MAXVEC + 3*MAXCPU
CALL FMSCMG (CMD, L_VGLOB, MAXVEC*NUMEQ)
C
C Accumulator vector identifiers:
C (Lists which global equations are currently held)
CALL FMSIMG (IMD, L_IDGLOB, MAXVEC)
C
C Process status (private for each process)
CALL FMSIMG (IMD, L_ISTAT, MAXCPU)
C
C Store scalar arguments in shared memory
C =======================================
C This is necessary if children are forked or spawned.
C This is not necessary (but doesn't hurt) if children are threads.
CALL FMSIMG (IMD, L_NXFACE, 5 + 25)
L_NAVAIL = L_NXFACE + 1
L_NFACES = L_NAVAIL + 1
L_NEDGES = L_NFACES + 1
L_MAXVEC = L_NEDGES + 1
L_LUA = L_MAXVEC + 1
IMD(L_NXFACE) = 0
IMD(L_NAVAIL) = MAXVEC
IMD(L_NFACES) = NFACES
IMD(L_NEDGES) = NEDGES
IMD(L_MAXVEC) = MAXVEC
C
C Initialize values:
C ==================
C Initialize accumulators to available:
DO I=1,MAXVEC
IMD(L_IDGLOB + I - 1) = 0
END DO
C
C Initialize NLEFT to MAXUSE:
DO I=1,NEDGES
IMD(L_NLEFT-1+I) = IMD(L_MAXUSE-1+I)
END DO
C
C Initialize addresses of private data:
LL_ISTAT = L_ISTAT
C
C Leave room for local vectors:
CALL FMSIGT ('MAXMD', MAXMD)
LENC16 = 2*(3*NUMEQ)
LENC16 = LENC16 + 2048
MD_NEW = MAXMD - MAXCPU*LENC16
CALL FMSIST ('MAXMD', MD_NEW)
C
C Initialize FMSROW:
C Use all remaining memory for buffers:
C NOTE: This is a common set of buffers for all processes.
C Therefore FMSROW is initialized once by parent before starting
C children.
CALL FMSROW (-1, CMD(L_VGLOB), LUA)
C
C Reset MAXMD to allow space for arrays allocated by children:
CALL FMSIST('MAXMD', MAXMD)
C
C Copy LUA(25) to shared memory for children:
C Note: If children are threads, you may use LUA directly.
DO I = 1,25
IMD(L_LUA+I-1) = LUA(I)
END DO
C
C Generate matrix elements in parallel
C ====================================
C Loop over children processes:
DO ICPU = 2,MAXCPU
C Increment addresses of private data:
LL_ISTAT = LL_ISTAT + 1
C Place work for children in a queue:
CALL FMSPAR (13, FILL,
1 IMD(L_NXFACE ), IMD(L_NAVAIL ),
2 IMD(L_NFACES ), IMD(L_NEDGES ), IMD(L_MAXVEC ),
3 IMD(L_LUA ), IMD(L_IDEDGE ), IMD(L_IORDER ),
4 IMD(L_MAXUSE ), IMD(L_NLEFT ),
5 CMD(L_VGLOB ), IMD(L_IDGLOB ), IMD(LL_ISTAT ))
END DO
C
C Start the children running:
IF(MAXCPU .GT. 1) CALL FMSRUN
C
C Do parent's part:
LL_ISTAT = L_ISTAT
CALL FILL(
1 IMD(L_NXFACE ), IMD(L_NAVAIL ),
2 IMD(L_NFACES ), IMD(L_NEDGES ), IMD(L_MAXVEC ),
3 LUA , IMD(L_IDEDGE ), IMD(L_IORDER ),
4 IMD(L_MAXUSE ), IMD(L_NLEFT ),
5 CMD(L_VGLOB ), IMD(L_IDGLOB ), IMD(LL_ISTAT ))
C
C Wait for the children to complete:
IF(MAXCPU .GT. 1) CALL FMSYNC
C
C Check the status of the processes:
WRITE(6,2002)
NTOTAL = 0
LL_ISTAT = L_ISTAT
DO I=1,MAXCPU
WRITE(6,2004) I, IMD(LL_ISTAT)
IF(IMD(LL_ISTAT) .LT. 0) STOP'Error during fill'
NTOTAL = NTOTAL + IMD(LL_ISTAT)
LL_ISTAT = LL_ISTAT + 1
END DO
WRITE(6,2003) NTOTAL
IF(NTOTAL .NE. NFACES) STOP'Not all faces processed'
C
C End FMSROW:
CALL FMSROW (NUMEQ+1, CMD(L_VGLOB), LUA)
C
C Free storage used for fill
C ==========================
CALL FMSIMR (IMD, L_NXFACE, 5 + 25)
CALL FMSIMR (IMD, L_ISTAT, MAXCPU)
CALL FMSIMR (IMD, L_IDGLOB, MAXVEC)
CALL FMSCMR (CMD, L_VGLOB, MAXVEC*NUMEQ)
CALL FMSIMR (IMD, L_NLEFT, NEDGES)
CALL FMSIMR (IMD, L_MAXUSE, NEDGES)
CALL FMSIMR (IMD, L_IORDER, NFACES)
CALL FMSIMR (IMD, L_IDEDGE, 3*MAXEQU)
PRINT *,'Ending fill'
C
C Write RHS vector
C ================
LENX = LUX(4)/2
CALL FMSCMG (CMD, L_RHS, LENX)
CMD(L_RHS) = CONE
DO I=2,NUMEQ
CMD(L_RHS+I-1) = CZERO
END DO
LOCD = 1
DO IVEC = 1,NUMVEC
CALL FMSWRT (LUX, LOCD, CMD(L_RHS), LUX(4))
LOCD = LOCD + LUX(4)
END DO
CALL FMSCMR (CMD, L_RHS, LENX)
C-----------------------------------------------------------------------
C (4) Perform matrix algebra:
C-----------------------------------------------------------------------
CALL CNDF (LUA, LUA, LUX, LUX, 0)
CALL CNDS (LUA, LUX, LUX, NUMVEC, 0)
C-----------------------------------------------------------------------
C (5) Read data from FMS files:
C-----------------------------------------------------------------------
CALL FMSCMG (CMD, L_X, LENX)
LOCD = 1
ERROR = 0.0D0
DO IVEC = 1,NUMVEC
CALL FMSRED (LUX, LOCD, CMD(L_X), LUX(4))
LOCD = LOCD + LUX(4)
DO I=1,NUMEQ
ETEST = ABS( CMD(L_X + I - 1) - ANSWER )
IF(ETEST .GT. ERROR) ERROR = ETEST
END DO
END DO
PRINT *,'MAXIMUM ERROR =', ERROR
CALL FMSCMR (CMD, L_X, LENX)
C-----------------------------------------------------------------------
C (6) End FMS:
C-----------------------------------------------------------------------
C Do the next problem.
CALL FMSCV (LUX)
CALL FMSCM (LUA)
GO TO 100
200 CONTINUE
CALL FMSPOP (MYNAME)
CALL FMSEND
C-----------------------------------------------------------------------
C FORMAT STATEMENTS
C-----------------------------------------------------------------------
2002 FORMAT (/
1 ' Number of faces computed by each process'/
2 ' PROCESS FACES'/
3 ' ======= =====')
2003 FORMAT (
1 ' -----'/
2 ' TOTAL', I12)
2004 FORMAT (I5,I13)
END
C=======================================================================
SUBROUTINE ORDERF
1 (NFACES, NEDGES, IDEDGE,
2 MAXVEC, IORDER, MAXUSE,
3 NLEFT, IUSE, IAGE)
C=======================================================================
C DESCRIPTION:
C
C This next section of code computes the order of processing
C the faces, IORDER(NFACES) that results in the minimum amount
C of storage required to accumulate the matrix elements.
C NOTE: For this problem, each matrix element has 2 contributions,
C one from each face. Before writing the matrix data to FMS,
C we must compute both contributions to make the term complete.
C We also note that the matrix terms are computed by integrals
C over the triangle faces. Therefore the order of computing the
C matrix terms is determined by the order of processing the
C faces.
C
C This approach uses a reasonable amount of memory to accumulate
C face integrals and then writes the parts of the matrix that
C become complete to FMS. The benefits are:
C * the integrals only need to be computed once
C * there is no rereading of data from disk to do the
C accumulation.
C * the algorithm easily runs in parallel.
C
C The general scheme is:
C COMPLEX*16 Z(3,3)
C COMPLEX*16 VGLOB (NUMEQ,MAXVEC)
C INTEGER IORDER(NFACES)
C Outer loop over faces:
C DO IFACE=1,NFACES
C Change to optimized processing order:
C MYFACE = IORDER(IFACE)
C DO JFACE=1,NFACES
C Compute the integral of MYFACE due to JFACE and store
C the results in VLOCAL.
C END DO
C The vectors in VLOCAL are now complete. They contain the
C contributions to MYFACE from all other faces.
C
C Before the vectors can be written to FMS, we must wait for
C their second contribution (each vector will appear twice,
C unless it is on an edge). We therefore store the vectors
C into VGLOB, until they become complete.
C
C Add VLOCAL TO VGLOB.
C IF(Any of the VLOCAL vectors are complete) THEN
C Write complete VGLOB vector to FMS
C Free VGLOB storage for another vector
C END IF
C END DO
C
C By properly picking the face processing order, we can reduce
C the number of VGLOB vectors to a reasonable level.
C The following code determins this level:
C
C FORMAL PARAMETERS:
C (R ) NFACES = Integer scalar
C Number of faces
C
C (R ) NEDGES = Integer scalar
C Number of edges
C
C (R ) IDEDGE(3,NFACES) = Integer array.
C Lists the edges associated with a face
C
C (RW) IORDER(NFACES) = Integer array.
C On input, specifies the list of starting faces.
C At least one face must be specified in IORDER(1).
C The value following the last face specified must be 0.
C On output, specifies the order for processing faces.
C
C ( W) MAXVEC = Integer array.
C Maximum number of active edges. This determines the
C number of vectors which must be held in memory to
C accumulate results. The objective of this subroutine
C is to select the IORDER array to minimize this value.
C
C ( W) MAXUSE(NEDGES) = Integer array.
C Maximum number of faces using each edge. For interior
C edges, the value will be 2. For edges on the boundary,
C the value will be 1.
C
C ( ) NLEFT(NEDGES) = Integer array.
C Scratch storage used to hold the number of times the
C edge must be found before it can be eliminated.
C This vector is initialized to MAXUSE.
C
C ( ) IUSE(NFACES) = Integer array.
C Scratch storage used to flag available faces.
C
C ( ) IAGE(NEDGES) = Integer array.
C Scratch storage used to hold the age of active edges.
C-----------------------------------------------------------------------
C Formal Parameters
C-----------------------------------------------------------------------
CHARACTER*6 MYNAME
PARAMETER (MYNAME='ORDERF')
INTEGER NFACES
INTEGER NEDGES
INTEGER IDEDGE(3,NFACES)
INTEGER IORDER(NFACES)
INTEGER MAXVEC
INTEGER MAXUSE(NEDGES)
INTEGER IUSE (NFACES)
INTEGER IAGE (NEDGES)
INTEGER NLEFT (NEDGES)
C-----------------------------------------------------------------------
C Local Variables
C-----------------------------------------------------------------------
INTEGER IFACE
INTEGER JFACE
INTEGER IEDGE
INTEGER MYEDGE
INTEGER MAXAGE
INTEGER NUMACT
C-----------------------------------------------------------------------
CALL FMSPSH (MYNAME)
C
C Initialize starting face:
C You should pick the face that is the furthest from the center.
C Here we will do something simple:
IORDER(1) = 1
IORDER(2) = 0
C
C ( 1) Initialize IUSE to unused:
DO IFACE = 1,NFACES
IUSE(IFACE) = 0
END DO
C
C ( 2) Find the maximum number of times each edge is used:
DO IEDGE = 1,NEDGES
MAXUSE(IEDGE) = 0
END DO
DO IFACE = 1,NFACES
DO I = 1,3
IEDGE = IDEDGE(I,IFACE)
MAXUSE(IEDGE) = MAXUSE(IEDGE) + 1
END DO
END DO
C
C ( 3) Initialize the number of times each edge must be found before
C it can be eliminated.
DO IEDGE = 1,NEDGES
NLEFT(IEDGE) = MAXUSE(IEDGE)
END DO
C
C ( 4) Initialize the age of each edge to inactive:
DO IEDGE = 1,NEDGES
IAGE(IEDGE) = 0
END DO
C
C ( 5) Find the face order:
NUMACT = 0
MAXVEC = 0
NSTART = 0
DO IFACE = 1,NFACES
IF(NSTART .EQ. 0) THEN
IF(IORDER(IFACE) .EQ. 0) THEN
C This is the first non-specified face:
NSTART = IFACE - 1
IF(NSTART .EQ. 0)
1 STOP 'ORDERF: No starting faces specified'
ELSE
C This face is specified:
GO TO 100
END IF
END IF
C
C Find the oldest edge.
MAXAGE = -1
MYEDGE = 0
DO IEDGE = 1,NEDGES
IF( IAGE(IEDGE) .GT. MAXAGE ) THEN
MYEDGE = IEDGE
MAXAGE = IAGE(IEDGE)
END IF
END DO
IF(MYEDGE .EQ. 0) STOP 'Oldest edge not found'
C
C Find a face which uses the oldest edge:
DO JFACE = 1,NFACES
IF( IUSE(JFACE) .EQ. 0 ) THEN
DO I= 1,3
IF(IDEDGE(I,JFACE) .EQ. MYEDGE) THEN
IORDER(IFACE) = JFACE
GO TO 100
END IF
END DO
END IF
END DO
STOP 'Face not found which uses oldest edge'
C
C The new face is IORDER(IFACE):
100 CONTINUE
JFACE = IORDER(IFACE)
C PRINT *,'ORDERF: IORDER(',IFACE,') =',JFACE
IUSE(JFACE) = 1
C
C Age all the active edges:
DO IEDGE = 1,NEDGES
IF(IAGE(IEDGE) .GT. 0) THEN
C Edge is active:
IAGE(IEDGE) = IAGE(IEDGE) + 1
END IF
END DO
C
C Start new edges:
DO I = 1,3
IEDGE = IDEDGE(I,JFACE)
IF(IAGE(IEDGE) .EQ. 0) THEN
C New active edge:
IAGE(IEDGE) = 1
NUMACT = NUMACT + 1
END IF
END DO
IF(NUMACT .GT. MAXVEC) MAXVEC = NUMACT
C Eliminate edges which are now complete:
DO I = 1,3
IEDGE = IDEDGE(I,JFACE)
NLEFT(IEDGE) = NLEFT(IEDGE) - 1
IF(NLEFT(IEDGE) .EQ. 0) THEN
C Edge is complete:
IAGE(IEDGE) = 0
NUMACT = NUMACT - 1
END IF
END DO
C
C Do the next face:
END DO
CALL FMSPOP (MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE FILL (
1 NXFACE, NAVAIL,
2 NFACES, NEDGES, MAXVEC,
3 LUA , IDEDGE, IORDER,
4 MAXUSE, NLEFT ,
5 VGLOB , IDGLOB, ISTAT)
C=======================================================================
C DESCRIPTION:
C This subroutine computes the matrix elements.
C It is designed to be run in parallel.
C
C FORMAL PARAMETERS:
C (RW) NXFACE = Integer scalar (shared, read, write)
C Outer face loop counter.
C
C (RW) NAVAIL = Integer scalar (shared, read, write)
C Number of global vectors available.
C (Initialized to MAXVEC)
C
C (R ) NFACES = Integer scalar (shared, read only)
C Number of faces
C
C (R ) NEDGES = Integer scalar (shared, read only)
C Number of edges
C
C (R ) MAXVEC = Integer scalar. (shared, read only)
C Number of vectors in memory
C
C (R ) LUA(25) = Integer array.
C Matrix file attributes.
C
C (R ) IDEDGE(3,NFACES) = Integer array. (shared, read only)
C Lists the edges associated with a face.
C
C (R ) IORDER(NFACES) = Integer array (shared, read only)
C Lists the order for processing faces
C
C (R ) MAXUSE(NEDGES) = Integer array (shared, read only)
C Lists the number of times each edge must be used.
C
C (RW) NLEFT(NEDGES) = Integer array (shared, read, write)
C Scratch storage, initialized to MAXUSE
C
C (RW) VGLOB(NEDGES,MAXVEC) = Complex*16 array (shared, r, w)
C Scratch storage for accumulating vectors
C
C (RW) IDGLOB(MAXVEC) = Integer array (shared, r, w)
C Scratch storage for vector numbers in V.
C
C ( W) ISTAT = Integer scalar (private, write)
C Return status.
C =-1, Error
C >=0, Number of faces processed.
C-----------------------------------------------------------------------
C Formal Parameters
C-----------------------------------------------------------------------
CHARACTER*4 MYNAME
PARAMETER (MYNAME='FILL')
INTEGER NXFACE
INTEGER NAVAIL
INTEGER NFACES
INTEGER NEDGES
INTEGER IDEDGE(3,NFACES)
INTEGER LUA(25)
INTEGER MAXVEC
INTEGER IORDER(NFACES)
INTEGER MAXUSE(NEDGES)
INTEGER NLEFT (NEDGES)
COMPLEX*16 VGLOB(NEDGES, MAXVEC)
INTEGER IDGLOB(MAXVEC)
INTEGER ISTAT
C VOLATILE NXFACE, NAVAIL
C-----------------------------------------------------------------------
C Local Variables
C-----------------------------------------------------------------------
INTEGER IFACE
INTEGER JFACE
INTEGER ITEMP
COMPLEX*16 CMD(0:1)
POINTER (CMD_PTR,CMD)
C-----------------------------------------------------------------------
CALL FMSPSH (MYNAME)
ISTAT = 0
CALL FMSIGT ('MYNODE', MYNODE)
CALL FMSIGT ('MEMPTR', CMD_PTR)
CALL FMSCMG (CMD, L_VLOCAL, 3*NEDGES)
C
C Get your starting face number:
C Get exclusive use:
CALL FMSONE
C NXFACE = NXFACE + 1
C MYFACE = NXFACE
MYFACE = INTINC(NXFACE)
CALL FMSALL
C NOTE: We saved the face number NXFACE in a local variable
C MYFACE, because NXFACE may be changed by another process
C as soon as we called FMSALL.
C
C Outer loop over faces:
100 CONTINUE
IF(MYFACE .GT. NFACES) THEN
C This process is done.
CALL FMSCMR (CMD, L_VLOCAL, 3*NEDGES)
CALL FMSPOP (MYNAME)
RETURN
END IF
IFACE = IORDER(MYFACE)
ISTAT = ISTAT + 1
C
C Inner loop over faces:
C NOTE: I may process these in any order, because I am doing
C all of them and I have room to store the results.
DO JFACE = 1,NFACES
C
C Compute the matrix relating IFACE to JFACE:
C There will need to be a lot more information passed
C through to this routine. The following is used to
C generate a test matrix only.
C
CALL FILL2 (IFACE, JFACE, IDEDGE, MAXUSE, NEDGES,
1 CMD(L_VLOCAL))
END DO
C
C The 3 vectors VLOCAL are now complete. They contain
C 1/2 of the terms for interior edges.
C We now need to accumulate these results into the global
C arrays VGLOB.
C
C Critical Section
C ================
CALL FMSONE
C
C This code writes to shared variables. Therefore it can
C be executed by only one process at a time.
C
C Transfer private vectors to global arrays:
L_V = L_VLOCAL - NEDGES
DO I = 1,3
L_V = L_V + NEDGES
IEQN = IDEDGE(I,IFACE)
C Look for vector IEQN in the accumulation list:
MYVEC = 0
DO IVEC = 1,MAXVEC
IF(IEQN .EQ. IDGLOB(IVEC)) THEN
C Add to existing vector:
MYVEC = IVEC
CALL CVADD (VGLOB(1,MYVEC), CMD(L_V), NEDGES)
GO TO 220
END IF
IF(IDGLOB(IVEC) .EQ. 0) THEN
C This slot is free:
IF(MYVEC .EQ. 0) MYVEC = IVEC
END IF
END DO
IF(MYVEC .EQ. 0) THEN
C Error:
WRITE(6,2000) MYNODE, IDEDGE(I,IFACE),
1 (II, IDGLOB(II),II=1,MAXVEC)
ISTAT = -1
CALL FMSCMR (CMD, L_VLOCAL, 3*NEDGES)
CALL FMSPOP (MYNAME)
RETURN
END IF
C New vector:
ITEMP = INTDEC(NAVAIL)
IDGLOB(MYVEC) = IEQN
CALL CVMOVE (VGLOB(1,MYVEC), CMD(L_V), NEDGES)
220 CONTINUE
NLEFT(IEQN) = NLEFT(IEQN) - 1
IF(NLEFT(IEQN) .EQ. 0) THEN
C Vector is ready to be written out:
CALL FMSROW (IEQN, VGLOB(1,MYVEC), LUA)
C Free the storage:
IDGLOB(MYVEC) = 0
ITEMP = INTINC(NAVAIL)
END IF
END DO
C
C Get your next face number:
C NXFACE = NXFACE + 1
C MYFACE = NXFACE
MYFACE = INTINC(NXFACE)
CALL FMSALL
C
C End of Critical Section
C =======================
C
C Do the next face:
GO TO 100
2000 FORMAT (/
1 ' FILL: No global vectors available for process',I5/
2 ' Trying to store vector',I10/
3 ' Vector Contents'/
4 ' ====== ========'/
5 (I13,I13))
END
C=======================================================================
SUBROUTINE FILL2 (IFACE, JFACE, IDEDGE, MAXUSE, NUMEQ, VLOCAL)
C=======================================================================
C This is a dummy subroutine to generate test matrix elements.
INTEGER IFACE, JFACE, IDEDGE(3,*), NUMEQ, MAXUSE(NUMEQ)
REAL*8 TEMP, DTEMP
COMPLEX*16 ZIJ
COMPLEX*16 VLOCAL(NUMEQ,3)
CHARACTER*5 MYNAME
PARAMETER (MYNAME='FILL2')
CALL FMSPSH (MYNAME)
DO I=1,3
IGLOB = IDEDGE(I,IFACE)
TEMP = (1.0D0)/(DFLOAT(MAXUSE(IGLOB)))
DO J=1,3
JGLOB = IDEDGE(J,JFACE)
IF( IGLOB .EQ. JGLOB ) THEN
C This is a diagonal element:
IF ( IGLOB .EQ. 1) THEN
DTEMP = DFLOAT(NUMEQ)*TEMP
ZIJ = DCMPLX(DTEMP,DTEMP)
ELSE
ZIJ = DCMPLX(TEMP,TEMP)
END IF
ELSE IF( (IGLOB .EQ. 1) .OR. (JGLOB .EQ. 1) ) THEN
C Term is in first row or column:
ZIJ = DCMPLX(-TEMP,-TEMP)
ELSE
C Interior element:
ZIJ = DCMPLX(0.0D0,0.0D0)
END IF
VLOCAL(JGLOB,I) = ZIJ
END DO
END DO
CALL FMSPOP (MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE MGEN (NEQEST, NFACES, NEDGES, IDEDGE)
C=======================================================================
INTEGER NEQEST
INTEGER NFACES
INTEGER NEDGES
INTEGER IDEDGE(3,NEQEST)
REAL*8 RN
CHARACTER*4 MYNAME
PARAMETER (MYNAME='MGEN')
CALL FMSPSH (MYNAME)
C Generate a mesh of quadrilaterals that is
C approximately square. Each quadrilateral will
C be divided into four triangles. Each interior
C quadrilateral will add 6 equations (4 interior
C edges and 2 sides).
C Compute the number of quadrilaterals:
RN = NEQEST/6
C Compute the number of quadrilaterals on an edge:
RN = SQRT(RN)
N = RN
N2 = N + N
N3 = N2 + N
N4 = N3 + N
N6 = N4 + N2
M = (NEQEST + N2)/N6
C PRINT *,'MGEN: Mesh is ',N,' by ',M
NFACES = 0
NEDGES = 0
C Loop down the cylinder:
DO J = 1,M
C Loop around the cylinder:
DO I = 1,N
I2 = I + I
C Generate the 4 face numbers in the square:
IFACE = NFACES + I
JFACE = NFACES + N + I2 - 1
KFACE = JFACE + 1
LFACE = NFACES + N3 + I
C Generate the edge list for each face:
IDEDGE(1,IFACE) = NEDGES + I
IDEDGE(2,IFACE) = NEDGES + N + I2
IDEDGE(3,IFACE) = IDEDGE(2,IFACE) - 1
IDEDGE(1,JFACE) = NEDGES + N3 + I
IDEDGE(2,JFACE) = IDEDGE(3,IFACE)
IDEDGE(3,JFACE) = NEDGES + N4 + I2 - 1
IDEDGE(1,KFACE) = IDEDGE(2,IFACE)
IDEDGE(2,KFACE) = IDEDGE(1,JFACE) + 1
IDEDGE(3,KFACE) = IDEDGE(3,JFACE) + 1
IDEDGE(1,LFACE) = IDEDGE(3,JFACE)
IDEDGE(2,LFACE) = IDEDGE(3,KFACE)
IDEDGE(3,LFACE) = NEDGES + N6 + I
END DO
C Wrap for a cylinder:
IDEDGE(2,KFACE) = IDEDGE(2,KFACE) - N
NFACES = NFACES + N4
NEDGES = NEDGES + N6
END DO
C Leave cylinder open:
NEDGES = NEDGES + N
C
C WRITE(6,5000) (I,(IDEDGE(J,I), J=1,3), I=1,NFACES)
CALL FMSPOP (MYNAME)
RETURN
C5000 FORMAT(/
C 1 ' FACE EDGE 1 EDGE 2 EDGE 3'/
C 2 ' ==== ====== ====== ======'/
C 3 (I5,3I9))
END
C=======================================================================
SUBROUTINE CVMOVE (A, B, N)
C=======================================================================
COMPLEX*16 A(N), B(N)
DO I = 1,N
A(I) = B(I)
END DO
RETURN
END
C=======================================================================
SUBROUTINE CVADD (A, B, N)
C=======================================================================
COMPLEX*16 A(N), B(N)
DO I = 1,N
A(I) = A(I) + B(I)
END DO
RETURN
END
C=======================================================================
INTEGER FUNCTION INTINC (I)
C=======================================================================
C This function increments a volatile shared variable. It is
C placed in a subroutine to prevent some compilers from storing
C the value in a register and not updating it.
INTEGER I
I = I + 1
INTINC = I
RETURN
END
C=======================================================================
INTEGER FUNCTION INTDEC (I)
C=======================================================================
C This function decrements a volatile shared variable. It is
C placed in a subroutine to prevent some compilers from storing
C the value in a register and not updating it.
INTEGER I
I = I - 1
INTDEC = I
RETURN
END
C=======================================================================
INTEGER FUNCTION ASK_I(STRING)
C=======================================================================
CHARACTER* (*) STRING
WRITE(6,2000) STRING
READ (5,*) ASK_I
RETURN
2000 FORMAT (1X,A,'>')
END