C ----- mark.f ---- Program to recreate a history from the Markov matrix. C compile: g77 -g mark.f -o mark C usage: mark mark.output C where mark.in is the input matrix file and mark.output is the generated C sequence. C Copyright (C) 1999 SAE Fatigue Design and Evaluation Committee C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of the C license, or (at your option) any later version. C C This program is distributed in the hope tha it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General PUblic License for more details. C C You should have received a copy of the GNU General PUblic License C along with this program; if not, write to the Free Software C Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA C Try also their web site: http://www.gnu.org/copyleft/gpl.html C--------------------------------------------------------------------------- c C This program was originally written by Jay Fash at the Univ. of Illinois for C testing some of the FD&E committee histories. It was altered somewhat by C F.A.Conle at Ford Mo. to run on the Cyber and other machines, again for C testing some of the committee's specimens. C For an explanation of the theory see the paper by C Schuetz, W. and Hueck, M., "Generating the Falstaff Load History by C Digital Mini-Computers," 8th ICAF Symp., Lausanne, June 1975. C The earliest paper on the subject found so far appears to be: C Scherschel, R., "Description of a Peak-Valley Counting Program," C Ford Mo. Co. Fatigue Workshop VIII, SRL Lab, Dec. 1969. C History generation and test results are given in: C Fash, J.W., F.A.Conle, G.L.Minter, "Analysis of Irregular Loading C Histories for the SAE Biaxial Fatigue Program," Chapter 3, C Multiaxial Fatigue, AE-14 SAE, ISBN 0-89883-780-4, 1989 pp.33-59 C Other works are available. C-------------------------------------------------------------------------- C C On the Cyber there is a short pseudo-random number generator C that simulates what happens in this routine: MARKOV,ID=PV C Arrays used: C IQ(32,32) Markovian matrix of transitions (input file) C IOUT(32,32) matrix collected during regerneration: check on IQ C IT1(32) sum of row elements C IM1(32),IM2(32) next largest binary multiple C IA1(32),IA2(32) multiply constant for random number C IB1(32),IB2(32) addition constant for random number C ID1(32),ID2(32) diff between total no of random nos needed (IT1) C and total no generated (IM1) C IR1(32),IR2(32) random nos (R0) C C ICOUNT total no of transitions in history C ISTOP no. of transitions genrated C ISGN +1 if going to peak C -1 if going to valley C R0 current random no. C R1 revised random no. C R index for next level C I0 previous level C I1 current level C C RLDR MARKRG FORD.LB? C implicit real*8 (a-h,o-z) LOGICAL DEBUG Character*80 ICOM,argv integer*2 mbuf(100000) DIMENSION IOUT(32,32) DIMENSION IQ(32,32), & IM1(32),ID1(32),IA1(32),IB1(32),IR1(32), & IM2(32),ID2(32),IA2(32),IB2(32),IR2(32), & IT1(32),IT2(32) C C COMMON PT,PT1,IO,I1,ISGN,ICNT,ILINE,ICOUNT,R,ISTOP C COMMON/A/IM2,ID2,IA1,IB2,IR2 C COMMON/B/IM1,ID1,IA1,IB1,IR1 DEBUG=.FALSE. NTASK=1 1 WRITE(6,25) 25 FORMAT("# ======= MARKOV REGEN Program Run=======") C get the filter value if(iargc().eq.0)then write(6,26) 26 format('# mark: No filter value specified, assuming 0') IFILT=0 go to 29 endif if(iargc().gt.1)then write(0,*)"# mark ERROR: too many arguments in calling line" write(0,*)"# Usage: mark ifilt sequence_file" write(0,*)"# where ifilt is the filter size in bins: 0,1,2...32" write(0,*)"# Stopping..." stop endif c Get one arg call getarg(1,argv) read(argv,*)xIFILT IFILT=IFIX(xIFILT) 29 continue write(6,106)IFILT 106 format("#IFILT= ",I3) C Get data from file READ(5,60)ICOM write(6,61)ICOM 60 FORMAT(A80) 61 FORMAT("#",A80) C C Now read data in DO 74 I=1,32 READ(5,*)(IQ(I,J),J=1,16) 74 CONTINUE DO 76 I=1,32 READ(5,*)(IQ(I,J),J=17,32) 76 CONTINUE C data check C check sums of rows and cols ID=0 DO 180 I=1,32 ITOT1=0 ITOT2=0 DO 170 J=I,32 ITOT1=ITOT1+IQ(I,J) 170 CONTINUE DO 171 K=I,32 ITOT2=ITOT2+IQ(K,I) 171 CONTINUE ID=ID+ITOT1+ITOT2 IF(ITOT1.EQ.ITOT2)GO TO 180 C WRITE(0,165)I,ITOT1,ITOT2 WRITE(6,165)I,ITOT1,ITOT2 165 FORMAT("# TRANSITION TO AND FROM LEVEL",I4/ &"ARE NOT EQUAL, THEY ARE",I6,"AND",I6) write(0,*)"continuing anyway..." 180 CONTINUE C get total number of transition in history 190 NPTS=ID WRITE(0,195)NPTS WRITE(6,195)NPTS 195 FORMAT("# TOTAL NUMBER of pts IN HISTORY IS ",I6) C C The next section of code is from Fash's calc. s/r C It calculates params from the matrix and prints them C ==================================================== C ZERO LISTS DO 105 I=1,32 C "1" refs to trans from V to P (above main diagonal) IM1(I)=0 ID1(I)=0 IA1(I)=0 IB1(I)=0 IR1(I)=0 IT1(I)=0 C "2" refs to trans from P to V (below main diagonal) IM2(I)=0 ID2(I)=0 IA2(I)=0 IB2(I)=0 IR2(I)=0 IT2(I)=0 C zero the output matrix DO 105 J=1,32 IOUT(I,J)=0 105 CONTINUE C C Calculate sums of the rows (Note: main diag. all 0) DO 110 I=1,32 C Sum row to right of diag DO 115 J=I,32 115 IT1(I)=IT1(I)+IQ(I,J) DO 109 K=1,I 109 IT2(I)=IT2(I)+IQ(I,K) 110 CONTINUE C C Get partial values for initial random nos. C This is used to randomize the seed R0 C Later it will be added to ID1,2,... C But R0 will never get larger than "M" the period. DO 117 I=1,16 DO 117 J=I,16 117 IR1(I)=IR1(I)+IQ(I,J) DO 119 I=17,32 DO 119 J=17,I 119 IR2(I)=IR2(I)+IQ(I,J) C C Calc. the values of constants used to determine the next random no. AL102=DLOG10(2.0D0) DO 151 I=1,32 IF(IT1(I).EQ.0) GO TO 120 IM1(I)=1+IDINT(DLOG10(0.0000000001+(DBLE(IT1(I))))/AL102) IM1(I)=2**IM1(I) ID1(I)=IM1(I)-IT1(I) IR1(I)=ID1(I)+IR1(I) IF(IR1(I).EQ.IM1(I))IR1(I)=IR1(I)-1 IA1(I)=IM1(I)/2-3 IF(IA1(I).LT.5)IA1(I)=5 IB1(I)=IM1(I)/4-1 IF(IB1(I).LT.3)IB1(I)=3 C 120 IF(IT2(I).EQ.0) GO TO 151 IM2(I)=1+IDINT(DLOG10(0.0000000001+(DBLE(IT2(I))))/AL102) IM2(I)=2**IM2(I) ID2(I)=IM2(I)-IT2(I) IR2(I)=ID2(I)+IR2(I) IF(IR2(I).EQ.IM2(I))IR2(I)=IR2(I)-1 IA2(I)=IM2(I)/2-3 IF(IA2(I).LT.5)IA2(I)=5 IB2(I)=IM2(I)/4-1 IF(IB2(I).LT.3)IB2(I)=3 151 CONTINUE C C Print the stuff if(DEBUG)then WRITE(6,140) 140 FORMAT("# ARRAY PRINT") DO 141 I=1,32 141 WRITE(6,142)(IQ(I,J),J=1,16) 142 FORMAT("#",16I5) DO 144 I=1,32 144 WRITE(6,142)(IQ(I,J),J=17,32) C WRITE(6,147) C do you need to change the format anyway? 147 FORMAT("# IR2 IA2 IB2 IM2 ID2 IT2 ...") DO 150 I=1,32 150 WRITE(6,148)IR2(I),IA2(I),IB2(I),IM2(I),ID2(I),IT2(I), & IT1(I),IR1(I),IA1(I),IB1(I),IM1(I),ID1(I) 148 FORMAT("#",6(1X,I5),2X,6(1X,I5)) endif C total no of points ICOUNT=0 DO 174 I=1,32 DO 174 J=1,32 ICOUNT=ICOUNT+IQ(I,J) 174 CONTINUE WRITE(6,175)ICOUNT 175 FORMAT("# ICOUNT=",I6) C next section of code contains Fash's runhis s/r C===================================================== I0=16 I1=16 ISGN=-1 ISTOP=0 ICNT=0 ILINE=0 NSGN=-1 NTFL=0 NOX=16 C NSGN, NOX are used for filtering c IOUT(X,Y) should equal IQ(X,Y) when history gen is done 200 CONTINUE MBUF(1)=NOX NTFL=NTFL+1 write(6,*)NOX 204 ISGN=-ISGN ISTOP=ISTOP+1 IF(ISTOP.EQ.(ICOUNT+1)) GO TO 390 C C======================================== C This section is from random s/r IF(ISGN.LT.0) GO TO 240 C R0=IR1(I1) GO TO 210 205 R0=R1 210 R1=IA1(I1)*R0+IB1(I1) R1=R1-IM1(I1)*(DBLE(IDINT(R1/IM1(I1)))) R=R1-ID1(I1) IF(R.LT.0) GO TO 205 IR1(I1)=R1 GO TO 260 C 240 R0=IR2(I1) GO TO 250 245 R0=R1 250 R1=IA2(I1)*R0+IB2(I1) R1=R1-IM2(I1)*(DBLE(IDINT(R1/IM2(I1)))) R=R1-ID2(I1) IF(R.LT.0) GO TO 245 IR2(I1)=R1 260 CONTINUE C========== END OF "RANDOM" ==================== C========================= Begin of DOCALC================ IF(ISGN.LT.0) GO TO 350 C DO 310 J=I1,32 R=R-IQ(I1,J) IF(R.LT.0)GO TO 315 310 CONTINUE 315 I0=I1 I1=J GO TO 370 C 350 DO 355 J=1,I1 R=R-IQ(I1,I1-J+1) IF(R.LT.0)GO TO 360 355 CONTINUE 360 I0=I1 I1=I1-J+1 370 CONTINUE C=================== END OF S/R DOCALC ====================== C C C END OF "RUNHIS"CODE C NEXT LEVEL SELECTED IS I1 IF(ISTOP.EQ.1)WRITE(6,380)I0,I1 380 FORMAT("# THE 1ST TRANS. IS FROM ",I3, " TO ",I3) INOX=I1 IOUT(I0,I1)=IOUT(I0,I1)+1 C IF(IFILT.EQ.0)GO TO 3999 C C FILTER CODE NOT ENTERED. SEE LISTING IF DESIRED. C 3999 CONTINUE C NTFL COUNTS THE POINTS NTFL=NTFL+1 C INOX IS NOW NEXT REV. STORE IT WRITE(6,*)INOX MBUF(NTFL)=INOX C OLD CYBER PACKING CODE IS HERE,BUT NOT USED. C MAKE MBUFF BIG ENOUGH FOR WHOLE HIST INSTEAD GO TO 204 390 CONTINUE C END OF HISTORY C============ END OF "RUNHIS" S/R ====================== C Print the stuff if(DEBUG)then WRITE(6,440) 440 FORMAT("# Output ARRAY :") DO 441 I=1,32 441 WRITE(6,442)(IOUT(I,J),J=1,16) 442 FORMAT("#",16I5) write(6,*)" " DO 444 I=1,32 444 WRITE(6,442)(IOUT(I,J),J=17,32) endif STOP END