c ================================================== c = 三重対角行列 LU分解ルーチン TridiagLU = c = = c ================================================== c ver. 1.00 1997 7/27 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 TridiagLU(alpha, beta, n, eig, + alpha_eig, beta_L, beta_U1, beta_U2, ip, k, errcode) c === 引数宣言 real*8, dimension (0:n) :: alpha, beta integer n real*8 eig real*8, dimension (1:n) :: alpha_eig real*8, dimension (1:n) :: beta_L, beta_U1, beta_U2 integer, dimension (1:n) :: ip integer k, errcode c === その他の変数 integer i integer maxi real*8 dtemp, dmax, pivot real*8 EPSLU c == ランク落ち判定定数 EPSLU = 1.0d-15 c == 初期設定 do i=1, n alpha_eig(i) = alpha(i) - eig beta_L(i) = beta(i) beta_U1(i) = beta(i) beta_U2(i) = 0.0d0 enddo c == LU分解 do i=1, n-2 c == ピポット選択 maxi = i dmax = dabs(alpha_eig(i)) if (dmax .lt. dabs(beta_L(i+1))) then maxi = i + 1 dmax = dabs(beta_L(i+1)) endif ip(i) = maxi c == ピボット交換 if (maxi .ne. i) then dtemp = alpha_eig(i) alpha_eig(i) = beta_L(i+1) beta_L(i+1) = dtemp dtemp = beta_U1(i+1) beta_U1(i+1) = alpha_eig(i+1) alpha_eig(i+1) = dtemp beta_U2(i+2) = beta_U1(i+2) beta_U1(i+2) = 0.0d0 endif if (dmax .lt. EPSLU) then c == ランク落ちの処理 errcode = 4200 write(*,10) k, i 10 format(i5,"th eigenvector lank down as ",i5,"th row ") alpha_eig(i) = 0.0d0 beta_L(i+1) = 1.0d0 / EPSLU alpha_eig(i+1) = alpha_eig(i+1) - beta_L(i+1) * beta_U1(i+1) beta_U2(i) = 0.0d0 else pivot = beta_L(i+1) / alpha_eig(i) beta_L(i+1) = pivot alpha_eig(i+1) = alpha_eig(i+1) - pivot * beta_U1(i+1) beta_U1(i+2) = beta_U1(i+2) - pivot * beta_U2(i+2) endif enddo c == LU分解の終り c == 最終段の処理 c == ピポット選択 maxi = n-1 dmax = dabs(alpha_eig(n-1)) if (dmax .lt. dabs(beta_L(n))) then maxi = n dmax = dabs(beta_L(n)) endif ip(n-1) = maxi ip(n) = n c == ピボット交換 if (maxi .ne. n-1) then dtemp = alpha_eig(n-1) alpha_eig(n-1) = beta_L(n) beta_L(n) = dtemp dtemp = beta_U1(n) beta_U1(n) = alpha_eig(n) alpha_eig(n) = dtemp endif if (dmax .lt. EPSLU) then c == ランク落ちの処理 errcode = 4200 write(*,10) k, n-1 alpha_eig(n-1) = 0.0d0 beta_L(n) = 1.0d0 / EPSLU alpha_eig(n) = alpha_eig(n) - beta_L(n) * beta_U1(n) else pivot = beta_L(n) / alpha_eig(n-1) beta_L(n) = pivot alpha_eig(n) = alpha_eig(n) - pivot * beta_U1(n) endif return end c ================================================== c = 三重対角行列 LU分解ルーチン の終り = c ================================================== c ================================================== c = 三重対角行列 前進、後退代入ルーチン = c = TridiagSub = c ================================================== c ver. 1.00 1997 7/27 c c All Rights Reserved, Copyright (C) 1998, T.Katagiri subroutine TridiagSub(W, alpha_eig, beta_L, beta_U1, + beta_U2, ip, n, p, k) c === 引数宣言 real*8, dimension (1:n, 1:n/p+1) :: W real*8, dimension (1:n) :: alpha_eig real*8, dimension (1:n) :: beta_L, beta_U1, beta_U2 integer, dimension (1:n) :: ip integer n, p, k c === その他の変数 integer i integer ipt real*8 dtemp c === 前進代入 do i=1, n-1 c === ピポット交換の処理 if (i .ne. ip(i)) then ipt = ip(i) dtemp = W(ipt, k) W(ipt, k) = W(i, k) W(i, k) = dtemp endif W(i+1, k) = W(i+1, k) - W(i, k) * beta_L(i+1) enddo c === 後退代入 W(n, k) = W(n, k) / alpha_eig(n) W(n-1, k) = (W(n-1, k) - beta_U1(n) * W(n, k)) / alpha_eig(n-1) do i=n-2, 1, -1 W(i, k) = (W(i, k) - beta_U1(i+1) * W(i+1, k) + - beta_U2(i+2) * W(i+2, k)) / alpha_eig(i) enddo return end c ================================================== c = 三重対角行列 前進、後退代入ルーチンの終り = c ==================================================