C********************************************************************** C leapflogmseq.f C Parallel Interger Random Number Generator C C 1998 01/01 C Composed by Hiroaki Watanabec C ver 1.00 C C*********************************************************************** C All Rights Reserved, Copyright (C) 1998, H.Watanabe c/*- c * Copyright (c) 1998 Information Promotion Agency of Japan and H.Watanabe 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 H.Watanabe c * c * 4. The name "Information Promotion Agency" or "H.Watanabe" 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 */ C -------------------------------------------------- C -------------------------------------------------- C 並列動作するプロセッサ数のチェック PE数は2の巾乗のみ C -------------------------------------------------- C -------------------------------------------------- Subroutine Check_PE_Size(iNumPE, iErr) Integer iNumPE Integer iErr Integer i, iTmpPE iErr = 1300 iTmpPE = 1 Do i = 1, 10 If (iNumPE .EQ. iTmpPE) Then iErr = 0 goto 9999 End If iTmpPE = iTmpPE * 2 End Do 9999 Continue End Subroutine C -------------------------------------------------- C 円周率の値を利用して、初期値を設計する。 C -------------------------------------------------- Subroutine InitialSeedVector3(iSeedvec, iPos, iErr) Parameter (iP = 250) Parameter (iQ = 103) Parameter (IARRAYSIZE = 131072) Integer iSeedvec(iP) Integer i, n, iPos, iIndex, iLocation Integer iTmp, iUnit Integer iTmp1, iTmp2, iTmp3, iTmp4 Integer Mask1, Mask2, Mask3, Mask4 Integer, Allocatable, Dimension(:) :: iArray Data Mask1, Mask2, Mask3, Mask4 1 /Z'1fffffff', Z'7fffffff',Z'001fffff',Z'00000003'/ C 初期化済チェック用フラグ Integer iFlag, iIOS Save iFlag C -------------------------------------------------- C 初期化済(iFlag == 1)の際は何もしないで終了する If (iFlag .EQ. 1) Then goto 9999 End If C -------------------------------------------------- C 小数点以下の桁数 = iSeed は 10^6 以下の整数でなければならない。 C それ以外の値の場合は 1 を初期値とする If (iSeed .LT. 1) Then C 種が0以下ならデフォルトの初期値を使う。 iSeed = 1 Else If (iSeed .GT. 1000000) Then C 種が 10^6 以上ならデフォルトの初期値を使う。 iSeed = 1 End If C -------------------------------------------------- C Step 1. C はじめの 9 個を円周率の小数点以下 第 iSeed 位の値から C 構成する。 C このサブプログラムは, 円周率の値を用いて初期化する。 C 1. 1 配列の割り付け Allocate(iArray(1:IARRAYSIZE), STAT=iErr) If (iErr .NE. 0) Then iErr = 1400 goto 9999 End If C 1. 2 配列サイズ、入出力装置番号の設定 n = IARRAYSIZE iUnit = 2 C 1. 3 ファイルのオープン、読み込み、クローズ Open (iUnit, FILE='pi.dat', STATUS='OLD', IOSTAT=iIOS, ERR=500) Read (iUnit, FMT=2000, IOSTAT=iIOS, ERR=500) (iArray(i), i = 1, n) 500 Continue If (iIOS .EQ. 0) Then Goto 600 Else C ファイル pi.dat が無い場合のエラー処理 iErr = 1500 Goto 9999 End If 600 Continue Close (iUnit) C -------------------------------------------------- C 1. 4 初期値ベクトルの作成 C 元データは10進8桁で詰められているので、それを 2^31 以内の値に設定し直す C 一番先頭の値を設定。 iIndex = iPos / 8 + 1 iLocation = MOD(iPos - 1, 8) iTmp = MOD(iArray(iIndex), 10**(8 - iLocation)) C iTmp(上位部分) が 2^31/10 を超えていないなら次の上位 2 桁を利用する If (iTmp .LT. 21474836) Then iSeedvec(1) = iTmp * 100 + iArray(iIndex + 1) / 1000000 C iTmp(上位部分) が 2^31/100 を超えているなら次の上位 1 桁を利用する Else iSeedvec(1) = iTmp * 10 + iArray(iIndex + 1) / 100000000 End If C 2 番から 9番目までの値を設定。2ワードの元データを一組にして、 C 初期値配列を生成する iIndex = iIndex + 2 Do i = 2, 9 If (iArray(iIndex) .LT. 21474836) Then iSeedvec(i) = 1 iArray(iIndex) * 100 + iArray(iIndex + 1) / 1000000 iIndex = iIndex + 2 Else iSeedvec(i) = 1 iArray(iIndex) * 10 + iArray(iIndex + 1) / 10000000 iIndex = iIndex + 2 End If End Do C テンポラリ配列の開放 Deallocate(iArray) C -------------------------------------------------- C Step 2. C X(9) は下位2ビット(a(249), a(250))以外は、 C a(i) = a(i - 250) xor a(i - 103) で構成する。すべて 2^31 ビット C 正整数として考えてください。 C iTmp1: X(1) の 下位 29 ビット C iTmp2: X(6) の 下位 21 ビットを 8 ビット左シフト C iTmp3: X(5) の 上位 8 ビットを 23 ビット右シフト C iTmp4: X(9) の 下位 2 ビット C X(9) = LSHFT(iTmp1 xor (iTmp2 + iTmp3), 2) + iTmp4) iTmp1 = IAND(iSeedvec(1), Mask1) iTmp2 = ISHFT(IAND(iSeedvec(6), Mask3), 8) iTmp3 = ISHFT(IAND(iSeedvec(5), Mask2), -23) iTmp4 = IAND(iSeedvec(9), Mask4) iSeedvec(9) = IOR(ISHFT(IEOR(iTmp1, IOR(iTmp2, iTmp3)), 2), iTmp4) C -------------------------------------------------- C Step 3. C 残りの系列も a(i) = a(i - 250) xor a(i - 103) に基づいて C 構成する C iSeedvec(i) 'iTmp1' 'iTmp2' C X(i) = (L_SHFT(X(i - 8), 2) + R_SHFT(X(i - 9), 29)) C 'iTmp3' 'iTmp4' C xor (L_SHFT(X(i - 3), 10) + R_SHFT(X(i - 4), 21) mod 2^31) Do i = 10, iP iTmp1 = ISHFT(iSeedvec(i - 8), 2) iTmp2 = IOR(iTmp1, ISHFT(iSeedvec(i - 9), -29)) iTmp3 = ISHFT(iSeedvec(i - 3), 10) iTmp4 = IOR(iTmp3, ISHFT(iSeedvec(i - 4), -21)) iSeedvec(i) = IAND(Mask2, IEOR(iTmp2, iTmp4)) End Do iFlag = 1 iErr = 0 9999 Continue 2000 Format(10I8) End Subroutine C ------------------------ End of InitialSeedVector3 C -------------------------------------------------- C 線形合同法 X(i) = 16709 * X(i-1) mod 2^31-1 C の最上位ビットを使って1ビットずつ設定する C サブプログラム C -------------------------------------------------- Subroutine InitialSeedVector2(iSeedvec, iSeed, iErr) Parameter (iP = 250) Parameter (iQ = 103) Integer iSeed, iSeedTmp Integer iSeedvec(iP) Integer i, j Integer iTmp1, iTmp2, iTmp3, iTmp4, iTmp Integer Mask1, Mask2, Mask3, Mask4 Data Mask1, Mask2, Mask3, Mask4 1 /Z'1fffffff', Z'7fffffff',Z'001fffff',Z'00000003'/ Real*8 dA, dM, dInv_M, dTmpMult C 初期化済チェック用フラグ Integer iFlag Save iFlag C -------------------------------------------------- C 初期化済(iFlag == 1)の際は何もしないで終了する If (iFlag .EQ. 1) Then goto 9999 End If C -------------------------------------------------- C 乱数種の初期値 = iSeed は 1 以上の奇数でなければならない。 C 奇数ならば、そのまま初期値とし、偶数が入ったら、 C 31415927 を初期値とする If (iSeed .LT. 1) Then C 種が0以下ならデフォルトの初期値を使う。 iSeedvec(1) = 31415927 Goto 100 End If iTmp1 = MOD(iSeed, 2) If (iTmp1 .EQ. 0) Then C 種が偶数ならデフォルトの初期値を使う。 iSeedvec(1) = 31415927 Else iSeedvec(1) = iSeed End If 100 Continue C -------------------------------------------------- C Step 1. C はじめの 9 個は線形合同法( X(i) = 16709 * X(i-1) mod 2^31-1 ) C で生成。X(i) の上位1ビットごとに初期値列として生成 dA = 16807.0 dM = 2147483647.0 dInv_M = 1.0 / dM iTmp = iSeedvec(1) Do i = 2, 9 iSeedTmp = 0 Do j = 1, 31 dTmpMult = iTmp * dA iTmp = IDINT(dTmpMult - dM * IDINT(dTmpMult * dInv_M)) iSeedTmp = iSeedTmp * 2 + ISHFT(iTmp, -30) End Do iSeedvec(i) = iSeedTmp End Do C -------------------------------------------------- C Step 2. C X(9) は下位2ビット(a(249), a(250))以外は、 C a(i) = a(i - 250) xor a(i - 103) で構成する。すべて 2^31 ビット C 正整数として考えてください。 C iTmp1: X(1) の 下位 29 ビット C iTmp2: X(6) の 下位 21 ビットを 8 ビット左シフト C iTmp3: X(5) の 上位 8 ビットを 23 ビット右シフト C iTmp4: X(9) の 下位 2 ビット C X(9) = LSHFT(iTmp1 xor (iTmp2 + iTmp3), 2) + iTmp4) iTmp1 = IAND(iSeedvec(1), Mask1) iTmp2 = ISHFT(IAND(iSeedvec(6), Mask3), 8) iTmp3 = ISHFT(IAND(iSeedvec(5), Mask2), -23) iTmp4 = IAND(iSeedvec(9), Mask4) iSeedvec(9) = IOR(ISHFT(IEOR(iTmp1, IOR(iTmp2, iTmp3)), 2), iTmp4) C -------------------------------------------------- C Step 3. C 残りの系列も a(i) = a(i - 250) xor a(i - 103) に基づいて C 構成する C iSeedvec(i) 'iTmp1' 'iTmp2' C X(i) = (L_SHFT(X(i - 8), 2) + R_SHFT(X(i - 9), 29)) C 'iTmp3' 'iTmp4' C xor (L_SHFT(X(i - 3), 10) + R_SHFT(X(i - 4), 21) mod 2^31) Do i = 10, iP iTmp1 = ISHFT(iSeedvec(i - 8), 2) iTmp2 = IOR(iTmp1, ISHFT(iSeedvec(i - 9), -29)) iTmp3 = ISHFT(iSeedvec(i - 3), 10) iTmp4 = IOR(iTmp3, ISHFT(iSeedvec(i - 4), -21)) iSeedvec(i) = IAND(Mask2, IEOR(iTmp2, iTmp4)) End Do iFlag = 1 iErr = 0 9999 Continue End Subroutine C ------------------------ End of InitialSeedVector2 C -------------------------------------------------- C -------------------------------------------------- C C 伏見の方法により、初期値ベクトルを生成する C [Ref. 乱数/伏見 1989] C C -------------------------------------------------- C -------------------------------------------------- Subroutine InitialSeedVector(iSeedvec, iSeed, iErr) Parameter (iP = 250) Parameter (iQ = 103) Integer iSeed Integer iSeedvec(iP) Integer i Integer iTmp1, iTmp2, iTmp3, iTmp4 Integer Mask1, Mask2, Mask3, Mask4 Data Mask1, Mask2, Mask3, Mask4 1 /Z'1fffffff', Z'7fffffff',Z'001fffff',Z'00000003'/ C 初期化済チェック用フラグ Integer iFlag Save iFlag C -------------------------------------------------- C 初期化済(iFlag == 1)の際は何もしないで終了する If (iFlag .EQ. 1) Then goto 9999 End If C -------------------------------------------------- C 乱数種の初期値 = iSeed は 1 以上の奇数でなければならない。 C 奇数ならば、そのまま初期値とし、偶数が入ったら、 C 31415927 を初期値とする If (iSeed .LT. 1) Then C 種が0以下ならデフォルトの初期値を使う。 iSeedvec(1) = 31415927 Goto 100 End If iTmp1 = MOD(iSeed, 2) If (iTmp1 .EQ. 0) Then C 種が偶数ならデフォルトの初期値を使う。 iSeedvec(1) = 31415927 Else iSeedvec(1) = iSeed End If 100 Continue C -------------------------------------------------- C Step 1. C はじめの 9 個は線形合同法( X(i+1) = 69069 * X(i - 1) mod 2^32, C で生成。X(i) の上位31ビットを初期値列として生成 Do i = 2, 9 iSeedvec(i) = ISHFT(69069 * iSeedvec(i - 1), -1) End Do C -------------------------------------------------- C Step 2. C X(9) は下位2ビット(a(249), a(250))以外は、 C a(i) = a(i - 250) xor a(i - 103) で構成する。すべて 2^31 ビット C 正整数として考えてください。 C iTmp1: X(1) の 下位 29 ビット C iTmp2: X(6) の 下位 21 ビットを 8 ビット左シフト C iTmp3: X(5) の 上位 8 ビットを 23 ビット右シフト C iTmp4: X(9) の 下位 3 ビット C X(9) = LSHFT(iTmp1 xor (iTmp2 + iTmp3), 2) + iTmp4) iTmp1 = IAND(iSeedvec(1), Mask1) iTmp2 = ISHFT(IAND(iSeedvec(6), Mask3), 8) iTmp3 = ISHFT(IAND(iSeedvec(5), Mask2), -23) iTmp4 = IAND(iSeedvec(9), Mask4) iSeedvec(9) = IOR(ISHFT(IEOR(iTmp1, IOR(iTmp2, iTmp3)), 2), iTmp4) C -------------------------------------------------- C Step 3. C 残りの系列も a(i) = a(i - 250) xor a(i - 103) に基づいて C 構成する C iSeedvec(i) 'iTmp1' 'iTmp2' C X(i) = (L_SHFT(X(i - 8), 2) + R_SHFT(X(i - 9), 29)) C 'iTmp3' 'iTmp4' C xor (L_SHFT(X(i - 3), 10) + R_SHFT(X(i - 4), 21) mod 2^31) Do i = 10, iP iTmp1 = ISHFT(iSeedvec(i - 8), 2) iTmp2 = IOR(iTmp1, ISHFT(iSeedvec(i - 9), -29)) iTmp3 = ISHFT(iSeedvec(i - 3), 10) iTmp4 = IOR(iTmp3, ISHFT(iSeedvec(i - 4), -21)) iSeedvec(i) = IAND(Mask2, IEOR(iTmp2, iTmp4)) End Do iFlag = 1 iErr = 0 9999 Continue End Subroutine C ------------------------ End of InitialSeedVector. C -------------------------------------------------- C -------------------------------------------------- C C ベクトル化したM系列乱数生成ルーチン C C -------------------------------------------------- C -------------------------------------------------- Subroutine VectorMsequence(iArray, n, iSeedvec, iErr) C iP, iQ : 原始多項式に現れる p と q Parameter (iP = 250) Parameter (iQ = 103) Integer iArray(n) Integer iSeedvec(iP) Integer i, iPk, iQk Integer iOffset C -------------------------------------------------- C Step 1. C 初期値ベクトル iSeedvec を乱数列 iArray の先頭に格納 Do i = 1, iP iArray(i) = iSeedvec(i) End Do C -------------------------------------------------- C Step 2. C iPO 個の乱数を通常の方法で生成する C iP0 は通常、iP * 16 位で取るみたいです。[Ref. Random Number C Generator on a Vector Processo / Ito, Kanada (1990)] iPk = iP * 16 iQk = iQ * 16 Do i = iP + 1, iPk iArray(i) = IEOR(iArray(i - iP), iArray(i - iQ)) End Do C 次に生成される数のインデックス iOffset = i C -------------------------------------------------- C Step 3. C 最内ループ長を2倍ずつ伸ばしながら、乱数列を生成 10 Continue Do i = iOffset, iOffset + iQk - 1 iArray(i) = IEOR(iArray(i - iPk), iArray(i - iQk)) End Do iOffset = i If (n - i .GT. iQk) Then C -------------------------------------------------- C 残りの生成する数列がベクトル長より短くなったら、 C 最終のループ計算を実行 goto 20 End If If (i .GT. iPk * 2) Then C -------------------------------------------------- C 現在の生成数 i が 2iP_k を超えたら、ベクトル長 iQ_k を2倍に C 伸ばして生成を繰り返す iPk = iPk * 2 iQk = iQk * 2 End If goto 10 20 Continue C -------------------------------------------------- C Step 5. C 最後のループ時の処理 Do i = iOffset, n iArray(i) = IEOR(iArray(i - iPk), iArray(i - iQk)) End Do C -------------------------------------------------- C Step 6. C iArray(n) から先のM系列 iArray(n + 1 : n + iP) をiSeedvec に格納 Do i = 1, iQ iSeedvec(i) = IEOR(iArray(n - iP + i), iArray(n - iQ + i)) End Do Do i = iQ + 1, iP iSeedvec(i) = IEOR(iArray(n - iP + i), iSeedvec(i - iQ)) End Do iErr = 0 End Subroutine C -------------------------- End of VectorMsequence. C -------------------------------------------------- C -------------------------------------------------- C C 各PEが独立して乱数列を生成できるように、初期値ベクトルを C 構成するサブルーチン。ParallelMseq() から呼び出す。 C C -------------------------------------------------- C -------------------------------------------------- C Subroutine SetSeedLeapFlog(iSeedVector, iMyID, iNumPE, iFlag, 1 iErr) Integer iP, iQ, iErr, iMyID, iNumPE Parameter (iP = 250) Parameter (iQ = 103) Integer iFlag(2) Integer iSeedVector(iP) Integer, Allocatable, Dimension(:) :: iTmpVec Integer i C -------------------------------------------------- C 初期化部 C -------------------------------------------------- C -------------------------------------------------- C テンポラリ配列の割り付け。割り付けに失敗したら iErr を 1400 にして C サブルーチンを終了する Allocate(iTmpVec(1:iP * iNumPE), STAT=iErr) If (iErr .NE. 0) Then iErr = 1400 goto 9999 End If C -------------------------------------------------- C 初期値ベクトルを iP 個生成する 生成アルゴリズムは3通り。 If (iFlag(1) .EQ. 3) Then C 円周率の値を利用 Call InitialSeedVector3(iSeedVector, iFlag(2), iErr) Else If (iFlag(1) .EQ. 2) Then C 線形合同法 X(i) = 16807 X(i-1) mod 2^31 - 1 Call InitialSeedVector2(iSeedVector, iFlag(2), iErr) Else C 線形合同法 X(i) = 69069 X(i-1) mod 2^32 Call InitialSeedVector(iSeedVector, iFlag(2), iErr) End If C エラー処理 If (iErr .NE. 0) Then Goto 9999 End If C -------------------------------------------------- C テンポラリ配列に初期値ベクトルをコピー Do i = 1, iP iTmpVec(i) = iSeedVector(i) End Do C -------------------------------------------------- C Sequential に p * PE 台数分のM系列を生成し, テンポラリに C 格納 Do i = iP + 1, iP * iNumPE iTmpVec(i) = IEOR(iTmpVec(i - iQ), iTmpVec(i - iP)) End Do C -------------------------------------------------- C 自分のノードIDを起点とした index でPE台数分おきに C テンポラリの乱数列をとびとびに初期値ベクトルに格納 Do i = 1, iP iSeedVector(i) = iTmpVec(iMyID + 1 + (i - 1) * iNumPE) End Do C -------------------------------------------------- C テンポラリ配列の解放をして、終了 Deallocate(iTmpVec) 9999 Continue End Subroutine C -------------------------- End of SetSeedLeapFlog. C -------------------------------------------------- C -------------------------------------------------- C C M 系列を各PEで通信無しで独立に生成するサブルーチン C C -------------------------------------------------- C -------------------------------------------------- Subroutine ParallelMseq(iA, iN, iFlag, iErr) C iA : (1, 2^31 - 1) の非負整数値が格納される配列。 C iN : 配列 iA のサイズ。(iN >= 10000 を満たさないとエラーになる) C iFlag : 初期値生成アルゴリズム設定配列(長さ2). C 初期値生成法として選択できる設定法は以下の3つである. C 1) 乗算合同法 x_i = 69069 x_{i - 1} \pmod{2^{32}} C 2) 乗算合同法 $x_i = 16807 x_{i -1} \pmod{2^{31} - 1}$ C 3) 並列高精度多倍長演算ルーチンを利用して計算された円周率 C の値の一部を用いて設定する。 C - 要素 1 : 0 から 2 の整数値を与える. C 0 : 1) を利用 (デフォルト) C 1 : 2) を利用 C 2 : 3) を利用 C C - 要素 2 : 1)または2)を初期値の設定方法とした場合は、種として, C 正の奇数を与える. (デフォルトは 31415927) C 3) の場合は、利用したい部分の位置(小数点以下の桁数)を C 与える。(デフォルトは 1) C C iErr : エラーコード (正常終了は 0) include 'mppf.h' Integer iP, iQ Parameter (iP = 250) Parameter (iQ = 103) Integer iN Integer iA(iN) Integer iFlag(2) Integer iIsInitialized, iErr, iMyID, iNumPE Integer iSeedVector(iP) Save iSeedVector, iIsInitialized, iMyID, iNumPE C -------------------------------------------------- C 配列のサイズのチェックと初期化の必要の有無の調査 C -------------------------------------------------- If (iN .LT. 10000) Then C 配列のサイズは10000以上を仮定。実際は8000程度が下限。 C 条件が偽なら、エラーコードを返して、9999行に飛び終了。 iErr = 1100 goto 9999 End If If (iIsInitialized .EQ. 1) Then C 2回目以降のサブルーチンコールなら10行に飛び、M系列を生成 goto 10 End If C -------------------------------------------------- C 1. 初期化部 C MPIの初期化確認、各プロセッサに初期値を与える C -------------------------------------------------- C -------------------------------------------------- C 1.1 MPI が初期化されているか否かをチェック call MPP_INITIALIZED(iMPIFlag, iErr) If (iMPIFlag .NE. 1) Then iErr = 1200 goto 9999 End If C -------------------------------------------------- C 1.2 MPI 環境情報の取得(自ノードID、 総プロセッサ数) call MPP_COMM_RANK( MPP_COMM_WORLD, iMyID, iErr ) call MPP_COMM_SIZE( MPP_COMM_WORLD, iNumPE, iErr ) C 1.2.1 総プロセッサ数のチェック iNumPE は2の巾乗のみ call Check_PE_Size(iNumPE, iErr) If (iErr .NE. 0) Then C 正しくなければ、そのままのエラーコードを与えて終了 goto 9999 End If C -------------------------------------------------- C 1.3 初期値ベクトルの作成(はじめてcallされるときだけ) call SetSeedLeapFlog(iSeedVector, iMyID, iNumPE, iFlag, iErr) If (iErr .NE. 0) Then goto 9999 End If C -------------------------------------------------- C -------------------------------------------------- C 2. M 系列の生成 C -------------------------------------------------- 10 Continue call VectorMsequence(iA, iN, iSeedVector, iErr) C -------------------------------------------------- C 3. 終了処理 C 初期化済のフラグ(iIsInitialized)を設定し, 正常終了の C エラーコード(0)を iErr を返す。 C -------------------------------------------------- iIsInitialized = 1 iErr = 0 9999 Continue End Subroutine C --------------------------------------------------