C*********************************************************************C
C**   Floating-Point Multiple-Precision Arithmetic Routines         **C
C**    for Massively Parallel Processor HITACHI SR2201              **C
C**                                                                 **C
C**   All rights reserved, Copyright (C) 1998, Daisuke TAKAHASHI    **C
C**                                                                 **C
C*********************************************************************C
c/*-
c * Copyright (c) 1998 Information Promotion Agency of Japan and D.Takahashi
c * All rights reserved.
c *
c * Redistribution and non-commercial use in source and binary forms,
c * with or without modification, are permitted provided that the following
c * conditions are met:
c *
c * 1. Redistributions of source code must retain the above copyright
c *    notice, this list of conditions and the following disclaimer.
c *
c * 2. Redistributions in binary form must reproduce the above copyright
c *    notice, this list of conditions and the following disclaimer in the
c *    documentation and/or other materials provided with the distribution.
c *
c * 3. All advertising materials mentioning features or use of this
c *    software must display the following acknowledgement:
c *
c *    This product includes software developed by Information Promotion
c *    Agency of Japan and D.Takahashi
c *
c * 4. The name "Information Promotion Agency" or "D.Takahashi" should not be
c *    used to endorse or promote products derived from this software
c *    without specific prior written permission.
c *
c * 5. Any use of the source code or the binary in a commercial product,
c *    whether may it be the origial representation or in some modified form,
c *    is not permitted without specific prior written permission.
c *
c * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
c * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
c * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
c * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
c * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
c * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
c * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
c * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
c * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
c * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
c * SUCH DAMAGE.
c */
      SUBROUTINE MPINIT(ME2,NPU2)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /CONST/BASE/CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
C
      ME=ME2
      NPU=NPU2
C
      IDCML=4
      BASE=10.0D0**IDCML
      BASE2=BASE**2
      RETURN
      END
      SUBROUTINE MPCLEAR(IA,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPCLEAR'
        END IF
        IERR=3000
        RETURN
      END IF
C
      DO 10 I=1,NN+2
        IA(I)=0
   10 CONTINUE
      IERR=0
      RETURN
      END
      SUBROUTINE MPMOVE(IA,IB,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPMOVE'
        END IF
        IERR=3000
        RETURN
      END IF
C
      DO 10 I=1,NN+2
        IB(I)=IA(I)
   10 CONTINUE
      IERR=0
      RETURN
      END
      SUBROUTINE MPADD(IA,IB,IC,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IC(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPADD'
        END IF
        IERR=3000
        RETURN
      END IF
C
      ISIGNA=IA(1)
      ISIGNB=IB(1)
      CALL MPADD0(IA(2),IB(2),IC(2),ISIGNA,ISIGNB,ISIGNC,N)
      IC(1)=ISIGNC
      IERR=0
      RETURN
      END
      SUBROUTINE MPSUB(IA,IB,IC,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IC(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPSUB'
        END IF
        IERR=3000
        RETURN
      END IF
C
      ISIGNA=IA(1)
      ISIGNB=-IB(1)
      CALL MPADD0(IA(2),IB(2),IC(2),ISIGNA,ISIGNB,ISIGNC,N)
      IC(1)=ISIGNC
      IERR=0
      RETURN
      END
      SUBROUTINE MPADD0(IA,IB,IC,ISIGNA,ISIGNB,ISIGNC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IC(*)
      REAL*8 N
C
      IEXPA=IA(1)
      IEXPB=IB(1)
      IF (ISIGNA*ISIGNB .EQ. 1) THEN
        CALL MPADD1(IA,IB,IC,N)
        ISIGNC=ISIGNA
      ELSE
        IF (IEXPA .GT. IEXPB) THEN
          CALL MPSUB1(IA,IB,IC,N)
          ISIGNC=ISIGNA
        ELSE IF (IEXPA .EQ. IEXPB) THEN
          IFLAG=ICMP2(IA(2),IB(2),N)
          IF (IFLAG .GE. 0) THEN
            CALL MPSUB1(IA,IB,IC,N)
            ISIGNC=ISIGNA
          ELSE
            CALL MPSUB1(IB,IA,IC,N)
            ISIGNC=ISIGNB
          END IF
        ELSE
          CALL MPSUB1(IB,IA,IC,N)
          ISIGNC=ISIGNB
        END IF
      END IF
      RETURN
      END
      SUBROUTINE MPADD1(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(*),IC(*),IW(NDB)
      COMMON /NPROC/NPU
      COMMON /C0ARY2/IW
      REAL*8 ISHIFT,N
C
      NN=IDINT(N/DBLE(NPU))
C
      IEXPA=IA(1)
      IEXPB=IB(1)
      IF (IEXPA .GT. IEXPB) THEN
        ISHIFT=DBLE(IEXPA)-DBLE(IEXPB)
        IF (ISHIFT .GE. N*8.0D0) THEN
          DO 10 I=1,NN+1
            IC(I)=IA(I)
   10     CONTINUE
          RETURN
        ELSE
          DO 20 I=1,NN+1
            IW(I)=IB(I)
   20     CONTINUE
          CALL SHIFTR(IW(2),DINT(ISHIFT/8.0D0),N)
          DO 30 I=1,IDINT(DMOD(ISHIFT,8.0D0))
            CALL DIV10(IW(2),N)
   30     CONTINUE
          CALL ADD(IA(2),IW(2),IC(2),N)
        END IF
        IEXPC=IEXPA
      ELSE IF (IEXPA .EQ. IEXPB) THEN
        CALL ADD(IA(2),IB(2),IC(2),N)
        IEXPC=IEXPA
      ELSE
        ISHIFT=DBLE(IEXPB)-DBLE(IEXPA)
        IF (ISHIFT .GE. N*8.0D0) THEN
          DO 40 I=1,NN+1
            IC(I)=IB(I)
   40     CONTINUE
          RETURN
        ELSE
          DO 50 I=1,NN+1
            IW(I)=IA(I)
   50     CONTINUE
          CALL SHIFTR(IW(2),DINT(ISHIFT/8.0D0),N)
          DO 60 I=1,IDINT(DMOD(ISHIFT,8.0D0))
            CALL DIV10(IW(2),N)
   60     CONTINUE
          CALL ADD(IW(2),IB(2),IC(2),N)
        END IF
        IEXPC=IEXPB
      END IF
      CALL NORMS(IC(2),IEXPC,N)
      IC(1)=IEXPC
      RETURN
      END
      SUBROUTINE MPSUB1(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(*),IC(*),IW(NDB)
      COMMON /C0ARY2/IW
      COMMON /NPROC/NPU
      REAL*8 ISHIFT,N
C
      NN=IDINT(N/DBLE(NPU))
C
      IEXPA=IA(1)
      IEXPB=IB(1)
      IF (IEXPA .GT. IEXPB) THEN
        ISHIFT=DBLE(IEXPA)-DBLE(IEXPB)
        IF (ISHIFT .GE. N*8.0D0) THEN
          DO 10 I=1,NN+1
            IC(I)=IA(I)
   10     CONTINUE
          RETURN
        ELSE
          DO 20 I=1,NN+1
            IW(I)=IB(I)
   20     CONTINUE
          CALL SHIFTR(IW(2),DINT(ISHIFT/8.0D0),N)
          DO 30 I=1,IDINT(DMOD(ISHIFT,8.0D0))
            CALL DIV10(IW(2),N)
   30     CONTINUE
          CALL SUB(IA(2),IW(2),IC(2),N)
        END IF
        IEXPC=IEXPA
      ELSE IF (IEXPA .EQ. IEXPB) THEN
        CALL SUB(IA(2),IB(2),IC(2),N)
        IEXPC=IEXPA
      ELSE
        ISHIFT=IEXPB-IEXPA
        IF (ISHIFT .GE. N*8.0D0) THEN
          DO 40 I=1,NN+1
            IC(I)=IB(I)
   40     CONTINUE
          RETURN
        ELSE
          DO 50 I=1,NN+1
            IW(I)=IA(I)
   50     CONTINUE
          CALL SHIFTR(IW(2),DINT(ISHIFT/8.0D0),N)
          DO 60 I=1,IDINT(DMOD(ISHIFT,8.0D0))
            CALL DIV10(IW(2),N)
   60     CONTINUE
          CALL SUB(IW(2),IB(2),IC(2),N)
        END IF
        IEXPC=IEXPB
      END IF
      CALL NORMS(IC(2),IEXPC,N)
      IC(1)=IEXPC
      RETURN
      END
      SUBROUTINE MPMUL(IA,IB,IC,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IC(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPMUL'
        END IF
        IERR=3000
        RETURN
      END IF
C
      ISIGNA=IA(1)
      ISIGNB=IB(1)
      IEXPA=IA(2)
      IEXPB=IB(2)
      CALL MUL(IA(3),IB(3),IC(3),N)
      ISIGNC=ISIGNA*ISIGNB
      IEXPC=IEXPA+IEXPB
      CALL NORMS(IC(3),IEXPC,N)
      IC(1)=ISIGNC
      IC(2)=IEXPC
      IERR=0
      RETURN
      END
      SUBROUTINE MPINV(IA,IB,IW,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IW(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPINV'
        END IF
        IERR=3000
        RETURN
      END IF
C
      ISIGNA=IA(1)
      IEXPA=IA(2)
      CALL INV(IA(3),IB(3),IW(3),N,IP+3,IERR)
      ISIGNB=ISIGNA
      IEXPB=-IEXPA
      CALL NORMS(IB(3),IEXPB,N)
      IB(1)=ISIGNB
      IB(2)=IEXPB
      IERR=0
      RETURN
      END
      SUBROUTINE MPDIV(IA,IB,IC,IW,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(*),IC(*),IW(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 L,LCMP,N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPDIV'
        END IF
        IERR=3000
        RETURN
      END IF
C
      L=LCMP(IB,0,N)
      IF (L .EQ. N) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*) ' ERROR     THE QUOTIENT IS ZERO IN MPDIV'
        END IF
        IERR=3400
        RETURN
      END IF
C
      CALL MPINV(IB,IC,IW,N,IERR)
      CALL MPMUL(IA,IC,IC,N,IERR)
      IERR=0
      RETURN
      END
      SUBROUTINE MPSQRT(IA,IB,IV,IW,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IV(*),IW(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      IP=IDINT(DLOG10(N)/DLOG10(2.0D0))
C
      IF (2.0D0**IP .NE. N .OR. NN .LE. 8*NPU) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*)
     1    ' ERROR     2.0D0**IP .NE. N .OR. N .LE. 8*NPU*NPU IN MPSQRT'
        END IF
        IERR=3000
        RETURN
      END IF
C
      ISIGNA=IA(1)
      IF (ISIGNA .EQ. -1) THEN
        IF (ME .EQ. 0) THEN
          WRITE(6,*) ' ERROR     NEGATIVE SQUARE ROOT IN MPSQRT'
        END IF
        IERR=3200
        RETURN
      END IF
C
      IEXPA=IA(2)
      IF (MOD(IABS(IEXPA),2) .EQ. 1) THEN
        CALL MULS(IA(3),10,IV(3),N)
        ISIGNV=ISIGNA
        IEXPV=IEXPA-1
      ELSE
        DO 10 I=1,NN+2
          IV(I)=IA(I)
   10   CONTINUE
        ISIGNV=ISIGNA
        IEXPV=IEXPA
      END IF
      CALL SQR(IV(3),IB(3),IW(3),N,IP+3,IERR)
      ISIGNB=ISIGNV
      IEXPB=IEXPV/2
      CALL NORMS(IB(3),IEXPB,N)
      IB(1)=ISIGNB
      IB(2)=IEXPB
      IERR=0
      RETURN
      END
      SUBROUTINE NORMS(IX,IEXP,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION IX(*)
      COMMON /CONST2/BASE2
      REAL*8 L,LCMP,N
C
      IBASE2=IDINT(BASE2)
      IBASE21=IBASE2-1
C
      ITEMP=IX(1)
      CALL MPI_BCAST(ITEMP,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      IF (ITEMP .GE. IBASE2) THEN
        CALL DIV10(IX,N)
        IEXP=IEXP+1
        RETURN
      END IF
C
      L=LCMP(IX,0,N)
      IF (L .GE. 1) THEN
        CALL SHIFTL(IX,L,N)
        IEXP=IDINT(DBLE(IEXP)-L*8.0D0)
      END IF
C
      IF (L .NE. N) THEN
        ITEMP=IX(1)
        CALL MPI_BCAST(ITEMP,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
        ITEMP=IBASE21/ITEMP
        CALL MPI_BCAST(ITEMP,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
        IEXP10=IDINT(DLOG10(DBLE(ITEMP)))
        DO 10 I=1,IEXP10
          CALL MULS(IX,10,IX,N)
   10   CONTINUE
        IEXP=IEXP-IEXP10
      END IF
      RETURN
      END
      SUBROUTINE SHIFTR(IA,ISHIFT,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(NDB)
      INTEGER*4 ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY3/IB
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 ISHIFT,N
C
      NN=IDINT(N/DBLE(NPU))
C
      J=IDINT(DMOD(ISHIFT,DBLE(NPU)))
      CALL MPI_ISEND(IA,NN,MPI_INTEGER,MOD(ME+J,NPU),1,
     1               MPI_COMM_WORLD,IREQ1,IERR)
      CALL MPI_IRECV(IB,NN,MPI_INTEGER,MOD(NPU+ME-J,NPU),1,
     1               MPI_COMM_WORLD,IREQ2,IERR)
      CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
      CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
      IS=IDINT(ISHIFT/DBLE(NPU))
      JS=IS
      IF (ME .LT. J) JS=JS+1
      DO 10 I=1,JS
        IA(I)=0
   10 CONTINUE
      DO 20 I=JS+1,NN
        IA(I)=IB(I-JS)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE SHIFTL(IA,ISHIFT,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(NDB)
      INTEGER*4 ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY3/IB
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 ISHIFT,N
C
      NN=IDINT(N/DBLE(NPU))
C
      J=IDINT(DMOD(ISHIFT,DBLE(NPU)))
      CALL MPI_ISEND(IA,NN,MPI_INTEGER,MOD(NPU+ME-J,NPU),1,
     1               MPI_COMM_WORLD,IREQ1,IERR)
      CALL MPI_IRECV(IB,NN,MPI_INTEGER,MOD(ME+J,NPU),1,
     1               MPI_COMM_WORLD,IREQ2,IERR)
      CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
      CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
      IS=IDINT(ISHIFT/DBLE(NPU))
      JS=IS
      IF (ME .GE. NPU-J) JS=JS+1
      DO 10 I=1,NN-JS
        IA(I)=IB(I+JS)
   10 CONTINUE
      DO 20 I=NN-JS+1,NN
        IA(I)=0
   20 CONTINUE
      RETURN
      END
      SUBROUTINE RANDOM(IA,IB,IS,IT,IV,IW,IX,IY,IZ,N,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IS(*),IT(*),IV(*),IW(*),IX(*),IY(*),IZ(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      M=IDINT(DLOG10(N*8.0D0)/DLOG10(2.0D0))
C
      CALL MPCLEAR(IA,N,IERR)
      CALL MPCLEAR(IB,N,IERR)
      CALL MPCLEAR(IS,N,IERR)
      CALL MPCLEAR(IT,N,IERR)
      CALL MPCLEAR(IV,N,IERR)
      CALL MPCLEAR(IW,N,IERR)
      CALL MPCLEAR(IX,N,IERR)
      CALL MPCLEAR(IY,N,IERR)
      CALL MPCLEAR(IZ,N,IERR)
C
      IA(1)=1
      IA(2)=1
      IF (ME .EQ. 0) THEN
        IA(3)=10000000
      END IF
      IW(1)=1
      IW(2)=0
      IF (ME .EQ. 0) THEN
        IW(3)=50000000
      END IF
      CALL MPSQRT(IW,IB,IV,IZ,N,IERR)
      IT(1)=1
      IT(2)=0
      IF (ME .EQ. 0) THEN
        IT(3)=25000000
      END IF
      IS(1)=1
      IS(2)=0
      IF (ME .EQ. 0) THEN
        IS(3)=50000000
      END IF
      IX(1)=1
      IX(2)=1
      IF (ME .EQ. 0) THEN
        IX(3)=10000000
      END IF
C
      DO 10 I=1,M
        CALL MPMOVE(IA,IY,N,IERR)
        CALL MPADD(IA,IB,IA,N,IERR)
        CALL MPMUL(IA,IS,IA,N,IERR)
        CALL MPMUL(IB,IY,IW,N,IERR)
        CALL MPSQRT(IW,IB,IV,IZ,N,IERR)
        CALL MPSUB(IY,IA,IW,N,IERR)
        CALL MPMUL(IW,IW,IW,N,IERR)
        CALL MPMUL(IX,IW,IW,N,IERR)
        CALL MPSUB(IT,IW,IT,N,IERR)
        CALL MPADD(IX,IX,IX,N,IERR)
   10 CONTINUE
      CALL MPADD(IA,IB,IA,N,IERR)
      CALL MPMUL(IA,IA,IA,N,IERR)
      IS(1)=1
      IS(2)=0
      IF (ME .EQ. 0) THEN
        IS(3)=25000000
      END IF
      CALL MPMUL(IA,IS,IA,N,IERR)
      CALL MPDIV(IA,IT,IB,IZ,N,IERR)
      CALL MPMOVE(IB,IA,N,IERR)
      RETURN
      END
      SUBROUTINE OUTPUT(IA,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(NDB),IC(NDB)
      COMMON /IPROC/ME/NPROC/NPU
      COMMON /C0ARY2/IB/C0ARY3/IC
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
C
      CALL MOV(IA,IB,N)
      CALL ITRANSP(IB,IC,NN,NPU)
      CALL ITRANS(IC,IB,NN/NPU,NPU)
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
      IF (ME .EQ. NPU-1) THEN
        WRITE(6,*) (IB(I),I=MAX0(1,NN-9),NN)
      END IF
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
      RETURN
      END
      SUBROUTINE CLR(IA,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
C
      DO 10 I=1,NN
        IA(I)=0
   10 CONTINUE
      RETURN
      END
      SUBROUTINE MOV(IA,IB,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
C
      DO 10 I=1,NN
        IB(I)=IA(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE ADD(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IC(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
C
      DO 10 I=1,NN
        IC(I)=IA(I)+IB(I)
   10 CONTINUE
C
      CALL NORM0(IC,NN)
      RETURN
      END
      SUBROUTINE SUB(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IC(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
C
      DO 10 I=1,NN
        IC(I)=IA(I)-IB(I)
   10 CONTINUE
C
      CALL NORM1(IC,NN)
      RETURN
      END
      SUBROUTINE SUBSR(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IB(*),IC(*)
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      IBASE2=IDINT(BASE2)
      IBASE21=IBASE2-1
      NN=IDINT(N/DBLE(NPU))
C
      IF (ME .EQ. 0) THEN
        IS=2
        IC(1)=IA*IBASE2-IB(1)-1
      ELSE
        IS=1
      END IF
      IF (ME .EQ. NPU-1) THEN
        IE=NN-1
        IC(NN)=IBASE2-IB(NN)
      ELSE
        IE=NN
      END IF
      DO 10 I=IS,IE
        IC(I)=IBASE21-IB(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE MULS(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IC(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
C
      DO 10 I=1,NN
        IC(I)=IA(I)*IB
   10 CONTINUE
C
      CALL NORM0(IC,NN)
      RETURN
      END
      SUBROUTINE HALF(IA,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      NX=NN/NDA
C
      R=0.0D0
      IF (NX .LE. 1) THEN
        CALL HALFSUB(IA,NN,R)
      ELSE
        NNX=NN/NX
        DO 10 J=1,NX
          NJX=NNX*(J-1)
          CALL HALFSUB(IA(NJX+1),NNX,R)
   10   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE HALFSUB(IA,NN,R)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION B(NDA+1),C(NDA+1),D(NDA+1)
      DIMENSION IA(*)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY1/B,C,D
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      BASE2H=BASE2*0.5D0
      N=DBLE(NN)*DBLE(NPU)
C
      IF (ME .EQ. 0) THEN
        IS=2
        U=DBLE(IA(1))+R*BASE2
        V=DINT(U*0.5D0)
        D(1)=U-V*2.0D0
        IA(1)=IDINT(V)
      ELSE
        IS=1
      END IF
      DO 10 I=IS,NN
        U=DBLE(IA(I))
        V=DINT(U*0.5D0)
        D(I)=U-V*2.0D0
        IA(I)=IDINT(V)
   10 CONTINUE
      CALL MPRANGEC(1.0D0,N-1.0D0,IMIN,IMAX,ME,NPU)
      CALL MPI_ISEND(D,IMAX-IMIN+1,MPI_DOUBLE_PRECISION,
     1               MOD(ME+1,NPU),1,MPI_COMM_WORLD,IREQ1,IERR)
      CALL MPRANGEC(1.0D0,N-1.0D0,IMIN,IMAX,MOD(NPU+ME-1,NPU),NPU)
      CALL MPI_IRECV(C(IS),IMAX-IMIN+1,MPI_DOUBLE_PRECISION,
     1               MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,IREQ2,IERR)
      CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
      CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
      DO 20 I=IS,NN
        IA(I)=IDINT(DBLE(IA(I))+C(I)*BASE2H)
   20 CONTINUE
      R=D(NN)
      CALL MPI_BCAST(R,1,MPI_DOUBLE_PRECISION,NPU-1,
     1               MPI_COMM_WORLD,IERR)
      RETURN
      END
      SUBROUTINE DIV10(IA,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*)
      COMMON /NPROC/NPU
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      NX=NN/NDA
C
      R=0.0D0
      IF (NX .LE. 1) THEN
        CALL DIV10SUB(IA,NN,R)
      ELSE
        NNX=NN/NX
        DO 10 J=1,NX
          NJX=NNX*(J-1)
          CALL DIV10SUB(IA(NJX+1),NNX,R)
   10   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE DIV10SUB(IA,NN,R)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION B(NDA+1),C(NDA+1),D(NDA+1)
      DIMENSION IA(*)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY1/B,C,D
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 N
C
      BASE2D=BASE2/10.0D0
      N=DBLE(NN)*DBLE(NPU)
C
      IF (ME .EQ. 0) THEN
        IS=2
        U=DBLE(IA(1))+R*BASE2
        V=DINT(U/10.0D0)
        D(1)=U-V*10.0D0
        IA(1)=IDINT(V)
      ELSE
        IS=1
      END IF
      DO 10 I=IS,NN
        U=DBLE(IA(I))
        V=DINT(U/10.0D0)
        D(I)=U-V*10.0D0
        IA(I)=IDINT(V)
   10 CONTINUE
      CALL MPRANGEC(1.0D0,N-1.0D0,IMIN,IMAX,ME,NPU)
      CALL MPI_ISEND(D,IMAX-IMIN+1,MPI_DOUBLE_PRECISION,
     1               MOD(ME+1,NPU),1,MPI_COMM_WORLD,IREQ1,IERR)
      CALL MPRANGEC(1.0D0,N-1.0D0,IMIN,IMAX,MOD(NPU+ME-1,NPU),NPU)
      CALL MPI_IRECV(C(IS),IMAX-IMIN+1,MPI_DOUBLE_PRECISION,
     1               MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,IREQ2,IERR)
      CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
      CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
      DO 20 I=IS,NN
        IA(I)=IDINT(DBLE(IA(I))+C(I)*BASE2D)
   20 CONTINUE
      R=D(NN)
      CALL MPI_BCAST(R,1,MPI_DOUBLE_PRECISION,NPU-1,
     1               MPI_COMM_WORLD,IERR)
      RETURN
      END
      FUNCTION ICMP(IA,IB,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION IA(*),IB(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 ICMP,ICMP2,N
C
      NN=IDINT(N/DBLE(NPU))
C
      ICMP=N
      DO 10 I=NN,1,-1
        IF (IA(I) .NE. IB(I)) THEN
          ICMP=DBLE(ME)+DBLE(I-1)*DBLE(NPU)
        END IF
   10 CONTINUE
      CALL MPI_ALLREDUCE(ICMP,ICMP2,1,MPI_DOUBLE_PRECISION,MPI_MIN,
     1                   MPI_COMM_WORLD,IERR)
      ICMP=ICMP2
      RETURN
      END
      FUNCTION ICMP2(IA,IB,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IA(*),IB(*),IW(NDB)
      COMMON /IPROC/ME
      REAL*8 N
C
      CALL SUB(IA,IB,IW,N)
      IF (ME .EQ. 0) THEN
        IF (IW(1) .GE. 0) THEN
          ICMP2=1
        ELSE
          ICMP2=-1
        END IF
      END IF
      CALL MPI_BCAST(ICMP2,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      RETURN
      END
      FUNCTION LCMP(IA,IB,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION IA(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 LCMP,LCMP2,N
C
      NN=IDINT(N/DBLE(NPU))
C
      LCMP=N
      DO 10 I=NN,1,-1
        IF (IA(I) .NE. IB) THEN
          LCMP=DBLE(ME)+DBLE(I-1)*DBLE(NPU)
        END IF
   10 CONTINUE
      CALL MPI_ALLREDUCE(LCMP,LCMP2,1,MPI_DOUBLE_PRECISION,MPI_MIN,
     1                   MPI_COMM_WORLD,IERR)
      LCMP=LCMP2
      RETURN
      END
      FUNCTION DITOR(IA)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION IA(*),IM(2)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME
C
      DITOR=0.0D0
      DO 10 J=1,2
        IF (ME .EQ. J-1) THEN
          CALL MPI_SEND(IA,1,MPI_INTEGER,0,1,MPI_COMM_WORLD,IERR)
        END IF
   10 CONTINUE
      IF (ME .EQ. 0) THEN
        DO 20 J=1,2
          CALL MPI_RECV(IM(J),1,MPI_INTEGER,J-1,1,MPI_COMM_WORLD,
     1                  ISTATUS,IERR)
   20   CONTINUE
      END IF
      CALL MPI_BCAST(IM,2,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      DO 30 J=1,2
        DITOR=DITOR+DBLE(IM(J))/(BASE2**J)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE RTOI(R,IA)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IM(2)
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME
C
      DO 10 J=1,2
        IM(J)=IDINT(R*(BASE2**J))
        R=R-DBLE(IM(J))/(BASE2**J)
   10 CONTINUE
      DO 20 J=1,2
        IF (ME .EQ. J-1) THEN
          IA(1)=IM(J)
        END IF
   20 CONTINUE
      RETURN
      END
      FUNCTION ILOG2(IA)
C
      ILOG2=0
      ITEMP=IA
   10 IF (ITEMP .GE. 2) THEN
        ITEMP=ITEMP/2
        ILOG2=ILOG2+1
        GO TO 10
      END IF
      RETURN
      END
      SUBROUTINE INV(IA,IB,IW,N,M,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IW(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 ICMP,L,N,NN
C
      CALL CLR(IW,N/2)
      R=1.0D0/DITOR(IA)
      CALL RTOI(R,IW)
      NN=DBLE(NPU)*DBLE(NPU)
C
      L=0.0D0
      DO 10 I=1,M-2
        CALL MUL(IW,IA,IB,NN)
        CALL SUBSR(2,IB,IB,NN)
        CALL MUL(IW,IB,IB,NN)
        L=ICMP(IB,IW,NN/2)
        IF (L .GE. N/4-1) GO TO 20
        IF (2.0D0**(I+1) .LT. DBLE(NPU)*DBLE(NPU) .OR. NN .EQ. N/2) THEN
          CALL MOV(IB,IW,NN/2)
        ELSE
          CALL MOV(IB,IW,NN)
          NN=NN*2.0D0
        END IF
   10 CONTINUE
C
      IF (ME .EQ. 0) THEN
        WRITE(6,*) ' ERROR     NOT CONVERGENT IN MPINV'
      END IF
      IERR=3300
      RETURN
C
   20 NN=N
      CALL MOV(IB,IW,NN/2)
      CALL MUL(IW,IA,IB,NN)
      CALL SUBSR(2,IB,IB,NN)
      CALL MUL(IW,IB,IB,NN)
      IERR=0
      RETURN
      END
      SUBROUTINE SQR(IA,IB,IW,N,M,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(*),IB(*),IW(*)
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 ICMP,L,N,NN
C
      CALL CLR(IW,N/2)
      R=1.0D0/DSQRT(DITOR(IA))
      CALL RTOI(R,IW)
      NN=DBLE(NPU)*DBLE(NPU)
C
      L=0.0D0
      DO 10 I=1,M-2
        CALL MUL(IW,IW,IB,NN)
        CALL MUL(IA,IB,IB,NN)
        CALL SUBSR(3,IB,IB,NN)
        CALL MUL(IW,IB,IB,NN)
        CALL HALF(IB,NN)
        L=ICMP(IB,IW,NN/2)
        IF (L .GE. N/4-1) GO TO 20
        IF (2.0D0**(I+1) .LT. DBLE(NPU)*DBLE(NPU) .OR. NN .EQ. N/2) THEN
          CALL MOV(IB,IW,NN/2)
        ELSE
          CALL MOV(IB,IW,NN)
          NN=NN*2.0D0
        END IF
   10 CONTINUE
C
      IF (ME .EQ. 0) THEN
        WRITE(6,*) ' ERROR     NOT CONVERGENT IN MPSQRT'
      END IF
      IERR=3100
      RETURN
C
   20 NN=N
      CALL MOV(IB,IW,NN/2)
      CALL MUL(IW,IW,IB,NN)
      CALL MUL(IA,IB,IB,NN)
      CALL SUBSR(3,IB,IB,NN)
      CALL MUL(IW,IB,IB,NN)
      CALL HALF(IB,NN)
      CALL MUL(IA,IB,IB,NN)
      IERR=0
      RETURN
      END
      SUBROUTINE MUL(IA,IB,IC,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION X(NDA+1),Y(NDA+1),Z(NDA+1)
      DIMENSION IA(*),IB(*),IC(*)
      COMMON /C0ARY1/X,Y,Z
      COMMON /NPROC/NPU
      REAL*8 N,N2
C
      N2=N*4
      NN2=IDINT(N2/DBLE(NPU))
C
      CALL MOVIR(IA,X,NN2)
      CALL RFFT(X,Z,N2,NPU,1)
      CALL MOVIR(IB,Y,NN2)
      CALL RFFT(Y,Z,N2,NPU,1)
      CALL CONV0(X,Y,Z,NN2)
      CALL RFFT(Z,Y,N2,NPU,2)
      CALL ROUND(Z,NN2/2)
      CALL NORM2(Z,IC,NN2/2)
      RETURN
      END
      SUBROUTINE MOVIR(IA,X,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(NN/4,*)
      DIMENSION IA(*)
C
      CALL MOVIR0(IA,X,X(1,3),X(1,4),NN)
      RETURN
      END
      SUBROUTINE MOVIR0(IA,X,C,D,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION X(*),C(*),D(*)
      DIMENSION IA(*)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /CONST/BASE
      COMMON /IPROC/ME/NPROC/NPU
C
      BASEH=BASE*0.5D0
C
      DO 10 I=1,NN/4
        U=DBLE(IA(I))
        V=DINT(U/BASE)
        X(I*2-1)=V
        X(I*2)=U-V*BASE
   10 CONTINUE
      IF (ME .EQ. 0) THEN
        V=DINT(X(2)/BASEH)
        X(1)=X(1)+V
        X(2)=X(2)-V*BASE
      END IF
      CALL MPRANGEC(2.0D0,DBLE(NN/4)*DBLE(NPU),IMIN,IMAX,ME,NPU)
      DO 20 I=IMIN,IMAX
        D(I)=DINT(X(I*2-1)/BASEH)
        U=X(I*2-1)-D(I)*BASE
        V=DINT(X(I*2)/BASEH)
        X(I*2-1)=U+V
        X(I*2)=X(I*2)-V*BASE
   20 CONTINUE
      CALL MPRANGEC(2.0D0,DBLE(NN/4)*DBLE(NPU),IMIN,IMAX,ME,NPU)
      CALL MPI_ISEND(D(IMIN),IMAX-IMIN+1,MPI_DOUBLE_PRECISION,
     1               MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,IREQ1,IERR)
      CALL MPRANGEC(2.0D0,DBLE(NN/4)*DBLE(NPU),IMIN,IMAX,
     1              MOD(ME+1,NPU),NPU)
      CALL MPI_IRECV(C,IMAX-IMIN+1,MPI_DOUBLE_PRECISION,
     1               MOD(ME+1,NPU),1,MPI_COMM_WORLD,IREQ2,IERR)
      CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
      CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
      IF (ME .EQ. NPU-1) THEN
        C(NN/4)=0.0D0
      END IF
      DO 30 I=1,NN/4
        X(I*2)=X(I*2)+C(I)
   30 CONTINUE
      DO 40 I=NN/2+1,NN
        X(I)=0.0D0
   40 CONTINUE
      RETURN
      END
      SUBROUTINE CONV0(X,Y,Z,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(*),Y(*),Z(*)
      COMMON /IPROC/ME
C
      IF (ME .EQ. 0) THEN
        IS=3
        Z(1)=X(1)*Y(1)
        Z(2)=X(2)*Y(2)
      ELSE
        IS=1
      END IF
      DO 10 I=IS,NN-1,2
        U=X(I)
        V=Y(I)
        Z(I)=U*V-X(I+1)*Y(I+1)
        Z(I+1)=U*Y(I+1)+X(I+1)*V
   10 CONTINUE
      RETURN
      END
      SUBROUTINE ROUND(X,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION X(*)
      COMMON /IPROC/ME
C
      IFLAG=0
      QMAX=0.0D0
      DO 10 I=1,NN
        U=X(I)
        X(I)=DNINT(X(I))
        V=U-X(I)
        Q=DABS(V)
        IF (Q .GE. 0.1D0) IFLAG=1
        IF (Q .GT. QMAX) QMAX=Q
   10 CONTINUE
      CALL MPI_ALLREDUCE(QMAX,QMAX2,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                   MPI_COMM_WORLD,IERR)
      QMAX=QMAX2
      IF (IFLAG .EQ. 1) THEN
        IF (QMAX .LT. 0.5D0) THEN
          IF (ME .EQ. 0) THEN
            WRITE(6,*) ' WARNING   RESIDUAL=',QMAX,NN
          END IF
        ELSE
          IF (ME .EQ. 0) THEN
            WRITE(6,*) ' RESIDUAL=',QMAX,NN
          END IF
          STOP ' ERROR     NOT EXACT IN ROUND'
        END IF
      END IF
      RETURN
      END
      SUBROUTINE NORM0(IX,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IX(*)
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
C
      IBASE2=IDINT(BASE2)
C
      NX=NN/(NDA*2)
      IF (NX .LE. 1) THEN
        CALL NORM0SUB(IX,NN)
      ELSE
        NNX=NN/NX
        DO 10 J=NX,2,-1
          NJX=NNX*(J-1)
          CALL NORM0SUB(IX(NJX+1),NNX)
          IF (ME .EQ. 0) THEN
            ICARRY=IX(NJX+1)/IBASE2
            IX(NJX+1)=IX(NJX+1)-ICARRY*IBASE2
          END IF
          CALL MPI_BCAST(ICARRY,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
          IF (ME .EQ. NPU-1) THEN
            IX(NJX)=IX(NJX)+ICARRY
          END IF
   10   CONTINUE
        CALL NORM0SUB(IX,NNX)
      END IF
      RETURN
      END
      SUBROUTINE NORM0SUB(IX,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IX(*),IB(NDA*2+2),IC(NDA*2+2),ID(NDA*2+2)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY1/IB,IC,ID
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 IE,IE2,IS,IS2
C
      IBASE2=IDINT(BASE2)
      IBASE21=IBASE2-1
C
      IE=DBLE(NN)*DBLE(NPU)
   10 IFLAG=0
      IE2=1.0D0
      CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
      DO 20 I=IMIN,IMAX
        IF (IX(I) .GE. IBASE2) THEN
          IFLAG=1
          IE2=DBLE(ME)+DBLE(I-1)*DBLE(NPU)+1.0D0
        END IF
   20 CONTINUE
      CALL MPI_ALLREDUCE(IE2,IE,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                   MPI_COMM_WORLD,IERR)
C
      CALL MPI_ALLREDUCE(IFLAG,IFLAG2,1,MPI_INTEGER,MPI_MAX,
     1                   MPI_COMM_WORLD,IERR)
      IFLAG=IFLAG2
C
      IF (IFLAG .EQ. 1) THEN
        IXMAX=0
        CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
        DO 30 I=IMIN,IMAX
          IF (IX(I) .GT. IXMAX) IXMAX=IX(I)
   30   CONTINUE
        CALL MPI_ALLREDUCE(IXMAX,IXMAX2,1,MPI_INTEGER,MPI_MAX,
     1                     MPI_COMM_WORLD,IERR)
        IXMAX=IXMAX2
        IS=1.0D0
        CALL MPRANGEC(2.0D0,IE-1.0D0,IMIN,IMAX,ME,NPU)
        DO 40 I=IMIN,IMAX
          IF (IX(I) .NE. IBASE21) IS=DBLE(ME)+DBLE(I-1)*DBLE(NPU)+1.0D0
   40   CONTINUE
        CALL MPI_ALLREDUCE(IS,IS2,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                     MPI_COMM_WORLD,IERR)
        IS=IS2
        IF (IE-IS .GE. 2.0D0 .AND. IXMAX .EQ. IBASE2) THEN
          CALL MPRANGEC(IE,IE,IMIN,IMAX,ME,NPU)
          DO 50 I=IMIN,IMAX
            IX(I)=IX(I)-IBASE2
   50     CONTINUE
          CALL MPRANGEC(IS+1.0D0,IE-1.0D0,IMIN,IMAX,ME,NPU)
          DO 60 I=IMIN,IMAX
            IX(I)=0
   60     CONTINUE
          CALL MPRANGEC(IS,IS,IMIN,IMAX,ME,NPU)
          DO 70 I=IMIN,IMAX
            IX(I)=IX(I)+1
   70     CONTINUE
        ELSE
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
          DO 80 I=IMIN,IMAX
            ID(I)=IX(I)/IBASE2
   80     CONTINUE
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
          CALL MPI_ISEND(ID(IMIN),IMAX-IMIN+1,MPI_INTEGER,
     1                   MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,IREQ1,IERR)
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,MOD(ME+1,NPU),NPU)
          CALL MPI_IRECV(IC,IMAX-IMIN+1,MPI_INTEGER,
     1                   MOD(ME+1,NPU),1,MPI_COMM_WORLD,IREQ2,IERR)
          CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
          CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
          CALL MPRANGEC(IE,IE,IMIN,IMAX,ME,NPU)
          DO 90 I=IMIN,IMAX
            IC(I)=0
   90     CONTINUE
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
          DO 100 I=IMIN,IMAX
            IX(I)=IX(I)+IC(I)-ID(I)*IBASE2
  100     CONTINUE
          IF (ME .EQ. 0) THEN
            IX(1)=IX(1)+IC(1)
          END IF
        END IF
        GO TO 10
      END IF
      RETURN
      END
      SUBROUTINE NORM1(IX,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IX(*)
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
C
      IBASE2=IDINT(BASE2)
      IBASE21=IBASE2-1
C
      NX=NN/(NDA*2)
      IF (NX .LE. 1) THEN
        CALL NORM1SUB(IX,NN)
      ELSE
        NNX=NN/NX
        DO 10 J=NX,2,-1
          NJX=NNX*(J-1)
          CALL NORM1SUB(IX(NJX+1),NNX)
          IF (ME .EQ. 0) THEN
            IBORROW=(IBASE21-IX(NJX+1))/IBASE2
            IX(NJX+1)=IX(NJX+1)+IBORROW*IBASE2
          END IF
          CALL MPI_BCAST(IBORROW,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
          IF (ME .EQ. NPU-1) THEN
            IX(NJX)=IX(NJX)-IBORROW
          END IF
   10   CONTINUE
        CALL NORM1SUB(IX,NNX)
      END IF
      RETURN
      END
      SUBROUTINE NORM1SUB(IX,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION IX(*),IB(NDA*2+2),IC(NDA*2+2),ID(NDA*2+2)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY1/IB,IC,ID
      COMMON /CONST2/BASE2
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 IE,IE2,IS,IS2
C
      IBASE2=IDINT(BASE2)
      IBASE21=IBASE2-1
C
      IE=DBLE(NN)*DBLE(NPU)
   10 IFLAG=0
      IE2=1
      CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
      DO 20 I=IMIN,IMAX
        IF (IX(I) .LT. 0) THEN
          IFLAG=1
          IE2=DBLE(ME)+DBLE(I-1)*DBLE(NPU)+1.0D0
        END IF
   20 CONTINUE
      CALL MPI_ALLREDUCE(IE2,IE,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                   MPI_COMM_WORLD,IERR)
C
      CALL MPI_ALLREDUCE(IFLAG,IFLAG2,1,MPI_INTEGER,MPI_MAX,
     1                   MPI_COMM_WORLD,IERR)
      IFLAG=IFLAG2
C
      IF (IFLAG .EQ. 1) THEN
        IXMIN=IBASE2
        CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
        DO 30 I=IMIN,IMAX
          IF (IX(I) .LT. IXMIN) IXMIN=IX(I)
   30   CONTINUE
        CALL MPI_ALLREDUCE(IXMIN,IXMIN2,1,MPI_INTEGER,MPI_MIN,
     1                     MPI_COMM_WORLD,IERR)
        IXMIN=IXMIN2
        IS=1.0D0
        CALL MPRANGEC(2.0D0,IE-1.0D0,IMIN,IMAX,ME,NPU)
        DO 40 I=IMIN,IMAX
          IF (IX(I) .NE. 0) IS=DBLE(ME)+DBLE(I-1)*DBLE(NPU)+1.0D0
   40   CONTINUE
        CALL MPI_ALLREDUCE(IS,IS2,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                     MPI_COMM_WORLD,IERR)
        IS=IS2
        IF (IE-IS .GE. 2.0D0 .AND. IXMIN .EQ. -1) THEN
          CALL MPRANGEC(IE,IE,IMIN,IMAX,ME,NPU)
          DO 50 I=IMIN,IMAX
            IX(I)=IX(I)+IBASE2
   50     CONTINUE
          CALL MPRANGEC(IS+1.0D0,IE-1.0D0,IMIN,IMAX,ME,NPU)
          DO 60 I=IMIN,IMAX
            IX(I)=IBASE21
   60     CONTINUE
          CALL MPRANGEC(IS,IS,IMIN,IMAX,ME,NPU)
          DO 70 I=IMIN,IMAX
            IX(I)=IX(I)-1
   70     CONTINUE
        ELSE
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
          DO 80 I=IMIN,IMAX
            ID(I)=(IBASE21-IX(I))/IBASE2
   80     CONTINUE
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
          CALL MPI_ISEND(ID(IMIN),IMAX-IMIN+1,MPI_INTEGER,
     1                   MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,IREQ1,IERR)
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,MOD(ME+1,NPU),NPU)
          CALL MPI_IRECV(IC,IMAX-IMIN+1,MPI_INTEGER,
     1                   MOD(ME+1,NPU),1,MPI_COMM_WORLD,IREQ2,IERR)
          CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
          CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
          CALL MPRANGEC(IE,IE,IMIN,IMAX,ME,NPU)
          DO 90 I=IMIN,IMAX
            IC(I)=0
   90     CONTINUE
          CALL MPRANGEC(2.0D0,IE,IMIN,IMAX,ME,NPU)
          DO 100 I=IMIN,IMAX
            IX(I)=IX(I)-IC(I)+ID(I)*IBASE2
  100     CONTINUE
          IF (ME .EQ. 0) THEN
            IX(1)=IX(1)-IC(1)
          END IF
        END IF
        GO TO 10
      END IF
      RETURN
      END
      SUBROUTINE NORM2(X,IA,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION X(*),B(NDA+1),C(NDA+1),D(NDA+1)
      DIMENSION IA(*)
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /C0ARY1/B,C,D
      COMMON /CONST/BASE
      COMMON /IPROC/ME/NPROC/NPU
      REAL*8 IE,IE2,IS,IS2
C
      BASE1=BASE-1.0D0
C
      CALL CTRANSP(X,C,NN/2,NPU)
      CALL CTRANS(C,X,NN/2/NPU,NPU)
      X2=X(NN)
      DO 10 I=NN,2,-1
        X(I)=X(I-1)
   10 CONTINUE
      CALL MPI_SEND(X2,1,MPI_DOUBLE_PRECISION,
     1              MOD(ME+1,NPU),1,MPI_COMM_WORLD,IERR)
      CALL MPI_RECV(X3,1,MPI_DOUBLE_PRECISION,
     1              MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,ISTATUS,IERR)
      IF (ME .EQ. 0) THEN
        X(1)=0.0D0
      ELSE
        X(1)=X3
      END IF
C
      IFLAG2=0
      IE=DBLE(NN)*DBLE(NPU)
   20 IFLAG=0
      IF (IFLAG2 .EQ. 0) THEN
        XMIN=BASE
        CALL MPRANGEB(2.0D0,DBLE(NN)*DBLE(NPU),IMIN,IMAX,ME,NN)
        DO 30 I=IMIN,IMAX
          IF (X(I) .LT. XMIN) XMIN=X(I)
   30   CONTINUE
        CALL MPI_ALLREDUCE(XMIN,XMIN2,1,MPI_DOUBLE_PRECISION,MPI_MIN,
     1                     MPI_COMM_WORLD,IERR)
        XMIN=XMIN2
        IF (XMIN .LT. 0.0D0) IFLAG=1
        IF (XMIN .GE. -BASE1) THEN
          CALL MPRANGEB(2.0D0,DBLE(NN)*DBLE(NPU)-1.0D0,IMIN,IMAX,ME,NN)
          DO 40 I=IMIN,IMAX
            X(I)=X(I)+BASE1
   40     CONTINUE
          CALL MPRANGEB(DBLE(NN)*DBLE(NPU),DBLE(NN)*DBLE(NPU),
     1                  IMIN,IMAX,ME,NN)
          DO 50 I=IMIN,IMAX
            X(I)=X(I)+BASE
   50     CONTINUE
          IFLAG2=1
        END IF
      END IF
      IE2=1
      CALL MPRANGEB(2.0D0,IE,IMIN,IMAX,ME,NN)
      DO 60 I=IMIN,IMAX
        IF (X(I) .GE. BASE) THEN
          IFLAG=1
          IE2=DBLE(ME)*DBLE(NN)+DBLE(I)
        END IF
   60 CONTINUE
      IF (IFLAG2 .EQ. 1) THEN
        CALL MPI_ALLREDUCE(IE2,IE,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                     MPI_COMM_WORLD,IERR)
      END IF
C
      CALL MPI_ALLREDUCE(IFLAG,IFLAG3,1,MPI_INTEGER,MPI_MAX,
     1                   MPI_COMM_WORLD,IERR)
      IFLAG=IFLAG3
C
      IF (IFLAG .EQ. 1) THEN
        XMAX=0.0D0
        CALL MPRANGEB(2.0D0,IE,IMIN,IMAX,ME,NN)
        DO 70 I=IMIN,IMAX
          IF (X(I) .GT. XMAX) XMAX=X(I)
   70   CONTINUE
        CALL MPI_ALLREDUCE(XMAX,XMAX2,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                     MPI_COMM_WORLD,IERR)
        XMAX=XMAX2
        IS=1.0D0
        CALL MPRANGEB(2.0D0,IE-1.0D0,IMIN,IMAX,ME,NN)
        DO 80 I=IMIN,IMAX
          IF (X(I) .NE. BASE1) IS=DBLE(ME)*DBLE(NN)+DBLE(I)
   80   CONTINUE
        CALL MPI_ALLREDUCE(IS,IS2,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     1                     MPI_COMM_WORLD,IERR)
        IS=IS2
        IF (IE-IS .GE. 2.0D0 .AND. XMAX .EQ. BASE) THEN
          CALL MPRANGEB(IE,IE,IMIN,IMAX,ME,NN)
          DO 90 I=IMIN,IMAX
            X(I)=X(I)-BASE
   90     CONTINUE
          CALL MPRANGEB(IS+1.0D0,IE-1.0D0,IMIN,IMAX,ME,NN)
          DO 100 I=IMIN,IMAX
            X(I)=0.0D0
  100     CONTINUE
          CALL MPRANGEB(IS,IS,IMIN,IMAX,ME,NN)
          DO 110 I=IMIN,IMAX
            X(I)=X(I)+1.0D0
  110     CONTINUE
        ELSE
          W=X(1)/BASE
          IF (W .LT. 0.0D0) W=W-1.0D0
          CARRY=DINT(W)
          CALL MPI_SEND(CARRY,1,MPI_DOUBLE_PRECISION,
     1                  MOD(NPU+ME-1,NPU),1,MPI_COMM_WORLD,IERR)
          CALL MPI_RECV(CARRY2,1,MPI_DOUBLE_PRECISION,
     1                  MOD(ME+1,NPU),1,MPI_COMM_WORLD,ISTATUS,IERR)
          IF (ME .EQ. NPU-1) THEN
            C(NN)=0.0D0
          ELSE
            C(NN)=CARRY2
          END IF
          DO 120 I=NN,2,-1
            W=X(I)/BASE
            IF (W .LT. 0.0D0) W=W-1.0D0
            C(I-1)=DINT(W)
  120     CONTINUE
          DO 130 I=NN,2,-1
            X(I)=X(I)+C(I)-C(I-1)*BASE
  130     CONTINUE
          IF (ME .EQ. 0) THEN
            X(1)=X(1)+C(1)
          ELSE
            X(1)=X(1)+C(1)-CARRY*BASE
          END IF
        END IF
        GO TO 20
      END IF
C
      IF (IFLAG2 .EQ. 1 .AND. ME .EQ. 0) X(1)=X(1)-1.0D0
C
      CALL CTRANS(X,C,NPU,NN/2/NPU)
      CALL CTRANSP(C,X,NN/2,NPU)
C
      DO 140 I=1,NN/2
        IA(I)=IDINT(X(I*2-1)*BASE+X(I*2))
  140 CONTINUE
      RETURN
      END
      SUBROUTINE MPRANGEB(IS,IE,IMIN,IMAX,ME,NN)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 IE,IS
C
      IMIN=IDINT(DMIN1(DMAX1(1.0D0,IS-DBLE(ME)*DBLE(NN)),
     1                 DBLE(NN)+1.0D0))
      IMAX=IDINT(DMAX1(0.0D0,DMIN1(IE-DBLE(ME)*DBLE(NN),DBLE(NN))))
      RETURN
      END
      SUBROUTINE MPRANGEC(IS,IE,IMIN,IMAX,ME,NPU)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 IE,IS
C
      IMIN=IDINT(((IS-DBLE(ME)-1.0D0)+(DBLE(NPU)-1.0D0))/DBLE(NPU))+1
      IMAX=IDINT((IE-DBLE(ME)-1.0D0+DBLE(NPU))/DBLE(NPU))
      RETURN
      END
      SUBROUTINE RFFT(A,B,N,NPU,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA=1048576,NDB=NDA/4)
      DIMENSION A(2,*),B(2,*),WW(NDA+1),WWW(NDA+1)
      DIMENSION IP(3)
      COMMON /IPROC/ME
      REAL*8 N
C
      NN=IDINT(N/DBLE(NPU))
      N2=IDINT(N/2.0D0)
      IP(1)=ILOG2(N2)
      IP(2)=0
      IP(3)=0
C
      IF (ME .EQ. 0) THEN
        IS=2
      ELSE
        IS=1
      END IF
      IF (IOPT .EQ. 1) THEN
        CALL FFT(A,B,WW,WWW,N2,IP,ME,NPU,1,IERR)
      END IF
      CALL REVERSE(A(1,NN/4+1),B,NN/4,NPU)
      CALL SETTBLR(A(1,NN/4+1),NN,NPU,IOPT)
      IF (IOPT .EQ. 1) THEN
*VOPTION INDEP
        DO 10 I=IS,NN/4
          H1R=A(1,I)+B(1,I)
          H1I=A(2,I)-B(2,I)
          H2R=A(1,I)-B(1,I)
          H2I=A(2,I)+B(2,I)
          WR=A(1,I+NN/4)
          WI=A(2,I+NN/4)
          TEMP=H2I*WI-H2R*WR
          H2R=H2R*WI+H2I*WR
          H2I=TEMP
          A(1,I)=0.5D0*(H1R+H2R)
          A(2,I)=-0.5D0*(H1I+H2I)
          B(1,I)=0.5D0*(H1R-H2R)
          B(2,I)=0.5D0*(H1I-H2I)
   10   CONTINUE
      ELSE
*VOPTION INDEP
        DO 20 I=IS,NN/4
          H1R=A(1,I)+B(1,I)
          H1I=B(2,I)-A(2,I)
          H2R=B(1,I)-A(1,I)
          H2I=A(2,I)+B(2,I)
          WR=A(1,I+NN/4)
          WI=A(2,I+NN/4)
          TEMP=H2I*WI-H2R*WR
          H2R=H2R*WI+H2I*WR
          H2I=TEMP
          A(1,I)=0.5D0*(H1R+H2R)
          A(2,I)=0.5D0*(H1I+H2I)
          B(1,I)=0.5D0*(H1R-H2R)
          B(2,I)=0.5D0*(H2I-H1I)
   20   CONTINUE
      END IF
      CALL REVERSE(B,A(1,NN/4+1),NN/4,NPU)
      IF (IOPT .EQ. 1) THEN
        IF (ME .EQ. 0) THEN
          TEMP=A(1,1)-A(2,1)
          A(1,1)=A(1,1)+A(2,1)
          A(2,1)=TEMP
        END IF
      ELSE
        IF (ME .EQ. 0) THEN
          TEMP=0.5D0*(A(1,1)-A(2,1))
          A(1,1)=0.5D0*(A(1,1)+A(2,1))
          A(2,1)=TEMP
        END IF
        CALL FFT(A,B,WW,WWW,N2,IP,ME,NPU,2,IERR)
      END IF
      RETURN
      END
      SUBROUTINE SETTBLR(W,NN,NPU,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION W(2,*)
      COMMON /IPROC/ME
C
      PI=4.0D0*DATAN(1.0D0)
      PX=PI/(DBLE(NN)*DBLE(NPU)/2.0D0)
      IF (IOPT .EQ. 1) PX=-PX
      DO 10 I=1,NN/4
        W(1,I)=DCOS(PX*(DBLE(I-1)*DBLE(NPU)+DBLE(ME)))
        W(2,I)=DSIN(PX*(DBLE(I-1)*DBLE(NPU)+DBLE(ME)))
   10 CONTINUE
      RETURN
      END
      SUBROUTINE REVERSE(A,B,NN,NPU)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      COMPLEX*16 A(*),B(*),TEMP
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      COMMON /IPROC/ME
C
      IF (ME .EQ. 0) THEN
        B(1)=A(1)
        B(NN/2+1)=A(NN/2+1)
*VOPTION INDEP
        DO 10 I=2,NN/2
          B(I)=A(NN-I+2)
          B(NN-I+2)=A(I)
   10   CONTINUE
      ELSE IF (ME .EQ. NPU/2) THEN
*VOPTION INDEP
        DO 20 I=1,NN/2
          B(I)=A(NN-I+1)
          B(NN-I+1)=A(I)
   20   CONTINUE
      ELSE
*VOPTION INDEP
        DO 30 I=1,NN/2
          TEMP=A(I)
          A(I)=A(NN-I+1)
          A(NN-I+1)=TEMP
   30   CONTINUE
        CALL MPI_ISEND(A,NN,MPI_DOUBLE_COMPLEX,
     1                 NPU-ME,1,MPI_COMM_WORLD,IREQ1,IERR)
        CALL MPI_IRECV(B,NN,MPI_DOUBLE_COMPLEX,
     1                 NPU-ME,1,MPI_COMM_WORLD,IREQ2,IERR)
        CALL MPI_WAIT(IREQ1,ISTATUS,IERR)
        CALL MPI_WAIT(IREQ2,ISTATUS,IERR)
      END IF
      RETURN
      END
      SUBROUTINE ITRANS(IA,IB,N1,N2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IA(N1,*),IB(N2,*)
C
      DO 20 J=1,N1
        DO 10 I=1,N2
          IB(I,J)=IA(J,I)
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ITRANSP(IA,IB,NN,NPU)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      DIMENSION IA(*),IB(*)
C
      NN2=NN/NPU
      CALL MPI_ALLTOALL(IA,NN2,MPI_INTEGER,IB,NN2,MPI_INTEGER,
     1                  MPI_COMM_WORLD,IERR)
      RETURN
      END
      SUBROUTINE CTRANS(A,B,N1,N2)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16 A(N1,*),B(N2,*)
C
      DO 20 J=1,N1
        DO 10 I=1,N2
          B(I,J)=A(J,I)
   10   CONTINUE
   20 CONTINUE
      RETURN
      END