c     ==================================================
c     =  固有値誤差解析ルーチン ChkEigErr              =        
c     ==================================================
c     ver. 1.00   1997 7/27
c
c     === 入出力仕様
c
c     [入力変数]
c     xi(0: n)                 :  固有値存在範囲の左端
c     n                        :  行列 A の大域的な次元数 
c     p                        :  総 PE 台数
c     cid                      :  一次元ラベル付時の自PE番号
c     MPP_DOUBLE_PRECISION,    :  MPI用
c     MPP_MAX, MPP_COMM_WORLD 
c         
c     [出力変数]
c         最大の相対誤差を画面に出力  
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 ChkEigErr(xi, n, p, cid, err)
      include  'mppf.h'

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

c     == general val.
      integer  i, k, np
      real*8   rs, tmpr1, max
      integer  ierr


      max = 0.0d0
      np = idceiling(dfloat(n)/dfloat(p))
      do i=1,n
        if ((i-1)/np .eq. cid) then
          k = n - i + 1
          rs = dfloat(2*k-1)/dfloat(2.0*n+1) * 3.14159265d0
          rs = 2.0d0*(1.0d0 - cos(rs))
          rs = 1.0d0 / rs

          tmpr1 = dabs(xi(i))
          tmpr1 = dabs(rs) - tmpr1
          tmpr1 = dabs(tmpr1) / dabs(rs)

          if (tmpr1 .gt. max) then 
            max = tmpr1
          endif

        endif
      enddo

      call MPP_REDUCE(max, tmpr1, 1, MPP_DOUBLE_PRECISION, MPP_MAX, 0, 
     &                MPP_COMM_WORLD, ierr)

      err = tmpr1
      
      return
      end
c     ==================================================
c     =  固有値誤差解析ルーチンの終り                  =        
c     ==================================================

c     ====================================================
c     =  固有ベクトル誤差解析ルーチン                    =
c     =  TestVec                                         =            
c     ====================================================
c     ver. 2.00   1997 8/4
c
c   All Rights Reserved, Copyright (C) 1998, T.Katagiri

      subroutine TestVec(W, n, nstart, nend, p, cid, err)
      include "mppf.h"
      common  /MPI2d/ MPP_COMM_X, MPP_COMM_Y, ISTATUS

c     === 引数宣言
      real*8,   dimension (1:n, 1:n/p+1) :: W
      integer   n
      integer   nstart, nend
      integer   p, cid
      real*8    err

c     == general val.
      integer  i, j, k
      integer  l, m, np
      integer  llocali, rlocali
      real*8   tmpr1
      integer  ierr
      real*8   VTV, VTVall
      real*8,   allocatable, dimension(:,:):: W_buf


      np = idceiling(dfloat(n)/dfloat(p))
      allocate(W_buf(1:n,1:np))

c     === 固有ベクトル直交検定      
      do l=0, p-1 
        if (cid .ne. l) then
c         === 担当外 
          call MPP_SEND(W(1, 1), n*np, MPP_DOUBLE_PRECISION, 
     +                    l, l, MPP_COMM_WORLD, ierr)
        else 
c         === cid = l が今の担当
          VTVall = 0.0d0
          do m=0, p-1
            if (m .ne. cid) then
c             === 受信PEの所有範囲
              llocali = m * np + 1
              rlocali = (m+1) * np
              if (llocali .gt. n) then
                llocali = n
                rlocali = n - 1
              endif
              if (rlocali .gt. n) then
                rlocali = n
              endif
  
              call MPP_RECV(W_buf(1, 1), n*np, MPP_DOUBLE_PRECISION, 
     +              m, l, MPP_COMM_WORLD, ISTATUS, ierr)

              do i=nstart, nend
                do j=llocali, rlocali
                  VTV = 0.0d0
                  do k=1, n 
                    VTV=VTV+W(k,i-nstart+1)*W_buf(k,j-llocali+1)
                  enddo

                  if (i .eq. j) then
                    VTVall = VTVall + (1.0d0 - VTV) * (1.0d0 - VTV)
                  else
                    VTVall = VTVall + VTV*VTV
                  endif

  
                enddo
              enddo      
c             === end of i loop 
            else
c           === 担当PE内にデータがある場合
              do i=nstart, nend
                do j=nstart, nend
                  VTV = 0.0d0
                  do k=1, n 
                    VTV = VTV + W(k, i-nstart+1) * W(k, j-nstart+1)
                  enddo

                  if (i .eq. j) then
                    VTVall = VTVall + (1.0d0 - VTV) * (1.0d0 - VTV)
                  else
                    VTVall = VTVall + VTV*VTV
                  endif

                enddo
              enddo  
            endif 
c           === PE内にデータがあるか？
          enddo
c         === end of m loop
        endif
c       === end of if (cid .ne. l) then 
      enddo
c     === end of l

      deallocate(W_buf)
      call MPP_REDUCE(VTVall, tmpr1, 1, MPP_DOUBLE_PRECISION, MPP_SUM,
     &                0, MPP_COMM_WORLD, ierr)
      tmpr1 = dsqrt(tmpr1)
      err = tmpr1

      return
      end
c     ===================================================
c     = 固有値固有ベクトル誤差解析ルーチン              =        
c     ===================================================


c     ====================================================
c     =  対称密行列用 残差，固有ベクトル精度検定(簡略版) =
c     =  TestSymResumeShort                              =            
c     ====================================================
c     ver. 1.01   1997 10/31 
c
c   All Rights Reserved, Copyright (C) 1998, T.Katagiri


      subroutine TestSymResumeShort(W, xi, n, nstart, nend, p, cid, 
     +          orgA, err)
      include "mppf.h"

c     === 引数宣言
      real*8,   dimension (1:n, 1:n/p+1) :: W
      real*8,   dimension (0:n) :: xi
      integer   n
      integer   nstart, nend
      integer   p, cid
      real*8,   dimension (1:n, 1:n) :: orgA
      real*8    err

c     == general val.
      integer  i, j, k
      integer  np
      integer  ierr
      real*8   max_err_norm, err_norm, norm
      real*8   dtemp
      real*8,  allocatable, dimension(:):: ri 
      real*8,  allocatable, dimension(:,:):: W_buf 
      integer  local_i


      ierr = cid

c     === 元の配列  
      np = idceiling(dfloat(n)/dfloat(p))    
      allocate(W_buf(1:n,1:np))
      allocate(ri(1:n))


c     === 残差ノルムの最大値
      max_err_norm = 0.0d0
      do i=nstart, nend
        local_i = i - nstart + 1
        err_norm = 0.0d0
        do j=1, n
          norm = 0.0d0
          do k=1, n
            norm = norm + orgA(k, j) * W(k, local_i)
          enddo
          dtemp = norm - xi(i) * W(j, local_i)
          ri(j) = dtemp
          err_norm = err_norm + dtemp * dtemp
        enddo
        err_norm = dsqrt(err_norm)

        if (err_norm .gt. max_err_norm) then
          max_err_norm = err_norm
        endif
      enddo
      call MPP_REDUCE(max_err_norm, dtemp, 1, MPP_DOUBLE_PRECISION, 
     &           MPP_MAX, 0, MPP_COMM_WORLD, ierr)
      err = dtemp

c     === 行列の領域解放
      deallocate(ri)
      deallocate(W_buf)

      return
      end
c     =========================================================
c     =  対称密行列 残差，固有ベクトル精度検定(簡略版) の終り =
c     =========================================================