!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004. All rights reserved. !--------------------------------------------------------------------- !=begin != Module Linlib ! ! * Developer: SUGIYAMA Ko-ichiro ! * Version: $Id: linlib.f90.g95,v 1.2.2.2 2005/11/14 01:55:12 kitamo Exp $ ! * Tag Name: $Name: arare3m-20051114 $ ! * Change History: ! !== Overview ! !線形計算ライブラリ SSL II と LAPACK の g95 用ラッパー ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! !=end module linlib implicit none private !=begin !== Public Interface public linlib_init !初期化ルーチン public linsolv !実 3 項行列の連立 1 次方程式(倍精度) integer :: LN = 10 !配列のサイズ. 値はダミー character(1) :: LLib !用いるライブラリのフラグ real(8), allocatable :: A(:) ! 係数行列 real(8), allocatable :: B(:) ! 係数行列 real(8), allocatable :: C(:) ! 係数行列 real(8), allocatable :: B2(:) ! 係数行列(LU 分解後) real(8), allocatable :: ss_F1BasicZ(:,:) ! 係数行列の計算に用いる配列 real(8), allocatable :: sf_F2BasicZ(:,:) ! 係数行列の計算に用いる配列 real(8), allocatable :: IPIV(:) ! save LLib, LN, A, B, C, B2, ss_F1BasicZ, sf_F2BasicZ, IPIV !=end contains !=begin !== Procedure Interface ! !=== Initialize module and acquire NAMELIST ! !利用する線形計算パッケージを NAMELIST から設定し, !配列の大きさを決定する ! subroutine linlib_init(cfgfile, N) !=== Dependency use dc_trace, only : BeginSub, EndSub use dc_message, only: MessageNotify use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, & & RegXMax, RegXMin, RegZMin, RegZMax, DelZ use timeset, only: DelTimeShort use basicset, only: ss_VelSoundBasicZ, ss_CpBasicZ, & & ss_DensBasicZ, ss_PotTempBasicZ use arareset, only: beta use ssl2_linear, only: ssl2_ltx_init use lapack_linear, only: lapack_dgtsv_init, lapack_dgttrf use average, only : sf_avr_ss !=== Input integer, intent(in) :: N !配列サイズ character(*), intent(in) :: cfgfile !設定ファイル (NAMELIST) !=== NAMELIST NAMELIST /linlib/ LLib !=end open (10, FILE=cfgfile) read(10, NML=linlib) close(10) call BeginSub("linlib_init", & & fmt="%c", & & c1="Initialize linear algebra library.") !--- 配列サイズを保管 LN = N !--- 配列の割り付け allocate( & & A(RegZMin+1:RegZMax), & & B(RegZMin+2:RegZMax), & & C(RegZMin+1:RegZMax-1), & & B2(RegZMin+3:RegZMax), & & ss_F1BasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & sf_F2BasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & IPIV(RegZMin+1:RegZMax) ) !--- 用いるライブラリの初期化 if (LLib == 's') then call MessageNotify("Message", "linlib.f90", & & "Can't Call SSL II in the calculation with g95") ! call ssl2_ltx_init(LN) elseif (LLib == 'l') then call MessageNotify("Message", "linlib.f90", & & "Call LAPACK in the calculation") call lapack_dgtsv_init(LN) else call MessageNotify("Error", "linlib.f90", & & "Unknown LLib ", c1=LLib) end if !--- 係数行列の計算 ss_F1BasicZ = ( ss_VelSoundBasicZ ** 2.0d0 ) & & / (ss_CpBasicZ * ss_DensBasicZ *(ss_PotTempBasicZ ** 2.0d0)) sf_F2BasicZ = sf_avr_ss( & & ss_CpBasicZ * ss_DensBasicZ * (ss_PotTempBasicZ ** 2.0d0) ) A(RegZMin+2: RegZMax-1) = & & 1.0d0 + (beta ** 2.0d0) & & * ss_F1BasicZ(RegXMax, RegZMin+2: RegZMax-1) & & * ( sf_F2BasicZ(RegXMax, RegZMin+2: RegZMax-1) & & + sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-2) ) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0) A(RegZMin+1) = 1.0d0 & & + (beta ** 2.0d0) & & * ss_F1BasicZ(RegXMax, RegZMin+1) & & * sf_F2BasicZ(RegXMax, RegZMin+1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0) A(RegZMax) = 1.0d0 & & + (beta ** 2.0d0) & & * ss_F1BasicZ(RegXMax, RegZMax) & & * sf_F2BasicZ(RegXMax, RegZMax-1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0) B(RegZMin+2: RegZMax) = & & - ( beta ** 2.0d0 ) & & * ss_F1BasicZ(RegXMax, RegZMin+1: RegZMax-1) & & * sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 ) C(RegZMin+1: RegZMax-1) = & & - ( beta ** 2.0d0 ) & & * ss_F1BasicZ(RegXMax, RegZMin+2: RegZMax) & & * sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 ) !--- 係数行列を LU 分解 (LAPACK) if (LLib == 'l') then call lapack_dgttrf(A, B, C, B2, IPIV) endif call EndSub("linlib_init") end subroutine linlib_init !=begin !== Procedure Interface ! !=== Initialize module ! !利用する線形計算パッケージに配列を渡す ! subroutine LinSolv(D) !=== Dependency use dc_trace, only : BeginSub, EndSub use dc_message, only: MessageNotify ! use ssl2_linear, only: ssl2_ltx use lapack_linear, only: lapack_dgttrs !=== In/Out real(8), intent(inout) :: D(LN) !=end !=== Work real(8) :: AA(LN) real(8) :: BB(LN-1) real(8) :: CC(LN-1) call BeginSub("linsolv", & & fmt="%c", & & c1="call linear algebra library.") !--- 連立 1 次方程式を解く if (LLib == 's') then ! call ssl2_ltx(A, B, C, D) call MessageNotify("Error", "linlib", & & "unknown LLib for g95") elseif (LLib == 'l') then !=== 属性が inout なので, 内部変数を使うことに. AA = A; BB = B; CC = C call lapack_dgttrs(AA, BB, CC, D, B2, IPIV) else call MessageNotify("Error", "linlib", & & "unknown LLib") end if call EndSub("linsolv") end subroutine LinSolv end module linlib