************************************************************************* * * Changes in BLUPF90 and related programs * ************************************************************************* ----------------------------------------------------------------- Changes from Nov 2003 are inversely chronological (newest on top) ----------------------------------------------------------------- Sept 25, 2007 Too long fields that used readln, nums and nums2 generated error messages because bounds were not checked. Increased input line length from 200 to 400 charaters. Corrected statistics in renumf90 on covariables. Done by I. Aguillar July 3, 2007 Change of May 24 reverted due to slow hashing with sparse matrices. May 25, 2007 Two additional operations added to FSPAK and FSPAK90: mutiplication by a factor and solutions with the factor as LHS; programming by Juan Pablo Sanchez. May 24, 2007 Storing elements in sparse_ija formt when the matrix was relatively dense resulted in long computing time and large memory requirements. Prime numbers in hashvr() were changed for bigger numbers. Modifications per Andres Legarra. April 9, 2007 In search.f90, sorting large alpha fields did not work with some compilers, and positioning for the last element returned the element next to the last. February 26, 2007 Sorting in search.f90 returned incorrect results with multiple sorting fields. Feb 2, 2007 fspak-metis modified to remove an error message with a compiler. Aug 1, 2006 Renumf90 generated cryptic messages when UPG were specified in_pedigrees and PED_DEPTH was not zero. July 3, 2006 Gibbs1 and later programs returned incorrect results under IBM AIX. For a quick fix, rename file gibbsibm.f90 as gibbs.f90 in libs. May 11, 2006 RENUMF90 stops with an error message if UPG is specified and the position of year of birth is missing Feb 2, 2006 CBLUP90IOD modified for improved stability; all thresholds are now estimated; previously last threshold was set to 0. Dec 10, 2005. RENUMF90 now accepts a COMBINE keyword that allows a combination of fields to be treated as a new field. Sept 8, 2005 A new version of fspak.f is available to work with METIS, a rapid ordering software. For details, read Readme.fspak-metis in the libs directory. Feb 8, 2005 Subroutine iounf is now overloaded to read/write integer and real*8 vectors. Feb 7, 2004 Renumf90 modified for PED_DEPTH 3 and flexible output format allowing for > 9,999,999 levels; also memory leak eliminated. blup90iod modified to correctly print levels greater than 9,999,999. Jan 5, 2005 Added version of Postgibbsf90 for SUSE64; does not seem to work under SUSE32. December 21, 2004 Renumf90 ocassionally produced an error message "pedigrees repeated". October 28, 2004 Solutions for covariables missing for some traits were not zero with some compilers. Subroutine find_addresses were replaced in all programs except MRF90. May 20, 2004 Presence of weight variable caused residual variance to approach zero in all gibbs sampling programs. April 13 Overloaded assignments (e.g., xx_ija=xx) caused memory leak and subsequently "out of memory" errors with some compilers including Lahey. This was corrected by changing intent(out) to intent(inout) in subroutines implementing those assignments. February 20, 2004 Memory leak in gibbs2f90 and gibbs3f90 eliminated. Weights in airemlf90 now function correctly. Jan 24, 2004 Blup90iod was modified to implement Lidauer et al acceleration. Speed-up depends on p = number of effects x number of traits; speedup is small for small p but could be 5 times for p=100. Dec 3 Conversion from hash to ija fromats caused deallocation error messages mainly on HP and IBM/AIX workstations. Nullifying a pointer in sparse.f90 fixed the problem. Dec 2 Added to SPARSEM overloaded routines for copying matrix structures of the same kind. Overloaded subroutine RESET now checks whether pointers are associated before deallocation. If they are not associated, there is no deallocation. Nov 14 AIREMLRES added to the distribution ----------------------------------------------------------------- Changes below are in chronological order (newest on bottom) ----------------------------------------------------------------- May 12/99 GIBBSF90 worked incorrectly when some variances were not estimated and were set to 0. Fixed random number generators in rndgen.f90 so that variances of 0 do not return NaN. April-May 1999 Added matrix-scalar and matrix-vector multiplication to SPARSEM. Added printing of rectangular matrices in DENSEOP. Added a function diag that changes matrix to vector and vice versa in DENSEOP. Added subroutine pos_def that checks whether a matrix is semipositive- definite, and if not, corrects it; in DENSEOP. Added checking for semi-positive dispersion matrices in BLUPS1.F90; also additional checks in REMLF90. June-July, 1999 Fixed denseop.f90 so that it would work with the Fujitsu compiler on Sun workstations. Added AIREMLF90 by Shogo Tsuruta that makes REMLF90 converge 10-20 times faster; the current version does not work yet with correlated effects. Refined GIBBSF90 so that it seems to work with all problems except very small examples. September, 1999 Fixed a bug in GIBBSF90 where covariances set to zero in initial covariances matrices were still estimated. Module PCG added that computes solutions by Preconditioned Conjugate Gradient with diagonal and non-diagonal preconditioners. Made PCG the deafult solver in BLUPF90. . The convergence is now faster, by a factor of 3-50 depending on model. The whole package no longer works with real*4 precision, when rh=r4 in kind.f90. Fixes are planned. October 1999 GIBBSF90 no longer updates covariances that are 0 originally. Same correction needs to be applied to REMLF90i when covariances are 0 but all corresponding variances are not zero. November 1999 Fixed crashes with IOUNF ocuuring with certain compilers where preopened units had to be closed before they could be opened by IOUNF. November 22, 1999 (log)determinant in FSPAK is now computed correctly with matrices not of full rank. November 26, 1999 In DENSEM, the Cholesky subroutines now create lower diagonal matrices rather than symmetric matrices. Subroutine rndgen.f in directory GIBBS modified to account for this. January 26, 2000 Fixed crashes with addm() to large SPARSE_HASHM matrices February 1, 2000 On HP, pointer arrays may not be nullified on creation by the system, causing crashes. With two identical models except that the animal effect was first in the first model and was last in the second model, remlf90 work fine with the first one and crashed with the second one. Partially fixed by inserting init(ainv(i)) before call zero(ainv(i),n) in one but not all places. Nullifying pointers on creation is a de-facto but not the real standard. February 8, 2000 Subroutine gen_invwishart destroyed the variance parameter. February 24, 2000 Made subroutines reading and writing the parameter file part of module textop. Now, all programs using these subroutines need to have use textop next to other "use" statements. The module includes a new procedure nums2 that tokenizes a line into numbers and strings separated by spaces. As opposed to the previous version nums(), it can read real numbers. Also included is a subroutine that can read optional Parameters from the parameter file. The optional parameters are at the end of the parameter file and have the form: OPTION name str1 str2 ... Description of these subroutines can be found as comments to those subroutines in libs/blups1.f90 March 21 The likelihood in REMLF90 is now computed correctly and takes into account not only G0 but also A so that comparisons of different models are possible. It is not clear whether the likelihood constants are correct in models with unknown parent groups and models with different random effects for each trait. April 13 Likelihood in REMLF90 was incorrect with some combination of missing trait. Variance components were incorrect when weights were used. Nested regression did not work if their number was < number of effects or larger than any effect or trait column. April 24 Took the gibbs.f90 module out of directory gibbs and moved to directory libs. Unit prob.f90 used for calculation of probabilities/cdf in threshold models was made a module and also moved to directory libs. The two programs were arranged then logically, where all statistical routines were moved to prob.f90 and all gibss sampling manipulation routines were retained in gibbs.f90. Both programs are now compiled into library stat.a. Appropriate changes were made to makefiles and sources programs in directories gibbs and thr*. The probability routines were updated with some for truncated normal distributions useful for implementation of gibbs sampling in threshold univariate and bivariate models. May 2 Files prob.f90 and rndgen.f are now combined as prob.f90. Uniform random number generator has been improved. It is also overloaded to operate for integer arguments and varying intervals. Generator for truncated distribution has been generalized. Argument "seed" has been removed from all functions. Seed can be changed using a new function. The manual for module stat.a is Readme in directory gibbs. May 18 Residuals in REMLF90 became zero when trait values were very large. Possibly, message about R being not positive-definite was printed. One of the fixes involved changing the operational zero in module denseop from 10d-6 to 10d-10. If residuals in R are estimated as 0 in round 2, and observations have very large values, scale the observations so that they are smaller. June 12 Compilation of programs under Service Pack 4 of Absoft Linux compiler version 6.1 (or 6.2?) and older SGI generated an error message due to optimization bugs. This has been corrected by a slight change in the store_solutions subroutine. June 14 Under sire model, the hash matrix would enlarge when it is less than 10% full, resulting in very large memory requirements. September 1 quad_form for SPARSE_HASH format generated accesses to 0 indices, which caused crashes under DEC computers. In prob.f90, gen_uniform did not sample randomly from all 100 elements, and sometimes indices of 0 were generated causing random crashes. In fspak90.f90, second routine has been added an external attribute, which is required by some compilers. In Gibbs1f90, df_random was not zeroed which caused strange message under DEC for some data files. September 13 New Gibbs sampling program available in gibbsfast; It creates only single-trait matrices in memory and is therefore much faster and memory efficient than gibbsf90 for multiple-tarit and correlated-effect models. October 18 Number of missing records computed incorrectly in gibbs1f90. January 11,2001 Assignment of hash to ija variable caused changes to the hash variable declared as intent(in), resulting in error messages with some compilers. During conversion, the hash variable is no longer modified. Added option to addm(a,i,j,x,'f'), where x is the hash matrix, to store it in full. Conversion of this matrix to the ija format is also in full. January 22,2001 Added optional printing to a file in printmat in DENSEOP January 24, 2001 Functions gen_wishart and gen_invwishart had nonstandard dispersion argument: standarized and inverted. Now, the definittion of these functions is compatible with literature. To obtain samples with quadratic form ee and degrees of freedom df, use: r=gen_invwishart(finverse_s(ee),df) Applications using these functions: gibbsf90 and gibbs1f90, have been modified accordingly. Gibbs1f90 can store every n-th sample in file gibbs_samples. Degrees of freedom for the iverted wishart generators in the gibbs programs were slightly off in models with correlated effects March 8, 2001 Gibbs1f90 did not handle correctly split fixed effects Added gibbs2f90, which samples correlated effects jointly. This results in much faster mixing in models with correlated random effects such as maternal and random-regression. March 12, 2001 IN REML and Gibbs programs, variance-covariances matrices with zeroed rows and columns were found non semi-positive definite. This was dues to rounding errors in a routine which calculates eigenvalues. Now, a matrix is declared as non semi-positive definite when eigenvalues < -1e-15 rather than < 0. March 22, 2001 REMLF90 was very slow when the number of effects was high, especially with random regression models. March 26, 2001 Small variances in REMLF90 were unnecessarily corrected for negative eigenvalues; decreased 1d-4 in pos_def to 1d-10. POS_DEF corrected to ignore zero rows/columns combinations. May 11, 2001 POS-DEF did not compile with some compilers. Added keyword "recursive" in front of it. May 29, 2001 Gibbs family of programs did not work with Fujitsu compiler on Sun. This compiler does not initialize dynamic structures by default. Added explicit initializations to elements of ainv(). June 10, 2001 With some compilers (incl. Lahey and DEC), pedigree files were not read or were overwritten if RANDOM_GROUP started with efefct 5 or 6. This was caused by compiler bug/feature where unit 5 and 6 could not be open for file access. File allocation has been changed so that units 5 and 6 are no longer used. See blups1, the bottom of model.f90, and blupf90.f90 for details. August 20,2001 (Changes made by Shogo Tsuruta) In DENSEOP, cholesky decomposition, inversion and related subroutines/functions were made faster by replacing dot_product with do loops. Speedup on matrices of 18x18 is over 2 times. Gibbs2f90 speeded up by calculating the quadratic forms with IJA instead of HASHM matrices August 23, 2001 PRINTMAT printed text to standard output even if un was present. August 30,2000 BLUPF90 modified to correct non-positive definite matrices of variance components; negative eigenvalues are replaced by 1d-10. November 7, 2001 Corrected formats for variances so that negative covariances are not written as *******. Decemeber 3, 2001 Modified POS_DEF subroutine so that the minimum eigenvalue is calculated relative to the maximum eigenvalue, and all eigenvalues below a limit are corrected; see denseop manual for reference. As a result, many programs and particularly AIREMLF90 and GIBBS*f90 are more stable when (co)variance matrices are close to singularity. Makefile for Absoft verions 7.0 and later changed for higher optimization level and f95 comatibility. January 10, 2002 BLUPF90 with the default PCG solver produced inaccurate solutions when variances were large. February 22, 2002 Parameter files for BLUPF90 and BLUP90IOD can now accept directives to change default convergence criterion and teh maximum number of rounds; additionally, parameter file for BLUPF90 can accept directives for solving methodology and relaxation factor for SOR. March 11 Gibbsf90 now writes samples to a file so that its output can be analyzed by postgibbsf90. May 16 cblup90reml now works more reliably for a single-trait threshold model. function diag() now works also with single precision variables June 14 AIREMLF90 now works with maternal and random regression models (fixes by Tom Druet, INRA). June 21 cblup90reml works faster with single-trait threshold models. Sept 9 AIREMLF90 no longer generates function diag() now works also with single precision variables June 14 AIREMLF90 now works with maternal and random regression models (fixes by Tom Druet, INRA). June 21 cblup90reml works faster with single-trait threshold models. Sept 9 AIREMLF90 no longer generates a message "wrong description of nested covariable' Sept 20 Added Gibbs3f90, a version of gibbs2f90 with support for heterogeneous residual variances. Sept 24 In eigenvalue/eigenvector computations in denseop.f90, in rare instances nonzero columns were recognized as 0 (Thanks Jan-Thijs van Kaam). Oct 20 Nested covariables with nesting value of 0 do not work as missing value (Detected by Tom Druet). Fix not yet applied. A workaround would be to set covariables to 0 when nesting variable is also 0. Nov 14 Renumf90 printed some statistics as *** or NaN. ALso some correlatiosn reported as NaN. Remlf90 printed correlations for nonestimated covariances as NaN. Dec 9 In remlf90, some formats were too small; parameters for nesting were incorrect in some multiple-trait cases. Dec 12 RENUMF90 computed statistics using real numbers truncated to integers Dec 13 Nested covariables with nesting value of 0 should work correctly Feb 17, 2003 Pointer variables in sparse.f90 and fspak90.f90 nullified in declaration Added support of user-defined relationships to remlf90 March 24 Missing values in gibbs3f90 were not correctly predicted. (Thanks to Andres Legarra) April 16 kron function moved from gibbs2 and gibbs3 to libs/denseop.f90 June 16 S. Tsuruta aded AIC to remlf90 and airemlf90 July 28 S. Tsuruta added options for convergence criteria to remlf90 July 31 Computations with dominance or user defined matrix may be inaccurate if these effects are not last. Corrected in remlf90 only. Changes include replacing neff with eff in add_g_dom, add_g_usr. Oct 15 Read buffer in renumf90 was increased from 200 to 300 characters. Data lines with greater length lead to long lines incorrectly read.