c********************************************************************** c Parallel Dense Real Symmetric Eigensolver c c ・Parallel Householder Tridiagonalization c ・Parallel Bisection Method c ・Parallel Inverse Iteratiom c ・Parallel Householder Inverse Transformation 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 */ subroutine PEigenSolver(A, xi, W, n, nx, ny, p, sw, cid, errcode) include 'mppf.h' common /MPI2d/ MPP_COMM_X, MPP_COMM_Y, ISTATUS real*8, dimension (0:nx-1,0:ny-1) :: A real*8, dimension (1:n, 1:n/p+1) :: W real*8, dimension (0:n) :: xi integer n, p, sw, errcode integer nstart, nend integer ncidx, ncidy, cid integer ncelx, ncely, nx, ny integer xlogp integer ISTATUS(MPP_STATUS_SIZE) c == matrix, vector val real*8, allocatable, dimension(:):: u_x, u_y real*8, allocatable, dimension(:):: x_k, y_k real*8, allocatable, dimension(:):: rbuf real*8, allocatable, dimension(:):: alpha real*8, allocatable, dimension(:):: beta, beta2 real*8, allocatable, dimension(:):: xs integer, allocatable, dimension(:):: nc real*8, allocatable, dimension(:):: buf, buf2 real*8 epsbi0, epsbi integer nstart, nend, NBI real*8, allocatable, dimension(:):: w_work real*8, allocatable, dimension(:):: alpha_eig real*8, allocatable, dimension(:):: beta_L real*8, allocatable, dimension(:):: beta_U1 real*8, allocatable, dimension(:):: beta_U2 integer, allocatable, dimension(:):: ip real*8 fact, gap0 integer gno real*8, allocatable, dimension(:,:):: U real*8, allocatable, dimension(:):: ALP integer xlogp real*8 dxlogp real*8 t1, t2, bt, t_all, t_total c === set errcode errcode = 0 c === Check PE number c ===================================================== dxlogp = dlog2(dsqrt(dfloat(p))) xlogp = idint(dxlogp) if ((dxlogp - dfloat(xlogp)) .ne. 0.0d0) then errcode = 4100 return endif if (p .eq. 1) then errcode = 4400 return endif c === Check matrix size c ===================================================== if (n .lt. p) then errcode = 4300 return endif c === MPI re-init. ncelx = dsqrt(dfloat(p)) ncely = ncelx xlogp = log2(dfloat(ncelx)) call cid2xy(cid, ncelx, ncidx, ncidy) c === MPI SPLIT c ===================================================== call MPP_COMM_SPLIT(MPP_COMM_WORLD, ncidy, ncidx, & MPP_COMM_X, ierr) call MPP_COMM_SPLIT(MPP_COMM_WORLD, ncidx, ncidy, & MPP_COMM_Y, ierr) allocate(U(1:ny,1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(ALP(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(rbuf(0:nx-1),stat=ierr) if (ierr .ne. 0) goto 10 allocate(u_x(0:ny-1),stat=ierr) if (ierr .ne. 0) goto 10 allocate(u_y(0:nx-1),stat=ierr) if (ierr .ne. 0) goto 10 allocate(x_k(0:nx-1),stat=ierr) if (ierr .ne. 0) goto 10 allocate(y_k(0:ny-1),stat=ierr) if (ierr .ne. 0) goto 10 goto 20 c ====== Memory overflow prcess 10 errcode = 1 return 20 continue call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() call PTred(A, u_x, u_y, x_k, y_k, rbuf, n, nx, ny, + ncelx, ncely, xlogp, cid, ncidx, ncidy, + U, ALP) 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 t_total = t_all if (cid .eq. 0) then write(*, 1000) t_all 1000 format(" Tridiagonalization Time = ",F10.3," [sec]") endif deallocate(u_x) deallocate(u_y) deallocate(x_k) deallocate(y_k) deallocate(rbuf) allocate(alpha(0:nx*ncelx+1),stat=ierr) if (ierr .ne. 0) goto 10 allocate(beta(0:nx*ncelx+1),stat=ierr) if (ierr .ne. 0) goto 10 allocate(beta2(0:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(xs(0:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(nc(0:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(alpha_eig(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(beta_L(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(beta_U1(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(beta_U2(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(w_work(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(ip(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(buf(1:n),stat=ierr) if (ierr .ne. 0) goto 10 allocate(buf2(1:n),stat=ierr) if (ierr .ne. 0) goto 10 call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() call RedistCyclic(A, buf, buf2, alpha, beta, + n, nx, ny, ncelx, ncely, ncidx, ncidy) 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 t_total = t_total + t_all if (cid .eq. 0) then write(*, 1001) t_all 1001 format(" Re-distribution Time = ",F10.3," [sec]") endif call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() NBI = 100 fact = 1.0d-3 call PBisect(alpha, beta, beta2, xi, xs, nc, nx, n, p, + ncelx, cid, epsbi0, epsbi, nstart, nend, NBI) 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 t_total = t_total + t_all if (cid .eq. 0) then write(*, 1002) t_all 1002 format(" Calculating Eigenvalues Time = ",F10.3," [sec]") write(*, 1010) t_total 1010 format(" All Eigenvalues Cal. Time = ",F10.3," [sec]") endif call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() call GatherAllEig(xi, xs, buf, buf2, n, p, cid) 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 if (cid .eq. 0) then write(*, 1003) t_all 1003 format(" Gathering Eigenvalues Time = ",F10.3," [sec]") endif call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() if (sw .eq. 1) then call PTriInvMGS(W, alpha, beta, alpha_eig, beta_L, + beta_U1, beta_U2, xi, xs, ip, w_work, n, nstart, nend, + p, epsbi, fact, cid, gap0, gno, errcode) endif if (sw .eq. 2) then call PTriInvCGS(W, alpha, beta, alpha_eig, beta_L, + beta_U1, beta_U2, xi, xs, ip, w_work, n, nstart, nend, + p, epsbi, fact, cid, gap0, gno, errcode) endif if (sw .le. 0) then print *, "Swich err!! sw=",sw endif if (sw .ge. 3) then print *, "Swich err!! sw=",sw endif 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 if (cid .eq. 0) then write(*, 1004) t_all 1004 format(" Calculating Eigenvectors Time = ",F10.3," [sec]") endif call MPP_BARRIER(MPP_COMM_WORLD, ierr) t1 = MPP_WTIME() call PHouseInvTrs(W, U, ALP, buf(1), buf2(1), n, nstart, nend, + ny, p, ncelx, ncidx, cid) 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 if (cid .eq. 0) then write(*, 1005) t_all 1005 format(" Householder Inverse Iteration Time = ",F10.3," [sec]") endif deallocate(alpha) deallocate(beta) deallocate(beta2) deallocate(xs) deallocate(nc) deallocate(alpha_eig) deallocate(beta_L) deallocate(beta_U1) deallocate(beta_U2) deallocate(w_work) deallocate(ip) deallocate(buf) deallocate(buf2) deallocate(U) deallocate(ALP) return end include 'PTred.f' include 'Redist.f' include 'PBisec.f' include 'PTriInvCGS.f' include 'PTriInvMGS.f' include 'TriDiagLU.f' include 'PHouseInvTrs.f'