JSPFS - single-trait REML program for animal models using sparse matrix solver by Ignacy Misztal University Of Illinois 1207 W Gregory Dr., Urbana, IL 61801, USA E-mail: ignacy@uiuc.edu, phone: (217)244-3164 last correction: 12/12/1994 last improvement - matrix need not be of full rank & inversion errors corrected Disclaimer This program was developed for my own research. It is believed to be reasonably free of serious errors, but few attempts have been made to make it resistant to errors in parameter and input files. Please use it at your own risk. Distribution The program is free for research applications. Please contact the author for information on commercial use. In any case, please acknowledge the use of this program in publications. Introduction JSPFS calculates variance components in mixed models using an accelerated EM algorithm and a sparse matrix solver. It supports fixed, crossclassified, covariate, and random effects, one of which could be the animal effect with unknown parent groups. Ideas used in the program are described in the Journal of Dairy Science (1990) 73:163-172. Additionally, the current version supports the true sparse matrix inversion, where only a small subset of inverse elements are actually computed. Therefore, this version is 50+ times faster. Parameters Data and relationship files should be in the same format as in program JAA. The ordering should be strict in a sense that all factors should be numbered from 1 consecutively, and that there should be no "holes" in the ordering. Parameter file is similar to that in JAA, with the exception that the regressions are specified in the factor definition lines as: (number_of_field) -1 (any_value) When run interactively, the program prints current estimates and Aitken extrapolated projections of converged variance ratios. Then it asks for the number of additional rounds. Aitken extrapolation The extrapolation should be performed only when the extrapolated ratios stay constant to 3-4 decimal digits for 2-3 rounds. The extrapolation may work in 2-6 rounds with a single random factor, but may not be possible until round 10-30 with more such factors. Restrictions The animal effect should be last on the effect list, with unknown parent groups being on top of that effect. If there no groups, one groups should be defined, which will have no effect on the results. Installation The program consists of several modules: 1. jspfs.for main program 2. jspfs1.for subroutines including sparse matrix package FSPAK 3. second.for computer dependent timing routine, described in FSPAK.MAN Before compiling, determine program limits coded in the first few lines in jspfs.for. Memory requirements of the program are model and data dependent, and increase as the number of effects and the number of relationships increase. Parameter optim determines whether the inversion will be for minimum memory (optim=1) or maximum speed (optim=0). The latter option takes 2-3 times more memory but executes 20-300% faster, with the difference being particularily large on the Cray supercomputer. Limitations of JSPFS JSPFS does not compute maternal effects, covariables of order higher than 1, and it does not consider inbreeding. Example The first example was using an older version of JSPFS. The results in the new version should be virtually identical. (small example; data due to Yang Da) meishan.animal.uiuc.edu_51% jspfs name of a file with parameters, con - if console: ydp number of groups in the relationship matrix? 2 initializing 0.940000 filled 26 out of 50200 elements in hash table 1 filled 33 out of 35140 elements in hash table 2 0.970000s before sort 1.10000 after sort 1.11000 before copying to fspak format 1.11000 before odrv 1.11000 reordering: flag &time = 0 1.14000 symb. fact.: iflag & smpesp & time 0 255972 1.14000 numerical factorization+solution: 1.14000 used 48 words out of 256020 Inverse:iflag,time 0 1.17000 Inverse used 122 out of 256020 words trace 1.17000 quadratic form 1.17000 vround 1 se= 2.6310 s( 2)= 1.3333 alfa= 1.9732 ex= 0.0000 how many iterations more?, extrapolate (1)? 3 0 trace 1.20000 quadratic form 1.20000 vround 2 se= 2.6292 s( 2)= 1.3453 alfa= 1.9543 ex= -1.9090 trace 1.22000 quadratic form 1.22000 vround 3 se= 2.6280 s( 2)= 1.3540 alfa= 1.9410 ex= -1.9090 how many iterations more?, extrapolate (1)? 2 1 trace 1.24000 quadratic form 1.24000 vround 4 se= 2.6250 s( 2)= 1.3750 alfa= 1.9090 ex= 0.0000 trace 1.26000 quadratic form 1.28000 vround 5 se= 2.6250 s( 2)= 1.3750 alfa= 1.9091 ex= -1.9091 how many iterations more?, extrapolate (1)? 0 0 solutions were stored in file solydp meishan.animal.uiuc.edu_52% jspfs Larger example with over 3,300 equations and 3 effects /home/ignacy/scr/jspfs cat hx2 hxrec2 hxrel2 hxpp2 3 7 .75 1 2 70 0 4 1299 1 4 2148 -1 /scratch/ignacy/jspfs jspfs name of a file with parameters, con - if console: hx2 Optimizing for memory read 2135 pedigree records 13 unknown parent groups filled 14354 out of 50000 elements in hash table 1 filled 13376 out of 40000 elements in hash table 2 1.16903s smpa used 24828 out of 70000 words ra used 13014 out of 40000 words before odrv 1.33849 read ordering information from file ordhx2 symbolic factorization time= 1.65549 used 13.74% memory numerical factorization time= 1.92976 used 19.89% memory rank= 3516 sparse inversion time= 2.74978 used 24.96% memory trace 2.76635 vround 1 se= 339.0740 s( 2)= 478.9960 alfa= 0.7079 ex= 0.0000 vround 1 se= 339.0740 s( 3)= 475.1993 alfa= 0.7135 ex= 0.0000 how many iterations more?, extrapolate (1)? 5 0 trace 4.06120 vround 2 se= 289.2350 s( 2)= 531.6018 alfa= 0.5441 ex= 0.3350 vround 2 se= 289.2350 s( 3)= 520.2158 alfa= 0.5560 ex= -0.3634 trace 5.35421 vround 3 se= 262.1219 s( 2)= 569.3130 alfa= 0.4604 ex= 0.3731 vround 3 se= 262.1219 s( 3)= 547.7828 alfa= 0.4785 ex= -0.4036 trace 6.64708 vround 4 se= 249.3165 s( 2)= 593.0851 alfa= 0.4204 ex= 0.3836 vround 4 se= 249.3165 s( 3)= 560.3889 alfa= 0.4449 ex= -0.4191 trace 7.93995 vround 5 se= 243.7313 s( 2)= 607.2960 alfa= 0.4013 ex= 0.3841 vround 5 se= 243.7313 s( 3)= 563.4492 alfa= 0.4326 ex= -0.4254 trace 9.23269 vround 6 se= 241.3690 s( 2)= 616.2249 alfa= 0.3917 ex= 0.3818 vround 6 se= 241.3690 s( 3)= 561.6829 alfa= 0.4297 ex= -0.4289 how many iterations more?, extrapolate (1)? 5 1 trace 10.52829 vround 7 se= 239.4589 s( 2)= 630.2066 alfa= 0.3800 ex= 0.0000 vround 7 se= 239.4589 s( 3)= 552.1171 alfa= 0.4337 ex= 0.0000 trace 11.8212 vround 8 se= 239.5478 s( 2)= 633.9509 alfa= 0.3779 ex= 0.3925 vround 8 se= 239.5478 s( 3)= 547.0750 alfa= 0.4379 ex= -0.4633 trace 13.1171 vround 9 se= 239.5692 s( 2)= 637.6162 alfa= 0.3757 ex= 0.5084 vround 9 se= 239.5692 s( 3)= 542.4233 alfa= 0.4417 ex= -0.4811 trace 14.4099 vround 10 se= 239.5641 s( 2)= 641.1354 alfa= 0.3737 ex= 0.3107 vround 10 se= 239.5641 s( 3)= 538.0729 alfa= 0.4452 ex= -0.4996 trace 15.7054 vround 11 se= 239.5489 s( 2)= 644.4866 alfa= 0.3717 ex= 0.3345 vround 11 se= 239.5489 s( 3)= 533.9795 alfa= 0.4486 ex= -0.5134 how many iterations more?, extrapolate (1)? 5 1 trace 17.0001 vround 12 se= 238.4931 s( 2)= 707.7000 alfa= 0.3370 ex= 0.0000 vround 12 se= 238.4931 s( 3)= 461.1466 alfa= 0.5172 ex= 0.0000 trace 18.2936 vround 13 se= 238.9002 s( 2)= 706.7047 alfa= 0.3380 ex= 0.3388 vround 13 se= 238.9002 s( 3)= 460.6543 alfa= 0.5186 ex= -0.5195 trace 19.5948 vround 14 se= 239.0652 s( 2)= 706.2504 alfa= 0.3385 ex= 0.3388 vround 14 se= 239.0652 s( 3)= 460.5202 alfa= 0.5191 ex= -0.5194 trace 20.8870 vround 15 se= 239.1323 s( 2)= 706.0159 alfa= 0.3387 ex= 0.3389 vround 15 se= 239.1323 s( 3)= 460.5264 alfa= 0.5193 ex= -0.5193 trace 22.1826 vround 16 se= 239.1597 s( 2)= 705.8726 alfa= 0.3388 ex= 0.3389 vround 16 se= 239.1597 s( 3)= 460.5861 alfa= 0.5193 ex= -0.5193 how many iterations more?, extrapolate (1)? 5 1 trace 23.4767 vround 17 se= 239.1836 s( 2)= 705.6603 alfa= 0.3390 ex= 0.0000 vround 17 se= 239.1836 s( 3)= 460.7435 alfa= 0.5191 ex= 0.0000 trace 24.7694 vround 18 se= 239.1812 s( 2)= 705.5949 alfa= 0.3390 ex= 0.3389 vround 18 se= 239.1812 s( 3)= 460.8320 alfa= 0.5190 ex= -0.5185 trace 26.0658 vround 19 se= 239.1805 s( 2)= 705.5297 alfa= 0.3390 ex= 0.3386 vround 19 se= 239.1805 s( 3)= 460.9137 alfa= 0.5189 ex= -0.5182 trace 27.3596 vround 20 se= 239.1804 s( 2)= 705.4665 alfa= 0.3390 ex= 0.3500 vround 20 se= 239.1804 s( 3)= 460.9902 alfa= 0.5188 ex= -0.5178 trace 28.6565 vround 21 se= 239.1805 s( 2)= 705.4062 alfa= 0.3391 ex= 0.3399 vround 21 se= 239.1805 s( 3)= 461.0623 alfa= 0.5188 ex= -0.5176 how many iterations more?, extrapolate (1)? 5 1 trace 29.9519 vround 22 se= 239.2325 s( 2)= 704.0868 alfa= 0.3398 ex= 0.0000 vround 22 se= 239.2325 s( 3)= 462.4384 alfa= 0.5173 ex= 0.0000 trace 31.2480 vround 23 se= 239.2054 s( 2)= 704.1521 alfa= 0.3397 ex= 0.3397 vround 23 se= 239.2054 s( 3)= 462.4709 alfa= 0.5172 ex= -0.5172 trace 32.5412 vround 24 se= 239.1944 s( 2)= 704.1822 alfa= 0.3397 ex= 0.3397 vround 24 se= 239.1944 s( 3)= 462.4798 alfa= 0.5172 ex= -0.5172 trace 33.8350 vround 25 se= 239.1900 s( 2)= 704.1979 alfa= 0.3397 ex= 0.3397 vround 25 se= 239.1900 s( 3)= 462.4793 alfa= 0.5172 ex= -0.5172 trace 35.1287 vround 26 se= 239.1881 s( 2)= 704.2074 alfa= 0.3397 ex= 0.3396 vround 26 se= 239.1881 s( 3)= 462.4753 alfa= 0.5172 ex= -0.5172 how many iterations more?, extrapolate (1)? 0 0 solutions were stored in file solhx2