c********************************************************************** c Parallel Dense Real Symmetric Eigensolver Test Routine c c 1997 10/30 c Composed by KATAGIRI Takahiro c c ver 1.00 c c Data Distribution Type : (Cyclic, Cyclic) 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 */ include 'mppf.f' include 'Houselib.f' include 'HouseCom.f' program main include 'mppf.h' real*8, allocatable, dimension(:):: eig real*8, allocatable, dimension(:,:):: W integer sw, show integer n, p integer np, nstart, nend character*1 c real*8, allocatable, dimension(:,:):: A real*8, allocatable, dimension(:,:):: orgA real*8 t1, t2, bt, t_all real*8 MPP_WTIME integer ncelx, ncely, nx, ny integer cid, ncidx, ncidy integer, allocatable, dimension(:):: LOC, OWNER real*8 dtemp real*8 EigErr, OrtErr, RemainErr integer i, j integer errcode c === start c === MPI Init. c ===================================================== call MPP_INIT( ierr ) call MPP_COMM_RANK( MPP_COMM_WORLD, cid, ierr ) call MPP_COMM_SIZE( MPP_COMM_WORLD, p, ierr ) c === Check Memory c ===================================================== n=1000 1 allocate(orgA(1:n,1:n),stat=ierr) if (ierr .ne. 0) goto 2 deallocate(orgA) n = n + 100 goto 1 2 n = n - 100 dtemp = float(n) * float(n) * 8.0d0 dtemp = dtemp / 1048576.0d0 if (cid .eq. 0) then print *, "" print *, "=== Parallel Eigensolver Check Program ===" print *, "Parallel Eigensolver on MPI Ver. 1.00 1997/10" print *, "Composed by KATAGIRI Takahiro" print *, "" write(*,3) dtemp 3 format(" Memory ", F7.2, " MByte OK.") print *, "" endif c === Setting Test Condition c ===================================================== if (cid .eq. 0) then write(6,*) "problem size >" read(5,*) n write(6,*) "Inverse Iteration" 10 write(6,*) "1)MGS 2)CGS >" read(5,*) sw if (sw .le. 0) goto 10 if (sw .ge. 3) goto 10 20 write(6, *) "Show Eigenvale&Eigenvector(y/n) >" read(5,*) c if (c .eq. "y") then show = 1 else if (c .eq. "Y") then show = 1 else show = 0 endif endif call MPP_BCAST(n,1,MPP_INTEGER,0,MPP_COMM_WORLD,ierr) call MPP_BCAST(sw,1,MPP_INTEGER,0,MPP_COMM_WORLD,ierr) call MPP_BCAST(show,1,MPP_INTEGER,0,MPP_COMM_WORLD,ierr) c === Set Parallel Control Val. c ===================================================== ncelx = dsqrt(dfloat(p)) ncely = ncelx call cid2xy(cid, ncelx, ncidx, ncidy) nx = n/ncelx+1 ny = n/ncely+1 c === Make Matrix c ===================================================== allocate(LOC(0:n-1),stat=ierr) if (ierr .ne. 0) goto 50 allocate(OWNER(0:n-1),stat=ierr) if (ierr .ne. 0) goto 50 allocate(A(0:nx-1,0:ny-1),stat=ierr) if (ierr .ne. 0) goto 50 goto 51 c === memory overflow 50 continue if (cid .eq. 0) then print *, "Memory Overflow!" print *, "Cannot allocate test matrix" endif stop 51 continue call MakeTestMat(A, nx, ny, n, LOC, OWNER, ncelx, ncely, + ncidx, ncidy) c === Call Eigensolver c ===================================================== allocate(eig(0:n)) allocate(W(1:n, 1:n/p+1)) call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() call PEigenSolver(A, eig, W, n, nx, ny, p, sw, cid, errcode) call MPP_BARRIER(MPP_COMM_WORLD, ierr) t2 = MPP_WTIME() t_all = t2 - t1 call MPP_REDUCE(t_all, bt, 1, MPP_DOUBLE_PRECISION, MPP_MAX, 0, & MPP_COMM_WORLD, ierr) t_all = bt c === Common Routine inline c =========================================================== if (errcode .eq. 0) then if (cid .eq. 0) print *, "正常終了" else if (cid .eq. 0) print *, "異常終了：エラーコード",errcode endif c === Check Eigenvalue Error c ===================================================== call ChkEigErr(eig, n, p, cid, EigErr) c === Check Orthogonality c ===================================================== np = idceiling(dfloat(n)/dfloat(p)) nstart = cid * np + 1 nend = (cid+1) * np if (nstart .gt. n) then nstart = n nend = n - 1 goto 100 endif if (nend .gt. n) then nend = n endif 100 continue call TestVec(W, n, nstart, nend, p, cid, OrtErr) c === Check Max Remainder c ===================================================== allocate(orgA(1:n,1:n)) call MakeMatSeq(orgA, n) call TestSymResumeShort(W, eig, n, nstart, nend, p, cid, + orgA, RemainErr) deallocate(orgA) deallocate(A) c === Test Result Output c ===================================================== if (cid .eq. 0) then print *, "=== Result " print *, "================================================ " print *, "n = ", n, "p = ", p print *, "Test Matrix: Frank " if (sw .eq. 1) print *, "Inverse Iteration : MGS " if (sw .eq. 2) print *, "Inverse Iteration : CGS " if (sw .eq. 3) print *, "Inverse Iteration : Hybrid CGS " print *, "" print *, "Eigenvalue Max Error = ", EigErr print *, "Orthogonality = ", OrtErr print *, "Max Remainder = ", RemainErr write(*, 1000) t_all 1000 format(" Total Calculation Time = ",F10.3," [sec]") print *, "" endif call MPP_BARRIER(MPP_COMM_WORLD, ierr) c === Show Eigenvalue & Eigenvector c ===================================================== if (show .eq. 1) then if (cid .eq. 0) print *, " No. | Eigenvalue | Eigenvector " do k=0, p if (k .eq. cid) then do j = nstart, nend write(*,1010) j, eig(j), (W(i,j-nstart+1),i=1,n) enddo endif call MPP_BARRIER(MPP_COMM_WORLD, ierr) enddo 1010 format(I6," |",F12.8," |",20F12.8) endif deallocate(eig) deallocate(W) c ===== MPI finazize call MPP_FINALIZE(ierr) stop end include 'MakeMat.f' include 'TestSym.f' include 'PEigenSolver.f'