c ================================================== c = 三重対角行列 全固有ベクトル探索ルーチン = c = 並列逆反復法 PTriInvCGS = c ================================================== c ver. 1.00 1997 8/20 直交化版プロトタイプ c 9/18 更新 c 初期ベクトル改良 c 9/19 テスト用プロトタイプ c 10/12 CGS版として認識 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 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) include 'mppf.h' common /MPI2d/ MPP_COMM_X, MPP_COMM_Y, ISTATUS c === 引数宣言 real*8, dimension (1:n, 1:n/p+1) :: W real*8, dimension (0:n) :: alpha real*8, dimension (0:n) :: beta real*8, dimension (1:n) :: alpha_eig real*8, dimension (1:n) :: beta_L, beta_U1, beta_U2 real*8, dimension (0:n) :: xi, xs integer, dimension (1:n) :: ip real*8, dimension (1:n) :: w_work integer n, nstart, nend integer p real*8 epsbi, fact integer cid real*8 gap0 integer gno, errcode c == general val. integer i, k real*8 eig, dtemp real*8 trinorm1, err_norm integer iter, ITER_MAX integer np real*8 norm integer local_k integer color, Gnum1, Gnum1PEno, present_k integer recvPEno, sendPEno c == 直交化スケジュール用 integer, allocatable, dimension(:):: GN integer, allocatable, dimension(:):: GNlocal integer, allocatable, dimension(:):: DoneW real*8, allocatable, dimension(:):: work_a real*4, allocatable, dimension(:):: tmpW c === 固有ベクトル担当範囲外のPEは終了 if (nstart .gt. nend) goto 120 c === 各PE所有ベクトルの上界 np = idceiling(dfloat(n)/dfloat(p)) c === 行列確保 c === 固有ベクトルグループ番号集合 allocate(GN(1:n+1)) c === 固有ベクトルグループ内番号集合 allocate(GNlocal(1:n+1)) c === 処理済フラグ 0:未処理 1:処理済 allocate(DoneW(1:np)) do i=1, np DoneW(i) = 0 enddo allocate(tmpW(1:n)) call RANDOM_SEED call RANDOM_NUMBER(tmpW) c === 直交化ワークエリア allocate(work_a(1:n)) c === 最大反復回数 ITER_MAX = 40 c === Peters, Wilkinsonの直交化範囲 10^-3 * ||T||1 trinorm1 = 0.0d0 do i=1, n trinorm1 = trinorm1 + dabs(alpha(i)) + dabs(beta(i)) enddo gap0 = fact * trinorm1 c === 初期化 k = 1 i = 1 j = 1 GN(1) = 1 GNlocal(1) = 1 GN(n+1) = 0 GNlocal(n+1) = 1 c === Peters , Wlikinson distance 10 continue dtemp = xi(k+1) - xi(k) if (dtemp .le. gap0) then j = j + 1 GN(k+1) = i GNlocal(k+1) = j else i = i + 1 j = 1 GN(k+1) = i GNlocal(k+1) = j endif k = k + 1 if (k .ne. n) goto 10 gno = GN(n) color = 1 c === 逆反復開始 1 continue c === 終了判定 k = nstart 1000 if (DoneW(k-nstart+1) .eq. 0) goto 1002 k = k + 1 if (k .eq. nend+1) goto 102 goto 1000 1002 continue c === 処理対象ベクトルの選択 k = nstart 1010 if ((GNlocal(k).eq.color).and.(DoneW(k-nstart+1).eq.0)) then goto 1012 endif k = k + 1 if (k .eq. nend+1) goto 1011 goto 1010 c === そのカラーが存在しないか全て処理済 1011 color = color + 1 k = nstart goto 1010 1012 continue c === 処理対象固有ベクトル決定 c === 局所位置設定 local_k = k - nstart + 1 do i=1, n W(i, local_k) = 1.0d0 - tmpW(i)*2.0d0 enddo c === 正規化 norm = 0.0d0 do i=1, n norm = norm + W(i, local_k) * W(i, local_k) enddo norm = dsqrt(norm) do i=1, n W(i, local_k) = W(i, local_k) / norm enddo c === 予想固有値 eig = (xi(k) + xs(k)) / 2.0d0 c === 三重対角行列のLU分解 call TridiagLU(alpha, beta, n, eig, + alpha_eig, beta_L, beta_U1, beta_U2, ip, k, errcode) c === 同一シフト点での逆反復 iter = 0 100 continue c === 前進代入、後退代入 call TridiagSub(W, alpha_eig, beta_L, beta_U1, + beta_U2, ip, n, p, local_k) c === Gramm-Schimidt 直交化 if (GNlocal(k) .ne. 1) then c === 直交化が必要 c === 同グループ内番号1番の位置 Gnum1 = k-1 1100 if (GNlocal(Gnum1) .eq. 1) goto 1102 Gnum1 = Gnum1 - 1 goto 1100 1102 continue c === 同グループ内番号1番の位置は自PE内にあるか？ Gnum1PEno = (Gnum1-1) / np if (Gnum1PEno .eq. cid) then c == 自PE内に有る do i=Gnum1, k-1 dtemp = 0.0d0 do j=1, n dtemp = dtemp + W(j,i-nstart+1) * W(j,local_k) enddo do j=1, n W(j,local_k) = W(j,local_k) - dtemp * W(j,i-nstart+1) enddo enddo else c == 自PE内に無い c == 同一グループの各プロセッサ群に転送 sendPEno = -1 i = Gnum1 2100 if ((i-1)/np .ne. sendPEno) then sendPEno = (i-1)/np call MPP_SEND(k, 1, MPP_INTEGER, sendPEno, 1, + MPP_COMM_WORLD, ierr) call MPP_SEND(W(1,local_k), n, MPP_DOUBLE_PRECISION, + sendPEno, 2, MPP_COMM_WORLD, ierr) endif i = i + 1 if (i .eq. nstart) goto 2102 goto 2100 2102 continue c == PE内データの直交化 do i=nstart, k-1 dtemp = 0.0d0 do j=1, n dtemp = dtemp + W(j,i-nstart+1) * W(j,local_k) enddo do j=1, n W(j,local_k) = W(j,local_k) - dtemp * W(j,i-nstart+1) enddo enddo c == 同一グループの各プロセッサ群から直交化済データ受信 recvPEno = -1 i = Gnum1 2200 if ((i-1)/np .ne. recvPEno) then recvPEno = (i-1)/np call MPP_RECV(w_work, n, MPP_DOUBLE_PRECISION, + recvPEno, 2, MPP_COMM_WORLD, ISTATUS, ierr) c == 部分的直交化 do j=1, n W(j,local_k) = W(j,local_k) + w_work(j) enddo endif i = i + 1 if (i .eq. nstart) goto 2202 goto 2200 2202 continue endif endif c === 並列版 Gramm-Schimidt 直交化 の終り c === 正規化 norm = 0.0d0 do i=1, n norm = norm + W(i, local_k) * W(i, local_k) enddo norm = dsqrt(norm) do i=1, n W(i, local_k) = W(i, local_k) / norm enddo c === レーリー商による固有値の改良 do i=1, n w_work(i) = alpha(i) * W(i, local_k) enddo do i=2, n w_work(i-1) = w_work(i-1) + beta(i) * W(i, local_k) enddo do i=2, n w_work(i) = w_work(i) + beta(i) * W(i-1, local_k) enddo eig = 0.0d0 do i=1, n eig = eig + W(i, local_k) * w_work(i) enddo c === 残差ノルムの計算 do i=1, n w_work(i) = alpha(i) * W(i, local_k) enddo do i=2, n w_work(i-1) = w_work(i-1) + beta(i) * W(i, local_k) enddo do i=2, n w_work(i) = w_work(i) + beta(i) * W(i-1, local_k) enddo do i=1, n w_work(i) = w_work(i) - eig * W(i, local_k) enddo err_norm = 0.0d0 do i=1, n err_norm = err_norm + w_work(i) * w_work(i) enddo err_norm = dsqrt(err_norm) c === 反復打ちきり判定 iter = iter + 1 if (iter .gt. ITER_MAX) then errcode = 4000 goto 101 endif if ((err_norm .gt. epsbi).or.(iter .le. 4)) goto 100 101 continue c === 固有値の格納と次の固有値探索へ xi(k) = eig DoneW(k-nstart+1) = 1 c **** 処理終了通告 c == このグループはプロセッサ間でまたいで所有されているか？ if ((GN(nstart-1) .eq. GN(k)).and.(GNlocal(k).ne.1)) then if ((k .eq. nend).or.(GNlocal(k+1).eq.1)) then if (k .eq. nend) then c === 現在処理しているグループ全体の終了 c == PE内処理終了通告 itemp = 0 endif if (GNlocal(k+1).eq.1) then c == 同一グループの各プロセッサ群にグループ処理終了通告 itemp = -1 endif c === 同一グループの各プロセッサ群に放送 sendPEno = -1 i = Gnum1 2300 if ((i-1)/np .ne. sendPEno) then sendPEno = (i-1)/np call MPP_SEND(itemp, 1, MPP_INTEGER, sendPEno, 1, + MPP_COMM_WORLD, ierr) endif i = i + 1 if (i .eq. nstart) goto 2302 goto 2300 2302 continue endif endif c **** 処理終了通告の終り goto 1 c === 最初に戻る 102 continue c === 逆反復の終了 c === 直交化に従事するPEの判定 c == プロセッサ間をまたぐグループが存在する if (GN(nend) .eq. GN(nend+1)) then c === 直交化すべき固有ベクトル番号 present_k = nend + 1 c === 同グループ内番号1番の位置 Gnum1 = nend 1400 if (GNlocal(Gnum1) .eq. 1) goto 1402 Gnum1 = Gnum1 - 1 goto 1400 1402 continue Gnum1PEno = (Gnum1-1) / np c === 直交化要求PEの番号 1300 recvPEno = (present_k-1) / np c === メッセージ受信待機 call MPP_RECV(itemp, 1, MPP_INTEGER, recvPEno, 1, + MPP_COMM_WORLD, ISTATUS, ierr) c === 終了通告か？ if (itemp .eq. -1) then c === 終了通告 goto 110 endif c === 今のPE内データが全て収束したか？ if (itemp .eq. 0) then present_k = present_k + 1 goto 1300 else c === 直交化処理 c == 現在処理中の固有ベクトル番号の変更 present_k = itemp c == データ受信 call MPP_RECV(work_a, n, MPP_DOUBLE_PRECISION, + recvPEno, 2, MPP_COMM_WORLD, ISTATUS, ierr) c == 直交化開始位置判定 k = nstart 1200 if (GN(nend) .eq. GN(k)) goto 1202 k = k + 1 goto 1200 1202 continue c == 直交化 do i=k, nend dtemp = 0.0d0 do j=1, n dtemp = dtemp + W(j,i-nstart+1) * work_a(j) enddo if (i .eq. k) then do j=1, n w_work(j) = - dtemp * W(j,i-nstart+1) enddo else do j=1, n w_work(j) = w_work(j) - dtemp * W(j,i-nstart+1) enddo endif enddo c == 直交化済ベクトルの転送 call MPP_SEND(w_work, n, MPP_DOUBLE_PRECISION, recvPEno, 2, + MPP_COMM_WORLD, ierr) endif c === 直交化処理の終り goto 1300 c === メッセージ受信待ち endif 110 continue c === 直交化に従事するPEの終了 c === 行列解放 deallocate(GN) deallocate(GNlocal) deallocate(DoneW) deallocate(work_a) deallocate(tmpW) 120 continue return end c ================================================== c = 三重対角行列 全固有ベクトル探索ルーチンの終り = c ==================================================