c     ==================================================
c     =    データ再分散ルーチン  RedistCyclic          =
c     =   (Cyclic,Cyclic) -> 全PE 三重対角行列所有     =
c     ==================================================
c     ver. 1.00   1997 7/26
c
c     === 入出力仕様
c
c     [入力変数]
c     A(0:nx-1, 0:ny-1)        : 三重対角行列
c     buf(0:ny-1),buf2(0:ny-1) : ワークエリア
c     alpha(0:n)               : 三重対角行列の対角要素エリア
c     beta(0:nx*ncelx+1)       : 三重対角行列の副対角要素エリア
c     n                        : 行列 A の大域的な次元数 
c     nx, ny                   : 行列 A の局所的な次元数
c     ncelx, ncely             : 二次元ラベル付時の各次元 x, y
c             におけるプロセッサ台数の総数 
c     cid                      : 一次元ラベル付における自プロセッサ番号
c     ncidx, ncidy             : 二次元ラベル付における自プロセッサ番号
c     MPP_DOUBLE_PRECISION,    : MPIのデータタイプ
c 
c     [出力変数]
c     alpha(0:n)               : 三重対角行列の対角要素
c     beta(0:nx*ncelx+1)       : 三重対角行列の副対角要素
c
c   All Rights Reserved, Copyright (C) 1998, T.Katagiri
c/*-
c * Copyright (c) 1998 Information Promotion Agency of Japan and T.Katagiri
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 T.Katagiri
c *
c * 4. The name "Information Promotion Agency" or "T.Katagiri" 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 RedistCyclic(A, buf, buf2, alpha, beta, 
     +       n, nx, ny, ncelx, ncely, ncidx, ncidy)
      include  'mppf.h'
      common  /MPI2d/ MPP_COMM_X, MPP_COMM_Y, ISTATUS

c     === 引数宣言
      real*8,   dimension (0:nx-1, 0:ny-1) :: A
      real*8,   dimension (0:ny-1) :: buf, buf2
      real*8,   dimension (0:n) :: alpha
      real*8,   dimension (0:nx*ncelx+1) :: beta
      integer   n, nx, ny
      integer   ncelx, ncely
      integer   ncidx, ncidy

c     == general val.
      integer   i, j
      
c     == MPI val
      integer   ierr 
      
c       === PE番号が対角成分なら X方向 に alphaを 1対全通信
        if (ncidx .eq. ncidy) then
c         === 対角成分パッケージ化
          do i=0, nx-1
             buf(i) = A(i, i)
          enddo

          call MPP_BCAST(buf, nx, MPP_DOUBLE_PRECISION,
     &                   ncidx, MPP_COMM_Y, ierr)
        else

          call MPP_BCAST(buf, nx, MPP_DOUBLE_PRECISION,
     &                   ncidx, MPP_COMM_Y, ierr)

        endif

c         ==== alphaの収納
        do i=0, nx-1
          alpha(i*ncely+1+ncidx) = buf(i)
        enddo

c       === Y方向に alpha を 全対全通信
        do j=0, ncelx-1
          if (ncidx .eq. j) then
            call MPP_BCAST(buf, nx, MPP_DOUBLE_PRECISION,
     &                     j, MPP_COMM_X, ierr)
c           === alphaの収納 
*voption vec, indep            
              do i=0, nx-1
                alpha(i*ncely+1+j) = buf(i)
              enddo
          else
            call MPP_BCAST(buf2, nx, MPP_DOUBLE_PRECISION,
     &                     j, MPP_COMM_X, ierr)
c           === alphaの収納          
*voption vec, indep    
            do i=0, nx-1
              alpha(i*ncely+1+j) = buf2(i)
            enddo
          endif
        enddo


        beta(1) = 0.0d0
c       === PE番号が (対角+1) mod xcely なら X方向 に betaを 1対全通信
        ipe = mod(ncidx+1, ncely)

        if (ncidy .eq. ipe) then
c         === 対角成分パッケージ化
          if (ncidy .ne. 0) then
            do i=0, nx-1
              buf(i) = A(i, i)
            enddo
          else
            do i=0, nx-2
              buf(i) = A(i, i+1)
            enddo
          endif
 
          call MPP_BCAST(buf, nx, MPP_DOUBLE_PRECISION,
     &                   ipe, MPP_COMM_Y, ierr)
        else
          call MPP_BCAST(buf, nx, MPP_DOUBLE_PRECISION,
     &                   ipe, MPP_COMM_Y, ierr)
        endif
      
  
c         ==== beta の収納
        do i=0, nx-1
          beta(i*ncely+2+ncidx) = buf(i)
        enddo


c       === Y方向に beta を 全対全通信
        do j=0, ncelx-1
          if (ncidx .eq. j) then
            call MPP_BCAST(buf, nx, MPP_DOUBLE_PRECISION,
     &                     j, MPP_COMM_X, ierr)
c           === betaの収納 
*voption vec, indep            
              do i=0, nx-1
                beta(i*ncely+2+j) = buf(i)
              enddo
          else
            call MPP_BCAST(buf2, nx, MPP_DOUBLE_PRECISION,
     &                     j, MPP_COMM_X, ierr)
c           === betaの収納 
*voption vec, indep    
            do i=0, nx-1
              beta(i*ncely+2+j) = buf2(i)
            enddo
           endif
         enddo


      return
      end
c     ==================================================
c     =         データ再分散ルーチンの終り             =
c     ==================================================


c     ==================================================
c     =   固有値収集ルーチン                           =
c     =   GatherAllEig                                 =
c     ==================================================
c     ver. 1.00   1997 8/6
c   All Rights Reserved, Copyright (C) 1998, T.Katagiri


      subroutine GatherAllEig(xi, xs, buf, buf2, n, p, cid)
      include  'mppf.h'
      common  /MPI2d/ MPP_COMM_X, MPP_COMM_Y, ISTATUS

c     === 引数宣言
      real*8,   dimension (0:n) :: xi, xs
      real*8,   dimension (1:n) :: buf, buf2
      integer   n, p, cid

c     === 一般変数
      integer   i, j
      integer   np, llocali, rlocali          
      integer   ierr


c     ==== 動作チェック用 固有値収集
                        
      np = idceiling(dfloat(n)/dfloat(p))

      if (cid .ne. 0) then

          call MPP_SEND(xi(1), n, MPP_DOUBLE_PRECISION, 
     +                0, 0, MPP_COMM_WORLD, ierr)
          call MPP_SEND(xs(1), n, MPP_DOUBLE_PRECISION, 
     +                0, 1, MPP_COMM_WORLD, ierr)
      else 
        do i=1, p-1

          llocali = i * np + 1
          rlocali = (i+1) * np
          if (llocali.gt. n) then
            llocali = n
            rlocali = n - 1
          endif
          if (rlocali .gt. n) then
            rlocali = n
          endif

          call MPP_RECV(buf(1), n, MPP_DOUBLE_PRECISION, 
     +              i, 0, MPP_COMM_WORLD, istatus, ierr)
          call MPP_RECV(buf2(1), n, MPP_DOUBLE_PRECISION, 
     +              i, 1, MPP_COMM_WORLD, istatus, ierr)
          
          do j=llocali, rlocali
            xi(j) = buf(j)
            xs(j) = buf2(j)
          enddo

        enddo
      endif

      call MPP_BCAST(xi, n+1, MPP_DOUBLE_PRECISION, 
     +                0, MPP_COMM_WORLD, ierr)
      call MPP_BCAST(xs, n+1, MPP_DOUBLE_PRECISION, 
     +                0, MPP_COMM_WORLD, ierr)

      return
      end

c     ==================================================
c     =   固有値収集ルーチンの終り                     =
c     ==================================================