c ================================================== c = 全固有値探索ルーチン = c = 並列二分法 PBisec = c ================================================== c ver. 1.00 1997 7/26 c c === 入出力仕様 c [入力変数] c alpha(0:n) : 三重対角行列の対角要素 c beta(0:nx*ncelx+1) : 三重対角行列の副対角要素 c beta2(0:n) : 三重対角行列の対角要素の二乗エリア c xi(0:n), xs(0:n) : 固有値存在範囲の右端と左端エリア c nc(0:n) : 密集固有値数エリア c n : 行列 A の大域的な次元数 c nx : 行列 A の局所的な次元数 c p : 総プロセッサ台数 c ncelx : 二次元ラベル付におけるx次元のプロセッサ総数 c cid : 一次元ラベル付における自プロセッサ番号 c c [出力変数] c beta2(0:n) : 三重対角行列の対角要素の二乗 c xi(0:n), xs(0:n) : 固有値存在範囲の右端と左端 c nc(0:n) : 密集固有値数 c epsbi0 : マシンε c epsbi : 固有値反復打ちきり値 c nstart : 計算担当の固有値番号(開始) c nend : 計算担当の固有値番号(終了) 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 PBisect(alpha, beta, beta2, xi, xs, nc, + nx, n, p, ncelx, cid, epsbi0, epsbi, nstart, nend, NBI) c === 引数宣言 real*8, dimension (0:n) :: alpha real*8, dimension (0:nx*ncelx+1) :: beta real*8, dimension (0:n) :: beta2 real*8, dimension (0:n) :: xi, xs integer, dimension (0:n) :: nc integer nx, n, p, ncelx, cid real*8 epsbi0, epsbi integer nstart, nend integer NBI c == general val. integer i, k real*8 max, eps, w, bisec_x, g real*8 tmpr1 integer nct, nb, np c === 計算機ε epsbi0 = 0.5e0 101 eps = epsbi0 + 1.0d0 if (eps .gt. 1.0e0) then epsbi0 = epsbi0 * 0.5d0 eps = epsbi0 + 1.0d0 goto 101 endif epsbi0 = epsbi0 * 2.0d0 c === betaの二乗 *voption vec, indep do i=1, n beta2(i) = beta(i) * beta(i) enddo c === max = || A ||1 max = dabs(beta(1)) + dabs(alpha(1)) do i=2, n-1 tmpr1 = dabs(beta(i)) + dabs(alpha(i)) + dabs(beta(i+1)) if (tmpr1 .gt. max) then max = tmpr1 endif enddo tmpr1 = dabs(beta(n)) + dabs(alpha(n)) if (tmpr1 .gt. max) max = tmpr1 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 152 endif if (nend .gt. n) then nend = n endif do i=nstart, nend xi(i) = -max xs(i) = max nc(i) = n enddo c === 反復打ち切り値 xs(nstart-1) = 0.0d0 epsbi = max * epsbi0 eps = 3.52e-15 * max c === カウント開始 k = nstart 150 if ((nc(k) .gt. k).or.((k .ne. 1).and.(xi(k) .eq. xs(k-1)))) then w = xs(k) - xi(k) if ((w .ge. epsbi).or.((w .ge. eps).and.(nc(k)-k .gt. 1))) then bisec_x = (xi(k)+xs(k)) / 2.0d0 c === スツルムカウント g = 1.0d0 nct = 0 do i=1, n if (g .ne. 0.0d0) then g = alpha(i) - bisec_x - beta2(i)/g else g = alpha(i) - bisec_x - dabs(beta(i))/epsbi0 endif if (g .lt. 0.0d0) nct = nct + 1 enddo c === スツルムカウントの終り if (nct .le. nend) then c === 固有値数 nct が nend より小さいが等しい (全て担当範囲内) do i=k, nct xs(i) = bisec_x nc(i) = nct enddo do i=nct+1,nend if (bisec_x .ge. xi(i)) then xi(i) = bisec_x endif enddo else c === 固有値数 nct は担当範囲外も含む do i=k, nend xs(i) = bisec_x nc(i) = nct enddo endif goto 150 endif else c === 孤立した点の処理 do nb=1, NBI if ((xs(k)-xi(k)) .gt. epsbi) then bisec_x = (xi(k)+xs(k)) / 2.0d0 c === スツルムカウント g = 1.0d0 nct = 0 do i=1, n if (g .ne. 0.0d0) then g = alpha(i) - bisec_x - beta2(i)/g else g = alpha(i) - bisec_x - dabs(beta(i))/epsbi0 endif if (g .lt. 0.0d0) nct = nct + 1 enddo c === スツルムカウントの終り if (nct .eq. k) then xs(k) = bisec_x else xi(k) = bisec_x endif else goto 151 endif enddo c === 孤立点処理の終り endif 151 continue k = nc(k) + 1 if (k .le. nend) goto 150 152 continue return end c ================================================== c = 全固有値探索ルーチンの終り = c ==================================================