C======================================================================= PROGRAM mtdfARM C============================================================== C C PURPOSE : To incorporate Anonymous Relationship Matrix (ARM), such C asmarker based relationship C C STRATEGY C 1) Row and column number in A corresponds to ID C 2) Dynamic determining size of A by reading it twice. The C first is to determine the size and the second one C is to load it to a array. C 3) Then call inverse sub and save C C INPUT C C - UNIT "333" : file with external A, containing col: I J Value C C - UNIT "5" : TERMINAL OR BATCH CONTROL FILE : INTERACTIVE C INPUT IS REQUIRED, SPECIFYING THE MINIMUM AND C MAXIMUM ORIGINAL ANIMAL ID TO BE ENCOUNTERED C C OUTPUT C - UNIT "11" : NUMBER OF ANIMALS, INA, AND THEN VECTOR C OF ANIMAL IDS to 'MTDF11' C FORMATTED new code, original ID C C - UNIT "44" : .5*DET|A|, THEN NON-ZERO ELEMENTS OF NRM INVERSE, C WRITTEN OUT ELEMENT BY ELEMENT, IN ORDER OF C DIAGONALS FOLLOWED BY UNSORTED LOWER TRIANGLE C TO 'MTDF44' UNFORMATTED (BINARY) C C - UNIT "66" : PAPER TRAIL - A LOG FILE, MTDF56, NOT MTDF66 C WHICH IS RESERVED FOR MTDFPREP AND NOT C MTDF76 WHICH IS RESERVED FOR MTDFRUN C C C PROGRAM HISTORY: C WRITTEN: April 23 2006, Zhiwu Zhang C Revised: May 6 2006, Zhiwu Zhang (MTDFREML in "Y" shape) C Revised: July 10 2006, Zhiwu Zhang (Defining zero) C C C C C----------------------------------------------------------------------- C C PARAMETERS DETERMINING PROGRAM DIMENSIONS PARAMETER (MAXAN=5000) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*60 FILENAM CHARACTER*25 OLD,UNKN,FORMA,UNFOR,NEW INTEGER*4 end_of_file,rowMax,ic,ICOLL,IROWL,icc INTEGER*4 NRANK,N,IOPT,chooseZero real(8) COEF,DET,ZERO real(8), dimension(:), allocatable :: A real(8), dimension(:), allocatable :: V real(8), dimension(:), allocatable :: W integer(4), dimension(:), allocatable :: IFLAG C C----------------------------------------------------------------------- C Initialization C----------------------------------------------------------------------- ZERO=1.0D-8 IUN5 = 5 IUN6 = 6 IUN66 = 66 IUN11 = 11 IUN44 = 44 IUN333=333 FILENAM = 'MTDF56' NEW = 'NEW' OLD = 'OLD' UNKN = 'UNKNOWN' UNFOR = 'UNFORMATTED' FORMA = 'FORMATTED' OPEN(IUN66,FILE=FILENAM,FORM=FORMA,STATUS=UNKN) C C----------------------------------------------------------------------- C Collecting information C----------------------------------------------------------------------- WRITE(IUN66,*) ' FILE NAME FOR RELATIONSHIP FILE ?' WRITE(IUN6,*) ' FILE NAME FOR RELATIONSHIP FILE ?' READ(IUN5,509) FILENAM write(IUN6,*) FILENAM write(IUN66,*) FILENAM 509 FORMAT(A) C Define threshold of zero 801 FORMAT('By default, element below following is treated as zero') 802 FORMAT('Do you want change the default?') 803 FORMAT('No - type 0') 804 FORMAT('yes - type 1') 805 FORMAT('Please enter a new threshold: ') WRITE(IUN6,801) write(IUN6,*)zero WRITE(IUN6,802) WRITE(IUN6,803) WRITE(IUN6,804) C WRITE(IUN66,801) write(IUN66,*)zero WRITE(IUN66,802) WRITE(IUN66,803) WRITE(IUN66,804) C READ(IUN5,*) chooseZero if(chooseZero.ne.0) then WRITE(IUN6,805) READ(IUN5,*) Zero WRITE(IUN6,*)"Zero is defined below: ",Zero WRITE(IUN66,*)"Zero is defined below: ",Zero endif C CALL FCONCT(IUN333,FILENAM,'FORMATTED','OLD') OPEN(IUN333,FILE=FILENAM,FORM=FORMA,STATUS=OLD) CALL CTIME(IH,IM,IS,ID,66) CALL CTIME(IH,IM,IS,ID,6) WRITE(IUN66,*) FILENAM 901 FORMAT(/80(1H+)/) 902 FORMAT( 1X,'PROGRAM "MTDFNRM" - Importing A for "MTFRUN" and * recode IDs for "MTDFPREP"') 903 FORMAT(1X,'New function added in April, 2006') 904 FORMAT(1X,A,T50,'=',I15) 905 FORMAT(1X,A,T50,'= ',A) 906 FORMAT(1X,A,T50,' ',I15) 907 FORMAT(1X,A,T50,'=',F15.8) CALL BLANK(IUN66,1) CALL TIMDATE(IUN66) CALL BLANK(6,4) WRITE(IUN6,901) WRITE(IUN66,901) WRITE(IUN6,902) WRITE(IUN66,902) WRITE(IUN6,903) WRITE(IUN66,903) WRITE(IUN6,901) WRITE(IUN66,901) CALL BLANK(6,2) FILENAM = 'MTDF11' OPEN(IUN11,FILE=FILENAM,FORM='FORMATTED',STATUS='UNKNOWN') FILENAM = 'MTDF44' OPEN(IUN44,FILE=FILENAM,FORM='UNFORMATTED',STATUS='UNKNOWN') CLOSE(IUN44,STATUS='DELETE') OPEN(IUN44,FILE=FILENAM,FORM='UNFORMATTED',STATUS='NEW') C C rowMax=0 C----------------------------------------------------------------------- C Load A matrix C----------------------------------------------------------------------- C First read to get rowmax write(IUN6,*) "Import A matrix" IC=0 read(IUN333,*,iostat=end_of_file)ICOLL,IROWL,COEF do while(end_of_file >= 0) IC=IC+1 C write(IUN6,*)IC,ICOLL,IROWL,COEF if(ICOLL>rowMax)rowMax=ICOLL read(IUN333,*,iostat=end_of_file)ICOLL,IROWL,COEF enddo REWIND(IUN333) write(IUN66,*) "Order of A: ",rowMax write(IUN6,*) "Order of A: ",rowMax write(IUN6,*) "Allocating memory" allocate(a(rowMax*(rowMax+1)/2)) allocate(V(rowMax)) allocate(W(rowMax)) allocate(IFLAG(rowMax)) C Initial A to zeros C A=0 do i=1,rowMax*(rowMax+1)/2 A(i)=0 enddo read(IUN333,*,iostat=end_of_file)ICOLL,IROWL,COEF do while(end_of_file >= 0) if(abs(COEF)>zero) then A(IHMSSF(ICOLL,IROWL,rowMax))=COEF endif read(IUN333,*,iostat=end_of_file)ICOLL,IROWL,COEF enddo CLOSE(IUN333,STATUS='KEEP') CALL CTIME(IHH,IMM,ISS,IDD,66) CALL ETIME(IHH-IH,IMM-IM,ISS-IS,IDD-ID,66) CALL CTIME(IHH,IMM,ISS,IDD,6) CALL ETIME(IHH-IH,IMM-IM,ISS-IS,IDD-ID,6) C----------------------------------------------------------------------- C Get inverse C----------------------------------------------------------------------- write(IUN6,*) "Inverting A" C V=? C W=? C DET=? C ZERO=1.0D-8 C IFLAG=? C NRANK=? N=rowMax IOPT=1 C call DKMWHF(A,V,W,DET,ZERO,IFLAG,NRANK,N,IOPT) write(IUN6,*) "Exporting A inverse" DETA=DET MAXGH=0 DETAH = .5D0*DETA C WRITE(IUN6,*)DETAH,MAXGH WRITE(IUN66,*)"LOG DETERMINANT OF NRM: ",DETA WRITE(IUN6,*)"LOG DETERMINANT OF NRM: ",DETA WRITE(IUN66,*)"Rank of A: ",NRANK WRITE(IUN6,*)"Rank of A: ",NRANK if(NRANK