******************************* Documentation on BLUPF90 distribution ******************************* ******************************* Readme ******************************* =========================================================================== Beta versions of SPARSEM, BLUPF90, REMLF90 and related programs. Ignacy Misztal and collaborators, University of Georgia 05/11/99 - 12/03/99 =========================================================================== This distribution contains files necessary to compile the sparse matrix module SPARSEM, the solution program BLUPF90, the variance component program REMLF90, and related programs. The distribution has been prepared for various Unix systems and requires a Fortran 90 compiler. The distribution contains the following directories: notes - Files specific for various Unix systems and installation details libs - libraries DENSEM, SPARSEM, and general BLUP subroutines blup - blupf90 program reml - remlf90 program gibbs - gibbsf90 program thr - bivariate threshold model program cblup90 thr1 - as above with thresholds computed + multiple linear traits cblup90thr thr2 - as above plus quasi-REML cblup90reml bin - space for executable programs course - simple programs from I.Misztal's course that help understand the design of blupf90 and other programs examples to example7 - Various example data sets Version of blupf90 with much smaller memory requirements is available separately. The documentation is available partly separately as PDF files, and partly in specific directories as Readme files. The programs are still in very active development. See file Changes on recent changes. ================ Installation ================ Please read file Installation for details on installing on specific systems. =============== Bugs =============== Many bug fixes are documented in file Changes. The programs have been tested and used in many practical applications but they still contain errors. Please use them at your risk. Bug reports are gratefully appreciated as well as sufggestions on improvements or extensions. ======== Use ======== The programs are free for research use but please acknowledge them in publications. For commercial use please contact Ignacy Misztal at the address below. FSPAK is a part of SPARSEM and related programs, and a separate contribution is requested for its commercial use. Please see details on anonymous FTP at num.ads.uga.edu in directory fspak. Ignacy Misztal University of Georgia ignacy@uga.edu http://nce.ads.uga.edu/~ignacy ******************************* Installation ******************************* ========================================================================= Installation of SPARSEM, BLUPF90, REMLF90 and related programs. Ignacy Misztal, University of Georgia 05/11/99 ========================================================================= =========================== Compilation in Unix systems =========================== 1. Set file "Makeinit", which contains compiler and timing options for various Unix systems. This can be done by typing (one time only): make xx where xx is a combination of system/compiler. Systems tested include: combination system ------------------------------ sun Sun Ultra 1/ Sun F90 compiler sun-fujitsu Sun Ultra 1 / Fujitsu compiler dec DEC Alpha hp HP sgi SGI rs6000 RS/6000 linux-absoft Linux with Absoft f90 compiler If your system is not included, try make generic If this does not work, make necessary corrections to files Makeinit.generic second.f.generic and type: make generic 2. To compile all programs, type make all The executables should be in the "bins" directory. Please test them on examples included with the documentation on specific programs. ============================== Other make options ============================== To eliminate "core" files from all subdirectories, type: make clean To remove all executables, object files and modules, type make cleanall To recompile files in one subdirectory only, type make in that subdirectory. =================================== Debugging =================================== If compilation with a debugging code is needed, change options optf90 and optf77 in Makeinit to debugging options (-g is most popular), erase all objects by make cleanall and recompile by make all If debugging is for Fortran 90 programs, then only optf90 needs to be changed. If only a program in a subdirectory is causing problems, don't issue make cleanall but instead type make in that subdirectory. =================================== Different compiler / system =================================== Please look at the main Makefile and construct options for a specific compiler. The meaning of specific variables is: f90 - name of f90 compiler optf90 - options for Fortran 90 programs and subroutines optf77 - options for Fortrabn 77 subroutines optdir - compiler option for specifying location of modules libf77 - option for linking extra libraries, if needed. Each installation needs a timning function second. If a specific function is not known, a "dummy function" can be included. ======================== Fortran 90 compilers ======================== The programs have been developed initially on Sun Ultra with the Sun f90 compiler and recently under Linux with the Absoft compiler. Successful installations were on SGI, RS/6000, HP, and DEC computers. In almost all of these installation, numerous compiler bugs have been found and workarounds were required. This package does not work, try upgrading the compiler, or run it without the optimization. Please read specific comments in directory notes. =========================================== Compilation under DOS or Windows =========================================== The programs have been successfully compiled under DOS using the Lahey 3.5 compiler using the following procedure: 1. Copy all library program and one application program to a subdirectory. 2. Compile using automake (command AM). The command: make win can help with part 1 by combining many library programs into larger units in directory mswin. If the automake does not correctly detect modules that need to be compiled first, compilations of these modules needs to be done separately using the command line. For compilers other than ABSOFT, enclosed make files may be modified. This will include changinng names of object and executable files. On most WIN/DOS systems, object files end with .obj, and executable files end with .exe. =============== Bugs =============== The programs have been tested and used in many practical applications but they still contain errors. Please use them at your risk. I would greatly appreciate bug reports and information on improvements or extensions. ======== Use ======== The programs are free for research use but please acknowledge them in publications. For commercial use please contact the author. FSPAK is a part of SPARSEM and related programs, and a separate contribution is requested for its commercial use. Please see details on anonymous FTP at num.ads.uga.edu in directory fspak. Ignacy Misztal University of Georgia ignacy@uga.edu ******************************* libs/Readme ******************************* ***************************************** * Library programs and modules * Ignacy Misztal, University of Georgia * 12-03-1999 ***************************************** Programs in this directory are used in application programs. Most of them are organized into libraries. Some programs are documented in separate documents, and in some the documentation is inside the programs. To compile all units under Unix, type: make To remove all compiled programs and libraries,for example to be able to recompile all the programs from scratch, type: make clean Descriptions ============= Library sparsem.a ----------------- Functionality: Contains sparse matrix module (module sparsem) with definitions of sparse matrix structures and basic operations on them, and sparse matrix factorization, inversion and multiplication routines (module sparseop). Both are described separately. Also included is an iterative solve by preconditioned-conjugate gradient (module pcg). Components: kind.f90 Precision definitions (module kind) sparse.f90 Main sparsem program (module sparsem) sparse2.f Helper subroutines for sparse.f90, incl. hash and sort fspak90.f90 F90 interface to FSPAK + multiplication routines fspak.o Sparse matrix factorization, inversion etc. by Perez-Enciso, Misztal and Elzo fspaksub.f Helper subroutines for fspak.f sparssub.f Low level sparse-matrix code by George and Liu pcg.f90 Iterative solvers by preconditioned-conjugate gradient for sparse matrix structure (module pcg) second.f Timing function, specific to each system Library denseop.a ----------------- Functionality: Contains the dense matrix module (module denseop), described in a separate document. Also contains module kinds, which holds definitions of numerical precisions that are used by most application and library programs. The module has been written by Tomasz Strabel Components: lapack90r.f90 Parts of LAPACK90 by Alan Miller kind.f90 Precision definitions (module kind) denseop.f90 The main body of the module Library blupsubs.a ------------------ Functionality: Helper subroutines for mixed-model programs. Components: blups1.f90 User interface for the BLUPF90 family of programs blups2.f Old Fortran 77 programs for splitting character lines into tokens and combining tokens into lines; also for determining residuals of the inverse of the parental dominance matrix. ginv.f old but fast generalized inverse for symmetric matrices Module model.f90 ---------------- Functionality: Contains definitions of all variables needed to describe a mixed model (module model); used by all mixed-model applications Use of libraries and modules ============================ To use a module, include it in the "use" statement in the application program. Compile the aplication program with an option specifying this directory as containing modules. Finally, link with the appropriate library. Example ------- Let program abc.f90 be in directory abc and assume that the libraries, already compiled, are in the directory libs. To use modules "model" and "sparsem", include the following in program abc.f90 program abc use model; use sparsem ... end program abc and compile f90 -P../libs abc.f90 ../sparselib.a ../libs/model.o The above line assumes that the compiler is called "f90" and that the include option in the compiler is "-P". The compilation us usually automated by make files. For examples, see Makefiles in directory blup or gibbs. ******************************* gibbs/Readme ******************************* Module Gibbs and Program GibbsF90 Ignacy Misztal, University of Georgia 04/29/99-06/07/99 Introduction ============ Module Gibbs is a collection of subroutines useful for implementing the Gibbs sampler predominantly in the BLUPF90 environment. It uses the features of Fortran 90 to simplify programnming and high-level optimization to reduce running time, with simplicity being as important as efficiency. To understand the module fully, please read the documentation on SPARSEM and on BLUPF90. Program GIBBSF90 is a conversion of BLUPF90 to Gibbs sampling using module Gibbs. Flat priors are used for (co)variances. The only statistics given are estimates of averages ond SD of variance components after a given number of rounds for burn-in. Luis Varona, IRTA, Spain, provided low level routines and answered many questions instrumental to getting GIBBSF90 programs running. Subroutine link_hash was inspired by Urs Schnyder, ETH, Switzerland. Module Gibbs ============ LINK_HASH ---------- call link_hash_ija(xx,xx_ija) xx - matrix in sparse_hashm format xx_ija - matrix in sparse_ija format Matrix in sparse_hashm form can be rapidly constructed but cannot easily be used. Matrix in sparse_ija form can easily be used but cannot be set up directly. Conversion from sparse_hashm to sparse_ija as provided in module SPARSEM is relatively slow. Subroutine link_hash provides a fast conversion from sparse_hashm to sparse_ija format when xx with the same nonzero structure needs to be converted repeatedly. When this routine is called the first time, a link is created that points to equivalent locations in xx and xx_ija. During subsequent calls, it is assumed that the nonzero structure of xx and xx_ija remains the same, and the conversion is done rapidly by using the link to update the numerical values only. Solve_iterm_block ----------------- call solve_iterm_block(xx_ija,xy,sol,i,j,diag,op) xx_ija - Matrix in sparse_ija format xy - (r8) vector of right hand sides sol - (r8) vector of solutions i - integer value of first equation to solve j - integer value of last equation to solve diag - (r8) matrix of dimension j-i+1 x j-i+1 op - a character variable containing either 'solve' or "update' When op='solve', the subroutine solves a block of equations from i to j, and puts the diagonal part of the matrix xx_ija into a dense matrix diag. When op='update', the right hand side is updated for the current block of solutions. The "update' option should always follow the "solve" option or the solutions will be incorrect. The need for the 'update' option arises from xx_ija being half-stored. Gen_normal ---------- x - gen_normal(mean,var,seed) mean - (r8) scalar or vector var - (r8) scalar or square matrix seed - an integer not equal to 0 x - (r8) scalar or square matrix Generates x=N(mean,V) when mean and var are scalars, or x=MVN(mean,V) when mean is a vector and V is a matrix. Gen_invwishart -------------- x=gen_invwishart(mean_q_form,df,seed2) mean_q_form - (r8) scalar or square matrix containing mean quadratic form df - an integer containing degrees of freedom seed - an integer not equal to 0 Generates samples from inverted chi square of inverted Wishart distributions. Example of Programming with module Gibbs ========================================= Program Gibbs use gibbs type (sparse_hashm)::xx type (sparse_ija)::sparse_ija .... call zerom(xx,neq) .... ! main loop do round=1,max xx%x(3,:) ! zero only numerical values to save setup time ... ! set up mixed model equations ... call link_hash_ija(xx,xx_ija) ! generate random samples for effects by blocks do i=1,neq,ntrait firsteq=i; lasteq=i+ntrait-1 call solve_iterm_block(xx_ija,xy,sol,firsteq,lasteq,diag,'solve') sol(firsteq:lasteq)=& gen_normal(sol(firsteq,lasteq,inverse(diag),seed1) call solve_iterm_block(xx_ija,xy,sol,firsteq,lasteq,diag,'update') enddo ... !calculate e'e and u'inv(a)u ... ................... ! Generate samples of (co)variances for random effects do i=1,nrandom df=df_random(i)-ntrait-1 !flat priors newG(i)=gen_invwishart(uAu(i)/df,df,seed2) enddo ! Generate samples of (co)variances for residual df=nrec-ntrait-1 newR=gen_invwishart(ee/df,df,seed2) enddo ... Bugs and Quirks =============== All random generators use the f90 uniform number generator called random_number. The quality of this generator is unknown. To change seeds, use the f90 procedure random_seed. Values of seeds in high level routines above are not used. Possible Speed Modifications for program GIBBSF90 ================================================= 1. Put data in memory or read data in binary, 2. Stop creating relationship matrices after 1st round and update the left hand side of MME directly using matrices AINV. Example for GIBBSF90 ==================== This is a three trait repeatability animal model with missing 2nd and third traits. The number of samples and the length of burn-in have been set to very low values for this example. In practice, the total number of samples needs to be 5,000 or larger, and the length of burn-in 1000-2000. Both numbers are model-dependent. See literature on Gibbs sampling for details on sampling. $cat exmr99sm # Three trait example from MTDFS, with missing traits DATAFILE renco99sm NUMBER_OF_TRAITS NUMBER_OF_EFFECTS OBSERVATION(S) 3 4 5 WEIGHT(S) EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT [EFFECT NESTED] 1 1 1 51 cross 2 2 2 2147 cross 2 2 2 2699 cross RANDOM_RESIDUAL VALUES 1825.258 909.822 842.392 909.822 2056.958 1309.826 842.392 1309.826 1716.444 RANDOM_GROUP RANDOM_TYPE diagonal FILES (CO)VARIANCES 2522.627 1901.389 1885.431 1901.389 1965.262 1808.874 1885.431 1808.874 1925.926 RANDOM_GROUP RANDOM_TYPE add_an_upg FILE renpe99s (CO)VARIANCES 613.747 -360.900 80.927 -360.900 475.232 319.516 80.927 319.516 666.600 $gibbsf90 name of parameter file?exmr99sm GIBBSF90 ver. 1.01 by I. Misztal Parameter file: exmr99sm Data file: renco99sm Number of Traits 3 Number of Effects 3 Position of Observations 3 4 5 Position of Weight (1) 0 Value of Missing Trait/Observation 0 EFFECTS # type position (2) levels [positions for nested] 1 cross-classified 1 1 1 51 2 cross-classified 2 2 2 2147 3 cross-classified 2 2 2 2699 Residual (co)variance Matrix 1825.3 909.82 842.39 909.82 2057.0 1309.8 842.39 1309.8 1716.4 Random Effect(s) 2 Type of Random Effect: diagonal trait effect (CO)VARIANCES 1 2 2522.6 1901.4 1885.4 2 2 1901.4 1965.3 1808.9 3 2 1885.4 1808.9 1925.9 Random Effect(s) 3 Type of Random Effect: additive animal with unknown parent groups Pedigree File: renpe99s trait effect (CO)VARIANCES 1 3 613.75 -360.90 80.927 2 3 -360.90 475.23 319.52 3 3 80.927 319.52 666.60 REMARKS (1) Weight position 0 means no weights utilized (2) Effect positions of 0 for some effects and traits means that such effects are missing for specified traits number of samples and length of burn-in 20 10 Data record length = 5 hash matrix increased from 5000 to 10000 hash matrix increased from 10000 to 20000 hash matrix increased from 20000 to 40000 hash matrix increased from 40000 to 80000 hash matrix increased from 80000 to 160000 hash matrix increased from 160000 to 320000 read 4539 records in 2.65000 s, 127095 nonzeroes hash matrix increased from 5000 to 10000 hash matrix increased from 10000 to 20000 read 2697 additive pedigrees finished peds in 3.64000 s, 191175 nonzeroes 0.296E+04 0.216E+04 0.222E+04 0.216E+04 0.224E+04 0.214E+04 0.222E+04 0.214E+04 0.236E+04 438. -247. 76.8 -247. 363. 272. 76.8 272. 585. 0.187E+04 942. 848. 942. 0.248E+04 0.182E+04 848. 0.182E+04 0.264E+04 round 1 0.307E+04 0.221E+04 0.226E+04 0.221E+04 0.227E+04 0.215E+04 0.226E+04 0.215E+04 0.235E+04 365. -200. 76.0 -200. 306. 234. 76.0 234. 532. 0.192E+04 954. 841. 954. 0.233E+04 0.162E+04 841. 0.162E+04 0.234E+04 round 2 0.291E+04 0.204E+04 0.209E+04 0.204E+04 0.208E+04 0.197E+04 0.209E+04 0.197E+04 0.218E+04 347. -193. 73.1 -193. 302. 233. 73.1 233. 527. 0.182E+04 902. 855. 902. 0.225E+04 0.155E+04 855. 0.155E+04 0.220E+04 .... .... round 19 0.307E+04 0.217E+04 0.218E+04 0.217E+04 0.202E+04 0.192E+04 0.218E+04 0.192E+04 0.210E+04 ave G 0.302E+04 0.214E+04 0.217E+04 0.214E+04 0.202E+04 0.193E+04 0.217E+04 0.193E+04 0.212E+04 SD G 71.6 44.5 37.2 44.5 38.6 35.7 37.2 35.7 38.4 281. -152. 76.2 -152. 252. 176. 76.2 176. 400. ave G 311. -169. 81.9 -169. 262. 179. 81.9 179. 419. SD G 27.3 12.4 11.1 12.4 6.56 8.89 11.1 8.89 12.6 0.184E+04 800. 731. 800. 0.196E+04 0.125E+04 731. 0.125E+04 0.167E+04 ave R 0.177E+04 824. 768. 824. 0.207E+04 0.134E+04 768. 0.134E+04 0.176E+04 SD R 42.1 37.5 22.5 37.5 58.7 50.8 22.5 50.8 45.4 round 20 solutions stored in file: "solutions" ******************************* thr/Readme ******************************* CBLUP90 - Bivariate threshold-linear model program CBLUP90 is the first approach to modifying blupf90 into a bivariate threshold-linear program. Models supported are all 2-trait models supported by blupf90, including those with missing traits. The modifications were done using the following sources: - Janns and Foulley paper in 1993 LPS, 33:183. - Dick Quaas'es notes, 1991. Currently there are many simplifying assumptions: - The first trait is categorical, second linear - Residual variance of the first trait on the underlying scale is 1 - Thresholds are not computed but programmed. To change thresholds, modify the program and recompile. Obtain approximate values of thresholds via univarriate threshold models. - No variances are calculated. Many optimizations are possible. For example, the CPU time can be cut a few times if the iterative solving routines use old solutions every round rather than by starting from 0. Goals: - compute thresholds, - add estimation of variance components, - create a large-data version for iteration on data. Directories thr1, thr2,... will contain future improvements to CBLUP90. Ignacy Misztal, 09/17/1998-04/30/99 Example ======== This example is from the paper by Foulley, Gianola and Thompson published in GSE 1993, 15:401. Parameter file exg2 -------------------- # Example of threshold model DATAFILE gse83dat NUMBER_OF_TRAITS NUMBER_OF_EFFECTS OBSERVATION(S) WEIGHT(S) EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT [EFFECT NESTED] 1 1 2 cross 3 3 2 cross 4 4 2 cross 2 2 6 cross RANDOM_RESIDUAL VALUES 1 1.2 1.2 25 RANDOM_GROUP RANDOM_TYPE diagonal FILE (CO)VARIANCES .05074 .098074 .098074 .9740 Data file gse83dat ------------------ 1 1 1 1 1 410 3280 1 1 1 1 1 375 3045 1 1 1 2 1 415 3548 1 1 2 2 1 400 3740 1 1 2 2 1 430 2850 1 1 2 2 1 420 3100 1 1 2 2 1 350 2700 2 1 1 2 1 460 3360 2 1 1 2 1 405 3740 2 1 2 2 1 390 3465 1 2 1 1 1 414 3465 1 2 1 1 2 430 3333 1 2 2 2 1 340 3465 1 2 2 1 2 470 3075 1 2 2 1 1 420 3150 2 2 2 1 1 445 3000 2 2 2 1 1 490 2730 1 3 1 1 1 416 3410 2 3 1 1 1 360 3000 2 3 1 2 1 427 3280 2 3 2 2 1 325 2660 2 3 2 2 1 444 3023 2 3 2 1 1 460 3200 1 4 2 1 2 470 3280 1 4 2 2 2 510 3440 1 4 2 2 1 390 2430 2 4 1 1 1 445 3333 1 5 1 1 1 405 3465 1 5 1 2 1 435 3045 1 5 2 1 1 425 3960 1 5 2 1 2 488 2850 1 5 2 1 1 385 3280 1 5 2 1 1 520 3465 1 5 2 2 1 480 3440 2 5 1 2 1 410 3570 2 5 1 1 2 505 3570 2 5 2 1 2 437 3178 2 5 2 1 2 510 3465 1 6 1 2 2 516 2900 1 6 1 1 2 453 2600 1 6 1 2 1 365 3570 1 6 2 1 1 505 3548 1 6 2 1 2 460 2730 1 6 2 1 1 450 3150 1 6 2 2 1 360 2565 2 6 1 2 1 435 3178 2 6 1 2 1 365 2900 Execution ---------- $cblup90 name of parameter file?exg2 CBLUP90 ver 1.02 by I. Misztal Parameter file: exg2 Data file: gse83dat Number of Traits 2 Number of Effects 4 Position of Observations 5 6 Position of Weight (1) 0 Value of Missing Trait/Observation 0 EFFECTS # type position (2) levels [positions for nested] 1 cross-classified 1 1 2 2 cross-classified 3 3 2 3 cross-classified 4 4 2 4 cross-classified 2 2 6 Residual (co)variance Matrix 1.0000 1.2000 1.2000 25.000 Random Effect(s) 4 Type of Random Effect: diagonal trait effect (CO)VARIANCES 1 4 0.50740E-01 0.98074E-01 2 4 0.98074E-01 .97400 REMARKS (1) Weight position 0 means no weights utilized (2) Effect positions of 0 for some effects and traits means that such effects are missing for specified traits Data record length = 6 original G .051 .098 .098 .974 inverted G 24.471 -2.464 -2.464 1.275 number of rounds? Bivariate Threshold Model, ver. 1.0 Hardwired for: Maximum Number of Categories: 5 Value of Threshold(s): .00 .53 .83 1.49 47 records finished peds in 3.000000E-02 s, 228 nonzeroes solution: -.94 415.92 -.96 422.34 -.02 -12.53 .00 .00 .34 32.12 -.00 .00 -.34 -3.48 -.20 -2.45 -.47 -4.85 .34 3.18 .42 5.09 .25 2.50 47 records finished peds in 5.000000E-02 s, 228 nonzeroes solution: -2.67 415.91 -2.70 422.33 .28 -12.53 .00 .00 1.20 32.11 .00 .00 -.30 -3.46 -.17 -2.46 -.46 -4.85 .31 3.14 .42 5.16 .20 2.47 47 records finished peds in 6.000000E-02 s, 228 nonzeroes solution: -4.17 415.92 -4.31 422.34 .51 -12.52 -.00 .00 2.52 32.10 .00 .00 -.34 -3.49 -.16 -2.45 -.48 -4.87 .33 3.16 .43 5.17 .21 2.48 47 records finished peds in 7.000000E-02 s, 228 nonzeroes solution: -4.96 415.92 -5.30 422.34 .44 -12.52 .00 .00 3.38 32.10 .00 -.00 -.35 -3.50 -.16 -2.45 -.49 -4.88 .36 3.18 .43 5.17 .21 2.48 47 records finished peds in 8.000000E-02 s, 228 nonzeroes solution: -5.09 415.92 -5.48 422.35 .40 -12.52 .00 .00 3.53 32.10 .00 -.00 -.35 -3.51 -.16 -2.45 -.49 -4.88 .37 3.19 .43 5.17 .21 2.48 solutions stored in file: "solutions" Solutions file --------------- trait/effect level solution 1 1 1 -5.0870 2 1 1 415.9205 1 1 2 -5.4778 2 1 2 422.3459 1 2 1 .4033 2 2 1 -12.5212 1 2 2 .0000 2 2 2 .0000 1 3 1 3.5290 2 3 1 32.1003 1 3 2 .0000 2 3 2 -.0000 1 4 1 -.3537 2 4 1 -3.5072 1 4 2 -.1615 2 4 2 -2.4520 1 4 3 -.4941 2 4 3 -4.8782 1 4 4 .3651 2 4 4 3.1897 1 4 5 .4313 2 4 5 5.1702 1 4 6 .2129 2 4 6 2.4776 ******************************* thr1/Readme ******************************* CBLUP90THR CBLUP90THR is a modification of CBLUP90 withh additional support for: - calculation of thresholds, - multiple linear traits, The modifications were based on the paper by Hoeschele, Tier and Graser published in 1995 JAS 73:1609, and were done by Benoit Auvray from the lab of Nicolas Gengler, FUSAGx, Gembloux, Belgium. Benoit Auvray's work on this program at the University of Georgia was supported by the GIFT EU committee. This program is still in the experimental phase. For assitance, please contact: Benoit Auvray, auvray.b@fsagx.ac.be or Ignacy Misztal ignacy@uga.edu Last modified 04/30/99 Example 1 ========== This example is a modified data set from Foulley and al. (198?). Parameter file ex1.par ----------------------- #Exemple of thershold model DATAFILE ex1.dat NUMBER_OF_TRAITS 3 NUMBER_OF_EFFECTS 4 OBSERVATION(S) 16 14 15 WEIGHT(S) EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT [EFFECT NESTED] 1 1 1 53 cross 5 5 5 2 cross 6 6 6 1 cross 10 0 10 1 cross RANDOM_RESIDUAL VALUES 1 15 3 15 726 29.8509 3 29.8509 22.078 RANDOM_GROUP 1 RANDOM_TYPE add_animal FILE ex1.ped (CO)VARIANCES 0.42857 -7.202 .6462 -7.202 484 15.198 .6462 15.198 3.9 Data file ex1.dat ------------------ 53 6 0 7 2 1 0 6 1 0 10 1 0 290.0 36.5 1 52 6 0 7 2 1 0 6 1 0 10 1 0 317.8 43.5 2 51 6 0 7 1 0 0 5 0 0 9 0 0 256.5 36.0 1 50 6 0 7 1 0 0 5 0 1 9 0 1 315.0 45.0 1 49 6 0 7 1 0 0 5 0 1 9 0 1 273.0 46.0 2 48 6 0 7 1 0 0 5 0 1 9 0 1 354.8 50.5 3 47 6 0 7 1 1 0 5 1 0 9 1 0 357.0 36.5 1 46 6 0 7 1 1 0 5 1 1 9 1 1 260.0 45.3 2 45 6 0 7 1 1 0 5 1 0 9 1 0 290.0 51.6 3 44 5 0 7 2 0 0 6 0 1 10 0 1 346.5 51.0 2 43 5 0 7 2 0 0 6 0 1 10 0 1 317.8 43.7 2 42 5 0 7 2 1 0 6 1 1 10 1 1 357.0 50.5 2 41 5 0 7 2 1 0 6 1 0 10 1 0 357.0 41.0 1 40 5 0 7 1 0 0 5 0 0 9 0 0 344.0 48.0 3 39 5 0 7 1 0 0 5 0 1 9 0 1 346.5 52.0 1 38 5 0 7 1 0 0 5 0 1 9 0 1 328.0 38.5 1 37 5 0 7 1 0 0 5 0 1 9 0 1 285.0 48.8 2 36 5 0 7 1 0 0 5 0 1 9 0 1 396.0 42.5 1 35 5 0 7 1 1 0 5 1 0 9 1 0 304.5 43.5 1 34 5 0 7 1 1 0 5 1 1 9 1 1 346.5 40.5 1 33 4 0 7 2 1 0 6 1 1 10 1 1 333.3 44.5 1 32 4 0 7 1 0 0 5 0 0 9 0 0 243.0 39.0 1 31 4 0 7 1 0 0 5 0 0 9 0 0 344.0 51.0 2 30 4 0 7 1 0 0 5 0 1 9 0 1 328.0 47.0 2 29 3 0 7 2 0 0 6 0 1 10 0 1 320.0 46.0 1 28 3 0 7 2 0 0 6 0 0 10 0 0 302.3 44.4 1 27 3 0 7 2 0 0 6 0 0 10 0 0 266.0 32.5 1 26 3 0 7 2 1 0 6 1 0 10 1 0 328.0 42.7 1 25 3 0 7 2 1 0 6 1 1 10 1 1 300.0 36.0 1 24 3 0 7 1 1 0 5 1 1 9 1 1 341.0 41.6 1 23 2 0 7 2 0 0 6 0 1 10 0 1 273.0 49.0 1 22 2 0 7 2 0 0 6 0 1 10 0 1 300.0 44.5 1 21 2 0 7 1 0 0 5 0 1 9 0 1 315.0 42.0 1 20 2 0 7 1 0 0 5 0 1 9 0 1 307.5 47.0 2 19 2 0 7 1 0 0 5 0 0 9 0 0 346.5 34.0 1 18 2 0 7 1 1 0 5 1 1 9 1 1 333.3 43.0 2 17 2 0 7 1 1 0 5 1 1 9 1 1 346.5 41.4 1 16 1 0 7 2 0 0 6 0 0 10 0 0 346.5 39.0 1 15 1 0 7 2 1 0 6 1 0 10 1 0 374.0 40.5 1 14 1 0 7 2 1 0 6 1 0 10 1 0 336.0 46.0 1 13 1 0 7 1 0 0 5 0 0 9 0 0 270.0 35.0 1 12 1 0 7 1 0 0 5 0 0 9 0 0 310.0 42.0 1 11 1 0 7 1 0 0 5 0 0 9 0 0 285.0 43.0 1 10 1 0 7 1 0 0 5 0 0 9 0 0 374.0 40.0 1 9 1 0 7 1 1 0 5 1 0 9 1 0 354.8 41.5 1 8 1 0 7 1 1 0 5 1 1 9 1 1 304.5 37.5 1 7 1 0 7 1 1 0 5 1 1 9 1 1 328.0 41.0 1 Pedigree file ex1.ped ---------------------- 53 6 0 52 6 0 51 6 0 50 6 0 49 6 0 48 6 0 47 6 0 46 6 0 45 6 0 44 5 0 43 5 0 42 5 0 41 5 0 40 5 0 39 5 0 38 5 0 37 5 0 36 5 0 35 5 0 34 5 0 33 4 0 32 4 0 31 4 0 30 4 0 29 3 0 28 3 0 27 3 0 26 3 0 25 3 0 24 3 0 23 2 0 22 2 0 21 2 0 20 2 0 19 2 0 18 2 0 17 2 0 16 1 0 15 1 0 14 1 0 13 1 0 12 1 0 11 1 0 10 1 0 9 1 0 8 1 0 7 1 0 6 0 0 5 0 0 4 0 0 3 0 0 2 0 0 1 0 0 Execution ---------- $ cblup90thr number of categories (numbered from 1)? 3 name of parameter file?ex1.par CBLUP90THR ver 1.0 by B. Auvray/I. Misztal Parameter file: ex1.par Data file: ex1.dat Number of Traits 3 Number of Effects 4 Position of Observations 16 14 15 Position of Weight (1) 0 Value of Missing Trait/Observation 0 EFFECTS # type position (2) levels [positions for nested] 1 cross-classified 1 1 1 53 2 cross-classified 5 5 5 2 3 cross-classified 6 6 6 1 4 cross-classified 10 0 10 1 Residual (co)variance Matrix 1.0000 15.000 3.0000 15.000 726.00 29.851 3.0000 29.851 22.078 Random Effect(s) 1 Type of Random Effect: additive animal Pedigree File: ex1.ped trait effect (CO)VARIANCES 1 1 .42857 -7.2020 .64620 2 1 -7.2020 484.00 15.198 3 1 .64620 15.198 3.9000 REMARKS (1) Weight position 0 means no weights utilized (2) Effect positions of 0 for some effects and traits means that such effects are missing for specified traits Data record length = 16 number of rounds? 10 G .429 -7.202 .646 -7.202 484.000 15.198 .646 15.198 3.900 Multivariate Threshold Model, ver. 1.0 47 records read 53 0 0 0 pedigrees finished peds in 8.000000E-02 s, 1982 nonzeroes solutions computed in 50 rounds of iteration reached convergence criterion 9.798506159770570E-006 round 1 thresholds: -.104577207169917 .583822238722925 G .429 -7.202 .646 -7.202 484.000 15.198 .646 15.198 3.900 47 records read 53 0 0 0 pedigrees finished peds in .210000 s, 1982 nonzeroes solutions computed in 50 rounds of iteration reached convergence criterion 5.320587274593940E-008 round 2 thresholds: 1.10275211935286 2.14983101505939 G .429 -7.202 .646 -7.202 484.000 15.198 .646 15.198 3.900 47 records read 53 0 0 0 pedigrees finished peds in .350000 s, 1982 nonzeroes solutions computed in 50 rounds of iteration reached convergence criterion 1.408987908274491E-009 round 3 thresholds: .830008028587783 2.06313685961232 G .429 -7.202 .646 -7.202 484.000 15.198 .646 15.198 3.900 47 records read 53 0 0 0 pedigrees finished peds in .900000 s, 1982 nonzeroes solutions computed in 3 rounds of iteration reached convergence criterion 9.345551199859900E-011 round 8 thresholds: .766802199974166 2.03990913282157 G .429 -7.202 .646 -7.202 484.000 15.198 .646 15.198 3.900 47 records read 53 0 0 0 pedigrees finished peds in .990000 s, 1982 nonzeroes solutions computed in 1 rounds of iteration reached convergence criterion 8.835166816357540E-011 round 9 thresholds: .766879992729754 2.03998703043051 G .429 -7.202 .646 -7.202 484.000 15.198 .646 15.198 3.900 47 records read 53 0 0 0 pedigrees finished peds in 1.07000 s, 1982 nonzeroes solutions computed in 1 rounds of iteration reached convergence criterion 7.838878061119310E-011 round 10 thresholds: .766975069550305 2.04007888555664 solutions stored in file: "solutions" Solutions file --------------- trait/effect level solution 1 1 1 -.2504 2 1 1 20.0153 3 1 1 1.2046 1 1 2 -.0526 2 1 2 -.2626 3 1 2 -.1812 1 1 3 -.0295 2 1 3 .2515 3 1 3 .2003 1 1 4 .0345 2 1 4 -4.0529 3 1 4 -.1752 1 1 5 -.2621 2 1 5 13.1412 3 1 5 -.1761 1 1 6 .5587 2 1 6 -29.0248 3 1 6 -.8723 1 1 7 -.0676 2 1 7 9.8071 3 1 7 .8779 1 1 8 .1467 2 1 8 2.9133 3 1 8 1.0492 1 1 9 -.3516 2 1 9 18.8020 3 1 9 .5976 1 1 10 -.5706 2 1 10 30.8871 3 1 10 .8030 1 1 11 .0514 2 1 11 .8861 3 1 11 .5591 1 1 12 -.1200 2 1 12 9.3260 3 1 12 .6388 1 1 13 .1703 2 1 13 -2.4967 3 1 13 .6424 1 1 14 -.2160 2 1 14 11.9884 3 1 14 .5417 1 1 15 -.3825 2 1 15 25.5904 3 1 15 1.0298 1 1 16 -.2882 2 1 16 22.3569 3 1 16 1.0875 1 1 17 -.2565 2 1 17 9.2179 3 1 17 -.0315 1 1 18 .1140 2 1 18 -9.9133 3 1 18 -.7147 1 1 19 -.3124 2 1 19 16.2794 3 1 19 .2737 1 1 20 .1555 2 1 20 -9.9122 3 1 20 -.4244 1 1 21 -.1277 2 1 21 4.1596 3 1 21 .0110 1 1 22 .0305 2 1 22 -1.2919 3 1 22 .1239 1 1 23 .1338 2 1 23 -9.8537 3 1 23 -.1442 1 1 24 -.2045 2 1 24 7.5507 3 1 24 .1436 1 1 25 .1967 2 1 25 -4.8225 3 1 25 .4896 1 1 26 -.0864 2 1 26 3.1099 3 1 26 .1641 1 1 27 .1465 2 1 27 -10.2159 3 1 27 -.1953 1 1 28 -.0333 2 1 28 -.1429 3 1 28 .0985 1 1 29 -.1517 2 1 29 5.6459 3 1 29 .2005 1 1 30 .0151 2 1 30 -3.3901 3 1 30 -.3494 1 1 31 -.2079 2 1 31 6.5915 3 1 31 -.1490 1 1 32 .4286 2 1 32 -20.5489 3 1 32 -.1020 1 1 33 -.1149 2 1 33 3.1669 3 1 33 -.0118 1 1 34 -.2446 2 1 34 13.6481 3 1 34 .1839 1 1 35 -.0536 2 1 35 -.3721 3 1 35 -.2832 1 1 36 -.7995 2 1 36 35.6717 3 1 36 .0614 1 1 37 .2780 2 1 37 -12.9555 3 1 37 -.3978 1 1 38 -.1725 2 1 38 13.4036 3 1 38 .4026 1 1 39 -.6090 2 1 39 23.4843 3 1 39 -.0153 1 1 40 .1661 2 1 40 -3.2380 3 1 40 -.3222 1 1 41 -.3029 2 1 41 17.3784 3 1 41 .2563 1 1 42 -.1337 2 1 42 5.7196 3 1 42 -.3011 1 1 43 .1955 2 1 43 -8.7404 3 1 43 -.6267 1 1 44 -.1585 2 1 44 7.9744 3 1 44 -.1909 1 1 45 .8580 2 1 45 -33.3675 3 1 45 -.1446 1 1 46 .7584 2 1 46 -40.1965 3 1 46 -1.0271 1 1 47 -.2063 2 1 47 3.8705 3 1 47 -.6002 1 1 48 .3451 2 1 48 -5.8407 3 1 48 .3324 1 1 49 .5498 2 1 49 -29.3260 3 1 49 -.8306 1 1 50 -.1209 2 1 50 -3.0774 3 1 50 -.6893 1 1 51 .4542 2 1 51 -23.9923 3 1 51 -.5663 1 1 52 .3452 2 1 52 -24.0415 3 1 52 -1.3160 1 1 53 .3692 2 1 53 -18.1925 3 1 53 -.3925 1 2 1 -1.2383 2 2 1 295.7527 3 2 1 31.1950 1 2 2 -1.3471 2 2 2 296.7625 3 2 2 33.2464 1 3 1 -.0342 2 3 1 15.6758 3 3 1 -.7914 1 4 1 -.2657 2 4 1 .0000 3 4 1 -1.8872 Example 2 ========== This example is a simulated data set (1000 observations,3 traits, first trait categorical (4 categories)). parameter file: ex2.par data file: ex2.dat pedigree file: ex2.ped Execution ---------- $ cblup90thr number of categories (numbered from 1)? 4 name of parameter file?ex2.par CBLUP90THR ver 1.0 by B. Auvray/I. Misztal Parameter file: ex2.par Data file: ex2.dat Number of Traits 3 Number of Effects 2 Position of Observations 3 4 5 Position of Weight (1) 0 Value of Missing Trait/Observation 0 EFFECTS # type position (2) levels [positions for nested] 1 cross-classified 2 2 2 50 2 cross-classified 1 1 1 1100 Residual (co)variance Matrix 1.0000 7.7460 7.7460 7.7460 1351.0 231.81 7.7460 231.81 1369.0 Random Effect(s) 2 Type of Random Effect: additive animal Pedigree File: ex2.ped trait effect (CO)VARIANCES 1 2 .66667 5.1640 5.1640 2 2 5.1640 1000.0 200.00 3 2 5.1640 200.00 1000.0 REMARKS (1) Weight position 0 means no weights utilized (2) Effect positions of 0 for some effects and traits means that such effects are missing for specified traits Data record length = 5 number of rounds? 10 G .667 5.164 5.164 5.164 1000.000 200.000 5.164 200.000 1000.000 Multivariate Threshold Model, ver. 1.0 hash matrix increased from 5000 to 10000 hash matrix increased from 10000 to 20000 hash matrix increased from 20000 to 40000 1000 records hash matrix increased from 40000 to 80000 hash matrix increased from 5000 to 10000 read 0 1100 pedigrees finished peds in 1.84000 s, 52002 nonzeroes solutions computed in 50 rounds of iteration reached convergence criterion 8.615375630517874E-008 round 1 thresholds: -.486352321325653 .225077497285302 .942505936509291 G .667 5.164 5.164 5.164 1000.000 200.000 5.164 200.000 1000.000 1000 records read 0 1100 pedigrees finished peds in 5.45000 s, 52002 nonzeroes solutions computed in 50 rounds of iteration reached convergence criterion 1.163042663602770E-009 round 2 thresholds: -.759987844887689 .204579230763401 1.12442986286466 G .667 5.164 5.164 5.164 1000.000 200.000 5.164 200.000 1000.000 1000 records read 0 1100 pedigrees finished peds in 9.02000 s, 52002 nonzeroes solutions computed in 41 rounds of iteration reached convergence criterion 9.651389369704020E-011 round 3 thresholds: -.872976604751947 .194634553181514 1.22120439357426 G .667 5.164 5.164 5.164 1000.000 200.000 5.164 200.000 1000.000 1000 records read 0 1100 pedigrees finished peds in 20.5900 s, 52002 nonzeroes solutions computed in 1 rounds of iteration reached convergence criterion 4.296803813044140E-011 round 8 thresholds: -.913602288524313 .192300843575052 1.25836992066762 G .667 5.164 5.164 5.164 1000.000 200.000 5.164 200.000 1000.000 1000 records read 0 1100 pedigrees finished peds in 22.4700 s, 52002 nonzeroes solutions computed in 1 rounds of iteration reached convergence criterion 3.242485327460530E-011 round 9 thresholds: -.913628221187042 .192350224059092 1.25846675157163 G .667 5.164 5.164 5.164 1000.000 200.000 5.164 200.000 1000.000 1000 records read 0 1100 pedigrees finished peds in 24.3300 s, 52002 nonzeroes solutions computed in 1 rounds of iteration reached convergence criterion 2.634563111772552E-011 round 10 thresholds: -.913622283751425 .192394816399360 1.25854033216819 solutions stored in file: "solutions" ******************************* Makefile ******************************* default: echo Makefile for building *F90 family of programs echo I. Misztal and collaborators @ University of Georgia echo echo Type make computer the first time to this Makefile is used to echo create options for specific computer/compilers. echo echo Combinations include sun, hp, rs6000, linux-absoft and generic echo echo To compile libraries type: make libs echo To compile everything, type: make all echo To remove all objects and executables, type make cleanall # initialization for native Sun compiler sun: echo f90=f90 >Makeinit echo optf90=-fast -fnonstop >>Makeinit echo optf77=-fast -fnonstop >>Makeinit echo optdir=-M >>Makeinit cp notes/second.f.unix libs/second.f # initialization for Fujitsu Sun compiler sun-fujitsu: echo f90=f90 >Makeinit echo optf90= -O3 -Am -X9 >>Makeinit echo optf77= -O3 -Am -X9 >>Makeinit echo optdir=-M >>Makeinit cp notes/second.f.unix libs/second.f linux-absoft: echo f90=f90 >Makeinit echo optf90=-O >>Makeinit echo optf77=-O >>Makeinit echo optdir=-p >>Makeinit echo libf77=-lU77 >>Makeinit cp notes/second.f.unix libs/second.f rs6000: echo f90=xlf90 >Makeinit echo optf90=-O -F:f90 -qcharlen=32767 >>Makeinit echo optf77=-O -qsource -qfixed >>Makeinit echo optdir=-I >>Makeinit cp notes/second.f.rs6000 libs/second.f dec: echo f90=f90 >Makeinit echo optf90=-O -fpe1 >>Makeinit echo optf77=-O -fpe1 >>Makeinit echo optdir="-module " >>Makeinit mv libs/fspak90.f90 libs/fspak90.f90.orig cp libs/fspak90.f90.dec libs/fspak90.f90 cp notes/second.f.unix libs/second.f # Should work for SGI and HP generic: echo f90=f90 >Makeinit echo optf90=-O >>Makeinit echo optf77=-O >>Makeinit echo optdir=-I >>Makeinit cp notes/second.f.null libs/second.f echo Note that the timing function always returns 0 libs: cd libs; make all: cd libs; make cd blup; make cd reml; make cd aireml; make cd gibbs; make cd thr; make cd thr1; make cd thr2; make clean: rm */core cleanall: rm */*.o */a.out bin/*90* */core */*.mod */*.a */*.M tar: tar czvf ../progsf90.tar.gz Readme Installation Makefile \ Changes libs/*.f* libs/Mak* blup/*.f90 */Read* */read*\ blup/Mak* reml/*.f90 reml/Mak* thr/*f90 thr/Mak* \ thr1/*f90 thr1/Mak* thr1/ex* thr2/*f90 thr2/Mak* thr2/ex* \ examples gibbs/*.f90 gibbs/*.f gibbs/Mak* gibbs/ex* notes course \ aireml/*f90 aireml/Makefile tarf: tar czvf ~/old/f90.for.gz.`date +%y%m%d` Mak* */Mak* */*.f* tarall: tar czvf ../allfsparse.o sparse2.o fspak90.o fspak.o fspaksub.o \ pcg.o sparssub.o second.o90.tar.gz * tarbin: tar czvf ../allf90.bin.tar.gz bin/remlf90 bin/blupf90 \ bin/cblup* bin/gibb* #generates combined programs so that they can be compiled under MS Windows win: cd libs; cat kind.f90 lapack90r.f90 denseop.f90 > ../mswin/denselib.f90 cd libs; cat sparse.f90 fspak90.f90 pcg.f90 iounf.f90 >../mswin/sparlib1.f90 cd libs; cat sparse2.f fspak.f fspaksub.f sparssub.f \ >../mswin/sparlib2.f cd libs; cat model.f90 blups1.f90 >../mswin/blups1.f90 cd libs; cat blups2.f ginv.f >../mswin/blups2.f cd blup; cat blupf90.f90 >../mswin/blupf90.f90 cd reml; cat remlf90.f90 >../mswin/remlf90.f90 cd notes ; cat second.f.null >../mswin/second.f ******************************* blup/Makefile ******************************* # Makefile for $(prog) include ../Makeinit prog=blupf90 dir= ../libs a.out: $(prog).o $(dir)/model.o $(dir)/blupsubs.a $(dir)/sparsem.a $(dir)/denseop.a $(f90) $(optf90) $(prog).o $(dir)/model.o $(dir)/blupsubs.a \ $(dir)/sparsem.a $(dir)/denseop.a $(libf77) cp a.out ../bin/$(prog) $(prog).o: $(prog).f90 $(dir)/model.o $(dir)/sparsem.a $(f90) $(optdir)$(dir) -c $(optf90) $(prog).f90 ******************************* libs/Makefile ******************************* # Makes libraries include ../Makeinit all: denseop.a sparsem.a model.o blupsubs.a iounf.o denseop.a: lapack90r.o kind.o denseop.o ar cr denseop.a kind.o denseop.o lapack90r.o sparsem.a: kind.o sparse.o sparse2.o fspak90.o fspak.o fspaksub.o \ pcg.o sparssub.o second.o ar cr sparsem.a kind.o sparse.o sparse2.o fspak90.o \ fspak.o fspaksub.o pcg.o sparssub.o second.o blupsubs.a: blups1.o blups2.o ginv.o ar cr blupsubs.a blups1.o blups2.o ginv.o kind.o: kind.f90 $(f90) -c $(optf90) kind.f90 model.o: model.f90 $(f90) -c $(optf90) model.f90 blups1.o: blups1.f90 $(f90) -c $(optf90) blups1.f90 blups2.o: blups2.f $(f90) -c $(optf77) blups2.f sparse.o: sparse.f90 $(f90) -c $(optf90) sparse.f90 sparse2.o: sparse2.f $(f90) -c $(optf77) sparse2.f fspak90.o: fspak90.f90 $(f90) -c $(optf90) fspak90.f90 fspak.o: fspak.f $(f90) -c $(optf77) fspak.f fspaksub.o: fspaksub.f $(f90) -c $(optf77) fspaksub.f ginv.o: ginv.f kind.f90 $(f90) -c $(optf77) ginv.f sparssub.o: sparssub.f $(f90) $(optf77) -c sparssub.f second.o: second.f $(f90) -c $(optf77) second.f denseop.o: denseop.f90 $(f90) -c $(optf90) denseop.f90 lapack90r.o: lapack90r.f90 $(f90) -c $(optf90) lapack90r.f90 pcg.o: pcg.f90 $(f90) -c $(optf90) pcg.f90 iounf.o: iounf.f90 $(f90) -c $(optf90) iounf.f90 clean: rm *.o *.a